zoukankan      html  css  js  c++  java
  • 数论及其应用——积性函数问题

      在学习快速幂的过程中,我们曾遇到过因子和函数σ(n),曾提及该函数是积性函数,不过当时并没有给出证明。在这篇文章中,我们将针对数论中的积性函数问题,讨论更多的模型。包括欧拉函数、和我们曾讨论过的gcd函数。

      首先我们给出一些定义

      定义1:定义在所有正整数上的函数成为算数函数。

      定义2:算术函数f如果满足对于任意两个互素的正整数m、n,均由f(mn) = f(m)f(n),就称其为积性函数。如果对于任意的m、n满足上述性质,则称其为完全积性函数。

      因子和函数σ(n):

      定义:对于整数n所有因子的和,我们记做σ(n)。

      下面我们来考虑这样一个问题,我们如何利用编程快速计算σ(n)呢?

      法一:首先我们最容易想到的就是暴力法,枚举n的每一个因子嘛,然后求和。

      法二:既然涉及到因子,我们很容易想到与整数的拆分密切相关的素数基本定理,n = ∏(p[i]^ei),我们可以构造一个母函数 f = ∏∑pi^j (j∈[0,ei]),便可不重不漏的构造出整数n的所有因子,然后求和即可。

      法三:我们还是从暴力枚举的角度来思考,如何更加优化简洁的完成这个枚举过程。其实非常类似桶排序和埃及筛法的思想,我们用一维数组记录a[n]来记录σ(n),我们设置参量i,然后设置j来遍历i的整数倍,表示i是j的一个因子,然后i这个元素"扔到"a[j]这个桶里面,即a[j] += i。
      这种自然描述是为了更好解释算法的逻辑细节,虽然显得有些繁琐,其伪码过程如下。
      for i
         for j  2i to ki
             a[j] += i
      我们通过一个简单的题目来完成代码实现。(Problem source : hdu 1215)
     

    Problem Description
    七夕节那天,月老来到数字王国,他在城门上贴了一张告示,并且和数字王国的人们说:"你们想知道你们的另一半是谁吗?那就按照告示上的方法去找吧!" 人们纷纷来到告示前,都想知道谁才是自己的另一半.告示如下:
    数字N的因子就是所有比N小又能被N整除的所有正整数,如12的因子有1,2,3,4,6. 你想知道你的另一半吗?


      我们这里利用法三,将数据的最大值考虑进去,参考代码如下。

    #include<cstdio>
    using namespace std;
    
    const int maxn = 500002;
    int a[maxn] = {0 , 0};
    
    void fun()
    {
        for(int i = 2;i < maxn;++i)
              a[i] = 1;
        for(int i = 2;i < maxn/2;++i)
              for(int j = i + i;j < maxn;j += i)
                 a[j] += i;
    }
    
    int main()
    {
         fun();
         int ncase;
         scanf("%d", &ncase);
         while(ncase--)
         {
              int n;
              scanf("%d",&n);
              printf("%d
    ",a[n]);
         }
    
         return 0;
    }

     
     

      欧拉函数φ(n):

      首先,该函数的含义表示不超过n且与n互素的正整数的个数。基于对该函数的含义,我们会得出如下的定理。

      定理1:如果p是素数,则φ(p) = p - 1。

      定理2:如果p是素数,则φ(p^a) = p^a - p^(a-1)。

      证明:不超过p^a的数有p^a 个,这其中我们找到与p^a是非互素关系的数——p、2p、3p、4p……p^(a-1)p,有p^(a-1)个,由此得证。

      定理3:设n =∏p^ai是正整数n的素数幂分解,那么有φ(n) = n ∏(1 - 1/pi)。

      证明:首先基于算数基本定理,对于任意大于1的正整数我们都可以写成素数幂分解的形式,然后再基于定理2,通过化简整理即可得到欧拉函数的通式。

      基于我们给出的这三条定理,我们通过一个题目来具体实现欧拉函数值的求解。(Problem source : pku 2407)

      

    Description

    Given n, a positive integer, how many positive integers less than n are relatively prime to n? Two integers a and b are relatively prime if there are no integers x > 1, y > 0, z > 0 such that a = xy and b = xz.

    Input

    There are several test cases. For each test case, standard input contains a line with n <= 1,000,000,000. A line containing 0 follows the last case.

    Output

    For each test case there should be single line of output answering the question posed above.

      题目大意:典型的计算欧拉函数值的问题。

      编程实现:知晓了的计算通式,我们只需要通过简单的编程技巧来实现即可。

      实现计算欧拉函数值的方法有很多,这里我们给出直接实现的方法。即遍历出n所有的素因子然后套用公式,优化技巧很类似与我们在《数论及其应用——素数问题》中探讨过的,找素因子只需穷举<=sqrt(n),即可。

      参考代码如下。

    #include<cstdlib>
    #include<iostream>
    using namespace std;
    
    int phi(int n)
    {
          int rea = n;
            for(int i = 2;i*i <=n;i++)
    
                 if(n%i == 0)
                 {
                     rea = rea - rea/i;
                     do
                        n /= i;
                     while(n%i == 0);
                 }
                 if(n > 1)
                      rea = rea - rea/n;
                 return rea;
    
    }
    
    int main()
    {
           int n;
           while(cin >> n && n)
           {
                      cout << phi(n) << endl;
           }
           return 0;
    }

      我们来看一道应用到欧拉函数的题目。(Problem soure : pku 1284)
      

    Description

    We say that integer x, 0 < x < p, is a primitive root modulo odd prime p if and only if the set { (xi mod p) | 1 <= i <= p-1 } is equal to { 1, ..., p-1 }. For example, the consecutive powers of 3 modulo 7 are 3, 2, 6, 4, 5, 1, and thus 3 is a primitive root modulo 7. Write a program which given any odd prime 3 <= p < 65536 outputs the number of primitive roots modulo p.

    Input

    Each line of the input contains an odd prime numbers p. Input is terminated by the end-of-file seperator.

    Output

    For each p, print a single number that gives the number of primitive roots in a single line.


      题目大意:给出了原根的概念,然后给出一个奇素数p,然你求解奇素数p有多少个原根。
      数理分析:这里我们引用数论中的一条现成结论——p是素数,则p有φ(p-1)个原根,基于这条结论,我们就能很好的解决这个问题了。由于时间原因,这里依然暂且折叠这条结论的证明过程,笔者会在后续的深入学习中补充。
      参考代码如下。

    #include<cstdlib>
    #include<iostream>
    using namespace std;
    int main()
    {
         int p;
         while(cin >> p)
         {
               p--;
               int re = p;
               for(int i = 2;i*i <= p;i++)
            if(p%i == 0)
            {
                 re = re  - re/i;
                 do
                     p /= i;
                 while(p%i == 0);
            }
            if(p > 1)
                 re = re - re/p;
            cout << re << endl;
         }
         return 0;
    }


      让我们来看一道有关欧拉函数φ(n)的应用。(Problem source : hdu 2588)

      Problem Description
      The greatest common divisor GCD(a,b) of two positive integers a and b,sometimes written (a,b),is the largest divisor common to a and b,For example,(1,2)=1,(12,18)=6.    (a,b) can be easily found by the Euclidean algorithm. Now Carp is considering a little more difficult problem:    Given integers N and M, how many integer X satisfies 1<=X<=N and (X,N)>=M.


      题目大意:给定n,m,求解满足x<=n并且gcd(x,n)>=m的x的数目。
      数理分析:我们从gcd(x,n) >= m这个条件入手。设gcd(x,n) = d,则有x = p*d , n = q*d,其中gcd(p,q) = 1。因此当符合要求gcd(x,n)一旦确定,我们通过寻找不大于q且与q互素的个数来确定x的个数,而这刚好是欧拉函数φ(q)能够做的事情。
      下面我们需要做的便是找到所有满足要求的gcd(x,n),这很简单,通过试除法遍历n所有的因子,然后判断其是否满足大于等于m即可。
      分析到这里,你是否有疑惑,对于gcd(x1,n) = d1 , gcd(x2,n) = d2,在两种情况中我们是否会构造出相同的x?即是否需要筛除重复出现的x?
      看来我们需要对我们算法的正确性做一个有效的证明。
      证明:假设x存在重复的情况
       因此有x = p1*d1 = p2*d2.
       而记n = q1*d1 = q2*d2.
      p1/p2 = q1/q2
      这与gcd(p1,q1) = gcd(p2,q2) = 1矛盾,因此假设不成立。
     
      简单的参考代码如下。

    #include<cstdlib>
    #include<cstdio>
    #include<iostream>
    using namespace std;
    
    int phi(int n)
    {
          int rea = n;
            for(int i = 2;i*i <=n;i++)
    
                 if(n%i == 0)
                 {
                     rea = rea - rea/i;
                     do
                        n /= i;
                     while(n%i == 0);
                 }
                 if(n > 1)
                      rea = rea - rea/n;
                 return rea;
    
    }
    
    int main()
    {
           int t;
           while(scanf("%d",&t)!=EOF)
            {
    
           while(t--)
           {
                int sum = 0;
                int n , m;
                scanf("%d%d",&n,&m);
                 for(int i = 1;i*i <= n;i++)
                 {
                       if(n%i == 0)
                       {
                             if(i >=m)
                                sum += phi(n/i);
                             if(n/i != i && (n/i)>= m)
                                sum += phi(i);
                       }
                 }
                 printf("%d
    ",sum);
           }
            }
    
           return 0;
    }

        基于我们之前对欧拉函数的初涉,这里我们再给出它与另一个积性函数——gcd有关的定理。

        定理4:对于定值n,如果d是n的一个约数,则符合gcd(i,n) = d 的i的个数有phi(n/d)个,这里的phi即是欧拉函数。

        同样由于时间原因,笔者会在以后的拓展中给出该定理的证明,这里暂时不给出证明。

        基于这个定理,我们来通过一个问题具体的应用它。(Problem source : pku 2480)

    Description

    Longge is good at mathematics and he likes to think about hard mathematical problems which will be solved by some graceful algorithms. Now a problem comes: Given an integer N(1 < N < 2^31),you are to calculate   ∑gcd(i, N) 1<=i <=N. 
    "Oh, I know, I know!" Longge shouts! But do you know? Please solve it.

    Input

    Input contain several test case. A number N per line.

    Output

    For each N, output ,∑gcd(i, N) 1<=i <=N, a line
     
      题目大意:计算∑gcd(i, N) 1<=i <=N。
      数理分析:基于N的数量级,我们遍历求然后加和显然是不现实的,那么我们就考虑用积性函数的一些性质来简化计算公式。
      首先,我们先来讨论gcd函数的积性,若有gcd(p,q) = 1 且n = p*q , 那么我们可以想象到gcd(i,n) = gcd(i,p*q) = gcd(i,p)*gcd(i,q)。这里如果把gcd函数看做一种运算符,就是在验证gcd运算和乘法运算是否符合结合律,如果符合,那么证明gcd(n)是一个积性函数。
      我们先来看这样一个命题:设i与p的公因子是p',i与q的公因子是q',这表明i既有因子p'又有q',p*q也是这样,显然,p'*q' 是i与p*q的公因子。
      那么我们现在让p' = gcd(i,p),q'=gcd(i,q),gcd(p,q) = 1,显然此时的p'*q'是p*q所有的公因子中最大的,而且由于p、q互素,p'和q'也一定互素,因此保证了p'*q'是最大公因子,上述定理得证。
      之所以给出这条性质,是为了像我们在前文中讨论因子和函数的时候,对gcd(i,n)计算式进行分解以方便编程来实现求值。
       我们再回到这道题目上来,我们看到给出的是gcd函数累加的形式,即∑gcd(i, N),因此我们这里还需要积性函数的另外一条性质——积性函数的和依然是积性函数。
      则记f(n) = ∑gcd(i, n),f(n)也是积性函数。
      此时我们再基于算术基本定理,即可以把n表示成p1^a1 * p2 ^ a2 * p3 ^ a3 ……的形式,带入到f(n)当中去,即可用积性函数的性质将其展开。
      f(n) = f(p1^a1) * f(p2 ^a2)*f(p3^a3) *……  ①
      基于定理4,我们得到某个因子的个数,基于n = pi^ai这种形式,其因子我们是可遍历的(pi,pi^2,……pi^ai).
      这样我们即可展开f(ai^pi) = f(pi^ai) =  Φ(pi^ai)+pi*Φ(pi^(ai-1))+pi^2*Φ(pi^(ai-2))+...+pi^(ai-1)* Φ(pi)+ pi^ai *Φ(1)
      带入欧拉函数φ化简后,我们得到f(pi^ai) = ai*(p^r - p^(r-1)) + p^r,再将此式子用到等式①当中,便可求解。
      透彻理解了以上的数理分析,编程实现就是简单的模拟。
      参考代码如下。
     
    #include<cstdlib>
    #include<iostream>
    #include<stdio.h>
    using namespace std;
    int main()
    {
         long long s , n;
         while(scanf("%I64d",&n) != -1)
         {
               s = 1;
               long long x  ,r;
               for(long long i = 2;i*i <= n;i++)
                if(n%i == 0)
               {
                   x = 1;r = 0;
                   do
                    {
                         n /= i;
                         x *= i;
                         r++;
                    }while(n%i == 0);
                    s *= (r+1)*x- r*x/i;
               }
                if(n > 1)
                      s *= (2*n-1);
                printf("%I64d
    ",s);
    
         }
         return 0;
    }

      我们再来看一道简单的有关因子和函数的问题。(Problem source : hdu 1999)
      

    Problem Description
    s(n)是正整数n的真因子之和,即小于n且整除n的因子和.例如s(12)=1+2+3+4+6=16.如果任何 数m,s(m)都不等于n,则称n为不可摸数.
     
    Input
    包含多组数据,首先输入T,表示有T组数据.每组数据1行给出n(2<=n<=1000)是整数。
     
    Output
    如果n是不可摸数,输出yes,否则输出no


      题目表述很简单,其实也不涉及很巧妙的数论证明,类似基于一个映射关系,原象对应着象,这里的原像是m,像则是s(m),现在给出一个值n,是象集合当中的元素,问你整数n是否有对应的原像。
      对于这个问题较为抽象的提炼,进行简单的模拟即可解决问题。
      值得注意的是,模拟过程首先我们要基于原象和象的一一对应关系,观察到象的数据范围——[1,1000],原像的数据范围应该是多少呢?
      简单分析一下可知,小于1000最大的素数是997,那么对于原象997 * 997 , 它对应的原象是997  + 1,这就已经很接近象的最大值了,由此我们不难看到,原象的范围在[1,1000000]是最严谨的。
      基于以上分析,参考代码如下。

    #include<cstdio>
    #include<cmath>
    #include<cstring>
    bool h[1001];
    int main()
    {
         int i , j , n , s , o;
    
         memset(h , false , sizeof(h));
         h[1] = true;
         for(n = 4;n < 1000000;n++)
         {
                s = 1;
                o = (int)sqrt(double(n));
                for(i = 2;i <= o;i++)
                {
                      if(n%i == 0)
                      {
                            s += i;
                      if(n/i != i)
                           s += n/i;
                      }
                      if(s > 1001)
                           break;
                }
               if(s < 1001) h[s] = true;
    
    
         }
         int t;
         scanf("%d",&t);
         while(t--)
         {
               scanf("%d",&n);
                 if(h[n]==true)  printf("no
    ");
                 else            printf("yes
    ");
         }
    }


     

      ——未完

  • 相关阅读:
    Qt QTimer定时器相关
    C#Datetime和long之间转换
    C# 把图片资源转成字节数组写入到数据库
    Qt QProcess启动和关闭外部程序
    Qt绘图
    有哪些十分惊艳的书值得推荐3
    Stack Overflow 推荐编程书单
    《编写可读代码的艺术》的读书笔记 (脑图)
    [apache spark]洞见纽约车辆事故|bluemix|apache spark
    [lean scala]|How to create a SBT project with Intellij IDEA
  • 原文地址:https://www.cnblogs.com/rhythmic/p/5266724.html
Copyright © 2011-2022 走看看