C语言名题精选百则——数字问题

DSQiu

浏览: 1409845 次
性别:
来自: 广州

最近访客更多访客>>

ymgjava

spike_wq

sdsunrunner

hackxboy

博主相关

博客

微博

相册

留言

关于我

博客专栏

: NGUI所见即所得
浏览量：196922

: Effective C# ...
浏览量：0

文章分类

社区版块

存档分类

博客分类：

C语言名题精选百则
问题解决

C语言名题精选百则——数字问题

尊重他人的劳动，支持原创

这篇博文，D.S.Qiu将对《C语言名题精选百则》第二章进行整理推出，本章一共16个问题，不光只是书上的名题，还会依据互联网的资源进行不断补充，加强（由于匆忙，暂时还只是罗列的过程，不过很快就有深度的发掘）。如果你有建议、批评或补充，请你不吝提出（email:gd.s.qiu@gmail.com，或者直接在本文末评论）。你的支持和鼓励，我将渐行渐远！

这篇博文整理的太艰难了，其中一度想放弃，个人觉得难度有点大，如果没有需求驱动的话，很难读下去，但最终为了保持文章的完整性，还是硬着头皮去截图整理。这16个问题主要经常的是对质数的求解，因式分解，整数自乘（快速求解），Fibonacci数（快速求解），二项式求解（快速求解），阶乘（快速求解），这些都是数论或者组合数学的基本内容，但要通过计算机基本求解的算法赴复杂度过高，文中给出了相应的快速算法，对有需求的读者来说，还是可以喜出望外的（也只是大概的理解了下），应该还有更多精彩的求解方法（日后定当补充）。终于写完了，没有你的支持和鼓励，很难继续下去，下一篇是《排列，组合与集合》。

问题2.1 Armstrong 数（ARMS1.C，ARMS2.C )

在三位的正整数中，例如abc,有一些可以满足a³+b³+c³ = abc的条件，也就是说，各个位数的立方和正好是该数的本身，这些数称为Armstrong数。试编写一个程序来求出

所有的三位Armstrong数。

【说明】

这是个常见的问题，许多课本、习题集都有。解决的基本方法就在于有没有办法从一个三位数中把它的3个位数取出（例如，把379分解成3、7与9),或者在已知3个位数时把该三位数算出来（例如，已知百位为5,十位为1,个位为9，则三位数为519)。

另一个变形是给出一个数，证明它是否为一个Armstromg数。这两者类似，由自己产生的数可以永远控制在三位，但是后者却一定要去测试它是不是三位数。

如何测试一个数是否为三位数呢？有的书中介绍这样一种方法（见程序UGLY.C), 使用对数log()以及指数exp()函数，不仅如此，还要用四舍五入（在C语言中只好加上 0.5再转型成整数）。为什么？原来只是在取出输入数据的某一个位数而已。一个简单的程序写得如此复杂，大多数初学者恐怕都很难看明白。在此给出一个提示，千万不要重蹈覆辙。

(1)如果要确认x是个五位数，就看x/100000是否为零。如果为零，那么x的位数不会超过五位。接着再检查x是否真的是五位，这就要看：

(x % 100000)/10000

是否为零，如果不为零，这个值就是x中从右向左数第五位数的值；如果是零，x就最多只能有四位而已。为什么？因为x%100000是把x用100000去除的余数，值在0〜99999 之间，正是x的右边五位，再用10000去除，那就是第五位的值了；如果x是五位数，那么这个位数就一定大于1,如程序2-1所示。

【程序】

程序2-1（UGLY.C）

#include  <stdio.h>
#include  <stdlib.h>
#include  <string.h>
#include  <math.h>

void main(void)
{
     int    number;           /* the input number         */
     int    pos, ten_pow;     /* # of digits and 10^x     */
     int    arm_no;           /* computed Armstrong Number*/
     int    temp, i;          /* working variables        */
     float  expon;            /* computed exponent        */
     char   line[100];        /* input buffer             */

     printf("\nPlease key in number : ");  /* get a number*/
     gets(line);
     number = atoi(line);

     pos = 0;                 /* count # of digits        */
     ten_pow = 1;             /* power of 10, 10^0=1      */
     do {                     /* start counting           */
          ten_pow *= 10;      /* one more digit...        */
          pos++;
     } while (ten_pow < number);

     if (ten_pow == number)   /* all 10^x are not Arm #   */
          printf("\n%d is not an Armstrong Number.", number);
     else {                   /* not 10^x                 */
          ten_pow = number;   /* now ten_pow is a temp var*/
          arm_no  = 0;        /* initial value for Arm. # */
          for (i = pos-1; i >= 1; i--) { /* scan from R->L*/
               expon = exp(i*log(10.0)); /* the ugly part */
               temp  = ten_pow / expon;  /* you do it...  */
               ten_pow %= (int)(expon + 0.5);
               arm_no  += exp(pos*log(temp)) + 0.5;
          }
          arm_no += exp(pos*log(ten_pow)) + 0.5;
          if (arm_no == number)
               printf("\n%d is an Armstrong Number.", number);
          else
               printf("\n%d is not an Armstrong Number.", number);
     }
}

(2)如果这两个条件都正确，则x是一个五位数。

(3)按照这个技巧，要找出x的第三位（自右往左算)，那就是

(x%1000）/100

(4)用下面的式子也可以求出x的第三位的值：

(x/100)%10

问题实现

/* ------------------------------------------------------ */
/* PROGRAM Armstrong Number Search :                      */
/*    This program searches all 3-digit Armstrong Numbers */
/* and thus numbers in the range of 0-99 are not counted. */
/*                                                        */
/* Copyright Ching-Kuang Shene               June/29/1989 */
/* ------------------------------------------------------ */

#include  <stdio.h>

void main(void)
{
     int  p, q, r;            /* three digits             */
     int  p_cube, q_cube, r_cube; /* their cubes          */
     int  p100, q10;          /* their position values    */
     int  number;             /* the computed number      */
     int  cube_sum;           /* the sum of the cubes     */
     int  count = 0;          /* counter                  */

     printf("\nArmstrong Number Search");
     printf("\n=======================");
     printf("\n\nCount  Number");
     printf(  "\n-----  ------");

     for (p = 1, p100 = 100; p <= 9; p++, p100+=100) {
          p_cube = p*p*p;
          for (q = q10 = 0; q <= 9; q++, q10+=10) {
               q_cube = q*q*q;
               for (r = 0; r <= 9; r++) {
                    r_cube   = r*r*r;
                    number   = p100 + q10 + r;
                    cube_sum = p_cube + q_cube + r_cube;
                    if (number == cube_sum) 
                         printf("\n%3d%9d", ++count, number);
               }
          }
     }
     printf("\n\nThere are %d 3-digit Armstrong Numbers.", count);
}

问题2.2数字谜（TRENTE.C )

在一些游戏与休闲的书刊中，经常会看到如下的数字谜

试编写一个程序进行破解（看下面的规则说明）。

【说明】

数字谜的规则是这样的：

(1)不同的英文字母表示不同的数字。因此，题目中最右边一行的T、Q、E就是不同的数。

(2)每一个数最左边一位不是0。所以上面的谜题中，V、C、T都不是0。

(3)上面的题目是加，那就表示VINGT这个五位数与CINQ这个四位数的两倍相加后，会得到TRENTE这个六位数，并且数字的安排如谜语所示。

在解这些数字谜时，很多人都会不知不觉地用最笨拙的方法，因此可能执行好几分钟, 甚至好几小时都得不到结果。就以上题为例，它一共有9个不同的字母，因此耗用10个数字符号中的9个，很多人就会写出下面的形式，如程序2-2所示。

【程序】

程序2-2

for(v=1;v<=9;v++)

for(I=0;I<=9;I++)

……

for(E=0;E<=9;E++)

{

VINGT=……;

CINQ=……;

TRENTE=……;

if(TRENTE == VINGT + CINQ +CINQ)

{

……

}

如果用这个程序，那么也许一两个小时也不一定能够算出结果，因为在最内层循环中的if 一共要执行将近10^9,近一亿次，所以这是最笨的方法。通常，用一点小学算术知识就可以节省大量的计算时间。换言之，可以用四则运算来找出某个字母值的范围（如偶数、奇数、5〜9、3或4等，而不是0〜9),范围越小，程序的工作就越少，执行起来就越快。

不妨从这一点出发，就本题而给出一些提示：

(1)因为3个位数相加，最大就是9+9+9=27,如果从这一个位数的右边一位会有进位，那么这个进位就最多是2，所以3个位数再加上进位的最大值是29，所以该位数和是 9，而进位是2。

(2)现在看V那位，加上进位之后得到TR,因为V最大是9，进位最大是2，所以 TR的最大值是11。因为T是结果的第一位，所以不能是0，因此只好是1;但结果前两位是TR, TR最大是11，现在T是1,于是R只能是0。注意，如果R是2〜9中任意一个，那就表示没有从R到T的进位，因而T就是0 了。所以从一点看来，两个值就固定了，即T=1， R=0。

(3)再来看V，V—定大于8。为什么？千位数的进位最大是2,要一个数加上2而有两位数，那么这个数不就至少是8吗？

(4)看个位数的E。它是由T加上两个Q得来的，已知T是1,而2Q—定是偶数, 所以E是个奇数。

就用这4点提示，T与R的循环就不必了； E的循环有3、5、7、9 (注意，T是1, E 就不能用1); V的循环只能用8与9。还没有决定的就剩下I、N、G、C与Q5个，所以循环就有7层，I、N、G与Q为0〜9; C为2〜9 (因为C是第一位，不能是0，也不能是1，因为T是1)，E是3、5、7、9； V是8与9，所以最内层的if就执行了 10^4x7x4x2次而已，于是用这个方法写成的程序所花的时间大约是前面所提到的：

10^4x7x4x2/10^9=0.056%

因此可以看到节省了多少时间。

思考一下，另外是否还有其他可以节省时间的方法呢?

问题实现

/* ------------------------------------------------------ */
/* PROGRAM Number Puzzle :                                */
/*    This is first solution of the number puzzle.  It is */
/* written under the guidelines in my book.  For each     */
/* digit, there corresponds a for loop and an if statement*/
/* to prevent duplicated digits.  Read it carefully.      */
/*                                                        */
/* Copyright Ching-Kuang Shene               June/30/1989 */
/* ------------------------------------------------------ */

#include <stdio.h>

void main(void)
{
     long  I, N, R=0, E, V, G, C, Q, T=1;
     long  IN, VINGT, CINQ, TRENTE;
     long  sum;

     printf("\nNumber Puzzle\n");
     printf("\n   VINGT");
     printf("\n    CINQ");
     printf("\n+)  CINQ");
     printf("\n--------");
     printf("\n  TRENTE\n");

     for (V=8; V<=9; V++)
       for (I=0; I<=9; I++)
         if (I != V && I!=T)
           for (N=0; N<=9; N++)
             if (N!=I && N!=V && N!=T) {
               IN = I*10 + N;
               for (G=0; G<=9; G++)
                 if (G!=N && G!=I && G!=V && G!=T && G!=R) 
                    for (C=2; C<=9; C++)
                      if (C!=G && C!=N && C!=I && C!= V && C!=T && C!=R)
                        for (Q=2; Q<=9; Q++)
                          if (Q!=C && Q!=G && Q!=N && Q!=I && Q!=V 
                                   && Q!=T && Q!=R)
                             for (E=3; E<=9; E+=2)
                               if (E!=Q && E!=C && E!=G && E!=N && E!=I
                                        && E!=V && E!=T && E!=R) {
                                 TRENTE=((((T*10+R)*10+E)*10+N)*10+T)*10+E;
                                 VINGT=((V*100+IN)*10+G)*10+T;
                                 CINQ =(C*100+IN)*10+Q;
                                 sum=VINGT+CINQ+CINQ;
                                 if (sum==TRENTE) {
                                   printf("\n\nThe answer is :\n");
                                   printf("\n%8ld", VINGT);
                                   printf("\n%8ld", CINQ);
                                   printf("\n+)%6ld", CINQ);
                                   printf("\n--------");
                                   printf("\n%8ld", TRENTE);
                                 }
                               }
             }
}

问题2.3求质数（PRIME1.C )

试编写一个程序，找出前N (如200)个质数。如果没有进一步要求，这不是难题。但在此希望从所知的、使用除法的方法中，用最快的办法来编写程序。

【说明】

可能最先想到的办法，就是让某个变量i从2变到N,然后检查它是不是质数，如果是就显示出来，如果不是，就检查下一个。这是正确的做法，但却没有注意到一个小细节, 因而使程序运行速度变慢。当然，2是质数，但所有2的倍数都不是质数，如果令i从2 变到N,不就很冤枉地测试了 4，6,8,10，…这些数？所以第一点提示是测试2，3，5，7，…，N，即 2与所有奇数就足够了。同理，3是质数，但6,9,12，…这些3的倍数却不是，因此，如果能够把2与3的倍数都跳过去而不测试，任意连续的6个数中，就只会测试两个而已。以6n，6n+1，6n+2，6n+3，6n+4，6n+5 为例，6n,6n+2，6n+4 是偶数，又 6n+3 是 3 的倍数，所以如果把2与3的倍数都不理会，要测试的数就只留下6n+1与6n+5而已，因而工作量只是前述想法的2/6=1/3,应该用这个方法去编写程序。

假设i是如上述的一个数（不是2与3的倍数)，如何测试i是个质数呢？按照定义，i 如果是质数，也就只有1与i可以除尽自己，所以可以用2,3,4,5,6，…，i-1去除i,如果都除不尽（余数不是0), i就是质数。观点也对，但却与上一点一样地笨拙。当i>2时，有哪 1个数i可以被i-1除尽的？没有，为什么？如果i不是质数，那么i=a*b,此地a与b既不是i又不是1;正因为a>1，a至少是2，因此b最多是i/2而已，去除i的数用不着是 2，3，4，…i-1；而用2，3，4,…,i/2就可以了。不但如此，因为i=a*b，a与b不能大于i^½，为什么呢？如果i不是质数,它的因子最大就是换言之，用2，3，4，5，…，i^½去除i就行了。

但是，用2,3，4，5，…，i^½去除 i 也是个浪费。就像是前一段所说的，2除不尽，2的倍数也除不尽；同理，3除不尽，3的倍数也除不尽……最理想的方法就是用质数去除i，因为在前一段的提示，i不是2与3的倍数，所以就用5，7，11，13，17，19，…这些比 i^½ 小的质数去除i即可。但问题是这些质数从何而来？这比较简单，可以准备一个数组prime[]，用来存放找出来的质数，一开始它应该有2、3与5。于是当i=5,7,11，13，17,19,23,25,29，…这些不是2与3 的倍数时，就用prime[]中小于 i^½ 的数去除 i 即可，如果都除不尽，i就是质数，把它放入 prime[]中，因此prime[]中的质数雜会越来越多，直到满足所要的个数为止。

不妨用这段说明来编写程序，不过应注意下面几点：

(1)不要处理2与3的倍数（见第一段)。

(2)用质数去除（见第二、三段）。

(3)用i^½可能会有问题，为什么？有什么方法可以解决？

【解答】

程序的结构本身并不复杂，只解释几个重点而已。

第一，看如何跳过2与3的倍数。在问题说明中，看过在6n，6n + 1，6n+2，611+3，6n+4 与6n+5中，只有6n+1与6n+5要做测试，看看是否是质数，看看这两个数的相对位置。上一组最后一个要测试的值，与这一组第一个要测试的值的距离是2; 而同一组中要测试的两个值之间的距离是4。所以从5起，依次加2、加4、加2、加4、… 就可以得到5，7，11，13，17，19，…这些要测试的值了。但是2+4=6,可以用一个变量（程序中是 gap),令它的值在开始时为2,下一个值就是6-gap (就是4),再下次6-gap不就是2了吗？这就是程序中gap=6-gap的意义。

第二，不能比较是否小于i^½，原因是对某些i而言，i^½可能并不精确。举个例子，如果i=121，在数学上i^½自然是11,但是计算机算出来的结果很可能是10.999999，于是去除i的数就是2，3，5，7而不含11，因此变成121是个质数了。解决之道很简单，不要用开方，用平方即可。例如，如果p*q<i，就用p去除，这就一定不会出现上面的问题。不但如此，平方的运算时间也一定会比开方有效率。

【程序说明】

在程序中，prime[]是用来存放找出来的质数的，一开始时就有了 2，3，5这3个最基本的质数（都小于6)。gap的初值是2，它的值依次是2，4，2，4，…，在解答中见过了。Count 计算到目前为止一共找到了多少个质数，2,3，5已经放入prime[],所以count的值是3。于是当count比MAXSIZE小的时候，就产生下一个要检查的数，存入may一be一prime,进入 for循环做检查工作；如果是个质数，那么for循环中的if就永远是假，is一prime的值就不曾被改变，因此就把它存入prime[]中，这就是程序的基本结构。

【问题实现】

/* ------------------------------------------------------ */
/* PROGRAM prime number generator :                       */
/*    This program generates first MAXSIZE-1 primes by    */
/* using the conventional division method, but multiples  */
/* of 2 and 3 are not tested.  Therefore it is faster.    */
/*                                                        */
/* Copyright Ching-Kuang Shene               July/02/1989 */
/* ------------------------------------------------------ */

#include  <stdio.h>

#define   MAXSIZE    100
#define   YES          1
#define   NO           0

void main(void)
{
     int  prime[MAXSIZE];     /* array to contains primes */
     int  gap = 2;            /* increasing gap = 2 and 4 */
     int  count = 3;          /* no. of primes            */
     int  may_be_prime = 5;   /* working  variable        */
     int  i, is_prime;

     prime[0] = 2;            /* Note that 2, 3 and 5 are */
     prime[1] = 3;            /* primes.                  */
     prime[2] = 5;
     while (count < MAXSIZE)  { /* loop until prime[] full*/
          may_be_prime += gap;  /* get next number        */
          gap           = 6 - gap; /* switch to next gap  */
          is_prime      = YES;  /* suppose it were a prime*/
          for (i = 2; prime[i]*prime[i] <= may_be_prime && is_prime; i++)
               if (may_be_prime % prime[i] == 0) /* NO    */
                    is_prime = NO; /* exit                */
          if (is_prime)       /* if it IS a prime...      */
               prime[count++] = may_be_prime;  /* save it */
     }

     printf("\nPrime Number Generation Program");
     printf("\n===============================\n");
     printf("\nFirst %d Prime Numbers are :\n", count);
     for (i = 0; i < count; i++) {
          if (i % 10 == 0) printf("\n");
          printf("%5d", prime[i]);
     }
}

【习题】

(1)笔者认为，使程序中跳过2,3，5的倍数也是可能的，请构想一个可以达成目标的做法。如果想出来了，请问它有多复杂，值不值得写到程序中？如果用这个办法，程序会快多少（当然会快的）？

(2)程序中把5放到prime[]$,请问这是不是多此一举？为什么？如果不这样做，有什么好建议？

(3)有人说，is_prime这个变量是多余的，其实也是如此，有没有办法修改这个程序，使它不用is_prime呢？程序会比较容易读吗？

(4)笔者曾在一本书中见到一个程序（BAD.C),它是用来求1〜100之间的质数用的，但该书的作者要求“尽可能地使循环的执行次数为最少”。请问，这位作者的要求达到没有？为什么？请列出所有的改良建议。

(5)请修改程序，让它可以求出1〜N (输入）的所有质数。请问，如果程序又改成读入M与N (M<N=,求出M〜N之间的所有质数时，程序的工作量可以少一些吗？为什么？

问题2.4筛法（SIEVE.C )

求质数是一个很普遍的问题，通常不外乎用数去除，到除不尽时，给定的数就是质数。但是早在2000年前人们就知道了一个不必用除法而出2〜N的质数的所有方法了。假设有一个很神奇的筛子，可以给出一个数，例如i,这个筛子有办法把i的所有倍数去掉。请用这个方法求出2〜N之间的所有质数。

这个方法称为Eratosthenes (人名）筛法（Sieve Mothod)。

【说明】

下面通过一个例子来加强对筛法的印象，求出2〜40之间的质数。首先，把2〜40这些数一字排开：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

2是质数，所以把2的倍数都筛掉:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

下一个数自然是质数3,所以把3的倍数都筛掉 •

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

下一个是5,把5的倍数筛掉:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

下一个是7,把7的倍数筛掉（其实在此之前都已经筛掉了）: f

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

再下来是11,比20/2大了，所以工作停止，没有被筛掉的就是质数，它们是2，7，11，13，17，19，23，29，31，37。

可以按照这一逻辑来编写程序，但是需注意下面几点：

(1)2是质数，所以2的倍数是一定会被删除的，所以在一开始时根本没有把2的倍数放到筛子中的必要，这就省下了一半的空间。

(2)如果要求2〜N之间的质数，做到N/2就可以停下来，因为大过N/2的数都不可能整除N。

(3)程序不可以使用乘法与除法，只能用加或减，以求加快速度。

请基于这3项要求来编制程序。

【解答】

用一个数组x[0]，x[1]，…来储存3，5，7，11，…这些奇数，因此x[i]中所储存的就是2i+3,或者是写成i+i+3。在开始的时候，把每一个元素都定成“没有筛掉”，然后再逐次把合数去掉。

如何把合数筛掉呢？办法其实很容易，把x[0]，x[1]，…一个接一个地查过去，如果x[i]没有被筛掉，那它一定是个质数（为什么？），但x[i]对应的是2i+3,因此把2i+3的所有倍数筛掉。2i+3的倍数在\[]数组中的什么地方呢？首先，知道2i+3的倍数可以写成(2i+3)j，此处j>l (如果j=l,就把自己筛掉了）；不过，因为2i+3是个奇数，而且x[]中所有的元素又是奇数，所以要筛掉的是2i+3的奇数倍数，换言之，即3(2i+3)，5(2i+3),7(2i+3)，…；这些奇数又是2n+l的形式，因此要删掉的其实是(2n+l)(2i+3)这些数，n=l，2，3,…。化简这个式子，得到：

(2n + 1)(2i + 3) = 2n(2i + 3) + 2i + 3= 2[n(2i + 3) + i] + 3

所以当知道x[i]所对应的数(2i+3)是个质数之后，所要筛掉的数在x[]中对应的位置是n(2i+3)+i，n=1，2，3，…。当 n=1时，这个式子是(2i+3)+i，n=2 时就是(2i+3)+i 加上(2i+3), n=3 时再加上2i+3;所以可以用(2i+3)+i作初值，筛掉x[]中那一个数，然后就每次加上2i+3, 筛掉对应的数等，一直到无数可筛为止。如果数组有N个元素，那么最后一个元素x[N-1] 就对应着2 • N+3,因此就可以求出2到2 • N+3之间的所有质数了。

在程序中sieve[]就对应着x[], MAXSIZE对应着n。一开始每一个元素存入KEPT，表示没有筛掉。

中间的大for循环就把sieve[]的元素一个接一个查过去，一旦查到sieve[i]是没被筛掉的，那么i+i+3就是个质数，把它存放在变量prime中，count是用来记录到目前为止找到了多少个质数，在开始时，2是质数，所以count为1。接下来的while循环中，k就是用来指出要筛掉prime的倍数，前面提过，k从prime+i起，每一次增加prime这个值，并且把sieve[k]筛掉。到sieVe[]中元素都处理完后，值为KEPT的就是没有被筛掉的元素；程序中剩余的部分把结果显示出来。

问题实现

/* ------------------------------------------------------ */
/* FUNCTION sieve :                                       */
/*    This program uses Eratosthenes Sieve method to      */
/* generate all primes between 2 and MAXSIZE*2+3.  This   */
/* is a very simple yet elegant method to generate primes */
/*                                                        */
/* Copyright Ching-Kuang Shene               July/01/1989 */
/* ------------------------------------------------------ */

#include  <stdio.h>

#define   MAXSIZE   200
#define   DELETED     1
#define   KEPT        0

void main(void)
{
     int  sieve[MAXSIZE+1];   /* the sieve array          */
     int  count = 1;          /* no. of primes counter    */
     int  prime;              /* a generated prime        */
     int  i, k;               /* working variable         */

     printf("\nEratosthenes Sieve Method");
     printf("\n=========================");
     printf("\n\nPrime Numbers between 2 and %d\n", MAXSIZE*2+3);

     for (i = 0; i <= MAXSIZE; i++) /* set all items to be*/
          sieve[i] = KEPT;    /* kept in the sieve        */

     for (i = 0; i <= MAXSIZE; i++) /* for each i, it     */
          if (sieve[i] == KEPT) {   /* corresponds to 2i+3*/
               prime = i + i + 3;   /* if it is not sieved*/
               count++;             /* prime=2i+3.        */
               for (k = prime + i; k <= MAXSIZE; k += prime)
                    sieve[k] = DELETED; /* screen multiple*/
          }

     printf("\n%6d", 2);      /* output part.             */
     for (i = 0, k = 2; i <= MAXSIZE; i++)
          if (sieve[i] == KEPT) {
               if (k > 10) {
                    printf("\n");
                    k = 1;
               }
               printf("%6d", 2*i+3);
               k++;
          }
     printf("\n\nThere are %d primes in total.", count);
}

【习题】

(1)为什么在程序运行过程中，若sieve[i]为KEPT, 2i+3就是一个质数呢？

(2)请分析一下程序一共查过了多少个元素。

提示：如果sieve[]有n个元素，在做到第i个时，要筛掉2i+3的倍数，这有多少个？当然这与i有关，当把所有与i对应的量都加起来之后，结果就做出来了。

(3)不处理2的倍数可以少处理一半的元素；因此有人建议，连3的倍数也不处理了，速度不就更快了吗？ 2与3的倍数是6个一组的：6n，6n+l,6n+2，6n+3,6n+4，6n+5; 6n,6n+2,6n+4是2的倍数，6n+3是3的倍数，所以6个数中可能会是质数的就只剩下6n+1， 6n+5而已，原来的2/6 = 1/3。换句话说，处理1/3的元素就行了。请把这个观点写成程序。

问题2.5线牲筛法（L_SIEVE.C )

在做筛法（Sieve)求质数的问题时，非质数的数据有很多是重复删除的。例如，如果有一个数是3x7x17x23,那么在删除3的倍数时会删到它，删除7、17与23的倍数时也都会去删它。看起来好像并没有多大关系，因为只是把一个值存入对应的位置而已，但是, 当要求质数的范围很大时，此类有很多质因子的合数越来越多，因此程序的效率自然就会大打折扣。基于筛法的观念，编写一个程序，使得在删除非质数时“绝对”不做重复的工作。

【说明】

解这个问题的诀窍是如何安排删除的次序，使得每一个非质数都只被删除一次。

中学时学过一个因式分解定理，它说任何一个非质（合）数都可以分解成质数的连乘积。例如，16 = 2^4，18 = 2 x 3^2 , 691488 = 2^5 x 3^2 x 7^4等。如果把因式分解中最小的质数写在最左边，有16 = 2^4, 18 = 2x9, 691488= 2^5x21609；换句话说，把合数n写成n = p^k*q，此时q当然是大于p的，因为p是因式分解中最小的质数。由于因式分解的惟一性，任何

一个合数n，写成n = p^k*q的方式也是惟一的。

由于q≥p的关系，因此在删除非质数时，如果已知P是质数，可以先删除P^2 ,P^3，P^4，…，接着找出比P大而没有被删除的数q，然后删除pq,p^2q,p^3q,…，一直到pq>n为止。试把这个概念编制成程序。

因为每一个非质数都只被删除一次，可想而知，这个程序的速度一定相当快。依据Gries与Misra的文章，线性的时间，也就是与n成正比的时间就足够了（此时要找出 2n的质数）。

【解答】

删除的技巧并不复杂，从2开始，先删除2^2，2^3，2^4，…，接着删除2·3,2^2·3，2^3·3，…(注意，3并没有被删除)，再删除…… 一般而言，当发现p是一个质数时，就先删除p^2，p^3，p^4，…这一串合数。正因为 p是个质数，p^2，p^3，p^4，…不可能被其他数整除，所以在此之前不可能被删除。接着，去找出比p大，但没有被删除的一个数，令为q;因为P 是质数，q是一个在此之前没有删除的数，因此p^i*q(i = 1，2，"*)在此之前也没有被删除过。为什么会这样呢？简单地说，如果p^i·q在此之前已经被删除了，那么依照上面所讲的删除方式，一定有一个质数a (a<p=存在使得a^kb正好就是pi_q,这样在删除a的倍数时就会把a^kb删掉的。不过，a^k*q=p^i*q, a<p,这是不可能的，因为a与p都是质数，彼此是互质 (不能整除对方）的，所以p^i就一定整除b才行，于是q = a^k*(b/p^i), q是a^k的倍数，在删除a^k的倍数时q早就被删除了才对，这与前面所说“q是比p大的第一个没有被删除的数” 相矛盾。因此p^i*q (i = l，2，…）都是第一次被删除，换句话说，删除的动作绝不重复。

问题实现

/* ------------------------------------------------------ */
/* PROGRAM linear sieve program of Gries and Misra :      */
/*    This program finds all prime numbers between 2 and  */
/* n, the input, by using Gries and Misra linear sieve    */
/* method.  This method does not use division, instead    */
/* multiplications are used.                              */
/*                                                        */
/* Copyright Ching-Kuang Shene               July/10/1989 */
/* ------------------------------------------------------ */

#include  <stdio.h>
#include  <stdlib.h>

#define   MAXSIZE    1000
#define   NEXT(x)    x = next[x]
#define   REMOVE(x)  { next[previous[x]] = next[x];          \
                       previous[next[x]] = previous[x];      \
                     }
#define   INITIAL(n) { unsigned long i;                      \
                       for (i = 2; i <= n; i++)              \
                            previous[i] = i-1, next[i] = i+1;\
                       previous[2] = next[n] = NULL;         \
                     }

void main(void)
{
     unsigned long  previous[MAXSIZE+1]; /* prev. pointer */
     unsigned long  next[MAXSIZE+1];     /* next pointer  */
     unsigned long  prime, fact, i, mult;
     unsigned long  n;
     unsigned long  count = 0;
     char           line[100], dummy;

     printf("\nLinear Sieve Program");
     printf("\n====================");
     printf("\n\nPrime Numbers between 2 to --> ");
     gets(line);
     n = strtoul(line, &dummy, 10);
  
     INITIAL(n);              /* initial the set          */
     for (prime = 2; prime*prime <= n; NEXT(prime))
          for (fact = prime; prime*fact <= n; NEXT(fact))
               for (mult = prime*fact; mult <= n; mult *= prime)
                    REMOVE(mult); /* remove multiples     */

     for (i = 2; i != NULL; NEXT(i)) { /* display result  */
          if (count % 8 == 0)  printf("\n");
          printf("%6ld", i);
          count++;
     }
     printf("\n\nThere are %ld Prime Numbers in Total", count);
}

【习题】

(1)这个程序花了一半的时间去删除2〜n之间的偶数，这是多余的，因为除了 2之外，所有偶数都不是质数。请修改此程序，让它不处理2之外的偶数，这就省下一半的时间（2〜n之间有一半元素是偶数)。提示：参考上一题的技巧。

(2)此处所用的内存量很明显地比上一题的传统筛法程序大许多，把它省下一半。首先不要previoiis[],只留下next[]。对于next[i]而言，如果i-1并没有被删除，i的上一个是i-1;但若i-1己经被删除了，用next[i-1]来储存原来的previous[i]。换句话说，i的上一个就是i-1与next[i-1]中值较小的那一个。请问，还有什么边界条件(Boundary Condition)要考虑，并且把这个概念写成程序。这个技巧在习题（1)中能不能适用呢？

(3)在这一题看看如何修改本问题而同时求出2〜n之间每一个数的因式分解，这是 Dexter Kozen提出来的方法。因为每一个合数r都可以写成r=p^i·q的形式，p是r的因式分解中最小的质数。当用本题求2〜n的质数时，一旦决定要删除r 了，就有p^i*q,可以修改程序算出p，i与q,再准备3个数组，即P[]、EMP[]与Q[],在P[r]、EXP[r]与Q[r]中分别储存P、i与q;储存的动作只会做一次。如果在程序执行之前，在？[]中存入一个奇怪的值（比如0),于是在程序执行完后，如果P[t]的值并未被改变，表示t不能写成p、q的形式，亦即不能分解，所以t是质数。如果原来P[t]的值已经被改变了，就可以去找Q[t],再用Q[t]作脚码去查P[Q[t]]而得到Q[t]的分解方式，因此一路追踪，就可以把t分解完成。

请把这个概念写成程序。

问题2.6因子分解(FACTOR.C )

编写一个程序，读入一个正整数，把它的所有质因子找出来。例如，如果输入的是72, 72 = 2^3x3^2,于是质因子就有2与3，如果输入的是181944, 181944=2^3x3^2x7x19^2,所以质因子为2、3、7与19。为了方便起见，2^3x3^2x7x19^2可以用2(3)3(2)7(1)19(2)作为输出形式，也就是说，如果分解开来有ab，输出时就是a(b)。

【说明】

传统的做法是把输入值（假设是n)用2去除，一直到除不尽为止。如果一共除了 i 次，就有2^i这一项，输出中就会出现2(i);接着再用3去除、5去除、7去除等，直到商数变成1为止。以181944为例，第一次用2除得到93972,再除一次是46896,第三次得到 23493,于是2就不能整除了。下来用3去除，第一次得到7831,第二次是2527，第三次就不能整除。对于2527而言，用7去除得到361，再用7就除不尽了，其次的11、13、15、 17,也都除不尽；但19可以，再用19去除得19;最后用19除，商为1，不能再除了，因此就得到181944 = 2^3x3^2x7x19^2的结果。

试用这个概念来编写程序。

【解答】

用一个数不断地去除另一个数的动作可以用for循环很漂亮地写出来，就是说要用k 去除work，一直到除不尽或被除数变成1为止。可以写成：

for(i=0;work%k==0 && work > 1; work /= k, i++)

变量i是用来计算整除了多少次的，work是要被除的数，至于k则是固定的除数，自然 i的值在原始时为0。每一次循环时，都用k去除work,看余数是不是0，而且被除数是否大于1;如果是，就进入循环主体，但是因为循环主体空无一物，这就相当于做到增加值的部分，把work被k除的商求出来，当作下一次除法的被除数，并且因为已经整除了一次，i的值也加上1。这就是用一个数反复地除另一个数，直到除不尽时的做法。了解这一点之后，令k=2，3，5，7，9，11,…这样反复地执行上面的循环，并且把k与i (前者是因子，i是指数）显示出来，一直到k大过work为止即可。注意，正整数被2除可以用右移1位完成，而余数部分正是最右边的那一位的值，这是下面FACT0R.C程序中第一个for循环的意义。

问题实现

/* ------------------------------------------------------ */
/* PROGRAM factorization of positive integers :           */
/*    Given un unsigned long integer, this program finds  */
/* all prime factor by using traditional division method. */
/*                                                        */
/* Copyright Ching-Kuang Shene               July/10/1989 */
/* ------------------------------------------------------ */

#include  <stdio.h>
#include  <stdlib.h>

#define   MAXSIZE  100
#define   SAVE_FACTOR(fact, exp) { if (exp > 0)                \
                                        factors[count] = fact, \
                                        exps[count++]  = i;    \
                                 }

void main(void)
{
     unsigned long  factors[MAXSIZE]; /* stores factors   */
     unsigned long  exps[MAXSIZE];    /* stores exps      */
     unsigned long  n, work;
     int            count = 0;        /* factor counter   */
     int            i, k;
     char           line[100], *dummy;

     printf("\nFactorization by Division Program");
     printf("\n=================================");
     printf("\n\nInput a positive integer --> ");
     gets(line);
     n = strtoul(line, &dummy, 10);

     for (i=0,work=n; (work & 0x01UL)==0 && work>1; work>>=1,i++)
          ;                   /* extract divisor 2        */
     SAVE_FACTOR(2, i);       /* save it and its exp.     */

     for (k = 3; k <= work; k += 2) { /* for k=3,5,7,9,.. */
          for (i = 0; work % k == 0 && work > 1; work /= k, i++)
               ;              /* extract divisor k        */
          SAVE_FACTOR(k, i);  /* save it and its exp.     */
     }

     printf("\n%ld = ", n);   /* display result.          */
     for (i = 0; i < count; i++)
          printf("%ld(%ld)", factors[i], exps[i]);
}

【习题】

(1)程序中都有work>l的测试，究竟这有没有必要？请说明理由。

(2)在第2个for循环中，分别用k=3，5，7，9，11，13，15，…去除，但事实上在用3去除时就已经去掉了所有3的倍数，所以9与15是不可能整除work的，做了多余的事。请修改这个程序，使得只会用质数去除，9，15，…等就不会出现了。

(3)请问，这个程序所输出的因子为什么都是质数？试证明它。

(4)在程序中，其实factors []与exps[]这两个用来存放因子与对应的指数的数组是不必要的，为什么？请修改程序，不用这两个数组，但效果完全相同。

问题2.7 数值自乘递归解（R_POWER.C ）

如果n与m是正整数，那么m^n就是把m连乘n次，这是一个效率很低的方法，请写一个计算效率高的程序，并且分析程序中一共用了多少个乘法，应该以n-1个乘法作为设计准则。

【说明】

这是一个典型的递归设计题目，应该注意以下几点：

(1)试用分而治之（Divide.and.Conquere )的策略。

(2)注意到x^4可以用x^2自乘的关系，由此可以大量地降低乘法数目。

(3)连乘n次要n-1个乘法，能做到只要2*

个乘法吗？

问题实现

unsigned long  recursive_power(unsigned long m, unsigned long n)
{
     unsigned long temp;

     if (n == 0)              /* m^0 = 1                  */
          return 1;
     else if (n & 0x01UL == 0) { /* if power is even then */
          temp = recursive_power(m, n >> 1); 
          return temp * temp; /* m^(2k) = (m^k)^2         */
     }
     else                     /* otherwise, m^n=m*m^(n-1) */
          return m * recursive_power(m, n-1);
}

问题2.8数值自乘非递归解（I_POWER.C)

继续求m^n问题（m与n是正整数)。前面的提示会得到一个递归的程序，请编制一个运算效率同样高的非递归的程序。

【说明】

或许读者有自己独特的看法，但在此提供一个简单的建议，可以釆用它来编写程序，当然也可以把它化简。建议是把指数n用二进制来看，比如若n=13，那么13(10)=1101(2)= 2^3 + 2^2 +2^0，1101(2)表示二进制数，所以求m^n时就相当于求m^(2^3 + 2^2 +2^0)=m^(2^3)·m^(2^2)·m^(2^0)；（输入不了公式抱歉），会发现二进制表示中对应那一位是1,在m^n中就有那一项。把这个观念编制成程序。

另外一个办法是可以把递归解法中每一个递归步骤的n提出来，看在什么时候用 (m^k)^2,什么时候用m(m^2k),然后用非递归方式写出来。了解了这些观点之后，编写这个程序就不难了。在编写完程序之后，还应该分析一下程序乘了多少次。

【问题实现】

unsigned long iterative_power(unsigned long m, unsigned long n)
{
     unsigned long  temp = 1;

     while (n > 0) {          /* if there exists 1 bits.. */
          if (n & 0x01UL == 1)/* the right most one ?     */
               temp *= m;     /* YES, the result get a 'm'*/
          m *= m;             /* anyway, compute m^k      */
          n >>= 1;            /* throw away this bit      */
     }
     return temp;
}

【习题】

(1)这个程序的乘法数目还可以再降低，其实保持完整的m，m^2，m^4，m^8，…是没多大必要的，请以此为出发点来改良程序。

(2)请用手算方式列出计算2^13，2^15，2^16的过程，是不是这个程序也会得到算2^15时要多几个乘法呢？请问，一般而言，当n是多少时会比较慢？

(3)如果用的硬件系统会检测出整数溢位（Overflow)，那么这个程序可能会有问题，因为求出m^45之后，程序还会多算m^64，这个m^64就可能溢位了。请修改这个程序，克服这

个问题。

问题2.9 Fibonacci数悲递归解(FIB_IT.C )

Fibonacci 数列 {fn} 的定义是:

fn=fn-1 +fn-2 n>2; fn=1, n=1,2;

不用递归的方法，也不用数组，编写一个函数，它接收n的值，返回fn。

【说明】

用递归方法算Fibonacci数列效率是很低的，要计算很多个重复的加法，这个题目要求不用递归，不用数组，把它求出来。不过应注意下面的事项：

(1)递归方式并非全然不好，但不能直接套用公式。

(2)因为当n>2时，fn = fn-1+fn-2 ,所以程序只保留fn-1与fn-2就可以算出fn。

传统的递归方法

unsigned long fibonacci(int n)

{

if (n<2)

return 1;

else

return fibonacci(n-1)+fibonacci(n-2);

}

递归方法多了很多重复的计算，每次算fibonacci(n-1)时都顺带计算了fibonacci(n-2)，但是后面算fibonacci(n-2)又给算了一次，为了避免这个弊端，用两个变量保存每次计算的fn-1和fn-2。

【问题实现】

unsigned long  fibonacci(int n)
{
     unsigned long  f0, f1, temp;
     int  i;

     if (n <= 2)
          return 1;
     else {
          for (f0 = f1 = 1, i = 3; i <= n; i++) {
               temp = f0 + f1; /* temp = f(n-2)+f(n-1)    */
               f0   = f1;      /* f0   = f(n-2)           */
               f1   = temp;    /* f1   = f(n-1)           */
          }
          return f1;
     }
}

问题 2.10 快速Fibonacci 数算法（FIB_MT.C )

继续讨论Fibonacci数列问题。在非递归的Fibonacci程序中，在算fn时最多不超过n-2 个加法，编写一个速度更快的程序，或许可以用乘法。如果每一个乘法用m单位时间，每一个加法用P单位时间，于是非递归的写法因为最多有n-2个加法，因此最多用(n-2)p时间。请问，改善的程序要用多少时间？假设只考虑加法与乘法而已。

【说明】

解决这个问题的技巧有不少，在此先提示一个很容易理解的方法。用矩阵来算，看下面的式子:

可以用这个观点来编写程序。

问题实现

C代码 




 
void matrix_power(unsigned long, unsigned long,   
                  unsigned long, unsigned long, int,  
                  unsigned long *, unsigned long *,   
                  unsigned long *, unsigned long *);  
  
  
unsigned long  fibonacci(int  n)  
{  
     unsigned long  a, b, c, d;  
  
     if (n <= 2)  
          return 1;  
     else {                   /* A^(n-2) of a matrix A    */  
          matrix_power(1UL, 1UL, 1UL, 0UL, n-2, &a, &b, &c, &d);  
          return a + b;  
     }  
}  
  
  
/* ------------------------------------------------------ */  
/* FUNCTION matrix_power :                                */  
/*    This function accepts a matrix X and computes       */  
/* X^(n-2) for n >= 3 using the power computation method, */  
/* where X is defined as follows                          */  
/*                                                        */  
/*     +-       -+   +-      -+(n-2)                      */  
/*     | aa   bb |   | a    b |                           */  
/*     |         | = |        |                           */  
/*     | cc   dd |   | c    d |                           */  
/*     +-       -+   +-      -+                           */  
/* ------------------------------------------------------ */  
  
void matrix_power(unsigned long a,   unsigned long b,   
                  unsigned long c,   unsigned long d, int n,  
                  unsigned long *aa, unsigned long *bb,  
                  unsigned long *cc, unsigned long *dd)  
{  
     unsigned long  xa, xb, xc, xd;  
  
     if (n == 1)  
          *aa = a, *bb = b, *cc = c, *dd = d;  
     else if (n & 0x01 == 1) { /* n odd: X^(2k+1)=X*X^(2k)*/  
          matrix_power(a, b, c, d, n-1, &xa, &xb, &xc, &xd);  
          *aa = a*xa + b*xc;  
          *bb = a*xb + b*xd;  
          *cc = c*xa + d*xc;  
          *dd = c*xb + d*xd;  
     }  
     else {                   /* n even: X^(2k)=(X^(2k))^2*/  
          matrix_power(a, b, c, d, n>>1, &xa, &xb, &xc, &xd);  
          *aa = xa*xa + xb*xc;  
          *bb = xa*xb + xb*xd;  
          *cc = xc*xa + xd*xc;  
          *dd = xc*xb + xd*xd;  
     }  
}  
  
  
/* ------------------------------------------------------ */  
  
#include  <stdio.h>  
#include  <stdlib.h>             /* for function atoi()   */  
  
void main(void)  
{  
     char  line[100];  
     int   n;  
  
     printf("\nO(log(n)) Fibonacci Number Computation\n");  
     printf("\nn please ---> ");  
     gets(line);  
     n = atoi(line);  
     printf("\nfib(%d) = %lu", n, fibonacci(n));  
}  

问题 2.11 扩充 Fibonacci 数（EX_FIB.C )

定义一组称为扩充的Fibonacci数如下，已知X与Y两个数，于是扩充Fibonacci数F为

【说明】

与许多求某个式子的和的程序一样，也不能马上下手编写程序。因为数学家用的式子与编写程序的式子不同，数学家讲求式子要能表达出原意，而编程的式子则要求好算，不仅如此，还要运算得快。如果用上面的式子，为了要用到f0与fn,可能要用一个数组把f0,f1,fn巧算出来后保存起来，再两两相乘、相加。这种方法固然不错，但当n比较大时，内存的要求就大了，所以限定不能用数组。因此，这里的提示是把式子换个面目，看有没有办法实现不用递归来算Fibonacci数的程序。

问题实现

C代码 




 
unsigned long ext_fibonacci(int n, int x, int y)  
{  
     unsigned long  f0, f1;  
     unsigned long  a0, a1;  
     unsigned long  temp1, temp2;  
     int  i;  
  
     if (n == 0)  
          return 1;  
     else if (n == 1)  
          return 2;  
     else {  
          for (f0 = f1 = 1, a0 = 1, a1 = 2, i = 2; i <= n; i++) {  
               temp1 = x*f1 + y*f0;  
               f0    = f1;  
               f1    = temp1;  
               temp2 = x*a0 + y*(a1 - f0) + f0 + f1;  
               a0    = a1;  
               a1    = temp2;  
          }  
          return a1;  
     }  
}  
  
  
/* ------------------------------------------------------ */  
  
#include  <stdio.h>  
#include  <stdlib.h>  
  
void main(void)  
{  
     char line[100];  
     int  n, x, y;  
  
     printf("\nExtended Fibonacci Number Computation\n");  
     printf("\nGiven x and y, define the Extended Fibonacci"  
            " Number as follows\n");  
     printf("\nXfib(0) = 1; Xfib(1) = 1;"  
            " Xfib(n) = x*Xfib(n-1) + y*Xfib(n-2)\n");  
     printf("\nNow given another integer n, compute the sum of");  
     printf("\nXfib(0)*Xfib(n) + Xfib(1)*Xfib(n-1) + ... +"  
            " Xfib(i)*Xfib(n-i) + ... ");  
     printf("\n                + Xfib(n-1)*Xfib(1) + Xfib(n)*Xfib(0)\n");  
  
     printf("\nx ---> ");  
     gets(line);  
     x = atoi(line);  
     printf("\ny ---> ");  
     gets(line);  
     y = atoi(line);  
     printf("\nn ---> ");  
     gets(line);  
     n = atoi(line);  
     printf("\n\nAnswer is %lu", ext_fibonacci(n, x, y));  
}  

问题2.12 二顶式系数加法解（CNR_ADD.C )

编写一个程序，只用加法，求出n中取r个的组合系数C(n，r);并且尽可能地使加法数目降低。

【说明】

有很多种办法可以编写这个程序，但是几乎所有的方法都与二项式系数C(n，r)有关，这个式子就是：

C(n，r)=C(n- 1，r)+C(n- 1，r-1)

于是马上就会有人提笔编写程序了，如程序2-3所示。

【程序】

程序2-3

int cnr(int n, int r)

{

if (n ==r || r == 0) return 1;

else

return cnr(n-1,r} + cnr(n-1,r-1);

}

当然这个程序的确只用了加法，而且也是正确的，但这可能是所有这一类程序中最差劲的一个，为什么呢？

【解答】

只要把用递归方法算C(n，r)时所经过的途径找出来（不过不去计较经过多少次)，那么答案就出来了。以求C(8，3)为例，递归的程序会最先从C(8，3)起，先是递归调用C(7，2)与 C(7，3),在C(7，2)会调用C(6，1)与C(6,2),在C(7，3)调用C(6,2), C(6，3)，如果把递归调用的路径仔细走一次，就会发现它正好经过图11-10中有画线段的部分。

对于每一个c(n，r)而言，如果它左上角与右上角的两个元素都算出来了，那么就可以算出C(n，r)。所以，在图11-10中可以从左上到右下一列一列地算，或者是从右上到左下一行一行地算；换句话说，可以算出C(2，1),C(3，2),C(4，3),再算C(3，1)，C(4，2)，C(5，3)，…；但也可以先算 C(2，1)，C(3，1),C(4，1),C(5，1)，C(6，1),再算 C(3，2)，C(4，2)，…，C(7，2)等。

此处打算用一列一列的方式。在开始时准备一个数组c[],其容量至少有r,它的内容全是1，代表图11-10中的C(0，0)，C(1，1),…,C(r,r)。接下来，对每一个c[i]而言，1≤r≤r, 令 c[i]为 c[i】+c[i-1],即为：C(2，1) = C(1，0) + C(1，1), C(3,2) = C(2,l) + C(2,2)。注意，在算C(2，l)时，它要用到C(1,0),这个值是1,目前在c[0];也要用到C(1，1),现在在c[1]，因此把c[0]与c[1]相加再放回c[1],于是c[1]中的值就变成C(1,1)了。

程序CNR_ADD.C就是用这个观点写成的，在一开始把一个足够大的数组c[]定成1，然后，就把计算的工作循环做n-r次，正好是图11-10中的第一列，第二列,…，第列的计算；对每一列而言，自第一个元素起把它与前一个元素相加。因此计算完后，c[r]中就是 C(n,r)的值。

问题实现

C代码 




 
#define  MAXSIZE  100  
  
unsigned long cnr(int n, int r)  
{  
     unsigned long c[MAXSIZE];  
     int           i, j;  
  
     for (i = 0; i <= r; i++) /* initial c[] to all 1's   */  
          c[i] = 1UL;  
  
     for (i = 1; i <= n-r; i++) /* compute C(n,r) by add. */  
          for (j = 1; j <= r; j++) /* but compute only    */  
               c[j] += c[j-1];  /* those needed.          */  
     return c[r];  
}  

【习题】

(1)请问CNR—ADD.C中一共用了多少个加法？是nr个吗？还是(n-l)(r-l)个？还是其他，请证明结果。

(2)请以C(8，3)为例，用手算追踪一下程序执行的动作，来证实的确了解此处所说的方法。

(3)这个程序有一个缺点，基本上可能会做一些重复的工作，这个潜在的问题就是r 可能会大于n/2。为什么r>n/2时就有重复的工作呢？请修改程序，把这个重复的工作消除。

(4)请改写程序，一行一行地算。请问一行一行地算与一列一列地算在加法个数上有没有差异？内存需求呢？请问哪一种算法比较好，为什么？

问题2.13快速二顶式系数算法（CNR_LOG.C )

看下面的式子：

如果使用的语言或编译程序（如Ada)有可以计算任意位整数加、减、乘、除的能力的话，能够从这个式子找出一个只需要与化

成正比的乘法次数求所有C(n，i)，i=0，l，…，n

的办法吗？下面假设C语言具有此能力，来编写程序。

【说明】

假设给了一个n,能否找出一种方法，用上式把所有C(n，i)，i=0，1，…，n都算出来，这是主题。为了方便起见，不妨假设所有运算都可以处理任意位数的整数，而不会造成溢位；其实有一些语言，如Ada,或是一些链接库都提供了这个能力，所以不必认为如此的假设不切实际。应用也还不少，包含有把71计算到一百万位，测定任何整数是不是质数以及做分解因式的程序，还有就是处理多项式；后者在计算几何学、机器人等应用上都有相当重要的地位。

提示：式子中有一个p，它的值是多少？这是要事先决定的。此外，能在大约

个乘法的限制下求出所有C(n，r)吗？此时n是输入，或许需要用到一些处理数元的技巧。

【问题实现】

/* --------------------------------------------------------- */
/* FUNCTION cnr :                                            */
/*    A special function to calculate all C(n,r)'s given n.  */
/* It employs a very special trick on bit operations.  NOTE  */
/* that depend on the number of bits for unsigned long type  */
/* this function might returns wrong answer.  The rule of    */
/* thumb is that if unsigned long has k bits then the max.   */
/* of n, the input, should be no more than following :       */
/*                                                           */
/*          -1 + sqrt(1 + 4*k)                               */
/*         --------------------                              */
/*                   2                                       */
/*                                                           */
/* Thus for k=32 as many microprocessor has, then the max.   */
/* of n is 5.  But this should not be viewd as a restriction */
/* because the same algorithm can be used in more advanced   */
/* non-numeric applications.                                 */
/*                                                           */
/* Copyright Ching-Kuang SHENE                   Feb/24/1989 */
/* --------------------------------------------------------- */

unsigned long  power(unsigned long, unsigned long);

void  cnr(int n, int answer[])
{
     unsigned long  x    = (1 << n) + 1;  /* 2^n + 1         */
     unsigned long  mask = (1 << n) - 1;  /* 2^n - 1, n 1's  */
     unsigned long  result;
     int      i;

     result = power(x, (unsigned long) n); /* (2^n+1)^n      */

     for (i = 0; i <= n; i++, result >>= n) /* retrieve data */
          answer[i] = result & mask;
}


/* ------------------------------------------------------ */
/* FUNCTION power :                                       */
/*    This is the function called 'iterative_power' in    */
/* file I_POWER.C of this book.                           */
/* ------------------------------------------------------ */

unsigned long power(unsigned long m, unsigned long n)
{
     unsigned long  temp = 1;

     while (n > 0) {          /* if there exists 1 bits.. */
          if (n & 0x01UL == 1)/* the right most one ?     */
               temp *= m;     /* YES, the result get a 'm'*/
          m *= m;             /* anyway, compute m^k      */
          n >>= 1;            /* throw away this bit      */
     }
     return temp;
}


/* ------------------------------------------------------ */

#include  <stdio.h>
#include  <stdlib.h>

#define   MAXSIZE  10

void main(void)
{
     int  answer[MAXSIZE];
     int  n, r;
     char line[100];

     printf("\nFast Combinatorial Coefficient Computation");
     printf("\n==========================================");
     printf("\n\nN ---> ");
     gets(line);
     n = atoi(line);

     cnr(n, answer);
     printf("\nAll Combinatorial Coefficients :\n");
     for (r = 0; r <= n; r++)
          printf("\nC(%d,%d) = %d", n, r, answer[r]);
}

问题2.14快速阶乘运算（FACTLOG2.C )

在问题2.13中己经看过用与

个乘法成正比的技巧来算出C(n，r)，请问有没有办法利用这个结果设计一个求n!的程序，它只用与(

)²成正比的乘法（所有有关硬件、语言的假设同问题2.13)。

【说明】

要做到只用

个乘法，老手们马上就会想到分而治之，完全正确。但是，如果分而治之不在合适的地方分，用错了方法去治，常常就会得不偿失。下面是作者在某篇文章中看到的做法，如程序2-5所示。

【程序】

程序2-5

unsigned long fact(int m, int n)

if (m==n)

return m

else

return fact(m, (m+n)/2)*fact((m+n)/2+1,n);

}

请问这个程序快了没？没有，不但如此，反而慢了。

因此这个方法用了 n-1个乘法，与传统的2,3，…，n —个接一个相乘所用的乘法数目样多，不但如此，还用了递归，因此不论n是多少都肯定会比传统的连乘法慢！所以这个程序的尝试并不成功。

但是，想法却很接近，只要把上面显而易见的分而治之的方法动动脑筋、换个角度，马上就可以得到乘法个数与（

）² 成正比的程序。另外，我们知道计算所有C(n，i)也不过用了

个乘法，因此能不能用某个C(n，i)来帮忙算n!呢？

问题实现

C代码 




 
unsigned long  power(unsigned long, unsigned long);  
  
void  cnr(int n, int answer[])  
{  
     unsigned long  x    = (1 << n) + 1;  /* 2^n + 1         */  
     unsigned long  mask = (1 << n) - 1;  /* 2^n - 1, n 1's  */  
     unsigned long  result;  
     int      i;  
  
     result = power(x, (unsigned long) n); /* (2^n+1)^n      */  
  
     for (i = 0; i <= n; i++, result >>= n) /* retrieve data */  
          answer[i] = result & mask;  
}  
  
  
/* ------------------------------------------------------ */  
/* FUNCTION power :                                       */  
/*    This is the function called 'iterative_power' in    */  
/* file I_POWER.C of this book.                           */  
/* ------------------------------------------------------ */  
  
unsigned long power(unsigned long m, unsigned long n)  
{  
     unsigned long  temp = 1;  
  
     while (n > 0) {          /* if there exists 1 bits.. */  
          if (n & 0x01UL == 1)/* the right most one ?     */  
               temp *= m;     /* YES, the result get a 'm'*/  
          m *= m;             /* anyway, compute m^k      */  
          n >>= 1;            /* throw away this bit      */  
     }  
     return temp;  
}  
  
  
/* ------------------------------------------------------ */  
  
#include  <stdio.h>  
#include  <stdlib.h>  
  
#define   MAXSIZE  10  
  
void main(void)  
{  
     int  answer[MAXSIZE];  
     int  n, r;  
     char line[100];  
  
     printf("\nFast Combinatorial Coefficient Computation");  
     printf("\n==========================================");  
     printf("\n\nN ---> ");  
     gets(line);  
     n = atoi(line);  
  
     cnr(n, answer);  
     printf("\nAll Combinatorial Coefficients :\n");  
     for (r = 0; r <= n; r++)  
          printf("\nC(%d,%d) = %d", n, r, answer[r]);  
}  

问题2.15更快的阶乘算法（FACTLOGC )

请仔细研究问题2.14,编写出一个只用与

成正比的乘法个数，计算n!的程序。

【说明】

问题2.14的解法是有进一步改良的余地的。用递归方法，通过下式计算:

一样，它也有一些被重复计算的地方，如果能够把它找出来，程序自然就可以做到乘法个数与lon2n成正比的地步。仔细做一做问题2.14的习题（1)与习题（2),必会有所收获。

问题实现

C代码 




 
unsigned long  mask;          /* mask for shifting bits     */  
unsigned cnrs[100];           /* storage for C(n,n/2)'s     */  
int      number;              /* index for cnrs[]           */  
int      p_size;              /* partition size in C(n,n/2) */  
  
unsigned long  factor(unsigned long);  
unsigned long  power(unsigned long, unsigned long);  
  
  
/* -------------------------------------------------------- */  
/* FUNCTION factorial :                                     */  
/*    Control routine of n! calculation.                    */  
/* -------------------------------------------------------- */  
  
unsigned long  factorial(unsigned long n)  
{  
     unsigned long  x = (1 << n) + 1;  
  
     number = 0;              /* set index to init. pos.    */  
     mask   = (1 << n) - 1;   /* set up mask for shifting   */  
     p_size = n;              /* set up partition size      */  
     (void) power(x, n);      /* now computing (2^n+1)^n    */  
     number = 0;              /* reset index to the start   */  
     return factor(n);        /* then compute n!            */  
}  
  
  
/* -------------------------------------------------------- */  
/* FUNCTION factor :                                        */  
/*    This is the working routine for n!.  The idea is      */  
/*                                                          */  
/*          +-  1                  if n = 1                 */  
/*          |                                               */  
/*    n! =  +   n * (n-1)!         if n is odd              */  
/*          |                                               */  
/*          +-  C(n,n/2)*[(n/2)]!  if n is even             */  
/*                                                          */  
/* Before factor() is called n_power() is used to compute   */  
/* (2^n+1)^n and during the course of this computation, it  */  
/* stores all relevant C(n,n/2) into array cnrs[].  Thus    */  
/* C(n,n/2)'s are available from array C(n,n/2) here.       */  
/* -------------------------------------------------------- */  
  
unsigned long  factor(unsigned long n)  
{  
     unsigned long  temp;  
  
     if (n == 1)  
          return 1;  
     else if (n & 0x01UL == 1)  
          return n * factor(n - 1);  
     else {  
          temp = factor(n >> 1);  
          return cnrs[number++] * temp * temp;  
     }  
}  
  
  
/* -------------------------------------------------------- */  
/* FUNCTION power :                                         */  
/*    The modified power function to compute (2^n+1)^n in   */  
/* recursive mode.  During its computation, n_power() will  */  
/* save C(n,n/2) into array cnrs[].                         */  
/* -------------------------------------------------------- */  
  
unsigned long  power(unsigned long n, unsigned long m)  
{  
     unsigned long  temp;  
  
     if (m == 1)  
          temp = n;  
     else if (m & 0x01UL != 0)  
          temp = n * power(n, m-1);  
     else {  
          temp = power(n, m >> 1);  
          temp *= temp;  
          cnrs[number++] = (temp >> ((m >> 1)*p_size)) & mask;  
     }  
     return temp;  
}  
  
  
/* -------------------------------------------------------- */  
  
#include  <stdio.h>  
#include  <stdlib.h>          /* for strtoul() function     */  
  
void main(void)  
{  
     char  line[100], *dummy_ptr;  
     unsigned long m, n;  
  
     printf("\nThe Fast N! Computation :");  
     printf("\n\nN for N! please --> ");  
     gets(line);  
     n = strtoul(line, &dummy_ptr, 10);  
     printf("\n%lu! = %lu", n, factorial(n));  
}  

问题2.16连续整数的固定和(GIVENSUM.C )

编写一个程序，读入一个正整数，把所有那些连续的和为给定的正整数的正整数找出来。例如，如果输入27,发现2〜7、8〜10、13与14的和是27，这就是解答；如果输入的是10000,应该有18〜142、297〜1328、388〜412、1998〜2002这4组。注意，不见得一定会有答案，譬如说4、16就无解；另外，排除只有一个数的情况，否则每一个输入就都至少有一个答案，就是它自己。

【说明】

任何人看到这个题目都会马上想到一个办法，把所有的和算出来，与输入比较。曾经看到过如下的一个解法，如程序2-6所示。

【程序】

程序 2-6 (NOT_GOOD.C)

#include <stdio.h>

#include <stdlib.h>

void main(void)

{

char line[100]; long sum, i, j, n;

gets(line); n = atol(line);

for (i = 1; i <= n/2; i++)

{

for (sum = i, j = i + 1; j <= n/2+1; j++) {

sum += j;

if (sum == n)

printf {" \n%ld - %ld", i, j);

}

它的做法是先固定一个i、sum变量，接着令j从i+1起变化，每次都把j的值加到sum 中，因此sum中的值就是这些连续整数的和。因此令i从1变到n/2 (n是给定的值)，而j从i+1变到n/2+l，如果有一个和与n相同，就显示i与j，表示n的值是i到j 这一串连续的正整数的和。为什么i要从1到n/2?很简单，如果i是n/2,下一个数就是 n/2+l,于是这两个（连续的）数的和n/2+(n/2+l)=n+l就大于n，所以i最多只能到n/2; 同理可以说明j不可以大过n/2+1。

这个程序当然是对的，但运行太慢了！用1_作为输入，它一共执行了 311.71秒，也就是5分钟多；但事实上，这个题目可以在不到一秒之内得出答案，而且当输入1000000 (100万）时，也不过用178.12秒（3分钟）左右而己，相比之下NOTLCOOD.C的效率实在太低了。

问题出在什么地方？加法次数太多了。在上面的程序中，i与j的关系永远满足 1≤i≤n , i + 1≤j≤n + 1, i<j;每一组i与j都会做一次加法，所以就一共做了大约n²/8个加法（这是个近似值)；当n=10000时，就大约是1250万个。

不过，一个好程序员应该研究是否有更好的方法存在，事实上就有一个，大约需要2n个加法就足够了，能想得出来吗？下面是几点有用的提示：

(1)如果在求和时是用 i+ (i + 1)+…+j表示，那么i≤n/2;这是上面提过的。

(2)如果某个和i+(i+ 1)+…+比给定的值大，那么把和变小，但仍然维持是一串连续整数的和时，拿掉j变成i + (i + 1)+…十(j-1),不如拿掉i变成(i + l) +…+j。为什么？因为j比i大，拿掉j，和就下降太快了，不如拿掉i,再慢慢降低（能用数学来证明吗？）

(3)如果和 i+ (i +1) +…+j比给定值小，加上j +1变成的i+…+j十(j-1);道理同前。

有了这几点，编程应该不会是件难事了。

问题实现

C代码 




 
#include <stdio.h>  
#include  <stdlib.h>  
  
void main(void)  
{  
     long  left, right;       /* left and right bound     */  
     long  sum;               /* computed sum             */  
     long  GIVEN;             /* the given number         */  
     int   count = 0;         /* solution counter         */  
     char  line[100];         /* input buffer             */  
  
     printf("\nConsecutive sum to a fixed given number");  
     printf("\n=======================================\n");  
     printf("\nYour number (> 0) please ---> ");  
     gets(line);  
     GIVEN = atol(line);  
  
     for (sum = 0, right = 1; sum < GIVEN; sum += right, right++)  
          ;  
  
     for (left = 1, right--; left <= GIVEN/2; )  
          if (sum > GIVEN)         /* if sum too large    */  
               sum -= (left++);    /* reduce from left    */  
          else {  
               if (sum == GIVEN) { /* if equal, print it  */  
                    printf("\n%ld = sum from %ld to %ld",   
                           GIVEN, left, right);  
                    count++;       /* and increase count  */  
               }                   /* if sum too small    */  
               sum += (++right);   /* increase from right */  
          }  
  
     if (count > 0)  
          printf("\n\nThere are %d solutions in total.", count);  
     else  
          printf("\n\nSorry, there is NO solution at all.");  
}  