立方体中的最大内接正方形

mathe

楼上矩阵最优解要求各行元素平方和相等，各列元素绝对值和为1，并且任意不同两行之间正交。可以看出，对于任意这样的矩阵，交换任意两行两列或改变任意行列所有数据的符号，还是满足条件的。利用这些性质，至少对于3*4阶矩阵，我们可以穷举所有情况，其中非零数不超过9个，反复利用上面性质变换，最后可以得出总共只会余下四种:
下面所有字母代表[0,1]区间中实数:
\(\begin{pmatrix}-a & d& 1-g& 1\\b& -e&g & 0\\1-a-b&1-d-e&0&0\end{pmatrix}\)
对应方程:
\(\begin{cases}a(1-a-b)=d(1-d-e)\\b(1-a-b)=e(1-d-e)\\ab+de=(1-g)g\\a^2+d^2+(1-g)^2+1=b^2+e^2+g^2=(1-a-b)^2+(1-d-e)^2\end{cases}\)
这个方程的关键是前两式相加，得a+b=c+d或a+b+c+d=1,其中第二中解有负数出现，而第一种情况得b是方程$16*b^4 + 8*b^3 - 23*b^2 + 14*b - 2$的根，选择$e=b=0.20490155350665129314330190383234265145$得$a=d=0.082734498297453867827044842101191253555,g=0.96486035057909231410839945902474434917$
对应边长的平方为1.0149247893784870917727027004682468871


\(\begin{pmatrix} -a&d&g&0\\ b&-e&1-g&0\\1-a-b&1-d-e&0&1\end{pmatrix}\)
对应方程
\(\begin{cases}a(1-a-b)=d(1-d-e)\\b(1-a-b)=e(1-d-e)\\ab+de=g(1-g)\\a^2+d^2+g^2=b^2+e^2+(1-g)^2=(1-a-b)^2+(1-d-e)^2+1\end{cases}\)
这个方程解法和上面方程类似，没有得出更好的解

\(\begin{pmatrix}-a&d&1-f&0\\b&1-d&0&h\\c&0&f&h-1\end{pmatrix}\)
对应方程
\(\begin{cases}ab=d(1-d)\\bc=h(1-h)\\ac=f(1-f)\\a^2+d^2+(1-f)^2=b^2+(1-d)^2+h^2=c^2+f^2+(1-h)^2\\a+b+c=1\end{cases}\)
这个解出来结果比较可怕，比如h需要满足方程:
$(h-1)*h*(h^2-h+1)*(2h^2-2h-1)(9h^2-9h+1)(96*h^4 - 176*h^3 + 92*h^2 + 6*h - 9)
   (4*h^6 - 12*h^5 + 16*h^4 - 12*h^3 + 6*h^2 - 2*h + 1)
  (1458*h^12 - 8748*h^11 + 18711*h^10 - 13365*h^9 - 13797*h^8 + 39150*h^7 - 29122*h^6 - 9564*h^5 + 25793*h^4 - 11355*h^3 + 555*h^2 + 284*h + 8)
  (4608000*h^14 - 32486400*h^13 + 106297280*h^12 - 208927040*h^11 + 273352416*h^10 - 252211768*h^9 + 175199545*h^8 - 104473284*h^7 + 54392414*h^6 - 14730172*h^5 - 1523377*h^4 + 436852*h^3 + 334918*h^2 + 46456*h + 3721)
  (165888*h^16 - 1327104*h^15 + 4783104*h^14 - 10167552*h^13 + 13055616*h^12 - 7273152*h^11 - 6556608*h^10 + 17191200*h^9 - 12065016*h^8 - 5695176*h^7 + 16083319*h^6 - 10125760*h^5 + 1059958*h^4 + 1398650*h^3 - 611586*h^2 + 150880*h + 51253)
(5308416*h^26 - 69009408*h^25 + 525975552*h^24 - 21776007168*h^23 + 279629595648*h^22 - 1755045642240*h^21 + 6774140980224*h^20 - 17688419204736*h^19 + 32024536790592*h^18 - 38324963108544*h^17 + 23050469214944*h^16 + 11097784110800*h^15 - 38178355003336*h^14 + 32961273735048*h^13 - 2566180746576*h^12 - 19901614525832*h^11 + 16553725604058*h^10 - 2217221289542*h^9 - 3810021021401*h^8 + 1602353883447*h^7 + 439842772479*h^6 - 325579847011*h^5 - 32978746973*h^4 + 41942923279*h^3 - 3484899190*h^2 - 547935318*h + 357216)
(7644119040000*h^39 - 199860269875200*h^38 + 2437351880785920*h^37 - 18759902729011200*h^36 + 103074607169667072*h^35 - 435377615689838592*h^34 + 1491209724313731072*h^33 - 4298851157566006272*h^32 + 10644024801160240128*h^31 - 22725803111402993472*h^30 + 41473159626822259712*h^29 - 63493882655040338784*h^28 + 78994789689497720864*h^27 - 74719500922192891184*h^26 + 42940885445432024752*h^25 + 9534775436809759632*h^24 - 61502974905383268024*h^23 + 89004913392991538292*h^22 - 80304104536362644692*h^21 + 43080843348906315238*h^20 + 836842045447274914*h^19 - 29487275391337400990*h^18 + 34209282358889258953*h^17 - 22554620937493320006*h^16 + 8738926791071167945*h^15 - 1443514589525468923*h^14 + 251119191467260523*h^13 - 920986793774588347*h^12 + 995856852080848484*h^11 - 531495875439132419*h^10 + 157583299487498045*h^9 - 23090925236014258*h^8 + 443993273326087*h^7 + 309966779646301*h^6 - 44102946274893*h^5 - 5733815275467*h^4 + 3295423854898*h^3 + 188330030823*h^2 + 5250435188*h + 110766728)$
最终只有合法解
$9h^2-9h+1=0=>d=f=h => b=1/3,a=1/3$也不是最优解


\(\begin{pmatrix}0&f&f&0\\1-f&1-f&f-1&1-f\\ -f&0&0&f\end{pmatrix}\)
对应解满足$f^2-4f+2=0,f=0.58578643762690495119831127579030192143$边长平方为$2f^2=0.68629150101523960958649020632241537145$不是最优解

mathe

根据20#的公式，我们知道对于m阶情况，我们需要构造$m\times n$阶矩阵，每行范数相等，任意两行正交。而每列绝对值和不超过1。
但是对于最优解，必然每列绝对值和等于1，而且矩阵中所有非零数字数目不超过${(m+2)(m-1)}/2+n$个。
首先由于每列绝对值和不超过1，所以每列数字平方和不超过1，所以$m\times n$个数字平方和不超过$n$,于是每行元素平方和不超过$n/m$,
我们得出最大内接m维超立方体棱长不超过$\sqrt(n/m)$,其中只有n正好是m的倍数的时候才能够取到等号，这时要求每列正好一个数绝对值是1。
而另外一方面，我们总是可以选择$m[n/m]$列，分成$[n/m]$组，每组$m$列取成一个$m\times m$单位阵，余下$n-m[n/m]$列全部为0，得到一个非最优解可以取到棱长为$\sqrt([n/m])$的m维超立方体。所以最优解棱长必然大于$\sqrt([n/m])$而不超过$\sqrt(n/m)$。
而对于$m=3$情况，我们应该可以计算出最优解
我们假设某个最优解有k列不止一个元素绝对值大于0，那么对于$m=3,k<=5$
那么对于这k列构成$3\times k$阶的子矩阵，它的任意两行还是正交，但是任意两行数字平方和的差值是整数。
而$k=0$对应n是3的倍数时的最优值，已经讨论过了，所以只要讨论$1<=k<=5$的情况。
第一种情况:
这$3\times k$阶的子矩阵有一行全部是0，那么余下两行每列都是两个绝对值和为1的数，相当于只有k个变量，额外约束是这两行数平方和是整数，而且这两行正交，
  k个变量只有3个约束。由于全局最优解必然取边界条件，对于$k>3$时必然还要取更多边界条件直到$k=3$
  而对于$k=1$和$k=2$,容易看出这种情况无解。所以这种情况只能$k=3$
于是对应矩阵\(\begin{pmatrix}a&b&c\\1-a&1-b&c-1\\0&0&0\end{pmatrix}\)
其中由于$3/4<=a^2+b^2+c^2+(1-a)^2+(1-b)^2+(c-1)^2<3$,而且$a^2+b^2+c^2>=1,(1-a)^2+(1-b)^2+(c-1)^2>=1$,而且它们都是整数，所以只能
\(\begin{cases}a^2+b^2+c^2=1\\(1-a)^2+(1-b)^2+(c-1)^2=1\\a(1-a)+b(1-b)=c(1-c)\end{cases}\)
前两式相减得出$a+b+c=3/2$,得出a,b满足方程$8x^2-8x+1=0$,即$a=0.14644660940672623779957781894757548036,b=0.85355339059327376220042218105242451964,c=1/2$
对应$n=3h+4$时，可以得到棱长平方为$h+1={n-1}/3$
第二种情况:
这$3\times k$阶的子矩阵三行都有非零数。其中如果某一行只有一个非零数，由于这一行和其它行正交，我们可以得出这个非零数所在列其它数都是0，与每列至少两个非零数矛盾。
所以每行至少两个非零数。另外每两行正交，所以三行线性无关，得出$k>=3$.
2.1 如果k=5, 由于非零数比列出最多多5，所以只能每列都正好两个数。
   其中如果一行5个非零数，那么另外两行就分别3个和2个非零数，对应矩阵
\(\begin{pmatrix}a&b&c&d&e\\1-a&1-b&c-1&0&0\\0&0&0&1-d&e-1\end{pmatrix}\)
由于5个变量其中只有4个约束条件（每两行平方和之差是整数两个条件，两行之间两两正交还是两个条件），所以还得有个变量取边界值，所以退化到$k<5$的情况
  另外两种模式4+4+2和4+3+3都容易利用正交性得出无解。
比如4+4+2，说明2个非零元那行都和另两行都只有一列同时有非零数，和它们正交矛盾。4+3+3得出3个非零元的两行只有一列同时有非零数，和它们正交矛盾。
所以k=5淘汰
2.2 如果k=4
  那么最多9个非零数(不超过k+5)，最少8个非零数(不少于2k), 对应三行包含的非零数数目有
    4+3+2
     对应矩阵\(\begin{pmatrix}a&b&c&d\\1-a&1-b&e&0\\0&0&1-c-e&d-1\end{pmatrix}\)，和第二第三行正交矛盾，淘汰
    4+2+2
       对应矩阵\(\begin{pmatrix}a&b&c&d\\1-a&b-1&0&0\\0&0&1-c&d-1\end{pmatrix}\)
    3+3+3
      对应矩阵\(\begin{pmatrix}a&b&c&0\\d&b-1&0&e\\1-a-d&0&c-1&e-1\end{pmatrix}\)
    3+3+2
      对应矩阵\(\begin{pmatrix}a&b&c&0\\1-a&b-1&0&e\\0&0&c-1&d-1\end{pmatrix}\)，其中第三行和前两行都无法正交，淘汰。
  四种模式。
2.3 如果k=3
   那么最多8个非零数(不超过k+5),最少6个非零数(不小于2k),对应三行包含的非零数数目有
   3+3+2
     对应矩阵\(\begin{pmatrix}a&b&c\\d&e&c-1\\1-a-d&b+e-1&0\end{pmatrix}\)
   3+2+2
     对应矩阵\(\begin{pmatrix}a&b&c\\d&0&c-1\\1-a-d&b-1&0\end{pmatrix}\)，其中第二第三行无法正交，淘汰
   2+2+2
     对应矩阵\(\begin{pmatrix}a&b&0\\1-a&0&d\\0&b-1&1-d\end{pmatrix}\),其中各行之间无法正交，淘汰
三种模式

mathe

2.2 4+2+2模式对应矩阵
\(\begin{pmatrix}a&b&c&d\\1-a&b-1&0&0\\0&0&1-c&d-1\end{pmatrix}\)
有$a(1-a)=b(1-b),c(1-c)=d(1-d)$,所以$a=b$或$a=1-b$,$c=d$或$c=1-d$
于是有
\(\begin{pmatrix}a&a&c&c\\1-a&a-1&0&0\\0&0&1-c&c-1\end{pmatrix}\)\(\begin{pmatrix}a&1-a&c&c\\1-a&-a&0&0\\0&0&1-c&c-1\end{pmatrix}\)\(\begin{pmatrix}a&1-a&c&1-c\\1-a&-a&0&0\\0&0&1-c&-c\end{pmatrix}\)
其中后面两行平方和相等或差1，同样它们和第一行的平方和也相差整数
第一个矩阵如果后面两行平方和相等，有$a=c$,所以分别得出
$a^2+2*a-1/2=0, a=0.22474487139158904909864203735294569598,4a^2=0.20204102886728760721086370117643443213,2(1-a)^2=1.20204102886728760721086370117643443213$
    对应n=3h+5, 立方体棱长平方$h+1.20204102886728760721086370117643443213=(n-2)/3+0.20204102886728760721086370117643443213$
$a^2+2a-1=0, a=0.41421356237309504880168872420969807857,4a^2=2(1-a)^2=0.68629150101523960958649020632241537145$
   对应n=3h+4,立方体棱长平方$h+0.68629150101523960958649020632241537145=(n-1)/3-0.31370849898476039041350979367758462856$不是最优情况 (n=1(mod 3))
$a^2+2a-3/2=0,a=0.58113883008418966599944677221635926686,4a^2=1.3508893593264826720044258222691258651,2(1-a)^2=0.3508893593264826720044258222691258651$
   对应n=3h+2+4,不是最优情况
$a^2 + 2*a - 2=0,a=0.73205080756887729352744634150587236694,4a^2=2.1435935394489816517804292679530210644,2(1-a)^2=0.1435935394489816517804292679530210644$
   对应n=3h+4+4, 立方体棱长平方为$h+2.1435935394489816517804292679530210644=(n-2)/3+0.1435935394489816517804292679530210644$,不如前面方案 (n=2(mod 3))
$a^2 + 2*a - 5/2=0,a=0.87082869338697069279187436615827465088,4a^2=3.0333704529042344576650050707338027930,2(1-a)62=0.0333704529042344576650050707338027930$
   对应n=3h+6+4, 立方体棱长平方为$h+3.0333704529042344576650050707338027930=(n-1)/3+0.0333704529042344576650050707338027930$ (n=1(mod 3))
第二个矩阵比较前两行平方和相差$2c^2$是整数，只能$2c^2=1$,得出无解，同样第三个矩阵比较前两行平方和要求$c^2+(1-c)^2$是整数同样不行

mathe

2.3的   3+3+2模式
     对应矩阵\(\begin{pmatrix}a&b&c\\d&e&c-1\\1-a-d&b+e-1&0\end{pmatrix}\)
对应方程组部分为
a(1-a-d)=b(1-b-e)
d(1-a-d)=e(1-b-e)
其中两式相加得到$a+d=b+e$或$a+d+b+e=1$
如果$a+d=b+e$,还可以得到$a=b,d=e$, 最后得到合法解
   $8*c^4 - 32*c^3 + 45*c^2 - 30*c + 1=0$
   c=0.035139649420907685891600540975255650831
   a=0.082734498297453867827044842101191253550,
   d=0.20490155350665129314330190383234265145,
   对应棱长平方为$(n-1)/3 + 0.014924789378487091772702700468246887087$ (n=1(mod 3))
$64*c^8 - 256*c^7 + 304*c^6 - 32*c^5 - 327*c^4 + 296*c^3 - 152*c^2 - 32*c + 16=0$
   c=0.27558461705654999784812328936037690179
   a=0.17987838803224237835524423692465721565
   d=0.55492418539619261753733901460781193790
  对应棱长平方为$(n-2)/3+0.14065935012036122417723850465213714931$ (n=2(mod 3))
  $16*c^8 - 64*c^7 + 64*c^6 + 24*c^5 - 96*c^4 + 24*c^3 + c^2 - 22*c + 3=0$
   c=0.13860489686232080596276176669196028601
   a=0.074380948371657708078488878468073906554
   d=0.80258172315533576161036783428400273187
$n=3h+3+4,  r^2=(n-1)/3 + 0.030276368395549006911606670186456882565$ (n=1(mod 3) 而且$n>=7$)
不过$a+d+b+e=1$情况的计算已经很复杂了，计算下来好像无解
而余下一个2.2的3+3+3模式矩阵对应情况会更加复杂 (如果这个方程解决了，那么m=3时就完全解决了）
除了这种情况，其余情况的结果
$n -=0(mod 3)$， 最大棱长平方 $n/3$
$n -=2(mod 3)$,  已知最大棱长平方${n-2}/3+0.20204102886728760721086370117643443213$
$n -=1(mod 3)$,  
  当$n<=4$时,已知最大棱长平方${n-1}/3+0.014924789378487091772702700468246887087$,
  而在$n=7$时，已知最大棱长平方为${n-1}/3+0.030276368395549006911606670186456882565$
  而在$n>=10$时，已知最大棱长平方${n-1}/3+0.0333704529042344576650050707338027930$

mathe

esol(0,0)
得到
9*c^2 - 9*c + 1
a=0.12732200375003505059847105521145396076;b=0.12732200375003505059847105521145396076;c=0.12732200375003505059847105521145396076;x=0.33333333333333333333333333333333333333;y=0.33333333333333333333333333333333333334;z=0.33333333333333333333333333333333333333
9*c^2 - 9*c + 1
a=0.87267799624996494940152894478879859054;b=0.87267799624996494940152894478796257612;c=0.87267799624996494940152894478854603924;x=0.33333333333333333333333333333202867014;y=0.33333333333333333333333333333431385034;z=0.33333333333333333333333333333365747953， 这个解不够好
$1327104*c^22 - 14598144*c^21 + 78520320*c^20 - 277198848*c^19 + 718562304*c^18 - 1445202432*c^17 + 2327525376*c^16 - 3056990976*c^15 + 3306544768*c^14 - 2954054848*c^13 + 2169816440*c^12 - 1288672208*c^11 + 592924988*c^10 - 185582420*c^9 + 14218336*c^8 + 27499858*c^7 - 22520955*c^6 + 10772277*c^5 - 3524433*c^4 + 663759*c^3 - 19311*c^2 - 11715*c + 1300$
$a=0.61645575611595019564366164451969012446;b=0.99964861889131459720628119890601018649;c=0.99914013369073816425105409912190416390;x=0.87543900054754931571057569585193373024;y=0.0028868720701670442017409435891338646630;z=0.12167412738228364008768336055893240510$
\(\begin{pmatrix}x&a&0&c-1\\y&a-1&b&0\\z&0&b-1&c\end{pmatrix}\)
    这组很好，对应$n -=1(mod 3)$, 棱长的平方$r^2=(n-1)/3+0.1464118822982497217046378348507596815$ 但是好像是增根？有条件检查好像通过不了,后面一个解也类似
11943936*c^30 - 214990848*c^29 + 1827090432*c^28 - 9753716736*c^27 + 36696038400*c^26 - 103401232128*c^25 + 226155106176*c^24 - 392322766464*c^23 + 547110063040*c^22 - 620405629664*c^21 + 584627278416*c^20 - 485442765424*c^19 + 400839008496*c^18 - 365668024160*c^17 + 346843699576*c^16 - 294398468288*c^15 + 201109515718*c^14 - 103025529162*c^13 + 35508759545*c^12 - 4751415874*c^11 - 3359452963*c^10 + 3163768959*c^9 - 1546309579*c^8 + 462530161*c^7 - 66873410*c^6 + 3478243*c^5 - 908424*c^4 - 214127*c^3 - 2115*c^2 - 2700*c - 200
a=0.71883970112472172592436784012526963004;b=0.27384239774228909248921926723362136821;c=0.68740721959507614742307169636665297282;x=0.028075493203415081208096501315243155710;y=0.67911134810728449916208563638479935889;z=0.29281315868930041962981786229995748540
增根

esol(0,-1)得到
1327104*c^24 - 15925248*c^23 + 91348992*c^22 - 335812608*c^21 + 889150464*c^20 - 1797152256*c^19 + 2860263936*c^18 - 3639095040*c^17 + 3706622592*c^16 - 2969565504*c^15 + 1755586808*c^14 - 577548952*c^13 - 195332444*c^12 + 488517392*c^11 - 449148924*c^10 + 281428234*c^9 - 126051203*c^8 + 33066423*c^7 + 4604354*c^6 - 11159397*c^5 + 6784002*c^4 - 2240179*c^3 + 367313*c^2 - 35036*c + 500
a=0.83075415866767342164081349566679015338;b=0.95646500394930519144717044967819810047;c=0.97650580758675700622895686025787971476;x=0.51496528479610405439835147108507665117;y=0.11146405329919368235620643373084308116;z=0.37357066190470226324544209518408026768
增根

esol(-1,-2)得到
1679616*c^28 - 20155392*c^27 + 109921536*c^26 - 357011712*c^25 + 746013888*c^24 - 952093440*c^23 + 386020224*c^22 + 1241707968*c^21 - 3242195168*c^20 + 4094308096*c^19 - 2678362336*c^18 - 536662928*c^17 + 3707328032*c^16 - 5104501336*c^15 + 4379168340*c^14 - 2497810860*c^13 + 717860116*c^12 + 277914444*c^11 - 507176717*c^10 + 333946194*c^9 - 106430119*c^8 - 11392092*c^7 + 26972507*c^6 - 9778522*c^5 + 103708*c^4 + 742752*c^3 - 100736*c^2 - 20824*c + 4600
a=0.022395238513895352908459101701502169208;b=0.16167701978121523069178099777660409280;c=0.91540024756583824365232716591630919668;x=0.11184581172946412415880434687636092276;y=0.19574887487748029835324145420547932654;z=0.69240531339305557748795419891815975070
对应$n -=1 (mod 3), n>=7$,对应边长平方为${n-1}/3+0.020168150421438454409348065456978559301$
esol(-1,-1);esol(0,-2);esol(0,-3)没有找到合法解

1. reduce_by_z(f1,f2)={
   
2.    local(f3,f4,d1,d2,r,h);f3=f1;f4=f2;
   
3.    d1=poldegree(f3,z);d2=poldegree(f4,z);h=polcoeff(f4,d2,z);
   
4.    while(d1>d2,
   
5.         f3=f3-polcoeff(f3,d1,z)/h *z^(d1-d2)*f4;
   
6.         d1=poldegree(f3,z)
   
7.    );
   
8.    if(d1 == d2,
   
9.         f3 = f3-polcoeff(f3,d1,z)/h*f4;
   
10.    );
    
11.    f3
    
12. }
    
13. 
    
14. reduce_by_b(f1,f2)={
    
15.    local(f3,f4,d1,d2,r,h,df);f3=f1;f4=f2;
    
16.    d1=poldegree(f3,b);d2=poldegree(f4,b);h=polcoeff(f4,d2,b);
    
17.    while(d1>d2,
    
18.         df=polcoeff(f3,d1,b)/h;
    
19.         f3=denominator(df)*f3-numerator(df) *b^(d1-d2)*f4;
    
20.         d1=poldegree(f3,b)
    
21.    );
    
22.    if(d1 == d2,
    
23.         df=polcoeff(f3,d1,b)/h;
    
24.         f3 = denominator(df)*f3-numerator(df)*f4;
    
25.    );
    
26.    f3
    
27. }
    
28. 
    
29. esol(N1,N2)={
    
30. local(f1,f2,f3,f4,f5,f6,r1,r2,r3,r4,r5, tmp,equ,rts);
    
31. local(x0,y0,z0,a0,b0,c0);
    
32. f1 = x+y+z-1;
    
33. f2 = x*y-a*(1-a);
    
34. f3 = y*z-b*(1-b);
    
35. f4 = z*x-c*(1-c);
    
36. f5 = x^2+a^2+(1-c)^2-y^2-(1-a)^2-b^2-N1;
    
37. f6 = x^2+a^2+(1-c)^2-z^2-(1-b)^2-c^2-N2;
    
38. if(poldegree(f5,a) !=1, print ("Error in r1"),
    
39.     r1 = - polcoeff(f5,0,a)/polcoeff(f5,1,a);
    
40.     f1 = subst(f1, a, r1);
    
41.     f2 = subst(f2, a, r1);
    
42.     f3 = subst(f3, a, r1);
    
43.     f4 = subst(f4, a, r1);
    
44.     f6 = subst(f6, a, r1);
    
45.     if(poldegree(f1, x) != 1, print("Error in r2");print(f1),
    
46.         r2 = - polcoeff(f1, 0, x)/polcoeff(f1,1,x);
    
47.         f2 = subst(f2, x, r2);
    
48.         f3 = subst(f3, x, r2);
    
49.         f4 = subst(f4, x, r2);
    
50.         f6 = subst(f6, x, r2);
    
51.         if(poldegree(f3, y)!=1, print("Error in r3");print(f3),
    
52.               r3 = - polcoeff(f3, 0, y)/polcoeff(f3, 1, y);
    
53.               f2 = subst(f2, y, r3)*4*z;
    
54.               f4 = subst(f4, y, r3);
    
55.               f6 = subst(f6, y, r3)*4*z^2;
    
56.               f2 = reduce_by_z(f2,f4);
    
57.               f6 = reduce_by_z(f6,f4);
    
58.               if(poldegree(f2,z)!=1, print("Error in r4"),
    
59.                  r4 = - polcoeff(f2, 0, z)/polcoeff(f2, 1, z);
    
60.                  f4 = numerator(subst(f4, z, r4));
    
61.                  f6 = numerator(subst(f6, z, r4))-f4;
    
62.                  while(poldegree(f6, b)>0,
    
63.                       tmp = reduce_by_b(f4,f6);
    
64.                       f4=f6;f6=tmp
    
65.                  );
    
66.                  f6=component(f6,1);
    
67.                  if(poldegree(f4, b)!=1, print("Error in r5"),
    
68.                  r5 = -polcoeff(f4,0,b)/polcoeff(f4,1,b);
    
69.                  equ = factor(f6);
    
70.                  for(s=1, length(equ[,1]),
    
71.                     rts=polroots(equ[s,1]);
    
72.                     for(t=1,length(rts),
    
73.                          if(0< real(rts[t]) && real(rts[t])<1 && abs(imag(rts[t]))<0.00001,
    
74.                              c0=real(rts[t]);
    
75.                              b0=subst(r5,c, c0);
    
76.                              z0=subst(r4,c, c0); z0=subst(z0, b, b0);
    
77.                              y0=subst(r3,c,c0);y0=subst(y0,b,b0);y0=subst(y0,z,z0);
    
78.                              x0=subst(r2,c,c0);x0=subst(x0,b,b0);x0=subst(x0,z,z0);x0=subst(x0,y,y0);
    
79.                              a0=subst(r1,c,c0);a0=subst(a0,b,b0);a0=subst(a0,z,z0);
    
80.                             a0=subst(a0,y,y0);
    
81.                             a0=subst(a0,x,x0);
    
82.                              if(0<b0&&b0<1&&0<z0&&z0<1&&0<y0&&y0<1&&0<x0&&x0<1&&0<a0&&a0<1,
    
83.           print(equ[s,1]);                       print("a="a0";b="b0";c="c0";x="x0";y="y0";z="z0";x^2+a^2+(1-c)^2="x0^2+a0^2+(1-c0)^2)
    
84.                              )
    
85.                          )
    
86.                     )
    
87.                  )
    
88.                  )
    
89.               )
    
90.         )
    
91.     )
    
92. )
    
93. }
    

复制代码

mathe

为了确认上面的结果，我又对上面部分方程采用密集采样点加牛顿迭代法进行数值解验算，结果一致，这说明24#得出的结果对于m=3已经是最优解了
比如下面代码数值验算了25#的方程:

1. #include <stdio.h>
   
2. #include <math.h>
   
3. 
   
4. #define MAXSOL   50
   
5. #define N 6
   
6. double solbuf[MAXSOL][N];
   
7. int solc;
   
8. #define A(i,j) a[(i)*N+(j)]
   
9. void solve(double *a, double *b)
   
10. {
    
11.     int i,j,k;
    
12.     for(i=0;i<N;i++){//in i column
    
13.        double maxv = -1.0;
    
14.        int index=-1;
    
15.        for(j=i;j<N;j++){
    
16.           if(fabs(A(i,j))>maxv){
    
17.               maxv=fabs(A(i,j));
    
18.               index = j;
    
19.           }
    
20.        }
    
21.        if(index != i){//swap row index and i
    
22.          double tmp;
    
23.          for(j=i;j<N;j++){
    
24.              tmp = A(i,j);
    
25.              A(i,j)=A(index,j);
    
26.              A(index,j)=tmp;
    
27.          }
    
28.          tmp = b[i];
    
29.          b[i]= b[index];
    
30.          b[index]=tmp;
    
31.        }
    
32.        //divid by A(i,i) in row i
    
33.        for(j=i+1;j<N;j++)A(i,j)/=A(i,i);
    
34.        b[i]/=A(i,i);A(i,i)=1.0;
    
35.        for(j=0;j<N;j++){
    
36.           if(j==i)continue;
    
37.           for(k=i+1;k<N;k++){
    
38.              A(j,k)-=A(i,k)*A(j,i);
    
39.           }
    
40.           b[j]-=b[i]*A(j,i);
    
41.           A(j,i)=0.0;
    
42.        }
    
43.     }
    
44. }
    
45. 
    
46. void init_b(double *b, double x0, double y0, double z0, double a0, double b0, double c0, int N1, int N2)
    
47. {
    
48.     b[0]=x0+y0+z0-1.0;
    
49.     b[1]=x0*y0-a0*(1-a0);
    
50.     b[2]=y0*z0-b0*(1-b0);
    
51.     b[3]=z0*x0-c0*(1-c0);
    
52.     b[4]=x0*x0+a0*a0+(1-c0)*(1-c0)-y0*y0-(1-a0)*(1-a0)-b0*b0-N1;
    
53.     b[5]=x0*x0+a0*a0+(1-c0)*(1-c0)-z0*z0-(1-b0)*(1-b0)-c0*c0-N2;
    
54. }
    
55. 
    
56. void init_a(double *a, double x0, double y0, double z0, double a0, double b0, double c0,int N1, int N2)
    
57. {
    
58.     A(0,0)=1.0;A(0,1)=1.0;A(0,2)=1.0;A(0,3)=0.0;A(0,4)=0.0;A(0,5)=0.0;
    
59.     A(1,0)=y0;A(1,1)=x0;A(1,2)=0.0;A(1,3)=-1.0+2*a0;A(1,4)=0.0;A(1,5)=0.0;
    
60.     A(2,0)=0.0;A(2,1)=z0;A(2,2)=y0;A(2,3)=0.0;A(2,4)=-1+2*b0;A(2,5)=0.0;
    
61.     A(3,0)=z0;A(3,1)=0.0;A(3,2)=x0;A(3,3)=0.0;A(3,4)=0.0;A(3,5)=-1+2*c0;
    
62.     A(4,0)=2*x0;A(4,1)=-2*y0;A(4,2)=0.0;A(4,3)=2.0;A(4,4)=-2*b0;A(4,5)=-2+2*c0;
    
63.     A(5,0)=2*x0;A(5,1)=0.0;A(5,2)=-2*z0;A(5,3)=2*a0;A(5,4)=2-2*b0;A(5,5)=-2.0;
    
64. }
    
65. 
    
66. #define MAXITER 20
    
67. #define MINERR  0.00001
    
68. #define DELT    (4*MINERR)
    
69. void pushs(double x0, double y0,double z0, double a0,double b0,double c0)
    
70. {
    
71.     int i;
    
72.     for(i=0;i<solc;i++){
    
73.        if(fabs(x0-solbuf[i][0])<DELT&&fabs(y0-solbuf[i][1])<DELT&&
    
74.           fabs(z0-solbuf[i][2])<DELT&&fabs(a0-solbuf[i][3])<DELT&&
    
75.           fabs(b0-solbuf[i][4])<DELT&&fabs(c0-solbuf[i][5])<DELT){
    
76.             break;//find a matched solution
    
77.        }
    
78.     }
    
79.     if(i>=solc){
    
80.        if(solc>=MAXSOL){
    
81.           fprintf(stderr, "buffer overflow\n");
    
82.           return;
    
83.        }
    
84.        solbuf[solc][0]=x0;solbuf[solc][1]=y0;solbuf[solc][2]=z0;
    
85.        solbuf[solc][3]=a0;solbuf[solc][4]=b0;solbuf[solc][5]=c0;
    
86.        printf("x=%f,y=%f,z=%f,a=%f,b=%f,c=%f\n",
    
87.             x0,y0,z0,a0,b0,c0);
    
88.        solc++;
    
89.     }
    
90. }
    
91. 
    
92. int main(int argc, char *argv[])
    
93. {
    
94.     int N1=atoi(argv[1]);
    
95.     int N2=atoi(argv[2]);
    
96.     int count;
    
97.     double a1,b1,c1;
    
98.     double a[N*N],b[N];
    
99.     for(a1=0.1;a1<=0.9;a1+=0.1)for(b1=0.1;b1<=0.9;b1+=0.1)for(c1=0.1;c1<=0.9;c1+=0.1){
    
100.         double a0=a1,b0=b1,c0=c1,x0,y0,z0;
     
101.         double r=sqrt(a0*(1-a0)*b0*(1-b0)*c0*(1-c0));
     
102.         z0=r/(a0*(1-a0));x0=r/(b0*(1-b0));y0=r/(c0*(1-c0));
     
103.         if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1){
     
104.           for(count=0;count<MAXITER;count++){
     
105.             init_a(a,x0,y0,z0,a0,b0,c0,N1,N2);
     
106.             init_b(b,x0,y0,z0,a0,b0,c0,N1,N2);
     
107.             solve(a,b);
     
108.             double maxerr=-1.0;
     
109.             int u;
     
110.             for(u=0;u<N;u++){
     
111.                if(fabs(b[u])>maxerr)maxerr=fabs(b[u]);
     
112.             }
     
113.             x0-=b[0];y0-=b[1];z0-=b[2];a0-=b[3];b0-=b[4];c0-=b[5];
     
114.             if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1&&
     
115.                 0<a0&&a0<1&&0<b0&&b0<1&&0<c0&&c0<1){
     
116.                  if(maxerr<MINERR){
     
117.                      
     
118.                      pushs(x0,y0,z0,a0,b0,c0);
     
119.                      break;
     
120.                  }
     
121.             }else{
     
122.                 break;
     
123.             }
     
124.           }
     
125.         }
     
126.     }
     
127. }
     

复制代码