数学研发论坛

 找回密码
 欢迎注册
查看: 332|回复: 10

[原创] mathematica精度问题

[复制链接]
发表于 2020-6-27 12:20:07 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
本帖最后由 wsc810 于 2020-6-27 12:42 编辑

以下是计算分解  $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. p[-1] = 0; p[0] = 1; p[i + 1] = a[i] p[i] + p[i - 1];
  11. q[-1] = 1; q[0] = 0; q[i + 1] = a[i] q[i] + q[i - 1];

  12. If[Q[i + 1] == Q[i] || P[i + 1] == P[i] ||
  13.    FractionalPart[Sqrt[Q[i + 1]]] == 0, Break[]];
  14. Print[{i, Q[i], P[i], a[i]}] i++];
  15. {i, Q[i], P[i], a[i]}
  16. {i + 1, Q[i + 1], P[i + 1], a[i + 1]}
  17. GCD[Q[0] p[i + 1] - P[0] q[i + 1] + Sqrt[Q[i + 1] Q[0]], d]
复制代码


毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-6-27 13:27:26 | 显示全部楼层
看到我的签名了吗?
没注释的程序,没人愿意看的,放弃吧
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-6-27 13:28:26 | 显示全部楼层
要懂得换位思考,要想别人帮助你,就加上你的程序注释,
要学会站在别人的观点上看问题,
懂了吗?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-6-27 16:57:50 | 显示全部楼层
在代码前部增加
{$MaxExtraPrecision = 10000}
试一试

点评

感谢dlpg070的解答,程序已经修改好了  发表于 2020-6-27 21:35
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-6-27 21:37:13 | 显示全部楼层
本帖最后由 wsc810 于 2020-6-27 21:42 编辑
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 控制最大迭代次数

谁能把上述代码改为动态形式,即不保留迭代过程中间的值
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-6-28 10:07:16 | 显示全部楼层
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位 都是素数
已经验证
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-7-4 13:23:08 | 显示全部楼层
wsc810 发表于 2020-6-27 21:37
剩下的代码和之前一样,这里用 i 控制最大迭代次数

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

把动态变量放在一起成为一个List,编写一个迭代函数实现\(List_i->List_{i+1}\),用NestWhile控制迭代终点。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-7-26 07:33:16 | 显示全部楼层
本帖最后由 wsc810 于 2020-7-26 07:42 编辑

利用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. d = 997331; pell = -1; P[0] = 0; Q[0] = 1;

  3. t[0] = (P[0] + Sqrt[d])/Q[0]
  4. a[0] = Round[t[0]];
  5. i = 0; While[
  6. i < 120000 && (t[i] != 1/(t[0] - a[0]) && P[i] != pell) || i == 1,
  7. P[i + 1] = Q[i] a[i] - P[i];
  8. Q[i + 1] = (d - P[i + 1]^2)/Q[i];
  9. t[i + 1] = (P[i + 1] + Sqrt[d])/Q[i + 1];
  10. a[i + 1] = Round[t[i + 1]];

  11. p[-1] = 0; p[0] = 1; p[i + 1] = Mod[a[i] p[i] + p[i - 1], d];
  12. q[-1] = 1; q[0] = 0; q[i + 1] = Mod[a[i] q[i] + q[i - 1], d];

  13. If[Q[i + 1] == Q[i] || P[i + 1] == P[i] ||
  14.    FractionalPart[Sqrt[Q[i + 1]]] == 0, Break[]];

  15. Print[{i, Q[i], P[i], a[i]}] i++];
  16. {i, Q[i], P[i], a[i]}
  17. {i + 1, Q[i + 1], P[i + 1], a[i + 1]}
  18. GCD[Q[0] p[i + 1] - P[0] q[i + 1] + (I)^(i + 1)*Sqrt[Q[i + 1] Q[0]],
  19.   d] // AbsoluteTiming
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2020-7-27 08:26:49 | 显示全部楼层
dlpg070 发表于 2020-6-28 10:07
查得因子分解结果
n = 2^512 + 1
n1 = 2424833

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

点评

回复错了对象,我没有分解---,我只是上网查询 !!!  发表于 2020-7-27 10:24
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2020-10-23 19:37 , Processed in 0.180966 second(s), 21 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表