欧拉函数+中国剩余定理

好久没写文章了。。

总结欧拉函数：

欧拉函数就是对于整数n，不超过n且与n互质的数的个数。

若mn为互质的正整数，那么f(mn)=f(m)*f(n);

因为任意正整数都可以唯一表示成如下形式:
                     k=p1^a1*p2^a2*……*pi^ai;(即分解质因数形式)
    可以推出:E(k)=(p1-1)(p2-1)……(pi-1)*(p1^(a1-1))(p2^(a2-1))……(pi^(ai-1))
               =k*(p1-1)(p2-1)……(pi-1)/(p1*p2*……pi);
               =k*(1-1/p1)*(1-1/p2)....(1-1/pk)
     ps：在程序中利用欧拉函数如下性质，可以快速求出欧拉函数的值(a为N的质因素)
若(N%a==0 && (N/a)%a==0) 则有:E(N)=E(N/a)*a;

若(N%a==0 && (N/a)%a!=0) 则有:E(N)=E(N/a)*(a-1);

要计算一个正整数n的欧拉函数的方法如下：
1. 将n表示成素数的乘积: n = p1 ^ k1 * p2 ^ k2 * ... * pn ^ kn（这里p1, p2, ..., pn是素数）
2. PHI(n) = (p1 ^ k1 - p1 ^ (k1 - 1)) * (p2 ^ k2 - p2 ^ (k2 - 1)) * ... *(pn ^ kn - pn ^ (kn - 1))
              = Mult { pi ^ ki - pi ^ (ki -1) }

证明过程如下：
1. 容易想到：当n为素数时，PHI(n) = n - 1。因为每个比n小的正整数都和n互素。当n为素数p的k次方时，PHI(n) = p ^ k - p ^ (k - 1)。因为在1到n之间的正整数只有p的倍数和n不互素，这样的数有(p ^ k / p)个。
2. 如果m和n互素，即GCD(m, n) = 1，那么PHI(m * n) = PHI(m) * PHI(n)。用中国剩余定理可以证明，证明的思路是建立这样一种一一对应的关系(a, b) <-> x，其中正整数a小于m并且gcd(a, m) = 1，正整数b小于n并且gcd(b, n) = 1，正整数x小于m*n并且gcd(m*n, x) = 1。证明过程如下：
    1)根据中国剩余定理，如果m和n互素，那么关于未知量x的方程组x % m = a, x % n = b(0 <= a < m, 0 <= b < n)，当0 <= x < m * n时存在并且仅存在一个解。容易证明，如果两个这样的方程组有相同的m, n但是a, b不同，那么他们的解x一定不同。
    2)首先用反正法证明：gcd(m, a) = 1且gcd(n, b) = 1是gcd(m*n, x) = 1的必要条件：假设gcd(a, m) = k > 1，由此可得：a = a' * k; m = m' * k => x = k' * m + a = k' * k * m' + k * a' = k * (k' * m' + a'); 所以gcd(x, m) = k > 1。同理可证，如果gcd(b, n) > 1， 那么gcd(x, n) > 1。所以x和m * n互素的必要条件是a和m互诉且b和n互素。
    3)接下来我们证明充分性：由x % m = a 可以得到x = k * m + a；由欧几里德算法求最大公约数的过程（就不证明了，呵呵，还得想）可以知道gcd(x, m) = gcd(m, a) = 1；同理可得，如果gcd(n, b) = 1那么gcd(x, n) = 1。接下来很容易得到：gcd(m*n, x) = 1。从而证明了充分性。
    4)上面三步的结论表明，数对(a, b)是可以和x建立起一一对应的关系的，所以有多少个不同的(a, b)，就有多少个不同的x。
3.将n分解成素数乘积后，显然对于任意的i, j(i != j)都满足 pi ^ ki和pj ^ kj是互素的，于是可以的到上面的公式。

跟据上面的公式，可以得到关于欧拉函数的递推关系：
假设素数p能整除n，那么
如果p还能整除n / p, PHI(n) = PHI(n / p) * p;
如果p不能整除n / p, PHI(n) = PHI(n / p) * (p - 1);

下面是两种求欧拉函数的不同编程方法： /*==================================================*\ |递推求欧拉函数phi(i) \*==================================================*/ for (i = 1; i <= maxn; i++) phi[i] = i; for (i = 2; i <= maxn; i += 2) phi[i] /= 2; for (i = 3; i <= maxn; i += 2) if(phi[i] == i) { for (j = i; j <= maxn; j += i) phi[j] = phi[j] / i * (i - 1); /*==================================================*\ |单独求欧拉函数phi(x) \*==================================================*/ unsigned euler(unsigned x) {// 就是公式 unsigned i, res=x; for (i = 2; i < (int)sqrt(x * 1.0) + 1; i++) if(x%i==0) { res = res / i * (i - 1); while (x % i == 0) x /= i; // 保证i一定是素数 } if (x > 1) res = res / x * (x - 1); return res; }

中国剩余定理：

对于n，已知n%m1=a1,n%m2=n1...一组同余方程，求n。

典经的、不同除数的同余式组解法

 

X≡R1  (mod  m₁ )    X≡2 (mod  7 )

X≡R2  (mod  m₂)     X≡5 (mod  9 )

X≡R3  (mod  m₃)     X≡1 (mod  5 )

   名词注释及计算步骤：

1 余数R：、R₁=2、R₂=5、R₃=1

2 模，亦即除数m：例中m₁=7、m₂=9、m₃=5

3 模的最小公倍数G：G=m₁*m₂*m₃，例中M=7*9*5=315

4 衍数（局部公倍数）y：Y₁=m₂m₃、Y₂=m₁m₃，Y₃=m₁m₂，例中Y₁=9*5=45、Y₂=7*5=35、Y₃=7*9=63

5 乘率C：这是解算中国剩余定理的关键，而计算“乘率”的方法，是秦九韶在《数书九章》一书中首次提出 的，称之为“大衍求一术”。“求一”就是使(衍数*乘率)除以模(除数)，而余数为1。即：

衍数Y*乘率C≡1  (mod m)，乘率C可以经过反算而得到。例中Y₁C₁≡1 (mod  7 )、

 Y₂C₂≡1 (mod  9 )、Y₃C₃≡1 (mod  5 )。

   计算C₁方法。由Y₁C₁≡1 (mod  7 )， →45C₁≡1 (mod  7 ) →(45C₁－1)  /  7=整数N，得C₁=5。因为45*5=225，225－1=224，224÷7=32，32是整数，合符要求。C₂、C₃之计算也相仿。乘率C之计算见下表：

这步可以用扩展欧几里得计算。

 

 

同余式 i	衍数Y	乘率C	余1	模m	检验  (Y*C-1)/m  = 整数
1	45	5	1	7	(45*C-1)/7 =N   (45*5-1)/7= 32
2	35	8	1	9	(35*C-1)/9 =N  (35*8-1)/9= 31
3	63	2	1	5	(63*C-1)/5 =N  (63*2-1)/5= 25

 

 

  

 

 

 

6  最终结果，X≡R₁Y₁C₁+R₂Y2C₂+R₃Y₃C₃ (mod G)

即X≡Σ余数*衍数*乘率 (mod G)，见下表计算：

 

 

i	余数R	衍数Y	乘率C	R*Y*C
1	2	45	5	450
2	5	35	8	1400
3	1	63	2	126
Σ	1976

 

 

 

 

 

 

 

X≡1976  (mod 315) 。1976 除去315的6倍后，剩下86，最终，

X≡86  (mod 315)

此为最小正整数解，通解为x+k×G。

代码：

int China()
{
    int M=1;
    for(int i=0;i<3;i++)
    {
        M*=m[i];
    }
    int x,y;
    int ans=0;
    for(int i=0;i<3;i++)
    {
        int mi=M/m[i];
        exgcd(mi,m[i],x,y);
        ans=(ans+mi*x*a[i])%M;
    }
    if(ans<=0)
        ans+=M;
    return ans;
}

void exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
        x=1,y=0;
        return ;
    }
    exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=(t-a/b*y);
}