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

zeroieme 发表于 2020-9-30 23:13:17

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

本帖最后由 zeroieme 于 2020-9-30 23:26 编辑

整数关系问题最早可被追溯到欧几里得时代……爬拉爬拉……不浪费版面了
其中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{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没啃完

zeroieme 发表于 2020-9-30 23:25:24

本帖最后由 zeroieme 于 2020-9-30 23:31 编辑

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

zeroieme 发表于 2020-10-1 00:18:58

各位仲秋快乐:loveliness:

wayne 发表于 2020-10-6 21:02:06

zeroieme 发表于 2020-9-30 23:25
防AB问题注解
之所以接触这个，是因为机器证明纯符号计算量越来越大还是我的电脑老了？并且发现单以代 ...既然声明了 防AB问题. 那我直接从A开始吧.哈哈哈.
如果先以数值计算出正数解，再用倒推回单未知数方程是否可行？后续问题接着是检验办法。
不太现实, 从数值值到代数表达,这个有如大海捞针,非封闭,可逆的操作, 而且也不具有唯一性.

zeroieme 发表于 2020-10-8 22:26:21

wayne 发表于 2020-10-6 21:02
既然声明了 防AB问题. 那我直接从A开始吧.哈哈哈.
不太现实, 从数值值到代数表达,这个有如大海捞针,非封 ...

小试验
求 \(\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
PSLQ :=
Block[
   {
    x,
    n = Length,
    \ = 2/Sqrt,
    A, B, H, D, Dinv, t, i, j, k, l, iter,
    \, \, \, \, r, R
    },
   (*Initialize*)
   x = N, prec];
   s = Sqrt - 1] &, x]];
   A = B = IdentityMatrix;
   H = Table[Which[
      i > j, (-x[]*x[])/(s[]*s[]),
      i == j, s[]/s[],
      i < j, 0
      ], {i, 1, n}, {j, 1, n - 1}];
   (* Reduce H *)
   t = HermiteReduce;
   D = First;
   Dinv = Inverse;
   (*Update*)
   H = Last; x = x.Dinv; A = D.A; B = B.Dinv;
   For[iter = 0, iter < $IterationLimit, ++iter,
    (* Step One *)
    r = MaxIndex[
      MapIndexed[\^First[#2] Abs[#1] &, Tr]];
    If = H[]; \ =
      H[]; \ = H[]; \ =
      Sqrt[\^2 + \^2]];
    R = IdentityMatrix; t = R[]; R[] = R[];
    R[] = 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, \/\,
         i == r && j == r + 1, -\/\,
         i == r + 1 && j == r, \/\,
         i == r + 1 && j == r + 1, \/\,
         i == j && j != r || i == j && j != r + 1, 1,
         True, 0],
      {i, 1, n - 1}, {j, 1, n - 1}]
   ];
    (* Step Three *)
    t = HermiteReduce;
    D = First;
    Dinv = Inverse;
    (*Update*)
    H = Last; x = x.Dinv; A = D.A; B = B.Dinv;
    (* Step Four *)
    If]]] <= 10^(-prec + 5), Break[]]
    ];(*Main Iteraton*)
   Return[]]]]
   ];
Recognize2 :=
PSLQ, 100], 100].Table;
MaxIndex :=
Block[{i = 1, j}, Do] > x[], i = j], {j, 2, Length}];
    i];
HermiteReduce := Block[
   {n = Length, i, j, k, q, H2 = H, D = IdentityMatrix]},
   Do[
    q = Round]/H[]];
    Do] -= q H2[], {k, 1, j}];
    Do] -= q D[], {k, 1, n}];
    , {i, 2, n}, {j, i - 1, 1, -1}];
   {D, H2}
   ];


开始探测表达式中r的次数
{{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//Factor//If[#[]===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 的次数。
Range//{#,#[]//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//FactorTermsList//Last//Sign[#/.r->10^4]#(*调整方程系数符号*)&//{{a,b,c},#}&]]@@#&/@#&//(*修补约掉的纯整数因子--开始*)FixedPoint],r,0]/Coefficient],r,0]]],1,If,1,Round],r,0]/Coefficient],r,0]]]}//Max//{X[],# X[]//Expand}&,{i,Length}]],#]&//{#[],FactorTermsList[#[]]//Sign[#[]]#&}&/@#&//Function[{X},Select]==1&]//({#[],Coefficient[#[],r,0]}//(\!\(
\*UnderoverscriptBox[\(\\), \(i = 0\), \(10\)]\(
\*SubscriptBox[\(t\), \(i\)]
\*SuperscriptBox[\(#[\(\)]\), \(i\)]\)\))==#[]&)&/@#&//Solve//First//(\!\(
\*UnderoverscriptBox[\(\\), \(i = 0\), \(10\)]\(
\*SubscriptBox[\(t\), \(i\)]
\*SuperscriptBox[\(n\), \(i\)]\)\))/.#&//Function[{tFunction},{#[],Round[(tFunction/.n->#[])/Coefficient[#[],r,0]]#[]//Expand}&/@X]](*修补约掉的纯整数因子--完成*)//Coefficient[#[],r,2]&/@#&//NestWhileList!={}&]&//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\}\)。

修改一下前面代码

Range//{#,#[]//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//FactorTermsList//Last//Sign[#/.r->10^4]#(*调整方程系数符号*)&//{{a,b,c},#}&]]@@#&/@#&//(*修补约掉的纯整数因子--开始*)FixedPoint],r,0]/Coefficient],r,0]]],1,If,1,Round],r,0]/Coefficient],r,0]]]}//Max//{X[],# X[]//Expand}&,{i,Length}]],#]&//{#[],FactorTermsList[#[]]//Sign[#[]]#&}&/@#&//Function[{X},Select]==1&]//({#[],Coefficient[#[],r,0]}//(\!\(
\*UnderoverscriptBox[\(\\), \(i = 0\), \(10\)]\(
\*SubscriptBox[\(t\), \(i\)]
\*SuperscriptBox[\(#[]\), \(i\)]\)\))==#[]&)&/@#&//Solve//First//(\!\(
\*UnderoverscriptBox[\(\\), \(i = 0\), \(10\)]\(
\*SubscriptBox[\(t\), \(i\)]
\*SuperscriptBox[\(n\), \(i\)]\)\))/.#&//Function[{tFunction},{#[],Round[(tFunction/.n->#[])/Coefficient[#[],r,0]]#[]//Expand}&/@X]](*修补约掉的纯整数因子--完成*)//{(\!\(
\*SubsuperscriptBox[\(s\), \(1\), \(4\)]\
\*SubscriptBox[\(t\), \(4, 1\)]\)+\!\(
\*SubsuperscriptBox[\(s\), \(1\), \(2\)]\
\*SubscriptBox[\(s\), \(2\)]\
\*SubscriptBox[\(t\), \(4, 2\)]\)+Subscript Subscript Subscript+\!\(
\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\
\*SubscriptBox[\(t\), \(4, 4\)]\))==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 Subscript Subscript Subscript+\!\(
\*SubsuperscriptBox[\(s\), \(2\), \(3\)]\
\*SubscriptBox[\(t\), \(6, 6\)]\)+\!\(
\*SubsuperscriptBox[\(s\), \(3\), \(2\)]\
\*SubscriptBox[\(t\), \(6, 7\)]\))==Coefficient[#[],r,0]}//.{Subscript->a+b+c,Subscript->a b+a c+b c,Subscript->a b c,a->#[],b->#[],c->#[]}&/@#&//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值不能完全消元，再补充几个数值。

{{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]//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 Subscript Subscript Subscript+\!\(
\*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 Subscript Subscript+\!\(
\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\
\*SubscriptBox[\(t\), \(4, 4\)]\))//.Flatten[{{Subscript->a+b+c,Subscript->a b+a c+b c,Subscript->a b c},{Subscript->-1,Subscript->4,Subscript->-8,Subscript->0,Subscript->0,Subscript->0,Subscript->53569616320232100/415680381919678231+(53569616320232100 Subscript)/415680381919678231,Subscript->-(12778274696348069/415680381919678231)-(12778274696348069 Subscript)/415680381919678231,Subscript->-(338451191942/556980108721)-(338451191942 Subscript)/556980108721,Subscript->53569616320232100/415680381919678231+(53569616320232100 Subscript)/415680381919678231}}]&//Solve]@@#&/@#&

\({{{Subscript->-1}},{{Subscript->-1}},{{Subscript->-1}},{{Subscript->-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 Subscript Subscript Subscript+\!\(
\*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 Subscript Subscript+\!\(
\*SubsuperscriptBox[\(s\), \(2\), \(2\)]\
\*SubscriptBox[\(t\), \(4, 4\)]\))r^2//.{Subscript->-1,Subscript->4,Subscript->-8,Subscript->0,Subscript->0,Subscript->0,Subscript->53569616320232100/415680381919678231+(53569616320232100 Subscript)/415680381919678231,Subscript->-(12778274696348069/415680381919678231)-(12778274696348069 Subscript)/415680381919678231,Subscript->-(338451191942/556980108721)-(338451191942 Subscript)/556980108721,Subscript->53569616320232100/415680381919678231+(53569616320232100 Subscript)/415680381919678231,Subscript->-1,Subscript->a+b+c,Subscript->a b+a c+b c,Subscript->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\)

查看完整版本: 整数关系探测算法，PSLQ 算法的缺陷？有什么新改进？