求n的阶乘的算法框图_欧拉计划381|Wilson定理+欧氏算法

关于欧拉计划的介绍请阅读上文：

Pika369：欧拉计划1,3|Python与Mathematica的比较​zhuanlan.zhihu.com

题目：

分析：

题目中定义了一种函数叫做

,而最终结果是要计算

其中p是从5到

这之间，所有的素数。

首先生成5到

之间的素数本身就已经会消耗很多的时间，这就要求计算单个S(p)的时间非常快才可以。

计算S(p)

1.如果直接计算?

如果直接计算出阶乘，然后再取模运算显然会非常慢，因为阶乘增长得极快。做一些简单的处理也是没有实质上的效果的，因为，当计算到10^5以下的S(p)的时候基本就需要消耗大量的时间了。

比如想到一边做模运算一边做乘法，最后直接得到(p-1)! mod p的结果避免计算阶乘的时候遇到大数乘法浪费时间，或者记忆计算的结果，这样当得到(p-1)! mod p结果的时候，实际上(p-2)! mod p等等结果就已经得到了。

等等以上的简单优化都不会有实质性的进展。

真正对这个问题最快的优化方式还是要来自于数学，Wilson定理。

2.Wilson定理对这个问题的优化

等价于 n 是素数

证明见Wiki：

https://en.wikipedia.org/wiki/Wilson's_theorem​en.wikipedia.org

上面的定理基本数学专业，非数学专业都能看。如果对代数了解较少，可以看初等证明，对抽象代数有了解可以看通过Sylow定理构造的证明。

引用Wilson定理这里的目的比较明显，就是看出来如果 p是素数，那么

.

这样我们似乎就可以不用计算阶乘也能得到结果。接下来我们试着推导(p-2)!的结果。

这里约定以下符号，

表示的是x做模p的运算以后，并且结果取

中的一个。也就是域

中的某个元素。并且有运算：

例子1：在模13的意义下，

于是：

用上面的符号表示，则由Wilson定理，

。假设

,则有

在

中消去律是可以允许的，因此

。(本质上是因为

这个方程式有唯一解的，当(a,n)=1的时候。此处

,

是模p以后的结果，p是素数，所以显然和p互素。)

那么这个No.381题就是一个纯粹的数学问题吗？我们是不是可以根本不用什么编程算法就解决呢？

不是的！在计算

的时候会遇到问题。

假设

因为

,所以有：

要解这个方程实际上就是要解

，这就相当于找在上面定义的

中的模的乘法意义下，

的逆元素，或者我们可以表示为

.

即

和上面同样的方式去推导，不难发现计算S(p)的主要任务就是计算这些元素的逆元素。

所以现在的关键就是计算

中：

计算这些元素的逆会用到欧氏算法求最大公因数过程中的一个副产物：

Bezout恒等式：
对于整数a,b,存在整数s,t使得

其中(a,b)表示整数a,b的最大公因数

之所以是欧氏算法的副产物，是因为一开始欧氏算法是求a,b的最大公因数的，但是这个过程也可以计算出s,t。计算最大公因数的时候，s,t是没有用处的，因此也不会输出，不过在计算元素在

中的逆的时候是非常有用的。

令b=p,(a,p)=1,则有：

也就是说，只要知道a对应的参数s就直接求得到了

的逆。

我们可以通过下面的算法得到s。

3.欧氏算法的一个变种(Extended Euclidean Algorithm)

先算一个例子，复习一下欧氏算法求最大公因数。

例子2：用欧氏算法求(13,97)

由此可知(2,97)=1

可是当其中一个数为素数的时候，欧氏算法计算最大公因数倒没什么用处。现在的问题是，例2中的过程，可以得到参数s吗？

例子3：在例子2的基础上，得到参数s,t

因此

实际上，得到参数s,t的过程和得到最大公因数的过程及其类似，是完全一样的递推式，只不过初值不相同。

假设，a>b

由第三行式子不难发现一个递推式，

，此处定义

这个递推式定义了

,并且当

的时候，

就是a,b的最大公因数。

在欧氏算法当中，最后一行式子，

中，如果通过递推消去

就会得到一个保留

组成的表达式，具有形式：

由欧氏算法中的递推我们可以得到关于s,t的递推：

其中由欧氏算法求最大公因数的第一行式子

，可以知道，

由上述递推，不妨假设

请注意这里,s,t的次序是参照a,b确定的，即a对应的参数为s,b对应的参数为t。而a,b的选择，约定的是较大的数字为a,较小的数字为b。

因此在本文体中，取

中元素的逆，p显然是更大的，于是

,因此求某个数x的逆，实际上是令

求

.

bezoutT[a_, b_] := Block[{A, B, r, t, rOld, tOld, q},
  (*t is the parameter corresponding to the smaller number in {a,b}*)
  
  If[a < b, A = b; B = a;, A = a; B = b;];
  (*initial values*)
  {rOld, r} = {A, B};
  {tOld, t} = {0, 1};
  (*When r =0 we have got the GCD of {a,b},which is rOld*)
  (*Thus the parameter we want is tOld*)
  While[r != 0, 
   q = Floor[rOld/r]; (*q is the quotient*)
   {rOld, r} = {r, rOld - q*r};   
   {tOld, t} = {t, tOld - q*t};
   ];
  tOld
  ]

不过，Mathematica 12自带的函数ExtendedGCD[]计算得更快，下面的最终版本的代码，用的是MMA12内置函数。

4.Mathematica12代码

S

最后的运行时间大致是，78s。还可以更快，比如ExtendedGCD当中实际上最大公因数，s,t三个都求出来了，而实际上求元素的逆只需要s一个输出。所以减少一些其他计算应该会更快。