数学研发论坛

 找回密码
 欢迎注册
查看: 181|回复: 6

[求助] pell程序怎么修改

[复制链接]
发表于 2019-4-24 17:49:37 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 kte 于 2019-4-24 20:03 编辑

  1. Clear[t, a, P]; d = 137*89  ; pell = -1; P[0] = 0; Q[0] = 1;
  2. t[0] = (P[0] + Sqrt[d])/Q[0]
  3. a[0] = IntegerPart[t[0]];
  4. i = 0; While[(t[i] != 1/(t[0] - a[0]) && P[i] != pell) || i == 1,
  5. P[i + 1] = Q[i] a[i] - P[i];
  6. Q[i + 1] = (d - P[i + 1]^2)/Q[i];
  7. t[i + 1] = (P[i + 1] + Sqrt[d])/Q[i + 1];
  8. a[i + 1] = IntegerPart[t[i + 1]];
  9. If[Q[i + 1] == Q[i] || P[i + 1] == P[i], Break[]];

  10. Print[{i, Q[i], P[i], a[i]}] i++];
  11. {i, Q[i], P[i], a[i]}
  12. {i + 1, Q[i + 1], P[i + 1], a[i + 1]}
复制代码


上面的是一般程序,

现在要求程序在Q[ i ]为平方数时跳出循环,并在模运算下计算p[n]的值,p[-1]=0;p[-2]=1

  1. p[n]=Mod[a[i]p[i-1]+p[i-2],d]  
复制代码


程序该如何修改

这个算法还是比较有意义的,由于$p_n^2-dq_n^2=(-1)^{n+1}Q_{n+1}$,当$Q_{n+1}$为完全平方数时, 两边对$d$ 取模

$p_n^2=(-1)^{n+1}Q_{n+1}^2  (mod d)$

则可以知道d的一个分解,或者可以将 d写成两个平方数的和

$(frac{p_n}{Q_{n+1}})^2=-1=i^2 ( mod d)$


$GCD[frac{p_n}{Q_{n+1}}+I,d]=a+bI$

$a^2+b^2=d$


补充内容 (2019-4-26 11:14):
程序中应该是判断Q[i+1]为平方数
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-4-26 11:12:36 | 显示全部楼层
怎么没人帮我
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-4-26 11:39:30 | 显示全部楼层
没注释
没缩进,
这样的代码谁看呀,
当然,还没钱!
本来mathematica软件调试功能就很垃圾,
你还让别人改你程序,怎么可能?
还是自己慢慢研究吧,
反正我特别讨厌调试mathematica代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-4-26 11:58:44 来自手机 | 显示全部楼层
前面的判断Q[i+1]为平方数时,跳出循环我会写了,现在仅仅要求写出mod那一行代码,该怎么写

点评

Mod函数  发表于 2019-4-26 13:08
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2019-4-26 13:21:20 | 显示全部楼层
本帖最后由 kte 于 2019-4-26 13:50 编辑

因子分解程序,修改好了


  1. Clear[t, a, P, p]; d = 29*45113; pell = -1; P[0] = 0; Q[0] = 1;
  2. t[0] = (P[0] + Sqrt[d])/Q[0]
  3. a[0] = IntegerPart[t[0]];
  4. p[0] = 1; p[-1] = 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[i + 1] = Mod[a[i] p[i] + p[i - 1], d];

  11. If[Q[i + 1] == Q[i] || P[i + 1] == P[i] ||
  12.    FractionalPart[Sqrt[Q[i + 1]]] == 0 && EvenQ[i + 1] == True, Break[]];


  13. Print[{i, Q[i], P[i], a[i], p[i + 1]}] i++]

  14. {i, Q[i], P[i], a[i], p[i + 1]}

  15. {i + 1, Q[i + 1], P[i + 1], a[i + 1], p[i + 1]}

  16. GCD[Sqrt[Q[i + 1]] - p[i + 1], d]

  17. GCD[Sqrt[Q[i + 1]] + p[i + 1], d]
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2019-4-26 16:16:27 | 显示全部楼层
mathematica是函数式语言,把它当C之类的命令式语言写是一种糟糕的写法。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2019-5-25 00:47 , Processed in 0.062892 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

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