fibonacci(10^14) % 1234567891011是多少? - 第6页 - 算法交流 - 数学研发论坛

nyy 发表于 2024-1-8 12:32:10

nyy 发表于 2024-1-8 12:16
用你的办法！

pari gp的代码

完全用mathematica的求解

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

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

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

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

921144120792 最终求解结果

nyy 发表于 2024-1-8 13:04:42

本帖最后由 nyy 于 2024-1-8 13:11 编辑

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

nyy 发表于 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 的解答，现在我 ...

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

看代码与输出结果

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

                  Copyright (C) 2000-2022 The PARI Group

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

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

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

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

parisize = 8000000, primelimit = 500000
(15:06) gp > Mod(1/(2*x+1),x^2-x-1)
%1 = Mod(2*x - 3, x^2 - x - 1)

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

nyy 发表于 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 的解答，现在我 ...

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

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

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

合并同类项之后
\

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

让三个系数都等于零，解三个未知数，得到三元一次方程组，求解得到
\[\{\{a\to 2,b\to -3,k\to 4\}\}\]

nyy 发表于 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的根），则

\=\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的代码：
n=10^14 /*n次幂*/
m=1234567891011 /*模*/
aa=Mod(1/x,x^2-5) /*求解分母上的根号5的逆*/
p1=Mod(Mod(1/2,m)+Mod(1/2,m)*x,x^2-5)^n /*求解正根的n次幂*/
p2=Mod(Mod(1/2,m)+Mod(-1/2,m)*x,x^2-5)^n /*求解负根的n次幂*/
out=aa*(p1-p2) /*得到最终结果*/

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

                  Copyright (C) 2000-2022 The PARI Group

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

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

parisize = 8000000, primelimit = 500000
(09:28) gp > n=10^14
%1 = 100000000000000
(09:28) gp > m=1234567891011
%2 = 1234567891011
(09:28) gp > aa=Mod(1/x,x^2-5)
%3 = Mod(1/5*x, x^2 - 5)
(09:28) gp > p1=Mod(Mod(1/2,m)+Mod(1/2,m)*x,x^2-5)^n
%4 = Mod(Mod(460572060396, 1234567891011)*x + Mod(973457049622, 1234567891011), x^2 - 5)
(09:28) gp > p2=Mod(Mod(1/2,m)+Mod(-1/2,m)*x,x^2-5)^n
%5 = Mod(Mod(773995830615, 1234567891011)*x + Mod(973457049622, 1234567891011), x^2 - 5)
(09:28) gp > out=aa*(p1-p2)
%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 发表于 2024-1-9 10:55:25

nyy 发表于 2024-1-8 12:16
用你的办法！

pari gp的代码

利用一下mathe的思路，

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

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

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

nyy 发表于 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代码如下：
Clear["Global`*"];(*Clear all variables*)
(*子函数,二次域w的模m的乘法.xx={x1,y1},表示x1+y1*w(w^2-w-1=0),yy同,m为整数*)
dmult:=Module[{x1=xx[],y1=xx[],x2=yy[],y2=yy[]},Mod[{x1*x2+y1*y2,x2*y1+x1*y2+y1*y2},m]]
n=10^14;(*n次方幂*)
m=1234567891011;(*模*)
n2=IntegerDigits;(*n次方,写成二进制表的形式*)
result=base={0,1};(*初始赋值,这个表示0+1*w=w,w^2-w-1=0*)
Do;(*平方一下,下标翻倍*)
If]==1,result=dmult](*如果是1,那么下标增加1*)
,{k,2,Length}];(*k就是从2开始的,没错!*)
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

nyy 发表于 2024-1-9 13:42:46

nyy 发表于 2024-1-8 09:20
还可以减少几行代码

利用pari计算如下：

n=10^14       /* n次幂 */
m=1234567891011 /* 模    */
mat= /* 递推矩阵 */
Mod(mat,m)^n /* 矩阵模幂 */

log文件如下（包含输入与输出结果）
Logging to pari.log
GPRC Done.

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

                  Copyright (C) 2000-2022 The PARI Group

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

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

parisize = 8000000, primelimit = 500000
(13:40) gp > n=10^14
%1 = 100000000000000
(13:41) gp > m=1234567891011
%2 = 1234567891011
(13:41) gp > mat=
%3 =

(13:41) gp > Mod(mat,m)^n
%4 =

nyy 发表于 2024-1-10 09:07:58

nyy 发表于 2024-1-9 12:10
我来全面地阐述一下你的思想，完全按照你的思路来。

w^2-w-1=0

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

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

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

                  Copyright (C) 2000-2022 The PARI Group

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

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

parisize = 8000000, primelimit = 500000
(09:01) gp > n=10^14
%1 = 100000000000000
(09:02) gp > m=1234567891011
%2 = 1234567891011
(09:02) gp > c1=Mod(1/(2*x-1),x^2-x-1)
%3 = Mod(2/5*x - 1/5, x^2 - x - 1)
(09:02) gp > p1=Mod(Mod(1,m)*x,x^2-x-1)^n
%4 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
(09:02) gp > c2=Mod(1/(2*x-1),x^2-x-1)
%5 = Mod(2/5*x - 1/5, x^2 - x - 1)
(09:02) gp > p2=Mod(Mod(1,m)*x,x^2-x-1)^n
%6 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)
(09:02) gp > cp1=c1*p1
%7 = Mod(Mod(636296398051, 1234567891011)*x + Mod(759707806876, 1234567891011), x^2 - x - 1)
(09:02) gp > cp2=c2*p2
%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) 这个就是我们需要的结果

nyy 发表于 2024-1-10 09:29:25

nyy 发表于 2024-1-10 09:07
用mathe的思路，来处理分母上的$\sqrt{5}$这个问题。

上pari的代码

还有更简单的。

n=10^14 /*n次幂*/
m=1234567891011 /*模*/
Mod(Mod(1,m)*x,x^2-x-1)^n

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

                  Copyright (C) 2000-2022 The PARI Group

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

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

parisize = 8000000, primelimit = 500000
(09:27) gp > n=10^14
%1 = 100000000000000
(09:27) gp > m=1234567891011
%2 = 1234567891011
(09:27) gp > Mod(Mod(1,m)*x,x^2-x-1)^n
%3 = Mod(Mod(921144120792, 1234567891011)*x + Mod(512884989226, 1234567891011), x^2 - x - 1)

上面的921144120792就是需要的结果。
原理是什么？
见
http://pari.math.u-bordeaux.fr/gp.html
这个页面下的
Mod(x,x^2-x-1)^100000==fibonacci(100000)*x+fibonacci(99999)

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

页: 1 2 3 4 5 [6] 7 8

数学研发论坛's Archiver