基本算法——数学

基本算法是解难题的基础，必须熟练掌握。这一讲将介绍跟数学密切相关的基本算法。

（１）素数

数学类的基本算法大多数属于初等数论范畴，相当大一部分与素数有直接关系，因此素数是一个很基本又很重要的内容。
　　我们先来看看怎么判断一个数是否素数。素数的定义为：如果一个数的正因子只有１和这个数本身，那么这个数就是素数。根据定义，我们立即能得到判断一个数Ｎ（大于１）是否素数的简单的算法：枚举２到Ｎ－１之间的整数，判断是否能整除Ｎ。该算法的 Pascal 代码。

program Prime_1;

var

  x:integer;

  function IsPrime(n:integer):boolean; {返回n是否素数,n>=2}

  var i:integer;

  begin

    for i:=2 to n-1 do                 {枚举因子i}

         if ( n mod i = 0 ) then        {i能整除n吗?}

         begin

           IsPrime:=false;              {n不是素数}

           exit;

         end;

    IsPrime:=true;                             {n是素数}

  end;

  begin

    readln(x);

         if IsPrime(x) then writeln(x,' is a prime!')

                  else writeln(x,' is not a prime!');

  end.

　　如果n很大，那么上面的程序就要运行比较长的一段时间，那么有没有更快一点的算法呢？回答是肯定的！因为如果n含有不为1和自身的因子，那么这些因子中必定有不大于sqrt(n)的（假设n有因子p，1<p<n，如果p<=sqrt(n)，那么p就不大于sqrt(n)，如果p>sqrt(n)，那么n/p也是n的因子，而且1<n/p<n，所以n/p不大于sqrt(n)）。于是我们可以改进上面的程序，得到另外一个 Pascal 程序。容易知道这个算法的时间复杂度为O(sqrt(n))。

（２）因式分解

因式分解的算法很简单，模拟手工分解的过程，我们得到分解n的算法：枚举所有不大于n的所有素数，判断这些素数能整除n多少次。判断2到n是否素数，总共要计算sqrt(2)+sqrt(3)+sqrt(4)…+sqrt(n)<=n*sqrt(n)次，因此算法的时间复杂度可以粗略地认为是O(n*sqrt(n))。事实上，我们有更好的算法。先看一个显而易见的结论：如果p是能整除n的所有大于1的数中最小的，那么p是n的一个素因子。基于这样一个结论，我们得到 Pascal 代码。

 Program YinShiFenJie;

var

  p:array[1..100] of integer;  {素因子}

  s:array[1..100] of integer;  {p[i]对应的指数}

  i,j,n,pnum:integer;

  {n为需要分解的数，n=(p[1]^s[1])*(p[2]^s[2])...*(p[pnum]^s[pnum])}

  procedure YSFJ(x:integer);

  var i:integer;

  begin

    pnum:=0;

    i:=2;

    while ( (x>1) and (i*i<=x) ) do

    begin

      if x mod i = 0 then    {当i是x的因子时}

      begin

        inc(pnum);           {找到一个新的素因子i}

        p[pnum]:=i;

        s[pnum]:=0;

             while ( x mod i = 0 ) do

             begin

               inc(s[pnum]);      {统计i的指数}

               x:=x div i;

             end;

          end;                   {这时，x已经不含小于I的因子了}

      inc(i);

    end;

    if (x>1) then            {表明x是一个素数}

    begin

                  inc(pnum);           {找到一个新的素因子x}

        p[pnum]:=x;

        s[pnum]:=1;

    end;

  end;

  begin

    readln(n);

    YSFJ(n);

    {以下为输出结果}

    write(n,' = ',p[1]);

    dec(s[1]);

    for i:=1 to pnum do

      for j:=1 to s[i] do

        write(' * ',p[i]);

    writeln;

  end.

（３）公因子的数量

问题描述：已知一个正整数N，问这个数有多少正公因子。
算法分析：最容易想到的算法是：枚举1..N，看看有多少个数能整除N，这种算法的复杂度为O( N )。可以优化一下：如果N有小于SQRT( N )的因子X，那么N必定有大于SQRT( N )的因子Y与X对应，而且XY=N。所以我们只需要枚举1..SQRT( N )的数即可，还要考虑N为完全平方数的特殊情况。程序：Pascal。

program factor_number;

  function NumberOfFactor(n:longint):longint;

  var s,i:Integer;

  begin

    s:=0;  {用于统计因子数目}

    for i:=1 to trunc(sqrt(n)) do

      if (n mod i=0) then inc(s);

    inc(s,s);

    if sqr(trunc(sqrt(n)))=n then dec(s);

    {上面这行判断N是否完全平方数}

    NumberOfFactor:=s;

  end;

  var x:longint;

  begin

    readln(x);

    while(x>0)do

    begin

      writeln(NumberOfFactor(x));

      readln(x);

    end;

  end.

上面这个算法的复杂度为O(sqrt(N))。其实我们可以利用因式分解的方法来做。假设我们已经分解N得到 N =(p[1]^s[1])*(p[2]^s[2])...*(p[pnum]^s[pnum])，其中p[i]为互不相同的素数，那么N的正因子的数量为（具体怎么推导请参考组合数学教材中的母函数一章）：(s[1]+1)*(s[2]+1)*…*(s[pnum]+1)。

（４）最大公因式

问题描述：已知两个正整数a和b，求这两个数的最大公因数GCD( a , b )。
(GCD是Greatest Common Divisor的缩写)
算法分析：不妨设a<=b，一种十分容易想到的算法是：枚举1到a的所有整数，在能同时整除a和b的数中取最大的。这个算法的时间复杂度为O(min(a,b))，当min(a,b)较大的时候程序要执行比较长的时间。我们可以利用初等数论中的一个定理：
GCD( a , b ) = GCD( a , b-a ) = GCD( a , b-2*a ) = GCD( a , b-3*a ) = … = GCD( a , b mod a )
关于这个定理的具体证明，请参考初等数论书（或者初中数学竞赛中的数论相关章节）。
下面给出利用这个定理来写的一个求最大公因式的程序，请读者仔细研究：Pascal。

Program GCD_1;

var

  a,b:integer;

  function GCD(a,b:integer):integer;

  begin

    if a=0 then GCD:=b

           else GCD:=GCD(b mod a , a);  {注意这里的顺序}

  end;

  begin

    readln(a,b);

    writeln(gcd(a,b));

  end.

此算法的时间复杂度为O(log(Max(a,b)))。

（５）最小公倍数

问题描述：已知两个正整数a和b，求这两个数的最小公倍数LCM ( a , b )。
(LCM是Least Common Multiply的缩写)
算法分析：直接利用公式：LCM ( a , b ) * GCD( a , b ) = a * b即可。

6）进制转换

我们平常计算都是用十进制数的，但是有的时候我们需要用到２进制数、16进制数等。一个k进制的数可以表示为：(As-1 As-2… A0)k = As-1 K^(s-1) + As-2 K^(s-2) + … + A0 K^(0) ，记为<1>式，其中0<=Ai<K（I=0,1,2..s-1）。对于一个已知的正整数n，如何得到n的K进制表示呢？换句话说，我们就是要求出As-1 As-2… A0来。具体的求解顺序是：先求出A0，然后是A1 ……，最后得到An-1。将<1>式等号两边同时取模k得：n mod K = a0。得到A0以后，(n-A0) div K = As-1 K^(s-2) + As-2 K^(s-3) + … + A1 K^(0)，用与求A0同样的方法可以得到A1，然后是A2……。Pascal 程序。

Program dec2k;  {进制转换，将n转换为k进制数}

const maxs = 100;

const stri:string = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ';

var

  n,k,s,i:integer;

  a:array[0..maxs-1] of integer;

  begin

    readln(n,k);    {从键盘读入n和k}

    s:=0;

    while( n>0 )do

    begin

      a[s] := n mod k;

      n := n div k;

      inc(s);

    end;

    for i:=s-1 downto 0 do

      write( stri[a[i]+1]);

  end.

运行这个程序，输入：155　16 就可以得到结果9B　（16进制的9B = 9*16+11=155）
　　怎么进行任意进制间的数的转换呢？已知一个数正整数n的P进制表示（As-1As-2…A0），求n的Q进制表示（Bl-1Bl-2…B0）。一种简单的方法是：根据P进制表示求出十进制的n，然后再将n转化为Q进制表示即可。
　　前面考虑的都是整数的问题，我们现在来看看怎么处理实有理数。由于实数跟整数的区别仅仅在于小数部分，所以现在只考虑实数r，r满足0 <= r <1的情况。定义r的K进制表示为：r =（0.A1 A2 A3 … As）K = A1 K^(-1) + A2 K^(-2) + A3 K^(-3) …. As K^(-s)。求解顺序为：A1、A2……As。解法：r K = A1 + A2 K^(-1) + A3 K^(-2) …. As K^(-(s-1)) = A1 + B。考察一下B的范围0<= B <= (K-1)K^(-1) +(K-1)K^(-2)… (K-1)K^(-(s-1)) = 1-K^(-s) < 1，也就是说B是 r K 的小数部分，A1是整数部分，于是 A1 = [r K]，[x]表示不大于x的最大整数。由于B = r-A1，所以用同样的方法分解 B 就可以得到 A2、A3……As。Pascal 程序。

Program dec2k2;  {进制转换，将正有理数 r 转换为 K 进制数}

const maxs = 100;

const stri:string = '0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ';

var

  k,s,i:integer;

  r:real;

  a:array[1..maxs] of integer;

  begin

    readln(r,k);    {从键盘读入n和k}

    s:=0;

    while( r>0 )do

    begin

      inc(s);

      r := r * k;

      a[s] := trunc(r);

      r := r - a[s];

    end;

    write('0.');

    for i:=1 to s do

      write( stri[a[i]+1]);

    writeln;

  end.