mathematica精度问题

wsc810

以下是计算分解  $2^512+1$  报错


IntegerPart::meprec: 当计算 IntegerPart[(98326482234939147319765934357793685660339022949204652705961168633314897265284+Sqrt[13407807929942597099574024998205846127479365820592393377723561443721764030073546976801874298166903427690031858186486050853753882811946569946433649006084097])/21411857147225534274333691936648159351360900761484521674541875264122802693051] 时，达到内部精确度极限值 MaxExtraPrecision = 50

用向下取整函数  Floor[ ]  还是一样

$P_0=\sqrt{d} (mod Q_0^2)$

分解整数所用的原理是用到如下pell 方程公式，迭代过程中$Q_i$为完全平方数终止迭代

$(Q_0^2p_i-P_0q_i)^2-dq_i^2=(-1)^{i+1}Q_0^2Q_i$

1. Clear[t, a, d, P, Q, i, p, q]; d = 2^512 + 1; pell = -1;
   
2. P[0] = 107450; Q[0] = 523^2;
   
3. t[0] = (P[0] + Sqrt[d])/Q[0]
   
4. a[0] = IntegerPart[t[0]];
   
5. i = 0; While[(t[i] != 1/(t[0] - a[0]) && P[i] != pell) || i == 1,
   
6. P[i + 1] = Q[i] a[i] - P[i];
   
7. Q[i + 1] = (d - P[i + 1]^2)/Q[i];
   
8. t[i + 1] = (P[i + 1] + Sqrt[d])/Q[i + 1];
   
9. a[i + 1] = IntegerPart [t[i + 1]];
   
10. 
    
11. p[-1] = 0; p[0] = 1; p[i + 1] = a[i] p[i] + p[i - 1];
    
12. q[-1] = 1; q[0] = 0; q[i + 1] = a[i] q[i] + q[i - 1];
    
13. 
    
14. If[Q[i + 1] == Q[i] || P[i + 1] == P[i] ||
    
15.    FractionalPart[Sqrt[Q[i + 1]]] == 0, Break[]];
    
16. Print[{i, Q[i], P[i], a[i]}] i++];
    
17. {i, Q[i], P[i], a[i]}
    
18. {i + 1, Q[i + 1], P[i + 1], a[i + 1]}
    
19. GCD[Q[0] p[i + 1] - P[0] q[i + 1] + Sqrt[Q[i + 1] Q[0]], d]
    

复制代码

mathematica

看到我的签名了吗？
没注释的程序，没人愿意看的，放弃吧

mathematica

要懂得换位思考，要想别人帮助你，就加上你的程序注释，
要学会站在别人的观点上看问题，
懂了吗？

dlpg070

在代码前部增加
{$MaxExtraPrecision = 10000}
试一试

wsc810

dlpg070 发表于 2020-6-27 16:57
在代码前部增加
{$MaxExtraPrecision = 10000}
试一试

1. Clear[t, a, d, P, Q, i, p, q]; $MaxExtraPrecision = 1000;
   
2. d = 3*(2^512 + 1); pell = -1; P[0] = 47654; Q[0] = 521^2;
   
3. t[0] = (P[0] + Sqrt[d])/Q[0]
   
4. a[0] = IntegerPart[t[0]];
   
5. i = 0; While[ i <= 1000 && (t[i] != 1/(t[0] - a[0]) && P[i] != pell) || i == 1,

复制代码

剩下的代码和之前一样，这里用 i 控制最大迭代次数

谁能把上述代码改为动态形式，即不保留迭代过程中间的值

dlpg070

wsc810 发表于 2020-6-27 21:37
剩下的代码和之前一样，这里用 i 控制最大迭代次数

谁能把上述代码改为动态形式，即不保留迭代 ...

查得因子分解结果
n = 2^512 + 1
n1 = 2424833
n2 = 7455602825647884208337395736200454918783366342657
n3 = 74164006262753080152478714190193747405994078109751902390582131614\
4415759504705008092818711693940737
n=n1*n2*n3
n1:7位, n2:49位, n3:99位 都是素数
已经验证

zeroieme

wsc810 发表于 2020-6-27 21:37
剩下的代码和之前一样，这里用 i 控制最大迭代次数

谁能把上述代码改为动态形式，即不保留迭代 ...

把动态变量放在一起成为一个List，编写一个迭代函数实现\(List_i->List_{i+1}\)，用NestWhile控制迭代终点。

wsc810

利用Round函数，将 Q，a 扩展为负数。计算$p_n,q_n$利用mod函数，减小对内存的占用

1. Clear[t, a, p, q];(*d=3*(2^512+1);pell=-1;P[0]=223787;Q[0]=521^2;*)
   
2. 
   
3. d = 997331; pell = -1; P[0] = 0; Q[0] = 1;
   
4. 
   
5. t[0] = (P[0] + Sqrt[d])/Q[0]
   
6. a[0] = Round[t[0]];
   
7. i = 0; While[
   
8. i < 120000 && (t[i] != 1/(t[0] - a[0]) && P[i] != pell) || i == 1,
   
9. P[i + 1] = Q[i] a[i] - P[i];
   
10. Q[i + 1] = (d - P[i + 1]^2)/Q[i];
    
11. t[i + 1] = (P[i + 1] + Sqrt[d])/Q[i + 1];
    
12. a[i + 1] = Round[t[i + 1]];
    
13. 
    
14. p[-1] = 0; p[0] = 1; p[i + 1] = Mod[a[i] p[i] + p[i - 1], d];
    
15. q[-1] = 1; q[0] = 0; q[i + 1] = Mod[a[i] q[i] + q[i - 1], d];
    
16. 
    
17. If[Q[i + 1] == Q[i] || P[i + 1] == P[i] ||
    
18.    FractionalPart[Sqrt[Q[i + 1]]] == 0, Break[]];
    
19. 
    
20. Print[{i, Q[i], P[i], a[i]}] i++];
    
21. {i, Q[i], P[i], a[i]}
    
22. {i + 1, Q[i + 1], P[i + 1], a[i + 1]}
    
23. GCD[Q[0] p[i + 1] - P[0] q[i + 1] + (I)^(i + 1)*Sqrt[Q[i + 1] Q[0]],
    
24.   d] // AbsoluteTiming
    

复制代码

mathematica

dlpg070 发表于 2020-6-28 10:07
查得因子分解结果
n = 2^512 + 1
n1 = 2424833

为什么要分解第9个费马数呢？直接网上查询不可以吗？