home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 135_01 / m.c < prev    next >
Encoding:
C/C++ Source or Header  |  1985-03-09  |  1.9 KB  |  84 lines
  1. /*
  2.     m.c---an implementation of Fermat's little theorem
  3.     as a practical test of primality for microcomputers.
  4.     In a nutshell, if N is a prime number and B is any
  5.     whole number, then B^N - B is a multiple of N, or
  6.     stated another way, mod(B^N-B,N) is congruent to 0 if N
  7.     is prime.  This does not guarentee that N is prime (the
  8.     existance of pseudoprimes and carmichael numbers precludes
  9.     that) but if carried out on a random selection of numbers
  10.     in the range of 2...N-1 passage of each test should allow
  11.     reasonable confidence of primality.  If we test say
  12.     100 numbers in this range then the chance that N is not
  13.     prime is 1/2^100, an exceedingly small number indeed!
  14.  
  15.     by Hugh S. Myers
  16.  
  17.     build:
  18.         n>cc m
  19.         n>clink m vli math
  20.  
  21.     3/13/84
  22.     4/3/84
  23. */
  24.  
  25. #include <bdscio.h>
  26. #define MAX 256
  27.  
  28. main()
  29. {
  30.     char N[MAX], B[MAX], r[MAX];
  31.  
  32.         printf("test for mod(B^N-B,N) === 0\n");
  33.     forever {
  34.         printf("?B ");
  35.         getline(B,MAX);
  36.         printf("?N ");
  37.         getline(N,MAX);
  38.     
  39.         strcpy(r,modp(B,N,N));
  40.         printf("mod(B^N,N) = %s\n",r);
  41.         if(isequal(r,B))
  42.             printf("N is possibly prime\n");
  43.         else {
  44.             strcpy(r,subl(r,B));
  45.             strcpy(r,modl(r,N));
  46.             if(iszero(r))
  47.                 printf("N is possibly prime\n");
  48.             else
  49.                 printf("N is not a prime\n");
  50.         }
  51.     }
  52. }
  53.  
  54. /*
  55.     char *modp(s1,s2,s3) return mod(s1^s2,s3) where s1^s2 is
  56.     based on Algorithm A of "Seminumerical Algorithms: The
  57.     Art of Computer Programming", Vol. 2, second edition, by
  58.     D.E. Knuth.
  59. */
  60. char *modp(s1,s2,s3)
  61. char *s1, *s2, *s3;
  62. {
  63.     char y[MAX], z[MAX], n[MAX], M[MAX];
  64.     int odd;
  65.  
  66.     strcpy(y,"1");
  67.     strcpy(z,s1);
  68.     strcpy(n,s2);
  69.     strcpy(M,s3);
  70.  
  71.     forever {
  72.         odd = (!iseven(n));
  73.         strcpy(n,sdivl(n,2));
  74.         if (odd) {
  75.             strcpy(y,mull(y,z));
  76.             strcpy(y,modl(y,M));
  77.             if (iszero(n)) 
  78.                 return y;
  79.         }
  80.         strcpy(z,mull(z,z));
  81.         strcpy(z,modl(z,M));
  82.     }
  83. }
  84.