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

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？欢迎注册

×

本帖最后由 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[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 没啃完

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

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

防AB问题注解

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

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

各位仲秋快乐

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[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的次数

{{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 的次数。

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\}$。

修改一下前面代码

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值不能完全消元，再补充几个数值。

{{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}}}$

最终结果

(\!\(
\*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$

账号		自动登录	找回密码
密码			欢迎注册

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

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

点评