OSX/libsecurity_cryptkit/lib/CurveParamDocs/tools.c

   1 /**************************************************************
   2  *
   3  *      tools.c
   4  *
   5  *      Number-theoretical algorithm implementations
   6  *
   7  *      Updates:
   8  *     30 Apr 99    REC  Modified init_tools type to void.
   9  *              3 Apr 98    REC  Creation
  10  *
  11  *
  12  *      c. 1998 Perfectly Scientific, Inc.
  13  *      All Rights Reserved.
  14  *
  15  *
  16  *************************************************************/
  17
  18 /* include files */
  19
  20 #include <stdio.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23 #include <time.h>
  24 #ifdef _WIN32
  25
  26 #include <process.h>
  27
  28 #endif
  29
  30 #include <string.h>
  31 #include "giants.h"
  32 #include "tools.h"
  33
  34 /* definitions */
  35
  36 #define STACK_COUNT 20
  37
  38 /* global variables */
  39
  40 int     pr[NUM_PRIMES];   /* External use allowed. */
  41 static giant tmp[STACK_COUNT];
  42 static int stack = 0;
  43 static giant popg();
  44 static void pushg();
  45
  46 /**************************************************************
  47  *
  48  *      Maintenance functions
  49  *
  50  **************************************************************/
  51
  52
  53 void
  54 init_tools(int shorts)
  55 {
  56         int j;
  57
  58         for(j = 0; j < STACK_COUNT; j++) {
  59                 tmp[j] = newgiant(shorts);
  60         }
  61         make_primes();  /* Create table of all primes < 2^16,
  62                                            to be used by other programs as array
  63                                            pr[0..NUM_PRIMES]. */
  64 }
  65
  66 static giant
  67 popg() {
  68         return(tmp[stack++]);
  69 }
  70
  71 static void
  72 pushg(int n) {
  73         stack -= n;
  74 }
  75
  76 /**************************************************************
  77  *
  78  *      Number-theoretical functions
  79  *
  80  **************************************************************/
  81
  82 int
  83 prime_literal(
  84         unsigned int    p
  85 )
  86 /* Primality test via small, literal sieve.
  87    After init, one should use primeq() instead.
  88  */
  89 {
  90         unsigned  int   j=3;
  91
  92         if ((p & 1)==0)
  93                 return ((p == 2)?1:0);
  94         if (j >= p)
  95                 return (1);
  96         while ((p%j)!=0)
  97         {
  98                 j += 2;
  99                 if (j*j > p)
 100                         return(1);
 101         }
 102         return(0);
 103 }
 104
 105 int
 106 primeq(
 107         unsigned int odd
 108 )
 109 /* Faster primality test, using preset array pr[] of primes.
 110    This test is valid for all unsigned, 32-bit integers odd.
 111  */
 112 {
 113         unsigned int p;
 114         unsigned int j;
 115
 116     if(odd < 2) return (0);
 117         if ((odd & 1)==0)
 118                 return ((odd == 2)?1:0);
 119         for (j=1; ;j++)
 120         {
 121                 p = pr[j];
 122                 if (p*p > odd)
 123                         return(1);
 124                 if (odd % p == 0)
 125                         return(0);
 126         }
 127 }
 128
 129 void
 130 make_primes()
 131 {   int k, npr;
 132         pr[0] = 2;
 133         for (k=0, npr=1;; k++)
 134         {
 135                 if (prime_literal(3+2*k))
 136                 {
 137                         pr[npr++] = 3+2*k;
 138                         if (npr >= NUM_PRIMES)
 139                                 break;
 140                 }
 141         }
 142 }
 143
 144 int
 145 prime_probable(giant p)
 146 /* Invoke Miller-Rabin test of given security depth.
 147    For MILLER_RABIN_DEPTH == 8, this is an ironclad primality
 148    test for suspected primes p < 3.4 x 10^{14}.
 149 */
 150 {
 151         giant t1 = popg(), t2 = popg(), t3 = popg();
 152     int j, ct, s;
 153
 154     if((p->n[0] & 1) == 0) {  /* Evenness test. */
 155                 pushg(3); return(0);
 156         }
 157     if(bitlen(p) < 32) {  /* Single-word case. */
 158                 pushg(3);
 159                 return(primeq((unsigned int)gtoi(p)));
 160         }
 161         itog(-1, t1);
 162         addg(p, t1);  /* t1 := p-1. */
 163         gtog(t1, t2);
 164         s = 1;
 165         gshiftright(1, t2);
 166         while(t2->n[0] & 1 == 0) {
 167                 gshiftright(1, t2);
 168                 ++s;
 169         }
 170         /* Now, p-1 = 2^s * t2. */
 171         for(j = 0; j < MILLER_RABIN_DEPTH; j++) {
 172                 itog(pr[j+1], t3);
 173                 powermodg(t3, t2, p);
 174                 ct = 1;
 175                 if(isone(t3)) continue;
 176                 if(gcompg(t3, t1) == 0) continue;
 177                 while((ct < s) && (gcompg(t1, t3) != 0)) {
 178                         squareg(t3); modg(p, t3);
 179                         if(isone(t3)) {
 180                                 goto composite;
 181                         }
 182                         ++ct;
 183                 }
 184                 if(gcompg(t1, t3) != 0) goto composite;
 185         }
 186         goto prime;
 187
 188 composite:
 189                 pushg(3); return(0);
 190 prime:  pushg(3); return(1);
 191 }
 192
 193 int
 194 jacobi_symbol(giant a, giant n)
 195 /* Standard Jacobi symbol (a/n).  Parameter n must be odd, positive. */
 196 {       int t = 1, u;
 197         giant t5 = popg(), t6 = popg(), t7 = popg();
 198
 199         gtog(a, t5); modg(n, t5);
 200         gtog(n, t6);
 201         while(!isZero(t5)) {
 202             u = (t6->n[0]) & 7;
 203                 while((t5->n[0] & 1) == 0) {
 204                         gshiftright(1, t5);
 205                         if((u==3) || (u==5)) t = -t;
 206                 }
 207                 gtog(t5, t7); gtog(t6, t5); gtog(t7, t6);
 208                 u = (t6->n[0]) & 3;
 209                 if(((t5->n[0] & 3) == 3) && ((u & 3) == 3)) t = -t;
 210                 modg(t6, t5);
 211         }
 212         if(isone(t6)) {
 213                 pushg(3);
 214                 return(t);
 215         }
 216         pushg(3);
 217         return(0);
 218 }
 219
 220 int
 221 pseudoq(giant a, giant p)
 222 /* Query whether a^(p-1) = 1 (mod p). */
 223 {
 224         int x;
 225         giant t1 = popg(), t2 = popg();
 226
 227         gtog(p, t1); itog(1, t2); subg(t2, t1);
 228         gtog(a, t2);
 229         powermodg(t2, t1, p);
 230         x = isone(t2);
 231         pushg(2);
 232         return(x);
 233 }
 234
 235 int
 236 pseudointq(int a, giant p)
 237 /* Query whether a^(p-1) = 1 (mod p). */
 238 {
 239         int x;
 240         giant t4 = popg();
 241
 242         itog(a, t4);
 243         x = pseudoq(t4, p);
 244         pushg(1);
 245         return(x);
 246 }
 247
 248
 249 void
 250 powFp2(giant a, giant b, giant w2, giant n, giant p)
 251 /* Perform powering in the field F_p^2:
 252    a + b w := (a + b w)^n (mod p), where parameter w2 is a quadratic
 253    nonresidue (formally equal to w^2).
 254  */
 255 {   int j;
 256     giant t6 = popg(), t7 = popg(), t8 = popg(), t9 = popg();
 257
 258         if(isZero(n)) {
 259                 itog(1,a);
 260                 itog(0,b);
 261                 pushg(4);
 262                 return;
 263         }
 264         gtog(a, t8); gtog(b, t9);
 265         for(j = bitlen(n)-2; j >= 0; j--) {
 266                 gtog(b, t6);
 267                 mulg(a, b); addg(b,b); modg(p, b);  /* b := 2 a b. */
 268                 squareg(t6); modg(p, t6);
 269                 mulg(w2, t6); modg(p, t6);
 270                 squareg(a); addg(t6, a); modg(p, a); /* a := a^2 + b^2 w2. */
 271                 if(bitval(n, j)) {
 272                         gtog(b, t6); mulg(t8, b); modg(p, b);
 273                         gtog(a, t7); mulg(t9, a); addg(a, b); modg(p, b);
 274                         mulg(t9, t6); modg(p, t6); mulg(w2, t6); modg(p, t6);
 275                         mulg(t8, a); addg(t6, a); modg(p, a);
 276                 }
 277         }
 278         pushg(4);
 279         return;
 280 }
 281
 282 int
 283 sqrtmod(giant p, giant x)
 284 /* If Sqrt[x] (mod p) exists, function returns 1, else 0.
 285    In either case x is modified, but if 1 is returned,
 286    x:= Sqrt[x] (mod p).
 287  */
 288 {   giant t0 = popg(), t1 = popg(), t2 = popg(), t3 = popg(),
 289                   t4 = popg();
 290
 291         modg(p, x);   /* Justify the argument. */
 292     gtog(x, t0);  /* Store x for eventual validity check on square root. */
 293     if((p->n[0] & 3) == 3) {  /* The case p = 3 (mod 4). */
 294                 gtog(p, t1);
 295                 iaddg(1, t1); gshiftright(2, t1);
 296                 powermodg(x, t1, p);
 297                 goto resolve;
 298         }
 299 /* Next, handle case p = 5 (mod 8). */
 300     if((p->n[0] & 7) == 5) {
 301                 gtog(p, t1); itog(1, t2);
 302                 subg(t2, t1); gshiftright(2, t1);
 303                 gtog(x, t2);
 304                 powermodg(t2, t1, p);  /* t2 := x^((p-1)/4) % p. */
 305                 iaddg(1, t1);
 306                 gshiftright(1, t1); /* t1 := (p+3)/8. */
 307                 if(isone(t2)) {
 308                         powermodg(x, t1, p);  /* x^((p+3)/8) is root. */
 309                         goto resolve;
 310                 } else {
 311                         itog(1, t2); subg(t2, t1);  /* t1 := (p-5)/8. */
 312                         gshiftleft(2,x);
 313                         powermodg(x, t1, p);
 314                         mulg(t0, x); addg(x, x); modg(p, x); /* 2x (4x)^((p-5)/8. */
 315                         goto resolve;
 316                 }
 317         }
 318
 319 /* Next, handle tougher case: p = 1 (mod 8). */
 320         itog(2, t1);
 321         while(1) {  /* Find appropriate nonresidue. */
 322                 gtog(t1, t2);
 323                 squareg(t2); subg(x, t2); modg(p, t2);
 324                 if(jacobi_symbol(t2, p) == -1) break;
 325                 iaddg(1, t1);
 326         }  /* t2 is now w^2 in F_p^2. */
 327     itog(1, t3);
 328     gtog(p, t4); iaddg(1, t4); gshiftright(1, t4);
 329         powFp2(t1, t3, t2, t4, p);
 330         gtog(t1, x);
 331
 332 resolve:
 333     gtog(x,t1); squareg(t1); modg(p, t1);
 334     if(gcompg(t0, t1) == 0) {
 335                 pushg(5);
 336                 return(1);  /* Success. */
 337         }
 338         pushg(5);
 339         return(0);  /* No square root. */
 340 }
 341
 342 void
 343 sqrtg(giant n)
 344 /* n:= Floor[Sqrt[n]]. */
 345 {   giant t5 = popg(), t6 = popg();
 346
 347         itog(1, t5); gshiftleft(1 + bitlen(n)/2, t5);
 348         while(1) {
 349                 gtog(n, t6);
 350                 divg(t5, t6);
 351                 addg(t5, t6); gshiftright(1, t6);
 352                 if(gcompg(t6, t5) >= 0) break;
 353                 gtog(t6, t5);
 354         }
 355     gtog(t5, n);
 356         pushg(2);
 357 }
 358
 359 int
 360 cornacchia4(giant n, int d, giant u, giant v)
 361 /* Seek representation 4n = u^2 + |d| v^2,
 362    for (negative) discriminant d and n > |D|/4.
 363    Parameter u := 0 and 0 is returned, if no representation is found;
 364    else 1 is returned and u, v properly set.
 365  */
 366 {   int r = n->n[0] & 7, sym;
 367         giant t1 = popg(), t2 = popg(), t3 = popg(), t4 = popg();
 368
 369         itog(d, t1);
 370         if((n->n[0]) & 7 == 1) {  /* n = 1 (mod 8). */
 371                 sym = jacobi_symbol(t1,n);
 372                 if(sym != 1) {
 373                         itog(0,u);
 374                         pushg(4);
 375                         return(0);
 376                 }
 377                 gtog(t1, t2);
 378                 sqrtmod(n, t2);  /* t2 := Sqrt[d] (mod n). */
 379     } else {  /* Avoid separate Jacobi/Legendre test. */
 380                 gtog(t1, t2);
 381                 if(sqrtmod(n, t2) == 0) {
 382                         itog(0, u);
 383                         pushg(4);
 384                         return(0);
 385                 }
 386         }
 387 /* t2 is now a valid square root of d (mod n). */
 388         gtog(t2, t3);
 389         subg(t1, t3); /* t3 := t2 - d. */
 390         if((t3->n[0] & 1) == 1) {
 391                 negg(t2);
 392                 addg(n, t2);
 393         }
 394         gtog(n, t3); addg(t3, t3);  /* t3 := 2n. */
 395         gtog(n, t4); gshiftleft(2, t4); sqrtg(t4); /* t4 = [Sqrt[4 p]]. */
 396         while(gcompg(t2, t4) > 0) {
 397                 gtog(t3, t1);
 398                 gtog(t2, t3);
 399                 gtog(t1, t2);
 400                 modg(t3, t2);
 401         }
 402         gtog(n, t4); gshiftleft(2, t4);
 403         gtog(t2, t3); squareg(t3);
 404         subg(t3, t4); /* t4 := 4n - t2^2. */
 405         gtog(t4, t3);
 406         itog(d, t1); absg(t1);
 407         modg(t1, t3);
 408         if(!isZero(t3)) {
 409                 itog(0,u);
 410                 pushg(4);
 411                 return(0);
 412         }
 413         divg(t1, t4);
 414         gtog(t4, t1);
 415         sqrtg(t4); /* t4 := [Sqrt[t4/Abs[d]]]. */
 416         gtog(t4, t3);
 417         squareg(t3);
 418         if(gcompg(t3, t1) != 0) {
 419                 itog(0, u);
 420                 pushg(4);
 421                 return(0);
 422         }
 423         gtog(t2, u);
 424         gtog(t4, v);
 425         pushg(4);
 426         return(1);
 427 }
 428
 429 /*
 430 rep[p_, d_] := Module[{t, x0, a, b, c},
 431                 If[JacobiSymbol[d,p] != 1, Return[{0,0}]];
 432                 x0 = sqrt[d, p];
 433                 If[Mod[x0-d,2] == 1, x0 = p-x0];
 434                 a = 2p; b = x0; c = sqrtint[4 p];
 435                 While[b > c, {a,b} = {b, Mod[a,b]}];
 436                 t = 4p - b^2;
 437                 If[Mod[t,Abs[d]] !=0, Return[{0,0}]];
 438                 v = t/Abs[d];
 439                 u = sqrtint[v];
 440                 If[u^2 != v, Return[{0,0}]];
 441                 Return[{b, u}]
 442                 ];
 443 */
 444
 445