#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"
#define STACK_COUNT 20
int pr[NUM_PRIMES];
static giant tmp[STACK_COUNT];
static int stack = 0;
static giant popg();
static void pushg();
void
init_tools(int shorts)
{
int j;
for(j = 0; j < STACK_COUNT; j++) {
tmp[j] = newgiant(shorts);
}
make_primes();
}
static giant
popg() {
return(tmp[stack++]);
}
static void
pushg(int n) {
stack -= n;
}
int
prime_literal(
unsigned int p
)
{
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
)
{
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)
{
giant t1 = popg(), t2 = popg(), t3 = popg();
int j, ct, s;
if((p->n[0] & 1) == 0) {
pushg(3); return(0);
}
if(bitlen(p) < 32) {
pushg(3);
return(primeq((unsigned int)gtoi(p)));
}
itog(-1, t1);
addg(p, t1);
gtog(t1, t2);
s = 1;
gshiftright(1, t2);
while(t2->n[0] & 1 == 0) {
gshiftright(1, t2);
++s;
}
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)
{ 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)
{
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)
{
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)
{ 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);
squareg(t6); modg(p, t6);
mulg(w2, t6); modg(p, t6);
squareg(a); addg(t6, a); modg(p, a);
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)
{ giant t0 = popg(), t1 = popg(), t2 = popg(), t3 = popg(),
t4 = popg();
modg(p, x);
gtog(x, t0);
if((p->n[0] & 3) == 3) {
gtog(p, t1);
iaddg(1, t1); gshiftright(2, t1);
powermodg(x, t1, p);
goto resolve;
}
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);
iaddg(1, t1);
gshiftright(1, t1);
if(isone(t2)) {
powermodg(x, t1, p);
goto resolve;
} else {
itog(1, t2); subg(t2, t1);
gshiftleft(2,x);
powermodg(x, t1, p);
mulg(t0, x); addg(x, x); modg(p, x);
goto resolve;
}
}
itog(2, t1);
while(1) {
gtog(t1, t2);
squareg(t2); subg(x, t2); modg(p, t2);
if(jacobi_symbol(t2, p) == -1) break;
iaddg(1, t1);
}
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);
}
pushg(5);
return(0);
}
void
sqrtg(giant 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)
{ 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) {
sym = jacobi_symbol(t1,n);
if(sym != 1) {
itog(0,u);
pushg(4);
return(0);
}
gtog(t1, t2);
sqrtmod(n, t2);
} else {
gtog(t1, t2);
if(sqrtmod(n, t2) == 0) {
itog(0, u);
pushg(4);
return(0);
}
}
gtog(t2, t3);
subg(t1, t3);
if((t3->n[0] & 1) == 1) {
negg(t2);
addg(n, t2);
}
gtog(n, t3); addg(t3, t3);
gtog(n, t4); gshiftleft(2, t4); sqrtg(t4);
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);
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);
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);
}