home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_200 / 247_01 / bnsmall.c < prev    next >
Text File  |  1989-04-19  |  4KB  |  183 lines
  1. /*
  2.  * MIRACL small number theoretic routines 
  3.  * bnsmall.c
  4.  */
  5.  
  6. #include <stdio.h>
  7. #include <math.h>
  8. #include "miracl.h"
  9.  
  10. /* access global variables */
  11.  
  12. extern int depth;     /* error tracing... */
  13. extern int trace[];   /* ....mechanism    */
  14. extern char *calloc();
  15.  
  16. int *gprime(maxp)
  17. int maxp;
  18. { /* generate all primes less than maxp */
  19.     int *primes;
  20.     double dn;
  21.     char *sv;
  22.     int pix,i,k,prime;
  23.     if (ERNUM) return NULL;
  24.     if (maxp<0)
  25.     { /* generate (-maxp) primes */
  26.         dn=(-maxp);
  27.         if (dn<20.0) dn=20.0;
  28.         maxp=(int)(dn*(log(dn*log(dn))-0.5)); /* Rossers upper bound */
  29.     }
  30.     depth++;
  31.     trace[depth]=70;
  32.     if (TRACER) track();
  33.     if (maxp>=TOOBIG)
  34.     {
  35.          berror(8);
  36.          depth--;
  37.          return NULL;
  38.     }
  39.     maxp=(maxp+1)/2;
  40.     sv=(char *)calloc(maxp,1);
  41.     if (sv==NULL)
  42.     {
  43.         berror(8);
  44.         depth--;
  45.         return NULL;
  46.     }
  47.     pix=1;
  48.     for (i=0;i<maxp;i++)
  49.         sv[i]=TRUE;
  50.     for (i=0;i<maxp;i++)
  51.     if (sv[i])
  52.     { /* using sieve of Eratosthenes */
  53.         prime=i+i+3;
  54.         for (k=i+prime;k<maxp;k+=prime)
  55.             sv[k]=FALSE;
  56.         pix++;
  57.     }
  58.     primes=(int *)calloc(pix+2,sizeof(int));
  59.     if (primes==NULL)
  60.     {
  61.         berror(8);
  62.         depth--;
  63.         return NULL;
  64.     }
  65.     primes[0]=2;
  66.     pix=1;
  67.     for (i=0;i<maxp;i++)
  68.         if (sv[i]) primes[pix++]=i+i+3;
  69.     primes[pix]=0;
  70.     free(sv);
  71.     depth--;
  72.     return primes;
  73. }
  74.  
  75. small smul(x,y,n)
  76. small x,y,n;
  77. { /* returns x*y mod n */
  78.     small r;
  79.     x%=n;
  80.     y%=n;
  81.     if (x<0) x+=n;
  82.     if (y<0) y+=n;
  83.     muldiv(x,y,(small)0,n,&r);
  84.     return r;
  85. }
  86.  
  87. small spmd(x,n,m)
  88. small x,n,m;
  89. { /*  returns x^n mod m  */
  90.     small r;
  91.     x%=m;
  92.     if (x<0) x+=m;
  93.     r=0;
  94.     if (x==0) return r;
  95.     r=1;
  96.     if (n==0) return r;
  97.     forever
  98.     {
  99.         if (n%2!=0) muldiv(r,x,(small)0,m,&r);
  100.         n/=2;
  101.         if (n==0) return r;
  102.         muldiv(x,x,(small)0,m,&x);
  103.     }
  104. }
  105.  
  106. small inverse(x,y)
  107. small x,y;
  108. { /* returns inverse of x mod y */
  109.     small r,s,q,t,p;
  110.     x%=y;
  111.     if (x<0) x+=y;
  112.     r=1;
  113.     s=0;
  114.     p=y;
  115.     while (p!=0)
  116.     { /* main euclidean loop */
  117.         q=x/p;
  118.         t=r-s*q;
  119.         r=s;
  120.         s=t;
  121.         t=x-p*q;
  122.         x=p;
  123.         p=t;
  124.     }
  125.     if (r<0) r+=y;
  126.     return r;
  127. }
  128.  
  129. small sqrmp(x,m)
  130. small x,m;
  131. { /* square root mod a small prime by Shanks method  *
  132.    * returns 0 if root does not exist or m not prime */
  133.     small q,z,y,v,w,t;
  134.     int i,r,e,n;
  135.     bool pp;
  136.     x%=m;
  137.     if (x==0) return 0;
  138.     if (x<0) x+=m;
  139.     if (x==1) return 1;
  140.     if (spmd(x,(m-1)/2,m)!=1) return 0;    /* Legendre symbol not 1   */
  141.     if (m%4==3) return spmd(x,(m+1)/4,m);  /* easy case for m=4.k+3   */
  142.     q=m-1;
  143.     e=0;
  144.     while (q%2==0) 
  145.     {
  146.         q/=2;
  147.         e++;
  148.     }
  149.     if (e==0) return 0;      /* even m */
  150.     for (r=2;;r++)
  151.     { /* find suitable z */
  152.         z=spmd((small)r,q,m);
  153.         if (z==1) continue;
  154.         t=z;
  155.         pp=FALSE;
  156.         for (i=1;i<e;i++)
  157.         { /* check for composite m */
  158.             if (t==(m-1)) pp=TRUE;
  159.             muldiv(t,t,(small)0,m,&t);
  160.             if (t==1 && !pp) return 0;
  161.         }
  162.         if (t==(m-1)) break;
  163.         if (!pp) return 0;   /* m is not prime */
  164.     }
  165.     y=z;
  166.     r=e;
  167.     v=spmd(x,(q+1)/2,m);
  168.     w=spmd(x,q,m);
  169.     while (w!=1)
  170.     {
  171.         t=w;
  172.         for (n=0;t!=1;n++) muldiv(t,t,(small)0,m,&t);
  173.         if (n>=r) return 0;
  174.         y=spmd(y,(small)(1<<(r-n-1)),m);
  175.         muldiv(v,y,(small)0,m,&v);
  176.         muldiv(y,y,(small)0,m,&y);
  177.         muldiv(w,y,(small)0,m,&w);
  178.         r=n;
  179.     }
  180.     return v;
  181. }
  182.  
  183.