楼上矩阵最优解要求各行元素平方和相等，各列元素绝对值和为1，并且任意不同两行之间正交。可以看出，对于任意这样的矩阵，交换任意两行两列或改变任意行列所有数据的符号，还是满足条件的。利用这些性质，至少对于3*4阶矩阵，我们可以穷举所有情况，其中非零数不超过9个，反复利用上面性质变换，最后可以得出总共只会余下四种:
下面所有字母代表区间中实数:
\(\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$不是最优解
根据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$列，分成$$组，每组$m$列取成一个$m\times m$单位阵，余下$n-m$列全部为0，得到一个非最优解可以取到棱长为$\sqrt()$的m维超立方体。所以最优解棱长必然大于$\sqrt()$而不超过$\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}\),其中各行之间无法正交，淘汰
三种模式
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$是整数同样不行
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$
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)没有找到合法解
reduce_by_z(f1,f2)={
local(f3,f4,d1,d2,r,h);f3=f1;f4=f2;
d1=poldegree(f3,z);d2=poldegree(f4,z);h=polcoeff(f4,d2,z);
while(d1>d2,
f3=f3-polcoeff(f3,d1,z)/h *z^(d1-d2)*f4;
d1=poldegree(f3,z)
);
if(d1 == d2,
f3 = f3-polcoeff(f3,d1,z)/h*f4;
);
f3
}
reduce_by_b(f1,f2)={
local(f3,f4,d1,d2,r,h,df);f3=f1;f4=f2;
d1=poldegree(f3,b);d2=poldegree(f4,b);h=polcoeff(f4,d2,b);
while(d1>d2,
df=polcoeff(f3,d1,b)/h;
f3=denominator(df)*f3-numerator(df) *b^(d1-d2)*f4;
d1=poldegree(f3,b)
);
if(d1 == d2,
df=polcoeff(f3,d1,b)/h;
f3 = denominator(df)*f3-numerator(df)*f4;
);
f3
}
esol(N1,N2)={
local(f1,f2,f3,f4,f5,f6,r1,r2,r3,r4,r5, tmp,equ,rts);
local(x0,y0,z0,a0,b0,c0);
f1 = x+y+z-1;
f2 = x*y-a*(1-a);
f3 = y*z-b*(1-b);
f4 = z*x-c*(1-c);
f5 = x^2+a^2+(1-c)^2-y^2-(1-a)^2-b^2-N1;
f6 = x^2+a^2+(1-c)^2-z^2-(1-b)^2-c^2-N2;
if(poldegree(f5,a) !=1, print ("Error in r1"),
r1 = - polcoeff(f5,0,a)/polcoeff(f5,1,a);
f1 = subst(f1, a, r1);
f2 = subst(f2, a, r1);
f3 = subst(f3, a, r1);
f4 = subst(f4, a, r1);
f6 = subst(f6, a, r1);
if(poldegree(f1, x) != 1, print("Error in r2");print(f1),
r2 = - polcoeff(f1, 0, x)/polcoeff(f1,1,x);
f2 = subst(f2, x, r2);
f3 = subst(f3, x, r2);
f4 = subst(f4, x, r2);
f6 = subst(f6, x, r2);
if(poldegree(f3, y)!=1, print("Error in r3");print(f3),
r3 = - polcoeff(f3, 0, y)/polcoeff(f3, 1, y);
f2 = subst(f2, y, r3)*4*z;
f4 = subst(f4, y, r3);
f6 = subst(f6, y, r3)*4*z^2;
f2 = reduce_by_z(f2,f4);
f6 = reduce_by_z(f6,f4);
if(poldegree(f2,z)!=1, print("Error in r4"),
r4 = - polcoeff(f2, 0, z)/polcoeff(f2, 1, z);
f4 = numerator(subst(f4, z, r4));
f6 = numerator(subst(f6, z, r4))-f4;
while(poldegree(f6, b)>0,
tmp = reduce_by_b(f4,f6);
f4=f6;f6=tmp
);
f6=component(f6,1);
if(poldegree(f4, b)!=1, print("Error in r5"),
r5 = -polcoeff(f4,0,b)/polcoeff(f4,1,b);
equ = factor(f6);
for(s=1, length(equ[,1]),
rts=polroots(equ);
for(t=1,length(rts),
if(0< real(rts) && real(rts)<1 && abs(imag(rts))<0.00001,
c0=real(rts);
b0=subst(r5,c, c0);
z0=subst(r4,c, c0); z0=subst(z0, b, b0);
y0=subst(r3,c,c0);y0=subst(y0,b,b0);y0=subst(y0,z,z0);
x0=subst(r2,c,c0);x0=subst(x0,b,b0);x0=subst(x0,z,z0);x0=subst(x0,y,y0);
a0=subst(r1,c,c0);a0=subst(a0,b,b0);a0=subst(a0,z,z0);
a0=subst(a0,y,y0);
a0=subst(a0,x,x0);
if(0<b0&&b0<1&&0<z0&&z0<1&&0<y0&&y0<1&&0<x0&&x0<1&&0<a0&&a0<1,
print(equ); 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)
)
)
)
)
)
)
)
)
)
}
为了确认上面的结果，我又对上面部分方程采用密集采样点加牛顿迭代法进行数值解验算，结果一致，这说明24#得出的结果对于m=3已经是最优解了
比如下面代码数值验算了25#的方程:
#include <stdio.h>
#include <math.h>
#define MAXSOL 50
#define N 6
double solbuf;
int solc;
#define A(i,j) a[(i)*N+(j)]
void solve(double *a, double *b)
{
int i,j,k;
for(i=0;i<N;i++){//in i column
double maxv = -1.0;
int index=-1;
for(j=i;j<N;j++){
if(fabs(A(i,j))>maxv){
maxv=fabs(A(i,j));
index = j;
}
}
if(index != i){//swap row index and i
double tmp;
for(j=i;j<N;j++){
tmp = A(i,j);
A(i,j)=A(index,j);
A(index,j)=tmp;
}
tmp = b;
b= b;
b=tmp;
}
//divid by A(i,i) in row i
for(j=i+1;j<N;j++)A(i,j)/=A(i,i);
b/=A(i,i);A(i,i)=1.0;
for(j=0;j<N;j++){
if(j==i)continue;
for(k=i+1;k<N;k++){
A(j,k)-=A(i,k)*A(j,i);
}
b-=b*A(j,i);
A(j,i)=0.0;
}
}
}
void init_b(double *b, double x0, double y0, double z0, double a0, double b0, double c0, int N1, int N2)
{
b=x0+y0+z0-1.0;
b=x0*y0-a0*(1-a0);
b=y0*z0-b0*(1-b0);
b=z0*x0-c0*(1-c0);
b=x0*x0+a0*a0+(1-c0)*(1-c0)-y0*y0-(1-a0)*(1-a0)-b0*b0-N1;
b=x0*x0+a0*a0+(1-c0)*(1-c0)-z0*z0-(1-b0)*(1-b0)-c0*c0-N2;
}
void init_a(double *a, double x0, double y0, double z0, double a0, double b0, double c0,int N1, int N2)
{
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;
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;
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;
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;
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;
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;
}
#define MAXITER 20
#define MINERR0.00001
#define DELT (4*MINERR)
void pushs(double x0, double y0,double z0, double a0,double b0,double c0)
{
int i;
for(i=0;i<solc;i++){
if(fabs(x0-solbuf)<DELT&&fabs(y0-solbuf)<DELT&&
fabs(z0-solbuf)<DELT&&fabs(a0-solbuf)<DELT&&
fabs(b0-solbuf)<DELT&&fabs(c0-solbuf)<DELT){
break;//find a matched solution
}
}
if(i>=solc){
if(solc>=MAXSOL){
fprintf(stderr, "buffer overflow\n");
return;
}
solbuf=x0;solbuf=y0;solbuf=z0;
solbuf=a0;solbuf=b0;solbuf=c0;
printf("x=%f,y=%f,z=%f,a=%f,b=%f,c=%f\n",
x0,y0,z0,a0,b0,c0);
solc++;
}
}
int main(int argc, char *argv[])
{
int N1=atoi(argv);
int N2=atoi(argv);
int count;
double a1,b1,c1;
double a,b;
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){
double a0=a1,b0=b1,c0=c1,x0,y0,z0;
double r=sqrt(a0*(1-a0)*b0*(1-b0)*c0*(1-c0));
z0=r/(a0*(1-a0));x0=r/(b0*(1-b0));y0=r/(c0*(1-c0));
if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1){
for(count=0;count<MAXITER;count++){
init_a(a,x0,y0,z0,a0,b0,c0,N1,N2);
init_b(b,x0,y0,z0,a0,b0,c0,N1,N2);
solve(a,b);
double maxerr=-1.0;
int u;
for(u=0;u<N;u++){
if(fabs(b)>maxerr)maxerr=fabs(b);
}
x0-=b;y0-=b;z0-=b;a0-=b;b0-=b;c0-=b;
if(0<x0&&x0<1&&0<y0&&y0<1&&0<z0&&z0<1&&
0<a0&&a0<1&&0<b0&&b0<1&&0<c0&&c0<1){
if(maxerr<MINERR){
pushs(x0,y0,z0,a0,b0,c0);
break;
}
}else{
break;
}
}
}
}
}