+/**************************************************************
+ *
+ * tools.c
+ *
+ * Number-theoretical algorithm implementations
+ *
+ * Updates:
+ * 30 Apr 99 REC Modified init_tools type to void.
+ * 3 Apr 98 REC Creation
+ *
+ *
+ * c. 1998 Perfectly Scientific, Inc.
+ * All Rights Reserved.
+ *
+ *
+ *************************************************************/
+
+/* include files */
+
+#include <stdio.h>
+#include <math.h>
+#include <stdlib.h>
+#include <time.h>
+#ifdef _WIN32
+
+#include <process.h>
+
+#endif
+
+#include <string.h>
+#include "giants.h"
+#include "tools.h"
+
+/* definitions */
+
+#define STACK_COUNT 20
+
+/* global variables */
+
+int pr[NUM_PRIMES]; /* External use allowed. */
+static giant tmp[STACK_COUNT];
+static int stack = 0;
+static giant popg();
+static void pushg();
+
+/**************************************************************
+ *
+ * Maintenance functions
+ *
+ **************************************************************/
+
+
+void
+init_tools(int shorts)
+{
+ int j;
+
+ for(j = 0; j < STACK_COUNT; j++) {
+ tmp[j] = newgiant(shorts);
+ }
+ make_primes(); /* Create table of all primes < 2^16,
+ to be used by other programs as array
+ pr[0..NUM_PRIMES]. */
+}
+
+static giant
+popg() {
+ return(tmp[stack++]);
+}
+
+static void
+pushg(int n) {
+ stack -= n;
+}
+
+/**************************************************************
+ *
+ * Number-theoretical functions
+ *
+ **************************************************************/
+
+int
+prime_literal(
+ unsigned int p
+)
+/* Primality test via small, literal sieve.
+ After init, one should use primeq() instead.
+ */
+{
+ unsigned int j=3;
+
+ if ((p & 1)==0)
+ return ((p == 2)?1:0);
+ if (j >= p)
+ return (1);
+ while ((p%j)!=0)
+ {
+ j += 2;
+ if (j*j > p)
+ return(1);
+ }
+ return(0);
+}
+
+int
+primeq(
+ unsigned int odd
+)
+/* Faster primality test, using preset array pr[] of primes.
+ This test is valid for all unsigned, 32-bit integers odd.
+ */
+{
+ unsigned int p;
+ unsigned int j;
+
+ if(odd < 2) return (0);
+ if ((odd & 1)==0)
+ return ((odd == 2)?1:0);
+ for (j=1; ;j++)
+ {
+ p = pr[j];
+ if (p*p > odd)
+ return(1);
+ if (odd % p == 0)
+ return(0);
+ }
+}
+
+void
+make_primes()
+{ int k, npr;
+ pr[0] = 2;
+ for (k=0, npr=1;; k++)
+ {
+ if (prime_literal(3+2*k))
+ {
+ pr[npr++] = 3+2*k;
+ if (npr >= NUM_PRIMES)
+ break;
+ }
+ }
+}
+
+int
+prime_probable(giant p)
+/* Invoke Miller-Rabin test of given security depth.
+ For MILLER_RABIN_DEPTH == 8, this is an ironclad primality
+ test for suspected primes p < 3.4 x 10^{14}.
+*/
+{
+ giant t1 = popg(), t2 = popg(), t3 = popg();
+ int j, ct, s;
+
+ if((p->n[0] & 1) == 0) { /* Evenness test. */
+ pushg(3); return(0);
+ }
+ if(bitlen(p) < 32) { /* Single-word case. */
+ pushg(3);
+ return(primeq((unsigned int)gtoi(p)));
+ }
+ itog(-1, t1);
+ addg(p, t1); /* t1 := p-1. */
+ gtog(t1, t2);
+ s = 1;
+ gshiftright(1, t2);
+ while(t2->n[0] & 1 == 0) {
+ gshiftright(1, t2);
+ ++s;
+ }
+ /* Now, p-1 = 2^s * t2. */
+ for(j = 0; j < MILLER_RABIN_DEPTH; j++) {
+ itog(pr[j+1], t3);
+ powermodg(t3, t2, p);
+ ct = 1;
+ if(isone(t3)) continue;
+ if(gcompg(t3, t1) == 0) continue;
+ while((ct < s) && (gcompg(t1, t3) != 0)) {
+ squareg(t3); modg(p, t3);
+ if(isone(t3)) {
+ goto composite;
+ }
+ ++ct;
+ }
+ if(gcompg(t1, t3) != 0) goto composite;
+ }
+ goto prime;
+
+composite:
+ pushg(3); return(0);
+prime: pushg(3); return(1);
+}
+
+int
+jacobi_symbol(giant a, giant n)
+/* Standard Jacobi symbol (a/n). Parameter n must be odd, positive. */
+{ int t = 1, u;
+ giant t5 = popg(), t6 = popg(), t7 = popg();
+
+ gtog(a, t5); modg(n, t5);
+ gtog(n, t6);
+ while(!isZero(t5)) {
+ u = (t6->n[0]) & 7;
+ while((t5->n[0] & 1) == 0) {
+ gshiftright(1, t5);
+ if((u==3) || (u==5)) t = -t;
+ }
+ gtog(t5, t7); gtog(t6, t5); gtog(t7, t6);
+ u = (t6->n[0]) & 3;
+ if(((t5->n[0] & 3) == 3) && ((u & 3) == 3)) t = -t;
+ modg(t6, t5);
+ }
+ if(isone(t6)) {
+ pushg(3);
+ return(t);
+ }
+ pushg(3);
+ return(0);
+}
+
+int
+pseudoq(giant a, giant p)
+/* Query whether a^(p-1) = 1 (mod p). */
+{
+ int x;
+ giant t1 = popg(), t2 = popg();
+
+ gtog(p, t1); itog(1, t2); subg(t2, t1);
+ gtog(a, t2);
+ powermodg(t2, t1, p);
+ x = isone(t2);
+ pushg(2);
+ return(x);
+}
+
+int
+pseudointq(int a, giant p)
+/* Query whether a^(p-1) = 1 (mod p). */
+{
+ int x;
+ giant t4 = popg();
+
+ itog(a, t4);
+ x = pseudoq(t4, p);
+ pushg(1);
+ return(x);
+}
+
+
+void
+powFp2(giant a, giant b, giant w2, giant n, giant p)
+/* Perform powering in the field F_p^2:
+ a + b w := (a + b w)^n (mod p), where parameter w2 is a quadratic
+ nonresidue (formally equal to w^2).
+ */
+{ int j;
+ giant t6 = popg(), t7 = popg(), t8 = popg(), t9 = popg();
+
+ if(isZero(n)) {
+ itog(1,a);
+ itog(0,b);
+ pushg(4);
+ return;
+ }
+ gtog(a, t8); gtog(b, t9);
+ for(j = bitlen(n)-2; j >= 0; j--) {
+ gtog(b, t6);
+ mulg(a, b); addg(b,b); modg(p, b); /* b := 2 a b. */
+ squareg(t6); modg(p, t6);
+ mulg(w2, t6); modg(p, t6);
+ squareg(a); addg(t6, a); modg(p, a); /* a := a^2 + b^2 w2. */
+ if(bitval(n, j)) {
+ gtog(b, t6); mulg(t8, b); modg(p, b);
+ gtog(a, t7); mulg(t9, a); addg(a, b); modg(p, b);
+ mulg(t9, t6); modg(p, t6); mulg(w2, t6); modg(p, t6);
+ mulg(t8, a); addg(t6, a); modg(p, a);
+ }
+ }
+ pushg(4);
+ return;
+}
+
+int
+sqrtmod(giant p, giant x)
+/* If Sqrt[x] (mod p) exists, function returns 1, else 0.
+ In either case x is modified, but if 1 is returned,
+ x:= Sqrt[x] (mod p).
+ */
+{ giant t0 = popg(), t1 = popg(), t2 = popg(), t3 = popg(),
+ t4 = popg();
+
+ modg(p, x); /* Justify the argument. */
+ gtog(x, t0); /* Store x for eventual validity check on square root. */
+ if((p->n[0] & 3) == 3) { /* The case p = 3 (mod 4). */
+ gtog(p, t1);
+ iaddg(1, t1); gshiftright(2, t1);
+ powermodg(x, t1, p);
+ goto resolve;
+ }
+/* Next, handle case p = 5 (mod 8). */
+ if((p->n[0] & 7) == 5) {
+ gtog(p, t1); itog(1, t2);
+ subg(t2, t1); gshiftright(2, t1);
+ gtog(x, t2);
+ powermodg(t2, t1, p); /* t2 := x^((p-1)/4) % p. */
+ iaddg(1, t1);
+ gshiftright(1, t1); /* t1 := (p+3)/8. */
+ if(isone(t2)) {
+ powermodg(x, t1, p); /* x^((p+3)/8) is root. */
+ goto resolve;
+ } else {
+ itog(1, t2); subg(t2, t1); /* t1 := (p-5)/8. */
+ gshiftleft(2,x);
+ powermodg(x, t1, p);
+ mulg(t0, x); addg(x, x); modg(p, x); /* 2x (4x)^((p-5)/8. */
+ goto resolve;
+ }
+ }
+
+/* Next, handle tougher case: p = 1 (mod 8). */
+ itog(2, t1);
+ while(1) { /* Find appropriate nonresidue. */
+ gtog(t1, t2);
+ squareg(t2); subg(x, t2); modg(p, t2);
+ if(jacobi_symbol(t2, p) == -1) break;
+ iaddg(1, t1);
+ } /* t2 is now w^2 in F_p^2. */
+ itog(1, t3);
+ gtog(p, t4); iaddg(1, t4); gshiftright(1, t4);
+ powFp2(t1, t3, t2, t4, p);
+ gtog(t1, x);
+
+resolve:
+ gtog(x,t1); squareg(t1); modg(p, t1);
+ if(gcompg(t0, t1) == 0) {
+ pushg(5);
+ return(1); /* Success. */
+ }
+ pushg(5);
+ return(0); /* No square root. */
+}
+
+void
+sqrtg(giant n)
+/* n:= Floor[Sqrt[n]]. */
+{ giant t5 = popg(), t6 = popg();
+
+ itog(1, t5); gshiftleft(1 + bitlen(n)/2, t5);
+ while(1) {
+ gtog(n, t6);
+ divg(t5, t6);
+ addg(t5, t6); gshiftright(1, t6);
+ if(gcompg(t6, t5) >= 0) break;
+ gtog(t6, t5);
+ }
+ gtog(t5, n);
+ pushg(2);
+}
+
+int
+cornacchia4(giant n, int d, giant u, giant v)
+/* Seek representation 4n = u^2 + |d| v^2,
+ for (negative) discriminant d and n > |D|/4.
+ Parameter u := 0 and 0 is returned, if no representation is found;
+ else 1 is returned and u, v properly set.
+ */
+{ int r = n->n[0] & 7, sym;
+ giant t1 = popg(), t2 = popg(), t3 = popg(), t4 = popg();
+
+ itog(d, t1);
+ if((n->n[0]) & 7 == 1) { /* n = 1 (mod 8). */
+ sym = jacobi_symbol(t1,n);
+ if(sym != 1) {
+ itog(0,u);
+ pushg(4);
+ return(0);
+ }
+ gtog(t1, t2);
+ sqrtmod(n, t2); /* t2 := Sqrt[d] (mod n). */
+ } else { /* Avoid separate Jacobi/Legendre test. */
+ gtog(t1, t2);
+ if(sqrtmod(n, t2) == 0) {
+ itog(0, u);
+ pushg(4);
+ return(0);
+ }
+ }
+/* t2 is now a valid square root of d (mod n). */
+ gtog(t2, t3);
+ subg(t1, t3); /* t3 := t2 - d. */
+ if((t3->n[0] & 1) == 1) {
+ negg(t2);
+ addg(n, t2);
+ }
+ gtog(n, t3); addg(t3, t3); /* t3 := 2n. */
+ gtog(n, t4); gshiftleft(2, t4); sqrtg(t4); /* t4 = [Sqrt[4 p]]. */
+ while(gcompg(t2, t4) > 0) {
+ gtog(t3, t1);
+ gtog(t2, t3);
+ gtog(t1, t2);
+ modg(t3, t2);
+ }
+ gtog(n, t4); gshiftleft(2, t4);
+ gtog(t2, t3); squareg(t3);
+ subg(t3, t4); /* t4 := 4n - t2^2. */
+ gtog(t4, t3);
+ itog(d, t1); absg(t1);
+ modg(t1, t3);
+ if(!isZero(t3)) {
+ itog(0,u);
+ pushg(4);
+ return(0);
+ }
+ divg(t1, t4);
+ gtog(t4, t1);
+ sqrtg(t4); /* t4 := [Sqrt[t4/Abs[d]]]. */
+ gtog(t4, t3);
+ squareg(t3);
+ if(gcompg(t3, t1) != 0) {
+ itog(0, u);
+ pushg(4);
+ return(0);
+ }
+ gtog(t2, u);
+ gtog(t4, v);
+ pushg(4);
+ return(1);
+}
+
+/*
+rep[p_, d_] := Module[{t, x0, a, b, c},
+ If[JacobiSymbol[d,p] != 1, Return[{0,0}]];
+ x0 = sqrt[d, p];
+ If[Mod[x0-d,2] == 1, x0 = p-x0];
+ a = 2p; b = x0; c = sqrtint[4 p];
+ While[b > c, {a,b} = {b, Mod[a,b]}];
+ t = 4p - b^2;
+ If[Mod[t,Abs[d]] !=0, Return[{0,0}]];
+ v = t/Abs[d];
+ u = sqrtint[v];
+ If[u^2 != v, Return[{0,0}]];
+ Return[{b, u}]
+ ];
+*/
+
+
projects / apple / security.git / blobdiff