/**************************************************************
 *
 *	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}]
		];
*/