整数关系探测算法，PSLQ 算法的缺陷？有什么新改进？

_only

整数关系问题最早可被追溯到欧几里得时代……爬拉爬拉……不浪费版面了

其中PSLQ算法是“20世纪十大算法之一”


https://mathworld.wolfram.com/IntegerRelation.html
https://mathworld.wolfram.com/PSLQAlgorithm.html

试验了一下Mathematica代码（来自 https://library.wolfram.com/infocenter/MathSource/4263/）
原实验\(\sqrt{5-\sqrt{2}}\)还好，改动为\(\sqrt{5-\sqrt{2}}+\sqrt{13}-\sqrt[3]{7}\)它的最小多项式是\(x^{24}-216 x^{22}+56 x^{21}+19812 x^{20}-6048 x^{19}-999956 x^{18}+295344 x^{17}+30220980 x^{16}-10738168 x^{15}-553263168 x^{14}+314972112 x^{13}+6032167718 x^{12}-5100616752 x^{11}-35472467208 x^{10}+24910873512 x^9+80640872964 x^8+161999430432 x^7-100710240772 x^6-1379769216672 x^5+826552295856 x^4+1397108782328 x^3+9646683219912 x^2+10193872879344 x+1668766870177\)

我实验一下并把方程次数提高为25、26。无法算出结果。

我搜索的结果，对于PSLQ的改进多是针对其实数/复数方面的。稳定性、效率上的改进没找到。

http://www.cas.cn/syky/201807/t20180703_4656888.shtml  没啃完

opukidefugowu

防AB问题注解
之所以接触这个，是因为机器证明纯符号计算量越来越大还是我的电脑老了？并且发现单以代数办法消元会引进虚数解、负数解等对实际几何问题无意义的增根。简单如\(\sqrt{x^2+y^2}=r\)用多项式表示就是\(x^2+y^2=r^2\)。就多一个负根。基本上多一组增根计算量就会翻一翻。
如果先以数值计算出正数解，再用倒推回单未知数方程是否可行？  后续问题接着是检验办法。

eqajoefiwae

各位仲秋愉快

opukidefugowu

zeroieme 发表于 2020-9-30 23:25
防AB问题注解
之所以接触这个，是因为机器证明纯符号计算量越来越大还是我的电脑老了？并且发现单以代 ...

既然声明了 防AB问题. 那我直接从A开始吧.哈哈哈.

如果先以数值计算出正数解，再用倒推回单未知数方程是否可行？  后续问题接着是检验办法。

不太现实, 从数值值到代数表达,这个有如大海捞针,非封闭,可逆的操作, 而且也不具有唯一性.

AbertVab

小试验
求 \(\triangle ABC\)，边长分别是a、b、c的外接圆半径r。
设点A在原点\({0,0}\)；点B在原点x轴正方向\({c,0}\)；点C坐标\({xC,yC}\)；外接圆心坐标\({xR,yR}\)
有以下方程组\(\left\{\text{xC}^2+\text{yC}^2=b^2,(\text{xC}-c)^2+\text{yC}^2=a^2,\text{xR}^2+\text{yR}^2=(\text{xR}-c)^2+\text{yR}^2=(\text{xR}-\text{xC})^2+(\text{yR}-\text{yC})^2=r^2,\text{yC}>0,r>0\right\}\)。

以下PSLQ代码来自 https://library.wolfram.com/infocenter/MathSource/4263  

1. PSLQ[inx_List, prec_] :=  Block[   {    x,    n = Length[inx],    \[Gamma] = 2/Sqrt[3],    A, B, H, D, Dinv, t, i, j, k, l, iter,    \[Alpha], \[Beta], \[Lambda], \[Delta], r, R    },   (*Initialize*)   x = N[inx /Sqrt[inx . inx], prec];   s = Sqrt[MapIndexed[Plus @@ Drop[x^2, First[#2] - 1] &, x]];   A = B = IdentityMatrix[n];   H = Table[Which[      i > j, (-x[[i]]*x[[j]])/(s[[j]]*s[[j + 1]]),      i == j, s[[i + 1]]/s[[i]],      i < j, 0      ], {i, 1, n}, {j, 1, n - 1}];   (* Reduce H *)   t = HermiteReduce[H];   D = First[t];   Dinv = Inverse[D];   (*Update*)   H = Last[t]; x = x.Dinv; A = D.A; B = B.Dinv;   For[iter = 0, iter < $IterationLimit, ++iter,    (* Step One *)    r = MaxIndex[      MapIndexed[\[Gamma]^First[#2] Abs[#1] &, Tr[H, List]]];    If[r < n - 1, \[Alpha] = H[[r, r]]; \[Beta] =       H[[r + 1, r]]; \[Lambda] = H[[r + 1, r + 1]]; \[Delta] =       Sqrt[\[Beta]^2 + \[Lambda]^2]];    R = IdentityMatrix[n]; t = R[[r]]; R[[r]] = R[[r + 1]];     R[[r + 1]] = t;    x = x.R; H = R.H; A = R.A; B = B.R;    (* Step Two *)    If[r < n - 1,     H = H.Table[        Which[         i == r && j == r, \[Beta]/\[Delta],         i == r && j == r + 1, -\[Lambda]/\[Delta],         i == r + 1 && j == r, \[Lambda]/\[Delta],         i == r + 1 && j == r + 1, \[Beta]/\[Delta],         i == j && j != r || i == j && j != r + 1, 1,         True, 0],        {i, 1, n - 1}, {j, 1, n - 1}]     ];    (* Step Three *)    t = HermiteReduce[H];    D = First[t];    Dinv = Inverse[D];    (*Update*)    H = Last[t]; x = x.Dinv; A = D.A; B = B.Dinv;    (* Step Four *)    If[Min[Abs[Union[x, Tr[H, List]]]] <= 10^(-prec + 5), Break[]]    ];(*Main Iteraton*)   Return[Transpose[B][[MaxIndex[-Abs[x]]]]]   ];Recognize2[n_, po_, v_] :=   PSLQ[N[Table[n^i, {i, 0, po}], 100], 100].Table[v^i, {i, 0, po}];MaxIndex[x_List] :=   Block[{i = 1, j}, Do[If[x[[j]] > x[[i]], i = j], {j, 2, Length[x]}];    i];HermiteReduce[H_] := Block[   {n = Length[H], i, j, k, q, H2 = H, D = IdentityMatrix[Length[H]]},   Do[    q = Round[H[[i, j]]/H[[j, j]]];    Do[H2[[i, k]] -= q H2[[j, k]], {k, 1, j}];    Do[D[[i, k]] -= q D[[j, k]], {k, 1, n}];    , {i, 2, n}, {j, i - 1, 1, -1}];   {D, H2}   ];

复制代码

开始探测表达式中r的次数

1. {{5,7,9},{11,13,15},{21,22,23},{33,34,35}}//Function[{a,b,c},{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}//NSolve[{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}/.{a->5,b->11,c->13},{xC,yC,xR,yR,r},WorkingPrecision->100]&//First//r/.#&//Function[{rApproximation},Table[Recognize2[rApproximation,i,r]//Factor//If[#[[0]]===Times,Select[#,((r/.NSolve[{#==0,r>0},{r},WorkingPrecision->50])-rApproximation=={0})&](*过滤多余因子*),#]&//Sign[#/.r->10^4]#(*调整方程系数符号*)&,{i,2,8}]//DeleteDuplicates//{{a,b,c},rApproximation,#}&]]@@#&/@#&

复制代码

结果\(
\left(
\begin{array}{ccc}
\{5,7,9\} & 4.522670168666454339702181004550936387173302562167300359625475653795206330830514336421150089720396342 & \left\{11 r^2-225\right\} \\
\{11,13,15\} & 7.701540462154053919411117443631364516316459094012752173067915657884708281454925563612634347927551431 & \left\{51 r^2-3025\right\} \\
\{21,22,23\} & 12.72816758217772681129554651634169209822120943578399772810145703249019261552293384091059863622668369 & \left\{160 r^2-25921\right\} \\
\{33,34,35\} & 19.64694897857340766262404934920350387306020375110037498847055496649207802752221427537178920302680188 & \left\{384 r^2-148225\right\} \\
\end{array}
\right)\)

可以推断 外接圆半径r符合 \(p r^2-q=0\)形式。

接着推断 \(p r^2-q=0\)中a 的次数。

1. Range[400,600]//{#,#[[1]]//PrimePi//#-1&//Prime,#[[-1]]//PrimePi//#+1&//Prime}&//Function[{A,b,c},Table[{a,b,c},{a,A}]]@@#&//Function[{a,b,c},{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}//NSolve[{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}/.{a->5,b->11,c->13},{xC,yC,xR,yR,r},WorkingPrecision->100]&//First//r/.#&//Function[{rApproximation},Recognize2[rApproximation,2,r]//FactorTermsList//Last//Sign[#/.r->10^4]#(*调整方程系数符号*)&//{{a,b,c},#}&]]@@#&/@#&//(*修补约掉的纯整数因子--开始*)FixedPoint[Function[{X},Table[{If[i==1,1,Round[Coefficient[X[[i-1,-1]],r,0]/Coefficient[X[[i,-1]],r,0]]],1,If[i==Length[X],1,Round[Coefficient[X[[i+1,-1]],r,0]/Coefficient[X[[i,-1]],r,0]]]}//Max//{X[[i,1]],# X[[i,-1]]//Expand}&,{i,Length[X]}]],#]&//{#[[1]],FactorTermsList[#[[2]]]//Sign[#[[1]]]#&}&/@#&//Function[{X},Select[X,#[[2,1]]==1&]//({#[[1,1]],Coefficient[#[[2,-1]],r,0]}//(\!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(10\)]\(\*SubscriptBox[\(t\), \(i\)] \*SuperscriptBox[\(#[\([1]\)]\), \(i\)]\)\))==#[[2]]&)&/@#&//Solve//First//(\!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(10\)]\(\*SubscriptBox[\(t\), \(i\)] \*SuperscriptBox[\(n\), \(i\)]\)\))/.#&//Function[{tFunction},{#[[1]],Round[(tFunction/.n->#[[1,1]])/Coefficient[#[[2,-1]],r,0]]#[[2,-1]]//Expand}&/@X]](*修补约掉的纯整数因子--完成*)//Coefficient[#[[2]],r,2]&/@#&//NestWhileList[Differences,#,Complement[#,{0}]!={}&]&//Length//#-2&

复制代码

结果为4


最后，根据纲量齐次性与三边对称性,设外接圆半径r方程为\(r^2 \left(s_1^4 t_{4,1}+s_2 s_1^2 t_{4,2}+s_3 s_1 t_{4,3}+s_2^2 t_{4,4}\right)+s_1^6 t_{6,1}+s_2 s_1^4 t_{6,2}+s_3 s_1^3 t_{6,3}+s_2^2 s_1^2 t_{6,4}+s_2 s_3 s_1 t_{6,5}+s_2^3 t_{6,6}+s_3^2 t_{6,7}=0\)其中\(\left\{s_1=a+b+c,s_2=a b+a c+b c,s_3=a b c\right\}\)。

修改一下前面代码

1. Range[400,600]//{#,#[[1]]//PrimePi//#-1&//Prime,#[[-1]]//PrimePi//#+1&//Prime}&//Function[{A,b,c},Table[{a,b,c},{a,A}]]@@#&//Function[{a,b,c},{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}//NSolve[{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}/.{a->5,b->11,c->13},{xC,yC,xR,yR,r},WorkingPrecision->100]&//First//r/.#&//Function[{rApproximation},Recognize2[rApproximation,2,r]//FactorTermsList//Last//Sign[#/.r->10^4]#(*调整方程系数符号*)&//{{a,b,c},#}&]]@@#&/@#&//(*修补约掉的纯整数因子--开始*)FixedPoint[Function[{X},Table[{If[i==1,1,Round[Coefficient[X[[i-1,-1]],r,0]/Coefficient[X[[i,-1]],r,0]]],1,If[i==Length[X],1,Round[Coefficient[X[[i+1,-1]],r,0]/Coefficient[X[[i,-1]],r,0]]]}//Max//{X[[i,1]],# X[[i,-1]]//Expand}&,{i,Length[X]}]],#]&//{#[[1]],FactorTermsList[#[[2]]]//Sign[#[[1]]]#&}&/@#&//Function[{X},Select[X,#[[2,1]]==1&]//({#[[1,1]],Coefficient[#[[2,-1]],r,0]}//(\!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(10\)]\(\*SubscriptBox[\(t\), \(i\)] \*SuperscriptBox[\(#[[1]]\), \(i\)]\)\))==#[[2]]&)&/@#&//Solve//First//(\!\(\*UnderoverscriptBox[\(\[Sum]\), \(i = 0\), \(10\)]\(\*SubscriptBox[\(t\), \(i\)] \*SuperscriptBox[\(n\), \(i\)]\)\))/.#&//Function[{tFunction},{#[[1]],Round[(tFunction/.n->#[[1,1]])/Coefficient[#[[2,-1]],r,0]]#[[2,-1]]//Expand}&/@X]](*修补约掉的纯整数因子--完成*)//{(\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(t\), \(4, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 2\)]\)+Subscript[s, 1] Subscript[s, 3] Subscript[t, 4,3]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 4\)]\))==Coefficient[#[[2]],r,2],(\!\(\*SubsuperscriptBox[\(s\), \(1\), \(6\)]\ \*SubscriptBox[\(t\), \(6, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 2\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(3\)]\ \*SubscriptBox[\(s\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 3\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 4\)]\)+Subscript[s, 1] Subscript[s, 2] Subscript[s, 3] Subscript[t, 6,5]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 6\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(3\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 7\)]\))==Coefficient[#[[2]],r,0]}//.{Subscript[s, 1]->a+b+c,Subscript[s, 2]->a b+a c+b c,Subscript[s, 3]->a b c,a->#[[1,1]],b->#[[1,2]],c->#[[1,3]]}&/@#&//Flatten//Solve

复制代码

\(\left\{\left\{t_{4,1}\to -1,t_{4,2}\to 4,t_{4,3}\to -8,t_{4,4}\to 0,t_{6,1}\to 0,t_{6,2}\to 0,t_{6,3}\to \frac{53569616320232100 t_{6,7}}{415680381919678231}+\frac{53569616320232100}{415680381919678231},t_{6,4}\to -\frac{12778274696348069 t_{6,7}}{415680381919678231}-\frac{12778274696348069}{415680381919678231},t_{6,5}\to -\frac{338451191942 t_{6,7}}{556980108721}-\frac{338451191942}{556980108721},t_{6,6}\to \frac{53569616320232100 t_{6,7}}{415680381919678231}+\frac{53569616320232100}{415680381919678231}\right\}\right\}\)
仅变动a值不能完全消元，再补充几个数值。

1. {{5,7,9},{11,13,15},{21,22,23},{33,34,35}}//Function[{a,b,c},{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}//NSolve[{xC^2+yC^2==b^2,(xC-c)^2+yC^2==a^2,xR^2+yR^2==(xR-c)^2+yR^2==(xR-xC)^2+(yR-yC)^2==r^2,yC>0,r>0}/.{a->5,b->11,c->13},{xC,yC,xR,yR,r},WorkingPrecision->100]&//First//r/.#&//Function[{rApproximation},Recognize2[rApproximation,2,r]]//Coefficient[#,r,0]/Coefficient[#,r,2]==(\!\(\*SubsuperscriptBox[\(s\), \(1\), \(6\)]\ \*SubscriptBox[\(t\), \(6, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 2\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(3\)]\ \*SubscriptBox[\(s\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 3\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 4\)]\)+Subscript[s, 1] Subscript[s, 2] Subscript[s, 3] Subscript[t, 6,5]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 6\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(3\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 7\)]\))/(\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(t\), \(4, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 2\)]\)+Subscript[s, 1] Subscript[s, 3] Subscript[t, 4,3]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 4\)]\))//.Flatten[{{Subscript[s, 1]->a+b+c,Subscript[s, 2]->a b+a c+b c,Subscript[s, 3]->a b c},{Subscript[t, 4,1]->-1,Subscript[t, 4,2]->4,Subscript[t, 4,3]->-8,Subscript[t, 4,4]->0,Subscript[t, 6,1]->0,Subscript[t, 6,2]->0,Subscript[t, 6,3]->53569616320232100/415680381919678231+(53569616320232100 Subscript[t, 6,7])/415680381919678231,Subscript[t, 6,4]->-(12778274696348069/415680381919678231)-(12778274696348069 Subscript[t, 6,7])/415680381919678231,Subscript[t, 6,5]->-(338451191942/556980108721)-(338451191942 Subscript[t, 6,7])/556980108721,Subscript[t, 6,6]->53569616320232100/415680381919678231+(53569616320232100 Subscript[t, 6,7])/415680381919678231}}]&//Solve]@@#&/@#&

复制代码

\({{{Subscript[t, 6,7]->-1}},{{Subscript[t, 6,7]->-1}},{{Subscript[t, 6,7]->-1}},{{Subscript[t, 6,7]->-1}}}\)

最终结果

1. (\!\(\*SubsuperscriptBox[\(s\), \(1\), \(6\)]\ \*SubscriptBox[\(t\), \(6, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 2\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(3\)]\ \*SubscriptBox[\(s\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 3\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 4\)]\)+Subscript[s, 1] Subscript[s, 2] Subscript[s, 3] Subscript[t, 6,5]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(3\)]\ \*SubscriptBox[\(t\), \(6, 6\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(3\), \(2\)]\ \*SubscriptBox[\(t\), \(6, 7\)]\))-(\!\(\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\ \*SubscriptBox[\(t\), \(4, 1\)]\)+\!\(\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\ \*SubscriptBox[\(s\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 2\)]\)+Subscript[s, 1] Subscript[s, 3] Subscript[t, 4,3]+\!\(\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\ \*SubscriptBox[\(t\), \(4, 4\)]\))r^2//.{Subscript[t, 4,1]->-1,Subscript[t, 4,2]->4,Subscript[t, 4,3]->-8,Subscript[t, 4,4]->0,Subscript[t, 6,1]->0,Subscript[t, 6,2]->0,Subscript[t, 6,3]->53569616320232100/415680381919678231+(53569616320232100 Subscript[t, 6,7])/415680381919678231,Subscript[t, 6,4]->-(12778274696348069/415680381919678231)-(12778274696348069 Subscript[t, 6,7])/415680381919678231,Subscript[t, 6,5]->-(338451191942/556980108721)-(338451191942 Subscript[t, 6,7])/556980108721,Subscript[t, 6,6]->53569616320232100/415680381919678231+(53569616320232100 Subscript[t, 6,7])/415680381919678231,Subscript[t, 6,7]->-1,Subscript[s, 1]->a+b+c,Subscript[s, 2]->a b+a c+b c,Subscript[s, 3]->a b c}//Collect[#,r,Factor]&

复制代码

\(r^2 (a-b-c) (a+b-c) (a-b+c) (a+b+c)-a^2 b^2 c^2\)