找回密码
 欢迎注册
楼主: gracias

[求助] fibonacci(10^14) % 1234567891011是多少?

[复制链接]
发表于 2024-1-8 12:32:10 | 显示全部楼层
nyy 发表于 2024-1-8 12:16
用你的办法!

pari gp的代码

完全用mathematica的求解

  1. Clear["Global`*"];(*Clear all variables*)
  2. n=10^14;(*n次方幂*)
  3. m=1234567891011;(*模*)
  4. f=PolynomialMod[x^n,x^2-x-1,Modulus->m](*x^n对x^2-x-1取模,并且系数对m取模*)
  5. ans=Solve[x^2-x-1==0,{x}](*求解递推方程的两个根,第二个根是正根*)
  6. aa=(f/.ans[[2]])-(f/.ans[[1]])(*代入斐波那契数列的通项公式的分子*)
  7. bb=Mod[aa/Sqrt[5],m]//Simplify(*再除以分母模一下,化简一下,得到结果*)
复制代码


求解结果
512884989226 + 921144120792 x    这个是x^(10^14)模x^2-x-1并且系数对m取模后的结果,至于为什么算那么快,我也不知道

{{x -> 1/2 (1 - Sqrt[5])}, {x -> 1/2 (1 + Sqrt[5])}} 这个是递推公式的两个根,第二个根是正的

-460572060396 (1 - Sqrt[5]) + 460572060396 (1 + Sqrt[5]) 两个根代入相减的结果

921144120792 最终求解结果
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-8 13:04:42 | 显示全部楼层
本帖最后由 nyy 于 2024-1-8 13:11 编辑

归纳一下,这个题,有5种办法!
第1种:硬算,目前来说计算机硬件达不到要求,目前来说计算不了!
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41744
第2种:利用矩阵模幂算法
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41737
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41742
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41750
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41792
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=44184
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=44188
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=44189
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98654
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98668
第3种:利用lucas数列递推模幂
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98653
第4种:利用根号5的二次域模幂
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98667
第5种:先x^n这个多项式对x^2-x-1模,得到a*x+b,此时这个a、b数值是非常的巨大,所以对m模一下,变小了,剩下的简单了
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=41755
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98670
https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98671



补充内容 (2024-1-9 14:44):
方法6:mathe的办法,利用w^2=w+1等量代换与模幂算法,把w^n化成a+b*w(mod m)的形式,见https://bbs.emath.ac.cn/forum.ph ... =4054&pid=98687
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-8 15:09:56 | 显示全部楼层
gxqcn 发表于 2012-2-21 17:20
$(2\omega+1)(2\omega-3)=4(\omega^2-\omega-1)+1=1$,我怎么就没想到呢?
非常感谢 mathe 的解答,现在我 ...

为什么“你怎么就没想到”?

看代码与输出结果

  1.                   GP/PARI CALCULATOR Version 2.15.4 (released)
  2.           amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
  3.            compiled: Jun 28 2023, gcc version 10-posix 20210110 (GCC)
  4.                             threading engine: single
  5.                  (readline v8.0 enabled, extended help enabled)

  6.                      Copyright (C) 2000-2022 The PARI Group

  7. PARI/GP is free software, covered by the GNU General Public License, and comes
  8. WITHOUT ANY WARRANTY WHATSOEVER.

  9. Type ? for help, \q to quit.
  10. Type ?18 for how to get moral (and possibly technical) support.

  11. parisize = 8000000, primelimit = 500000
  12. (15:04) gp > w=quadgen(5)
  13. %1 = w
  14. (15:05) gp > (2*w+1)^(-1)
  15. %2 = -3 + 2*w
  16. Logging to pari.log
  17. GPRC Done.


  18. parisize = 8000000, primelimit = 500000
  19. (15:05) gp > Mod(2*x+1,x^2-x-1)^-1
  20. %1 = Mod(2*x - 3, x^2 - x - 1)
  21. Logging to pari.log
  22. GPRC Done.


  23. parisize = 8000000, primelimit = 500000
  24. (15:06) gp > Mod(1/(2*x+1),x^2-x-1)
  25. %1 = Mod(2*x - 3, x^2 - x - 1)
复制代码


这个证明了mathe比你聪明那么一点点

点评

nyy
提供了三种办法  发表于 2024-1-8 15:10
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-9 08:54:55 | 显示全部楼层
gxqcn 发表于 2012-2-21 17:20
$(2\omega+1)(2\omega-3)=4(\omega^2-\omega-1)+1=1$,我怎么就没想到呢?
非常感谢 mathe 的解答,现在我 ...


我来手把手地教你如何想到!待定系数法就可以!没有任何神秘的地方!

直接看代码与代码中的注释,很容易看明白!
  1. Clear["Global`*"];(*Clear all variables*)
  2. (*待定系数法求解多项式的逆元*)
  3. f=(2*x+1)(a*x+b)-(k*(x^2-x-1)+1)(*两个多项式相减,应该恒等*)
  4. ff=Collect[f,x](*合并同类项*)
  5. aa=CoefficientList[ff,x](*提取多项式系数*)
  6. bb=Solve[aa==0,{a,b,k}](*多项式恒等,所以系数都等于零,解方程组*)
复制代码



求解结果:
\[(2 x+1) (a x+b)-k \left(x^2-x-1\right)-1\]

合并同类项之后
\[x (a+2 b+k)+x^2 (2 a-k)+b+k-1\]

提取系数
\[\{b+k-1,a+2 b+k,2 a-k\}\]

让三个系数都等于零,解三个未知数,得到三元一次方程组,求解得到
\[\{\{a\to 2,b\to -3,k\to 4\}\}\]
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-9 09:31:21 | 显示全部楼层
本帖最后由 nyy 于 2024-1-9 09:51 编辑
nyy 发表于 2024-1-9 08:54
我来手把手地教你如何想到!待定系数法就可以!没有任何神秘的地方!

直接看代码与代码中的注释,很容易 ...


再来换一种办法,
把通项公式中的\(\sqrt{5}\)替换成x(x是方程x^2-5的根),则

\[a[n]=\frac{\left(\frac{1}{2} \left(1+\sqrt{5}\right)\right)^n-\left(\frac{1}{2} \left(1-\sqrt{5}\right)\right)^n}{\sqrt{5}}=\frac{\left(\frac{1+x}{2}\right)^n-\left(\frac{1-x}{2}\right)^n}{x}\]

现在对这个多项式模x^2-5,在模的过程中同时对系数取模(模m=1234567891011 ),然后求解结果

上pari的代码:
  1. n=10^14 /*n次幂*/
  2. m=1234567891011 /*模*/
  3. aa=Mod(1/x,x^2-5) /*求解分母上的根号5的逆*/
  4. p1=Mod(Mod(1/2,m)+Mod(1/2,m)*x,x^2-5)^n /*求解正根的n次幂*/
  5. p2=Mod(Mod(1/2,m)+Mod(-1/2,m)*x,x^2-5)^n /*求解负根的n次幂*/
  6. out=aa*(p1-p2) /*得到最终结果*/
复制代码


上求解结果
  1.                   GP/PARI CALCULATOR Version 2.15.4 (released)
  2.           amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
  3.            compiled: Jun 28 2023, gcc version 10-posix 20210110 (GCC)
  4.                             threading engine: single
  5.                  (readline v8.0 enabled, extended help enabled)

  6.                      Copyright (C) 2000-2022 The PARI Group

  7. PARI/GP is free software, covered by the GNU General Public License, and comes
  8. WITHOUT ANY WARRANTY WHATSOEVER.

  9. Type ? for help, \q to quit.
  10. Type ?18 for how to get moral (and possibly technical) support.

  11. parisize = 8000000, primelimit = 500000
  12. (09:28) gp > n=10^14
  13. %1 = 100000000000000
  14. (09:28) gp > m=1234567891011
  15. %2 = 1234567891011
  16. (09:28) gp > aa=Mod(1/x,x^2-5)
  17. %3 = Mod(1/5*x, x^2 - 5)
  18. (09:28) gp > p1=Mod(Mod(1/2,m)+Mod(1/2,m)*x,x^2-5)^n
  19. %4 = Mod(Mod(460572060396, 1234567891011)*x + Mod(973457049622, 1234567891011), x^2 - 5)
  20. (09:28) gp > p2=Mod(Mod(1/2,m)+Mod(-1/2,m)*x,x^2-5)^n
  21. %5 = Mod(Mod(773995830615, 1234567891011)*x + Mod(973457049622, 1234567891011), x^2 - 5)
  22. (09:28) gp > out=aa*(p1-p2)
  23. %6 = Mod(Mod(921144120792, 1234567891011), x^2 - 5)
复制代码


看最后一行 Mod(Mod(921144120792, 1234567891011), x^2 - 5)
这里面的921144120792就是我们需要求解的结果!
不是太能说明白其中的原理,但是我感觉这样肯定没错!

注:其实分母上的根号5不用求逆,看p1-p2的求解结果
(09:49) gp > p1-p2
%7 = Mod(Mod(921144120792, 1234567891011)*x, x^2 - 5)
正好有一个x,再除以x,正好约掉,得到Mod(921144120792, 1234567891011),得到921144120792

点评

nyy
模m不能是2的倍数,因为p1=Mod(Mod(1/2,m)+Mod(1/2,m)*x,x^2-5)^n中的Mod(1/2,m)求解不了  发表于 2024-1-9 10:07
nyy
可以是5的倍数,但是只要不求Mod(1/x,x^2-5)这个就行,直接把x放在分母上,因为后来的p1-p2在分子上,上面有一个x,正好约掉  发表于 2024-1-9 10:07
nyy
这个模m也不能是5的倍数  发表于 2024-1-9 09:57
nyy
这个方法的缺点是模m不能是偶数  发表于 2024-1-9 09:45
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-9 10:55:25 | 显示全部楼层
nyy 发表于 2024-1-8 12:16
用你的办法!

pari gp的代码

利用一下mathe的思路,

上pari的代码
  1. n=10^14                           /* n次幂                                                        */
  2. m=1234567891011                   /* 模                                                           */
  3. p1=Mod(Mod(1,m)*x,x^2-x-1)^n      /* x=(sqrt(5)+1)/2,求解正根(=x)的n次幂                          */
  4. p2=Mod(Mod(1,m)*(-1/x),x^2-x-1)^n /* 求解负根(=-1/x)的n次幂                                       */
  5. out=(p1-p2)/Mod(2*x-1,x^2-x-1)    /* 再除以分母sqrt(5)=2*x-1,得到最终求解结果,应该是个a*x的多项式 */
复制代码


输出求解结果
  1. parisize = 8000000, primelimit = 500000
  2. (10:53) gp > n=10^14
  3. %1 = 100000000000000
  4. (10:53) gp > m=1234567891011
  5. %2 = 1234567891011
  6. (10:53) gp > p1=Mod(Mod(1,m)*x,x^2-x-1)^n
  7. %3 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
  8. (10:53) gp > p2=Mod(Mod(1,m)*(-1/x),x^2-x-1)^n
  9. %4 = Mod(Mod(313423770219, 1234567891011)*x + Mod(199461219007, 1234567891011), x^2 - x - 1)
  10. (10:53) gp > out=(p1-p2)/Mod(2*x-1,x^2-x-1)
  11. %5 = Mod(Mod(921144120792, 1234567891011), x^2 - x - 1)
复制代码


这样只需要借助pari,不需要借助mathematica就完成了求解。
%5 = Mod(Mod(921144120792, 1234567891011), x^2 - x - 1)这里面的921144120792就是最终求解结果

点评

nyy
这个思路的坏处是模是5的倍数不能用  发表于 2024-1-9 11:15
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-9 12:10:28 | 显示全部楼层
mathe 发表于 2012-2-21 08:48
其实写成$1/{sqrt{5}}(\omega^n+(-\omega)^{-n})$也可以用模幂运算
只是这时需要符号运算比较好。在模不是5 ...

我来全面地阐述一下你的思想,完全按照你的思路来。

w^2-w-1=0
先来定义一下乘法
(x1 + y1*w)*(x2 + y2*w)
=x1 x2 + w^2 y1 y2 + w (x2 y1 + x1 y2)
=x1 x2 + (1+w) y1 y2 + w (x2 y1 + x1 y2)
=x1 x2 + y1 y2 + w (x2 y1 + x1 y2 + y1 y2)
对于模幂运算来说,这个很重要。
下面计算w^n模m的结果,

mathematica代码如下:
  1. Clear["Global`*"];(*Clear all variables*)
  2. (*子函数,二次域w的模m的乘法.xx={x1,y1},表示x1+y1*w(w^2-w-1=0),yy同,m为整数*)
  3. dmult[xx_,yy_,m_]:=Module[{x1=xx[[1]],y1=xx[[2]],x2=yy[[1]],y2=yy[[2]]},Mod[{x1*x2+y1*y2,x2*y1+x1*y2+y1*y2},m]]
  4. n=10^14;(*n次方幂*)
  5. m=1234567891011;(*模*)
  6. n2=IntegerDigits[n,2];(*n次方,写成二进制表的形式*)
  7. result=base={0,1};(*初始赋值,这个表示0+1*w=w,w^2-w-1=0*)
  8. Do[result=dmult[result,result,m];(*平方一下,下标翻倍*)
  9.    If[n2[[k]]==1,result=dmult[result,base,m]](*如果是1,那么下标增加1*)
  10. ,{k,2,Length[n2]}];(*k就是从2开始的,没错!*)
  11. result (*w^n模m后的求解结果*)
复制代码


求解结果如下:
{512884989226, 921144120792}
这个结果表示512884989226+w*921144120792
因此
\[\frac{\text{w1}^n-\text{w2}^n}{\sqrt{5}}=\frac{(512884989226+921144120792 \text{w1})-(512884989226+921144120792 \text{w2})}{\sqrt{5}}=\frac{921144120792 (\text{w1}-\text{w2})}{\sqrt{5}}=921144120792\pmod{m}\]
其中假设w1>w2,且w1-w2=根号5
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-9 13:42:46 | 显示全部楼层
nyy 发表于 2024-1-8 09:20
还可以减少几行代码

利用pari计算如下:

  1. n=10^14         /* n次幂    */
  2. m=1234567891011 /* 模       */
  3. mat=[1,1;1,0]   /* 递推矩阵 */
  4. Mod(mat,m)^n    /* 矩阵模幂 */
复制代码


log文件如下(包含输入与输出结果)
  1. Logging to pari.log
  2. GPRC Done.

  3.                   GP/PARI CALCULATOR Version 2.15.4 (released)
  4.           amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
  5.            compiled: Jun 28 2023, gcc version 10-posix 20210110 (GCC)
  6.                             threading engine: single
  7.                  (readline v8.0 enabled, extended help enabled)

  8.                      Copyright (C) 2000-2022 The PARI Group

  9. PARI/GP is free software, covered by the GNU General Public License, and comes
  10. WITHOUT ANY WARRANTY WHATSOEVER.

  11. Type ? for help, \q to quit.
  12. Type ?18 for how to get moral (and possibly technical) support.

  13. parisize = 8000000, primelimit = 500000
  14. (13:40) gp > n=10^14
  15. %1 = 100000000000000
  16. (13:41) gp > m=1234567891011
  17. %2 = 1234567891011
  18. (13:41) gp > mat=[1,1;1,0]
  19. %3 =
  20. [1 1]

  21. [1 0]

  22. (13:41) gp > Mod(mat,m)^n
  23. %4 =
  24. [Mod(199461219007, 1234567891011) Mod(921144120792, 1234567891011)]

  25. [Mod(921144120792, 1234567891011) Mod(512884989226, 1234567891011)]

复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-10 09:07:58 | 显示全部楼层
nyy 发表于 2024-1-9 12:10
我来全面地阐述一下你的思想,完全按照你的思路来。

w^2-w-1=0

用mathe的思路,来处理分母上的\(\sqrt{5}\)这个问题。

上pari的代码
  1. n=10^14 /*n次幂*/
  2. m=1234567891011 /*模*/
  3. c1=Mod(1/(2*x-1),x^2-x-1) /*(1+sqrt(5))/2=x1,故1/sqrt(5)=1/(2*x1-1)*/
  4. p1=Mod(Mod(1,m)*x,x^2-x-1)^n /*计算x1^n mod x^2-x-1 */
  5. c2=Mod(1/(2*x-1),x^2-x-1) /*(1-sqrt(5))/2=x2,故-1/sqrt(5)=1/(2*x2-1)*/
  6. p2=Mod(Mod(1,m)*x,x^2-x-1)^n /*计算x2^n mod x^2-x-1 */
  7. cp1=c1*p1 /*计算通项公式的第1部分(正根)*/
  8. cp2=c2*p2 /*计算通项公式的第2部分(负根)*/
复制代码


输出结果log
  1.                   GP/PARI CALCULATOR Version 2.15.4 (released)
  2.           amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
  3.            compiled: Jun 28 2023, gcc version 10-posix 20210110 (GCC)
  4.                             threading engine: single
  5.                  (readline v8.0 enabled, extended help enabled)

  6.                      Copyright (C) 2000-2022 The PARI Group

  7. PARI/GP is free software, covered by the GNU General Public License, and comes
  8. WITHOUT ANY WARRANTY WHATSOEVER.

  9. Type ? for help, \q to quit.
  10. Type ?18 for how to get moral (and possibly technical) support.

  11. parisize = 8000000, primelimit = 500000
  12. (09:01) gp > n=10^14
  13. %1 = 100000000000000
  14. (09:02) gp > m=1234567891011
  15. %2 = 1234567891011
  16. (09:02) gp > c1=Mod(1/(2*x-1),x^2-x-1)
  17. %3 = Mod(2/5*x - 1/5, x^2 - x - 1)
  18. (09:02) gp > p1=Mod(Mod(1,m)*x,x^2-x-1)^n
  19. %4 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
  20. (09:02) gp > c2=Mod(1/(2*x-1),x^2-x-1)
  21. %5 = Mod(2/5*x - 1/5, x^2 - x - 1)
  22. (09:02) gp > p2=Mod(Mod(1,m)*x,x^2-x-1)^n
  23. %6 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
  24. (09:02) gp > cp1=c1*p1
  25. %7 = Mod(Mod(636296398051, 1234567891011)*x + Mod(759707806876, 1234567891011), x^2 - x - 1)
  26. (09:02) gp > cp2=c2*p2
  27. %8 = Mod(Mod(636296398051, 1234567891011)*x + Mod(759707806876, 1234567891011), x^2 - x - 1)
复制代码


从上面的结果可以看出cp1为(此处的x为x1,也就是正根)
Mod(Mod(636296398051, 1234567891011)*x + Mod(759707806876, 1234567891011), x^2 - x - 1)
cp2也是(此处的x为x2,也就是负根)
Mod(Mod(636296398051, 1234567891011)*x + Mod(759707806876, 1234567891011), x^2 - x - 1)
而根据韦达定理x1+x2=1(x1与x2是x^2-x-1=0的两个根)
Mod(636296398051, 1234567891011)*x1 + Mod(759707806876, 1234567891011)+Mod(636296398051, 1234567891011)*x2 + Mod(759707806876, 1234567891011)
=Mod(636296398051, 1234567891011)*(x1+x2) + Mod(759707806876, 1234567891011)*2
=Mod(636296398051, 1234567891011)*1 + Mod(759707806876, 1234567891011)*2
=Mod(921144120792, 1234567891011) 这个就是我们需要的结果
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2024-1-10 09:29:25 | 显示全部楼层
nyy 发表于 2024-1-10 09:07
用mathe的思路,来处理分母上的\(\sqrt{5}\)这个问题。

上pari的代码

还有更简单的。

  1. n=10^14 /*n次幂*/
  2. m=1234567891011 /*模*/
  3. Mod(Mod(1,m)*x,x^2-x-1)^n
复制代码


输出结果
  1.                   GP/PARI CALCULATOR Version 2.15.4 (released)
  2.           amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
  3.            compiled: Jun 28 2023, gcc version 10-posix 20210110 (GCC)
  4.                             threading engine: single
  5.                  (readline v8.0 enabled, extended help enabled)

  6.                      Copyright (C) 2000-2022 The PARI Group

  7. PARI/GP is free software, covered by the GNU General Public License, and comes
  8. WITHOUT ANY WARRANTY WHATSOEVER.

  9. Type ? for help, \q to quit.
  10. Type ?18 for how to get moral (and possibly technical) support.

  11. parisize = 8000000, primelimit = 500000
  12. (09:27) gp > n=10^14
  13. %1 = 100000000000000
  14. (09:27) gp > m=1234567891011
  15. %2 = 1234567891011
  16. (09:27) gp > Mod(Mod(1,m)*x,x^2-x-1)^n
  17. %3 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
复制代码


上面的921144120792就是需要的结果。
原理是什么?

http://pari.math.u-bordeaux.fr/gp.html
这个页面下的
  1. Mod(x,x^2-x-1)^100000==fibonacci(100000)*x+fibonacci(99999)
复制代码


这个网页是在线使用pari/gp

点评

nyy
如果w^2-w-1=0,则w^n=fibonacci(n-1) +fibonacci(n)*w,用数学归纳法很容易证明  发表于 2024-1-10 11:39
nyy
从上面可以知道Mod(x,x^2-x-1)^n==fibonacci(n)*x+fibonacci(n-1)  发表于 2024-1-10 09:49
nyy
这个在线使用pari/gp的网页真不错,还有很多例子!  发表于 2024-1-10 09:32
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-5-2 17:09 , Processed in 0.045175 second(s), 15 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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