《算法导论》——数论

1 . 欧几里得算法(递归法球两个数的最大公约数)

  算法比较简单就直接贴代码了:

    int gcd(int a , int b){          

       return b ==0 ? a : gcd(b , a%b);     

    }

  在这个算法的基础上可以该进,gcd(a*n , b*n)=n*gcd(a , b);这个定理给我们提供了思路:可不可以在算法中加如判断减少递归的次数呢?当然可以。

int gac(int m,int n){

   if(m==0 || n==0)   return (m==0) ? n : m;

  if(m<n)

    return gac(n,m);     //m<n时返回gac(n,m)的最大公约数

  else{

      if(m%2 !=1){

        if(n%2 !=1)     //两个数都是k的倍数时函数返回gac(m/k,n/k)的最大公约数

          return 2*gac(m/2,n/2);

        else

          return gac(m/2,n);   //两个数中若只有一个数是k的倍数则返回gac(m/k,n)的最大公约数

        }

       else{

        if(n%2!=1)     return gac(m,n/2);   //两个数中若只有一个数是k的倍数则返回gac(m,n/k)的最大公约数

        else     return gac(n,m-n);   //两个数都不是k的倍数则返回两个数中小的数和两个数的差的最大公约数

         }

    }

}

2 . 扩展形式的欧几里得算法,再求两个数的最大公约数的时候,还有可以利用的信息,假设有这样的一个方程:d=a*x+b*y ; 怎么求? 观察一下,不难发现可以令 d=gcd(a , b) , 由于gcd(a , b)=gcd(b , a%b ),令d’=gcd(b ,a%b),那 d’=gcd(b ,a%b)=b*x’+a%b*y’,由于a%b=a -(a/b)*b ,于是有等式 b*x’+(a -(a/b)*b)*y’=d=a*x+b*y成立,在合并同类项:d=a*y‘+b*(x’-(a/b)*y‘),选择x=y‘,y=x’-(a/b)*y‘;就是方程的一组解。

  struct equation{
    int d=0;
    int x=0;
    int y=0;
};

equation extend_euclid(int a , int b){    //利用欧几里得算法求解二元方程的一组解
    static  equation result;              //gcd(a , b)=gcd(b , a % b)
    if(b == 0){                           //
        result.d=a;
        result.x=1;
        result.y=0;
        return result;
    }else {
        equation temp=extend_euclid(b , a%b);
        result.d=temp.d;
        result.x=temp.y;
        result.y=temp.x - (a / b) * temp.y;
        return result;
    }
}

3 . 模线性方程:a*x=b(mod n)

  a.根据《算法导论》中数论的讨论中对这类方程的解法,第一步就是确定该方程有没有解,另d=gcd(a , n),当d整除b的时候该方程才会有解,并且其中一个特解 x0=x’(b/d) % n ;其中x’为gcd(a , n)=a*x+n*y的解;

  b.假设方程有解,并且x‘是该方程的一个解,那么该方程有d个解,并且都可以用x‘表示出来,x[i]=x‘+i*(n/d);其中i从0增长到d-1;

void modular_linear_equation_solver(int a ,int  b ,int  n){

  equation result=extend_euclid(a , n) ;

if(b % result.d == 0){

     int x=(result.x*(b/result.d)) % n;

     for(int i=0 ; i<result.d ; ++i)

      printf("%d ",(x+i*(n/result.d))%n);

   }

  else

    printf("no solutions\n");

}

4. 反复平方法求元素的幂: 

 int   repede_pow(int a , int b){
      int result=1;
      for(int i = 0 ; i < 64 ; ++i ){
          result *= result ;                      //反复平方
          if(b & 1 << i)                          //考虑到b为奇数的情况
              result *= a ;
      }
      return result ;
 }

5. 反复平方法幂的模运算:

long long int modular_exponentiation(long long a , long long b , long long n){   //反复平方法求元素幂

   long long result=1;

  long long c=0;

   long long temp = n ;

  int  bit_size = 0 ;

while(temp){

    ++bit_size ;

     temp >>= 1 ;

  }

for(int i = bit_size ; i >= 0 ; --i){

      c <<= 1 ;                                                               //通过c的成倍增长,达到反复平方

      result *= result ;

      result %= n ;

      if(n << i & 1){                                                         //考虑到奇数的情况

        c += 1 ;

         result *= a ;

         result %= n ;

      }

  }

  return result ;

}

6 . 好了下面说一下素数判定的问题

  a.暴力测试;

  b.miller_rabin测试法;

  暴力测试方法不用多说,说一下miller_rabin测试方法吧,基于前面的快速幂模运算,和乘积模运算,我们可以多次测试多次测试一个数是不是和数判定这个数是不是素数。

  可以编程实现:

#include<stdio.h>

long long multi(long long a , long long b , long long n){

    long long result=0;

     a = a % n;

    while(b){

      if(b & 1){

        result = result + a ;

        if(result >= n)                 result -= n;

      }

        a <<= 1;

       if(a > n)             a -= n;

       b >>= 1;

    }

    return result;

}

long long pow_mod(long long a , long long b , long long n){

    if(b == 0)

      return 1;

     if(b == 1)

       return a%n;

    long long answer=1;

    while(b){

       if(b & 1)             answer = multi(answer , a , n);             //answer = answer * a % n;

       b >>= 1;         //a = a * a % n;

      a = multi(a , a , n);

    }

    return answer;

}

bool witness(long long a , long long n){

     if(n == 2 || n == a)         return 0;

     if((n & 1) == 0)         return 1;

   long long u=n-1;

     while(!(u & 1))  u >>= 1;

     long long y,x=pow_mod(a , u , n);

     while(u != n-1){

        u <<= 1;

          y = multi(x , x , n);

          if(y == 1 && x != 1 && x != n-1)             return 1;

        x=y;

    }

    if(y != 1)         return 1;

    return 0;

}

bool miller_rabin(long long n){

     if(n < 2)         return 0;

    long long test[]={2 , 3 , 7};

    for(int i=0 ; i<3 ; ++i){

        if(witness(test[i] , n))             return 0;

    }

     return 1;

}

int main(){

    printf("%lld",pow_mod(7,560,561));

    printf("\n\n\n");

    long long date;

    while(scanf("%lld",&date)){

        if(miller_rabin(date))

          printf("%lldyes\n\n",date);

    }

}

c.还有一种就是打表的方法:

//写得太好了,每一个函数都是一个模版呀~~~~

#include<cstdio>

#include<cstdlib>

#define gcc 10007

#define MAX ((INT)1<<63)-1

typedef unsigned long long INT;

INT p[10]={2,3,5,7,11,13,17,19,23,29};

inline INT gcd(INT a,INT b) {

  INT m=1;

   if(!b)

     return a;

  while(m)     {

    m=a%b;

    a=b;

    b=m;

  }

  return a;

}

//计算a*b%n

inline INT multi_mod(INT a,INT b,INT mod) {

  INT sum=0;

  while(b)     {

    if(b&1)    sum=(sum+a)%mod;

    a=(a+a)%mod;

    b>>=1;

  }

  return sum;

}

//计算a^b%n;

inline INT quickmod(INT a,INT b,INT mod) {

    INT sum=1;

    while(b)     {

       if(b&1)    sum=multi_mod(sum,a,mod);

      a=multi_mod(a,a,mod);

      b>>=1;

     }

    return sum;

}

bool miller_rabin(INT n) {

   int i,j,k=0;

  INT u,m,buf;     //将n分解为m*2^k

  if(n==2)

    return true;

  if(n<2||!(n&1))

     return false;

   m=n-1;

  while(!(m&1))         k++,m>>=1;

   for(i=0;i<9;i++)     {

    if(p[i]>=n)

      return true;

     u=quickmod(p[i],m,n);

    if(u==1)

      continue;

     for(j=0;j<k;j++) {

      buf=multi_mod(u,u,n);

       if(buf==1&&u!=1&&u!=n-1)

         return false;

       u=buf;

    }         //如果p[i]^(n-1)%n!=1那么n为合数

    if(u-1)

     return false;

   }

  return true;

}

//寻找n的一个因子,该因子并不一定是最小的,所以下面要二分查找最小的那个因子

INT pollard(INT n) {

     INT i=1;

     INT x=rand()%(n-1)+1;

     INT y=x;

     INT k=2;

    INT d;

    do     {

       i++;

      d=gcd(n+y-x,n);

      if(d>1&&d<n)

      return d;

       if(i==k)

         y=x,k*=2;

      x=(multi_mod(x,x,n)+n-gcc)%n;

    }while(y!=x);

    return n;

}

INT pollard_min(INT n) {

     INT p,a,b=MAX;

     if(n==1)    return MAX;

    if(miller_rabin(n))    return n;

    p=pollard(n);

    a=pollard_min(p);//二分查找

    INT y=n/p;

    b=pollard_min(y);

    return a<b?a:b;

}

时间: 2024-11-04 07:51:28

《算法导论》——数论的相关文章

NKOJ1236 a^b (数论定理的应用)

          a^b 对于任意两个正整数a,b(0<=a,b<10000)计算a b各位数字的和的各位数字的和的各位数字的和的各位数字的和. Input 输入有多组数据,每组只有一行,包含两个正整数a,b.最后一组a=0,b=0表示输入结束,不需要处理. Output 对于每组输入数据,输出ab各位数字的和的各位数字的和的各位数字的和的各位数字的和. Sample Input 2 3 5 7 0 0 Sample Output 8 5 思路: 数论定理:任何数除以9的余数等于各位数的和除

CodeForces 396A 数论 组合数学

题目:http://codeforces.com/contest/396/problem/A 好久没做数论的东西了,一个获取素数的预处理跟素因子分解写错了,哭瞎了,呵呵, 首先ai最大值为10^9,n为500,最坏的情况 m最大值为500个10^9相乘,肯定不能获取m了,首选每一个ai肯定是m的一个因子,然后能分解就把ai给分解素因子,这样全部的ai都分解了  就能得到m的 所有素因子 以及 所有素因子的个数,题目求的 是n个因子的 不同序列的个数,所以每次 只能选出n个因子,这n个因子由素因子

HDU 4861 Couple doubi(数论)

HDU 4861 Couple doubi 题目链接 题意:给定k,p,有k个球,每个球的值为1^i+2^i+...+(p-1)^i (mod p) (1 <= i <= k),现在两人轮流取球,最后球的值总和大的人赢,问先手是否能赢 思路:先手不可能输,非赢即平,那么只要考虑每种球的值, 利用费马小定理或欧拉定理,很容易得到该函数的循环节为p - 1, 那么i如果为p - 1的倍数,即为循环节的位置,那么每个值都为1,总和为p - 1 如果i不在循环节的位置,任取一个原根g,根据原根的性质,

UVA 10548 - Find the Right Changes(数论)

UVA 10548 - Find the Right Changes 题目链接 题意:给定a,b,c,表示货物的价值,求由A货物和B货物组成C货物有几种方法,判断有无解和是否有无限多种 思路:扩展欧几里得求通解,去计算上限和下限就能判断 代码: #include <cstdio> #include <cstring> #include <cmath> #include <algorithm> using namespace std; const long l

hdu 4542 数论 + 约数个数相关 腾讯编程马拉松复赛

题目:http://acm.hdu.edu.cn/showproblem.php?pid=4542 小明系列故事--未知剩余系 Time Limit: 500/200 MS (Java/Others)    Memory Limit: 65535/32768 K (Java/Others) Total Submission(s): 889    Accepted Submission(s): 207 Problem Description "今有物不知其数,三三数之有二,五五数之有三,七七数之有

[施工中]良心数论.

/* Copyright: xjjppm Author: xjjppm Date: 08-08-17 11:36 Description: Number Theory */ #include <map> #include <cmath> #include <cstdio> #include <cstring> #include <algorithm> inline int input() { char c=getchar();int x=0,a=

信息学中的数论(一)

做oi题目的时候,遇到数论题会令我兴奋不已. 这一篇让我来聊一聊我学过的gcd,lcm,扩展欧几里得算法,逆元,组合数等. 这篇贴的代码都是未经过编译运行的,所以如果有错或有疑问请评论. 恩 那么什么是数论 和数学有关的非几何都是数论? 嘛,我也不知道定义,那么就草率地认为所有和数学有关的非计算几何知识都是数论吧. 我们先来聊一聊gcd. 这个东西,非常的有用. 它的名字叫最大公约数. 正常人都知道,有一个方法叫辗转相除法(证明略): int gcd(int a,int b) { if(!b)r

【数论Day3】进制问题 题解

数论进入第三天,进制问题是常用提醒,是数论的一个重要知识点,常考! 题面:http://www.cnblogs.com/ljc20020730/p/6935255.html 1.K进制数(Kbased.pas/c/cpp) 首先明确数据范围: [数据规模和约定] 对于40%的数据,a的长度不超过5. 对于100%的数据,a的长度不超过100000. 对于40%暴力枚举不多说,上代码: var t,i,k,tt:longint; a:qword; s:string; function pow(x,

Bzoj2219 数论之神

Time Limit: 3 Sec  Memory Limit: 259 MBSubmit: 954  Solved: 268 Description 在ACM_DIY群中,有一位叫做“傻崽”的同学由于在数论方面造诣很高,被称为数轮之神!对于任何数论问题,他都能瞬间秒杀!一天他在群里面问了一个神题: 对于给定的3个非负整数 A,B,K 求出满足 (1) X^A = B(mod 2*K + 1) (2) X 在范围[0, 2K] 内的X的个数!自然数论之神是可以瞬间秒杀此题的,那么你呢? Inpu

【bzoj4872】[Shoi2017]分手是祝愿 数论+期望dp

题目描述 Zeit und Raum trennen dich und mich. 时空将你我分开. B 君在玩一个游戏,这个游戏由 n 个灯和 n 个开关组成,给定这 n 个灯的初始状态,下标为从 1 到 n 的正整数.每个灯有两个状态亮和灭,我们用 1 来表示这个灯是亮的,用 0 表示这个灯是灭的,游戏的目标是使所有灯都灭掉.但是当操作第 i 个开关时,所有编号为 i 的约数(包括 1 和 i)的灯的状态都会被改变,即从亮变成灭,或者是从灭变成亮.B 君发现这个游戏很难,于是想到了这样的一个