-/**************************************************************
- *
- * giants.c
- *
- * Library for large-integer arithmetic.
- *
- * The large-gcd implementation is due to J. P. Buhler.
- * Special mod routines use ideas of R. McIntosh.
- * Contributions from G. Woltman, A. Powell acknowledged.
- *
- * Updates:
- * 18 Jul 99 REC Added routine fer_mod(), for use when Fermat
- giant itself is available.
- * 17 Jul 99 REC Fixed sign bug in fermatmod()
- * 17 Apr 99 REC Fixed various comment/line wraps
- * 25 Mar 99 REC G. Woltman/A. Powell fixes Karat. routines
- * 05 Mar 99 REC Moved invaux, binvaux giants to stack
- * 05 Mar 99 REC Moved gread/gwrite giants to stack
- * 05 Mar 99 REC No static on cur_den, cur_recip (A. Powell)
- * 28 Feb 99 REC Error detection added atop newgiant().
- * 27 Feb 99 REC Reduced wasted work in addsignal().
- * 27 Feb 99 REC Reduced wasted work in FFTmulg().
- * 19 Feb 99 REC Generalized iaddg() per R. Mcintosh.
- * 2 Feb 99 REC Fixed comments.
- * 6 Dec 98 REC Fixed yet another Karatsuba glitch.
- * 1 Dec 98 REC Fixed errant case of addg().
- * 28 Nov 98 REC Installed A. Powell's (correct) variant of
- Karatsuba multiply.
- * 15 May 98 REC Modified gwrite() to handle huge integers.
- * 13 May 98 REC Changed to static stack declarations
- * 11 May 98 REC Installed Karatsuba multiply, to handle
- * medregion 'twixt grammar- and FFT-multiply.
- * 1 May 98 JF gshifts now handle bits < 0 correctly.
- * 30 Apr 98 JF 68k assembler code removed,
- * stack giant size now invariant and based
- * on first call of newgiant(),
- * stack memory leaks fixed.
- * 29 Apr 98 JF function prototyes cleaned up,
- * GCD no longer uses personal stack,
- * optimized shifts for bits%16 == 0.
- * 27 Apr 98 JF scratch giants now replaced with stack
- * 20 Apr 98 JF grammarsquareg fixed for asize == 0.
- * scratch giants now static.
- * 29 Jan 98 JF Corrected out-of-range errors in
- * mersennemod and fermatmod.
- * 23 Dec 97 REC Sped up divide routines via split-shift.
- * 18 Nov 97 JF Improved mersennemod, fermatmod.
- * 9 Nov 97 JF Sped up grammarsquareg.
- * 20 May 97 RDW Fixed Win32 compiler warnings.
- * 18 May 97 REC Installed new, fast divide.
- * 17 May 97 REC Reduced memory for FFT multiply.
- * 26 Apr 97 REC Creation.
- *
- * c. 1997,1998 Perfectly Scientific, Inc.
- * All Rights Reserved.
- *
- **************************************************************/
-
-
-/* Include Files. */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "giants.h"
-
-
-/* Compiler options. */
-
-#ifdef _WIN32
-#pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
-#endif
-
-
-/* Global variables. */
-
-int error = 0;
-int mulmode = AUTO_MUL;
-int cur_prec = 0;
-int cur_shift = 0;
-static int cur_stack_size = 0;
-static int cur_stack_elem = 0;
-static int stack_glen = 0;
-static giant *stack;
-giant cur_den = NULL,
- cur_recip = NULL;
-int current_max_size = 0,
- cur_run = 0;
-double * sinCos=NULL;
-int checkFFTerror = 0;
-double maxFFTerror;
-static giant u0=NULL, u1=NULL, v0=NULL, v1=NULL;
-static double *z = NULL,
- *z2 = NULL;
-
-/* stack handling functions. */
-static giant popg(void);
-static void pushg(int);
-
-
-/* Private function prototypes. */
-
-int gerr(void);
-double gfloor(double);
-int radixdiv(int, int, giant);
-void columnwrite(FILE *, short *, char *, short, int);
-
-void normal_addg(giant, giant);
-void normal_subg(giant, giant);
-void reverse_subg(giant, giant);
-int binvaux(giant, giant);
-int invaux(giant, giant);
-int allzeros(int, int, giant);
-void auxmulg(giant a, giant b);
-void karatmulg(giant a, giant b);
-void karatsquareg(giant b);
-void grammarmulg(giant a, giant b);
-void grammarsquareg(giant b);
-
-int lpt(int, int *);
-void addsignal(giant, double *, int);
-void FFTsquareg(giant x);
-void FFTmulg(giant y, giant x);
-void scramble_real();
-void fft_real_to_hermitian(double *z, int n);
-void fftinv_hermitian_to_real(double *z, int n);
-void mul_hermitian(double *a, double *b, int n);
-void square_hermitian(double *b, int n);
-void giant_to_double(giant x, int sizex, double *z, int L);
-void gswap(giant *, giant *);
-void onestep(giant, giant, gmatrix);
-void mulvM(gmatrix, giant, giant);
-void mulmM(gmatrix, gmatrix);
-void writeM(gmatrix);
-static void punch(giant, gmatrix);
-static void dotproduct(giant, giant, giant, giant);
-void fix(giant *, giant *);
-void hgcd(int, giant, giant, gmatrix);
-void shgcd(int, int, gmatrix);
-
-
-
-/**************************************************************
- *
- * Functions
- *
- **************************************************************/
-
-
-/**************************************************************
- *
- * Initialization and utility functions
- *
- **************************************************************/
-
-double
-gfloor(
- double f
-)
-{
- return floor(f);
-}
-
-
-void
-init_sinCos(
- int n
-)
-{
- int j;
- double e = TWOPI/n;
-
- if (n<=cur_run)
- return;
- cur_run = n;
- if (sinCos)
- free(sinCos);
- sinCos = (double *)malloc(sizeof(double)*(1+(n>>2)));
- for (j=0;j<=(n>>2);j++)
- {
- sinCos[j] = sin(e*j);
- }
-}
-
-
-double
-s_sin(
- int n
-)
-{
- int seg = n/(cur_run>>2);
-
- switch (seg)
- {
- case 0: return(sinCos[n]);
- case 1: return(sinCos[(cur_run>>1)-n]);
- case 2: return(-sinCos[n-(cur_run>>1)]);
- case 3: return(-sinCos[cur_run-n]);
- }
- return 0;
-}
-
-
-double
-s_cos(
- int n
-)
-{
- int quart = (cur_run>>2);
-
- if (n < quart)
- return(s_sin(n+quart));
- return(-s_sin(n-quart));
-}
-
-
-int
-gerr(void)
-{
- return(error);
-}
-
-
-giant
-popg (
- void
-)
-{
- int i;
-
- if (current_max_size <= 0) current_max_size = MAX_SHORTS;
-
- if (cur_stack_size == 0) {
-/* Initialize the stack if we're just starting out.
- * Note that all stack giants will be whatever current_max_size is
- * when newgiant() is first called. */
- cur_stack_size = STACK_GROW;
- stack = (giant *) malloc (cur_stack_size * sizeof(giant));
- for(i = 0; i < STACK_GROW; i++)
- stack[i] = NULL;
- if (stack_glen == 0) stack_glen = current_max_size;
- } else if (cur_stack_elem >= cur_stack_size) {
-/* Expand the stack if we need to. */
- i = cur_stack_size;
- cur_stack_size += STACK_GROW;
- stack = (giant *) realloc (stack,cur_stack_size * sizeof(giant));
- for (; i < cur_stack_size; i++)
- stack[i] = NULL;
- } else if (cur_stack_elem < cur_stack_size - 2*STACK_GROW) {
-/* Prune the stack if it's too big. Disabled, so the stack can only expand */
- /* cur_stack_size -= STACK_GROW;
- for (i = cur_stack_size - STACK_GROW; i < cur_stack_size; i++)
- free(stack[i]);
- stack = (giant *) realloc (stack,cur_stack_size * sizeof(giant)); */
- }
-
-/* Malloc our giant. */
- if (stack[cur_stack_elem] == NULL)
- stack[cur_stack_elem] = malloc(stack_glen*sizeof(short)+sizeof(int));
- stack[cur_stack_elem]->sign = 0;
-
- return(stack[cur_stack_elem++]);
-}
-
-
-void
-pushg (
- int a
-)
-{
- if (a < 0) return;
- cur_stack_elem -= a;
- if (cur_stack_elem < 0) cur_stack_elem = 0;
-}
-
-
-giant
-newgiant(
- int numshorts
-)
-{
- int size;
- giant thegiant;
-
- if (numshorts > MAX_SHORTS) {
- fprintf(stderr, "Requested giant too big.\n");
- fflush(stderr);
- }
- if (numshorts<=0)
- numshorts = MAX_SHORTS;
- size = numshorts*sizeof(short)+sizeof(int);
- thegiant = (giant)malloc(size);
- thegiant->sign = 0;
-
- if (newmin(2*numshorts,MAX_SHORTS) > current_max_size)
- current_max_size = newmin(2*numshorts,MAX_SHORTS);
-
-/* If newgiant() is being called for the first time, set the
- * size of the stack giants. */
- if (stack_glen == 0) stack_glen = current_max_size;
-
- return(thegiant);
-}
-
-
-gmatrix
-newgmatrix(
- void
-)
-/* Allocates space for a gmatrix struct, but does not
- * allocate space for the giants. */
-{
- return((gmatrix) malloc (4*sizeof(giant)));
-}
-
-int
-bitlen(
- giant n
-)
-{
- int b = 16, c = 1<<15, w;
-
- if (isZero(n))
- return(0);
- w = n->n[abs(n->sign) - 1];
- while ((w&c) == 0)
- {
- b--;
- c >>= 1;
- }
- return (16*(abs(n->sign)-1) + b);
-}
-
-
-int
-bitval(
- giant n,
- int pos
-)
-{
- int i = abs(pos)>>4, c = 1 << (pos&15);
-
- return ((n->n[i]) & c);
-}
-
-
-int
-isone(
- giant g
-)
-{
- return((g->sign==1)&&(g->n[0]==1));
-}
-
-
-int isZero(
- giant thegiant
-)
-/* Returns TR if thegiant == 0. */
-{
- register int count;
- int length = abs(thegiant->sign);
- register unsigned short * numpointer = thegiant->n;
-
- if (length)
- {
- for(count = 0; count<length; ++count,++numpointer)
- {
- if (*numpointer != 0 )
- return(FA);
- }
- }
- return(TR);
-}
-
-
-void
-gtog(
- giant srcgiant,
- giant destgiant
-)
-/* destgiant becomes equal to srcgiant. */
-{
- int numbytes = sizeof(int) + abs(srcgiant->sign)*sizeof(short);
-
- memcpy((char *)destgiant,(char *)srcgiant,numbytes);
-}
-
-
-void
-itog(
- int i,
- giant g
-)
-/* The giant g becomes set to the integer value i. */
-{
- unsigned int j = abs(i);
-
- if (i==0)
- {
- g->sign = 0;
- g->n[0] = 0;
- return;
- }
- g->n[0] = (unsigned short)(j & 0xFFFF);
- j >>= 16;
- if (j)
- {
- g->n[1] = (unsigned short)j;
- g->sign = 2;
- }
- else
- {
- g->sign = 1;
- }
- if (i<0)
- g->sign = -(g->sign);
-}
-
-
-signed int
-gtoi(
- giant x
-)
-/* Calculate the value of an int-sized giant NOT exceeding 31 bits. */
-{
- register int size = abs(x->sign);
- register int sign = (x->sign < 0) ? -1 : 1;
-
- switch(size)
- {
- case 0:
- break;
- case 1:
- return sign * x->n[0];
- case 2:
- return sign * (x->n[0]+((x->n[1])<<16));
- default:
- fprintf(stderr,"Giant too large for gtoi\n");
- break;
- }
- return 0;
-}
-
-
-int
-gsign(
- giant g
-)
-/* Returns the sign of g. */
-{
- if (isZero(g))
- return(0);
- if (g->sign >0)
- return(1);
- return(-1);
-}
-
-
-#if 0
-int gcompg(a,b)
-/* Returns -1,0,1 if a<b, a=b, a>b, respectively. */
- giant a,b;
-{
- int size = abs(a->sign);
-
- if(isZero(a)) size = 0;
- if (size == 0) {
- if (isZero(b)) return(0); else return(-gsign(b));
- }
-
- if (b->sign == 0) return(gsign(a));
- if (gsign(a)!=gsign(b)) return(gsign(a));
- if (size>abs(b->sign)) return(gsign(a));
- if (size<abs(b->sign)) return(-gsign(a));
-
- do {
- --size;
- if (a->n[size] > b->n[size]) return(gsign(a));
- if (a->n[size] < b->n[size]) return(-gsign(a));
- } while(size>0);
-
- return(0);
-}
-#else
-
-int
-gcompg(
- giant a,
- giant b
-)
-/* Returns -1,0,1 if a<b, a=b, a>b, respectively. */
-{
- int sa = a->sign, j, sb = b->sign, va, vb, sgn;
-
- if(sa > sb)
- return(1);
- if(sa < sb)
- return(-1);
- if(sa < 0)
- {
- sa = -sa; /* Take absolute value of sa. */
- sgn = -1;
- }
- else
- {
- sgn = 1;
- }
- for(j = sa-1; j >= 0; j--)
- {
- va = a->n[j];
- vb = b->n[j];
- if (va > vb)
- return(sgn);
- if (va < vb)
- return(-sgn);
- }
- return(0);
-}
-#endif
-
-
-void
-setmulmode(
- int mode
-)
-{
- mulmode = mode;
-}
-
-
-/**************************************************************
- *
- * Private I/O Functions
- *
- **************************************************************/
-
-
-int
-radixdiv(
- int base,
- int divisor,
- giant thegiant
-)
-/* Divides giant of arbitrary base by divisor.
- * Returns remainder. Used by idivg and gread. */
-{
- int first = TR;
- int finalsize = abs(thegiant->sign);
- int j = finalsize-1;
- unsigned short *digitpointer=&thegiant->n[j];
- unsigned int num,rem=0;
-
- if (divisor == 0)
- {
- error = DIVIDEBYZERO;
- exit(error);
- }
-
- while (j>=0)
- {
- num=rem*base + *digitpointer;
- *digitpointer = (unsigned short)(num/divisor);
- if (first)
- {
- if (*digitpointer == 0)
- --finalsize;
- else
- first = FA;
- }
- rem = num % divisor;
- --digitpointer;
- --j;
- }
-
- if ((divisor<0) ^ (thegiant->sign<0))
- finalsize=-finalsize;
- thegiant->sign=finalsize;
- return(rem);
-}
-
-
-void
-columnwrite(
- FILE *filepointer,
- short *column,
- char *format,
- short arg,
- int newlines
-)
-/* Used by gwriteln. */
-{
- char outstring[10];
- short i;
-
- sprintf(outstring,format,arg);
- for (i=0; outstring[i]!=0; ++i)
- {
- if (newlines)
- {
- if (*column >= COLUMNWIDTH)
- {
- fputs("\\\n",filepointer);
- *column = 0;
- }
- }
- fputc(outstring[i],filepointer);
- ++*column;
- }
-}
-
-
-void
-gwrite(
- giant thegiant,
- FILE *filepointer,
- int newlines
-)
-/* Outputs thegiant to filepointer. Output is terminated by a newline. */
-{
- short column;
- unsigned int i;
- unsigned short *numpointer;
- giant garbagegiant, basetengrand;
-
- basetengrand = popg();
- garbagegiant = popg();
-
- if (isZero(thegiant))
- {
- fputs("0",filepointer);
- }
- else
- {
- numpointer = basetengrand->n;
- gtog(thegiant,garbagegiant);
-
- basetengrand->sign = 0;
- do
- {
- *numpointer = (unsigned short)idivg(10000,garbagegiant);
- ++numpointer;
- if (++basetengrand->sign >= current_max_size)
- {
- error = OVFLOW;
- exit(error);
- }
- } while (!isZero(garbagegiant));
-
- if (!error)
- {
- i = basetengrand->sign-1;
- column = 0;
- if (thegiant->sign<0 && basetengrand->n[i]!=0)
- columnwrite(filepointer,&column,"-",0, newlines);
- columnwrite(filepointer,&column,"%d",basetengrand->n[i],newlines);
- for( ; i>0; )
- {
- --i;
- columnwrite(filepointer,&column,"%04d",basetengrand->n[i],newlines);
- }
- }
- }
- pushg(2);
-}
-
-
-void
-gwriteln(
- giant theg,
- FILE *filepointer
-)
-{
- gwrite(theg, filepointer, 1);
- fputc('\n',filepointer);
-}
-
-
-void
-gread(
- giant theg,
- FILE *filepointer
-)
-{
- char currentchar;
- int isneg,size,backslash=FA,numdigits=0;
- unsigned short *numpointer;
- giant basetenthousand;
- static char *inbuf = NULL;
-
- basetenthousand = popg();
- if (inbuf == NULL)
- inbuf = (char*)malloc(MAX_DIGITS);
-
- currentchar = (char)fgetc(filepointer);
- if (currentchar=='-')
- {
- isneg=TR;
- }
- else
- {
- isneg=FA;
- if (currentchar!='+')
- ungetc(currentchar,filepointer);
- }
-
- do
- {
- currentchar = (char)fgetc(filepointer);
- if ((currentchar>='0') && (currentchar<='9'))
- {
- inbuf[numdigits]=currentchar;
- if(++numdigits==MAX_DIGITS)
- break;
- backslash=FA;
- }
- else
- {
- if (currentchar=='\\')
- backslash=TR;
- }
- } while(((currentchar!=' ') && (currentchar!='\n') &&
- (currentchar!='\t')) || (backslash) );
- if (numdigits)
- {
- size = 0;
- do
- {
- inbuf[numdigits] = 0;
- numdigits-=4;
- if (numdigits<0)
- numdigits=0;
- basetenthousand->n[size] = (unsigned short)strtol(&inbuf[numdigits],NULL,10);
- ++size;
- } while (numdigits>0);
-
- basetenthousand->sign = size;
- theg->sign = 0;
- numpointer = theg->n;
- do
- {
- *numpointer = (unsigned short)
- radixdiv(10000,1<<(8*sizeof(short)),basetenthousand);
- ++numpointer;
- if (++theg->sign >= current_max_size)
- {
- error = OVFLOW;
- exit(error);
- }
- } while (!isZero(basetenthousand));
-
- if (isneg)
- theg->sign = -theg->sign;
- }
- pushg(1);
-}
-
-
-
-/**************************************************************
- *
- * Private Math Functions
- *
- **************************************************************/
-
-
-void
-negg(
- giant g
-)
-/* g becomes -g. */
-{
- g->sign = -g->sign;
-}
-
-
-void
-absg(
- giant g
-)
-{
- /* g becomes the absolute value of g. */
- if (g->sign <0)
- g->sign=-g->sign;
-}
-
-
-void
-iaddg(
- int i,
- giant g
-)
-/* Giant g becomes g + (int)i. */
-{
- int w,j=0,carry = 0, size = abs(g->sign);
- giant tmp;
-
- if (isZero(g))
- {
- itog(i,g);
- }
- else if(g->sign < 0) {
- tmp = popg();
- itog(i, tmp);
- addg(tmp, g);
- pushg(1);
- return;
- }
- else
- {
- w = g->n[0]+i;
- do
- {
- g->n[j] = (unsigned short) (w & 65535L);
- carry = w >> 16;
- w = g->n[++j]+carry;
- } while ((carry!=0) && (j<size));
- }
- if (carry)
- {
- ++g->sign;
- g->n[size] = (unsigned short)carry;
- }
-}
-
-
-/* New subtract routines.
- The basic subtract "subg()" uses the following logic table:
-
- a b if(b > a) if(a > b)
-
- + + b := b - a b := -(a - b)
- - + b := b + (-a) N.A.
- + - N.A. b := -((-b) + a)
- - - b := (-a) - (-b) b := -((-b) - (-a))
-
- The basic addition routine "addg()" uses:
-
- a b if(b > -a) if(-a > b)
-
- + + b := b + a N.A.
- - + b := b - (-a) b := -((-a) - b)
- + - b := a - (-b) b := -((-b) - a)
- - - N.A. b := -((-b) + (-a))
-
- In this way, internal routines "normal_addg," "normal_subg,"
- and "reverse_subg;" each of which assumes non-negative
- operands and a non-negative result, are now used for greater
- efficiency.
- */
-
-void
-normal_addg(
- giant a,
- giant b
-)
-/* b := a + b, both a,b assumed non-negative. */
-{
- int carry = 0;
- int asize = a->sign, bsize = b->sign;
- long k;
- int j=0;
- unsigned short *aptr = a->n, *bptr = b->n;
-
- if (asize < bsize)
- {
- for (j=0; j<asize; j++)
- {
- k = *aptr++ + *bptr + carry;
- carry = 0;
- if (k >= 65536L)
- {
- k -= 65536L;
- ++carry;
- }
- *bptr++ = (unsigned short)k;
- }
- for (j=asize; j<bsize; j++)
- {
- k = *bptr + carry;
- carry = 0;
- if (k >= 65536L)
- {
- k -= 65536L;
- ++carry;
- }
- *bptr++ = (unsigned short)k;
- }
- }
- else
- {
- for (j=0; j<bsize; j++)
- {
- k = *aptr++ + *bptr + carry;
- carry = 0;
- if (k >= 65536L)
- {
- k -= 65536L;
- ++carry;
- }
- *bptr++ = (unsigned short)k;
- }
- for (j=bsize; j<asize; j++)
- {
- k = *aptr++ + carry;
- carry = 0;
- if (k >= 65536L)
- {
- k -= 65536L;
- ++carry;
- }
- *bptr++ = (unsigned short)k;
- }
- }
- if (carry)
- {
- *bptr = 1; ++j;
- }
- b->sign = j;
-}
-
-
-void
-normal_subg(
- giant a,
- giant b
-)
-/* b := b - a; requires b, a non-negative and b >= a. */
-{
- int j, size = b->sign;
- unsigned int k;
-
- if (a->sign == 0)
- return;
-
- k = 0;
- for (j=0; j<a->sign; ++j)
- {
- k += 0xffff - a->n[j] + b->n[j];
- b->n[j] = (unsigned short)(k & 0xffff);
- k >>= 16;
- }
- for (j=a->sign; j<size; ++j)
- {
- k += 0xffff + b->n[j];
- b->n[j] = (unsigned short)(k & 0xffff);
- k >>= 16;
- }
-
- if (b->n[0] == 0xffff)
- iaddg(1,b);
- else
- ++b->n[0];
-
- while ((size-- > 0) && (b->n[size]==0));
-
- b->sign = (b->n[size]==0) ? 0 : size+1;
-}
-
-
-void
-reverse_subg(
- giant a,
- giant b
-)
-/* b := a - b; requires b, a non-negative and a >= b. */
-{
- int j, size = a->sign;
- unsigned int k;
-
- k = 0;
- for (j=0; j<b->sign; ++j)
- {
- k += 0xffff - b->n[j] + a->n[j];
- b->n[j] = (unsigned short)(k & 0xffff);
- k >>= 16;
- }
- for (j=b->sign; j<size; ++j)
- {
- k += 0xffff + a->n[j];
- b->n[j] = (unsigned short)(k & 0xffff);
- k >>= 16;
- }
-
- b->sign = size; /* REC, 21 Apr 1996. */
- if (b->n[0] == 0xffff)
- iaddg(1,b);
- else
- ++b->n[0];
-
- while (!b->n[--size]);
-
- b->sign = size+1;
-}
-
-void
-addg(
- giant a,
- giant b
-)
-/* b := b + a, any signs any result. */
-{
- int asgn = a->sign, bsgn = b->sign;
-
- if (asgn == 0)
- return;
- if (bsgn == 0)
- {
- gtog(a,b);
- return;
- }
- if ((asgn < 0) == (bsgn < 0))
- {
- if (bsgn > 0)
- {
- normal_addg(a,b);
- return;
- }
- absg(b);
- if(a != b) absg(a);
- normal_addg(a,b);
- negg(b);
- if(a != b) negg(a);
- return;
- }
- if(bsgn > 0)
- {
- negg(a);
- if (gcompg(b,a) >= 0)
- {
- normal_subg(a,b);
- negg(a);
- return;
- }
- reverse_subg(a,b);
- negg(a);
- negg(b);
- return;
- }
- negg(b);
- if(gcompg(b,a) < 0)
- {
- reverse_subg(a,b);
- return;
- }
- normal_subg(a,b);
- negg(b);
- return;
-}
-
-void
-subg(
- giant a,
- giant b
-)
-/* b := b - a, any signs, any result. */
-{
- int asgn = a->sign, bsgn = b->sign;
-
- if (asgn == 0)
- return;
- if (bsgn == 0)
- {
- gtog(a,b);
- negg(b);
- return;
- }
- if ((asgn < 0) != (bsgn < 0))
- {
- if (bsgn > 0)
- {
- negg(a);
- normal_addg(a,b);
- negg(a);
- return;
- }
- negg(b);
- normal_addg(a,b);
- negg(b);
- return;
- }
- if (bsgn > 0)
- {
- if (gcompg(b,a) >= 0)
- {
- normal_subg(a,b);
- return;
- }
- reverse_subg(a,b);
- negg(b);
- return;
- }
- negg(a);
- negg(b);
- if (gcompg(b,a) >= 0)
- {
- normal_subg(a,b);
- negg(a);
- negg(b);
- return;
- }
- reverse_subg(a,b);
- negg(a);
- return;
-}
-
-
-int
-numtrailzeros(
- giant g
-)
-/* Returns the number of trailing zero bits in g. */
-{
- register int numshorts = abs(g->sign), j, bcount=0;
- register unsigned short gshort, c;
-
- for (j=0;j<numshorts;j++)
- {
- gshort = g->n[j];
- c = 1;
- for (bcount=0;bcount<16; bcount++)
- {
- if (c & gshort)
- break;
- c <<= 1;
- }
- if (bcount<16)
- break;
- }
- return(bcount + 16*j);
-}
-
-
-void
-bdivg(
- giant v,
- giant u
-)
-/* u becomes greatest power of two not exceeding u/v. */
-{
- int diff = bitlen(u) - bitlen(v);
- giant scratch7;
-
- if (diff<0)
- {
- itog(0,u);
- return;
- }
- scratch7 = popg();
- gtog(v, scratch7);
- gshiftleft(diff,scratch7);
- if (gcompg(u,scratch7) < 0)
- diff--;
- if (diff<0)
- {
- itog(0,u);
- pushg(1);
- return;
- }
- itog(1,u);
- gshiftleft(diff,u);
-
- pushg(1);
-}
-
-
-int
-binvaux(
- giant p,
- giant x
-)
-/* Binary inverse method. Returns zero if no inverse exists,
- * in which case x becomes GCD(x,p). */
-{
-
- giant scratch7, u0, u1, v0, v1;
-
- if (isone(x))
- return(1);
- u0 = popg();
- u1 = popg();
- v0 = popg();
- v1 = popg();
- itog(1, v0);
- gtog(x, v1);
- itog(0,x);
- gtog(p, u1);
-
- scratch7 = popg();
- while(!isZero(v1))
- {
- gtog(u1, u0);
- bdivg(v1, u0);
- gtog(x, scratch7);
- gtog(v0, x);
- mulg(u0, v0);
- subg(v0,scratch7);
- gtog(scratch7, v0);
-
- gtog(u1, scratch7);
- gtog(v1, u1);
- mulg(u0, v1);
- subg(v1,scratch7);
- gtog(scratch7, v1);
- }
-
- pushg(1);
-
- if (!isone(u1))
- {
- gtog(u1,x);
- if(x->sign<0) addg(p, x);
- pushg(4);
- return(0);
- }
- if(x->sign<0)
- addg(p, x);
- pushg(4);
- return(1);
-}
-
-
-int
-binvg(
- giant p,
- giant x
-)
-{
- modg(p, x);
- return(binvaux(p,x));
-}
-
-
-int
-invg(
- giant p,
- giant x
-)
-{
- modg(p, x);
- return(invaux(p,x));
-}
-
-int
-invaux(
- giant p,
- giant x
-)
-/* Returns zero if no inverse exists, in which case x becomes
- * GCD(x,p). */
-{
-
- giant scratch7, u0, u1, v0, v1;
-
- if ((x->sign==1)&&(x->n[0]==1))
- return(1);
-
- u0 = popg();
- u1 = popg();
- v0 = popg();
- v1 = popg();
-
- itog(1,u1);
- gtog(p, v0);
- gtog(x, v1);
- itog(0,x);
-
- scratch7 = popg();
- while (!isZero(v1))
- {
- gtog(v0, u0);
- divg(v1, u0);
- gtog(u0, scratch7);
- mulg(v1, scratch7);
- subg(v0, scratch7);
- negg(scratch7);
- gtog(v1, v0);
- gtog(scratch7, v1);
- gtog(u1, scratch7);
- mulg(u0, scratch7);
- subg(x, scratch7);
- negg(scratch7);
- gtog(u1,x);
- gtog(scratch7, u1);
- }
- pushg(1);
-
- if ((v0->sign!=1)||(v0->n[0]!=1))
- {
- gtog(v0,x);
- pushg(4);
- return(0);
- }
- if(x->sign<0)
- addg(p, x);
- pushg(4);
- return(1);
-}
-
-
-int
-mersenneinvg(
- int q,
- giant x
-)
-{
- int k;
- giant u0, u1, v1;
-
- u0 = popg();
- u1 = popg();
- v1 = popg();
-
- itog(1, u0);
- itog(0, u1);
- itog(1, v1);
- gshiftleft(q, v1);
- subg(u0, v1);
- mersennemod(q, x);
- while (1)
- {
- k = -1;
- if (isZero(x))
- {
- gtog(v1, x);
- pushg(3);
- return(0);
- }
- while (bitval(x, ++k) == 0);
-
- gshiftright(k, x);
- if (k)
- {
- gshiftleft(q-k, u0);
- mersennemod(q, u0);
- }
- if (isone(x))
- break;
- addg(u1, u0);
- mersennemod(q, u0);
- negg(u1);
- addg(u0, u1);
- mersennemod(q, u1);
- if (!gcompg(v1,x)) {
- pushg(3);
- return(0);
- }
- addg(v1, x);
- negg(v1);
- addg(x, v1);
- mersennemod(q, v1);
- }
- gtog(u0, x);
- mersennemod(q,x);
- pushg(3);
- return(1);
-}
-
-
-void
-cgcdg(
- giant a,
- giant v
-)
-/* Classical Euclid GCD. v becomes gcd(a, v). */
-{
- giant u, r;
-
- v->sign = abs(v->sign);
- if (isZero(a))
- return;
-
- u = popg();
- r = popg();
- gtog(a, u);
- u->sign = abs(u->sign);
- while (!isZero(v))
- {
- gtog(u, r);
- modg(v, r);
- gtog(v, u);
- gtog(r, v);
- }
- gtog(u,v);
- pushg(2);
-}
-
-
-void
-gcdg(
- giant x,
- giant y
-)
-{
- if (bitlen(y)<= GCDLIMIT)
- bgcdg(x,y);
- else
- ggcd(x,y);
-}
-
-void
-bgcdg(
- giant a,
- giant b
-)
-/* Binary form of the gcd. b becomes the gcd of a,b. */
-{
- int k = isZero(b), m = isZero(a);
- giant u, v, t;
-
- if (k || m)
- {
- if (m)
- {
- if (k)
- itog(1,b);
- return;
- }
- if (k)
- {
- if (m)
- itog(1,b);
- else
- gtog(a,b);
- return;
- }
- }
-
- u = popg();
- v = popg();
- t = popg();
-
- /* Now neither a nor b is zero. */
- gtog(a, u);
- u->sign = abs(a->sign);
- gtog(b, v);
- v->sign = abs(b->sign);
- k = numtrailzeros(u);
- m = numtrailzeros(v);
- if (k>m)
- k = m;
- gshiftright(k,u);
- gshiftright(k,v);
- if (u->n[0] & 1)
- {
- gtog(v, t);
- negg(t);
- }
- else
- {
- gtog(u,t);
- }
-
- while (!isZero(t))
- {
- m = numtrailzeros(t);
- gshiftright(m, t);
- if (t->sign > 0)
- {
- gtog(t, u);
- subg(v,t);
- }
- else
- {
- gtog(t, v);
- negg(v);
- addg(u,t);
- }
- }
- gtog(u,b);
- gshiftleft(k, b);
- pushg(3);
-}
-
-
-void
-powerg(
- int m,
- int n,
- giant g
-)
-/* g becomes m^n, NO mod performed. */
-{
- giant scratch2 = popg();
-
- itog(1, g);
- itog(m, scratch2);
- while (n)
- {
- if (n & 1)
- mulg(scratch2, g);
- n >>= 1;
- if (n)
- squareg(scratch2);
- }
- pushg(1);
-}
-
-#if 0
-void
-jtest(
- giant n
-)
-{
- if (n->sign)
- {
- if (n->n[n->sign-1] == 0)
- {
- fprintf(stderr,"%d %d tilt",n->sign, (int)(n->n[n->sign-1]));
- exit(7);
- }
- }
-}
-#endif
-
-
-void
-make_recip(
- giant d,
- giant r
-)
-/* r becomes the steady-state reciprocal
- * 2^(2b)/d, where b = bit-length of d-1. */
-{
- int b;
- giant tmp, tmp2;
-
- if (isZero(d) || (d->sign < 0))
- {
- exit(SIGN);
- }
- tmp = popg();
- tmp2 = popg();
- itog(1, r);
- subg(r, d);
- b = bitlen(d);
- addg(r, d);
- gshiftleft(b, r);
- gtog(r, tmp2);
- while (1)
- {
- gtog(r, tmp);
- squareg(tmp);
- gshiftright(b, tmp);
- mulg(d, tmp);
- gshiftright(b, tmp);
- addg(r, r);
- subg(tmp, r);
- if (gcompg(r, tmp2) <= 0)
- break;
- gtog(r, tmp2);
- }
- itog(1, tmp);
- gshiftleft(2*b, tmp);
- gtog(r, tmp2);
- mulg(d, tmp2);
- subg(tmp2, tmp);
- itog(1, tmp2);
- while (tmp->sign < 0)
- {
- subg(tmp2, r);
- addg(d, tmp);
- }
- pushg(2);
-}
-
-void
-divg_via_recip(
- giant d,
- giant r,
- giant n
-)
-/* n := n/d, where r is the precalculated
- * steady-state reciprocal of d. */
-{
- int s = 2*(bitlen(r)-1), sign = gsign(n);
- giant tmp, tmp2;
-
- if (isZero(d) || (d->sign < 0))
- {
- exit(SIGN);
- }
-
- tmp = popg();
- tmp2 = popg();
-
- n->sign = abs(n->sign);
- itog(0, tmp2);
- while (1)
- {
- gtog(n, tmp);
- mulg(r, tmp);
- gshiftright(s, tmp);
- addg(tmp, tmp2);
- mulg(d, tmp);
- subg(tmp, n);
- if (gcompg(n,d) >= 0)
- {
- subg(d,n);
- iaddg(1, tmp2);
- }
- if (gcompg(n,d) < 0)
- break;
- }
- gtog(tmp2, n);
- n->sign *= sign;
- pushg(2);
-}
-
-#if 1
-void
-modg_via_recip(
- giant d,
- giant r,
- giant n
-)
-/* This is the fastest mod of the present collection.
- * n := n % d, where r is the precalculated
- * steady-state reciprocal of d. */
-
-{
- int s = (bitlen(r)-1), sign = n->sign;
- giant tmp, tmp2;
-
- if (isZero(d) || (d->sign < 0))
- {
- exit(SIGN);
- }
-
- tmp = popg();
- tmp2 = popg();
-
- n->sign = abs(n->sign);
- while (1)
- {
- gtog(n, tmp); gshiftright(s-1, tmp);
- mulg(r, tmp);
- gshiftright(s+1, tmp);
- mulg(d, tmp);
- subg(tmp, n);
- if (gcompg(n,d) >= 0)
- subg(d,n);
- if (gcompg(n,d) < 0)
- break;
- }
- if (sign >= 0)
- goto done;
- if (isZero(n))
- goto done;
- negg(n);
- addg(d,n);
-done:
- pushg(2);
- return;
-}
-
-#else
-void
-modg_via_recip(
- giant d,
- giant r,
- giant n
-)
-{
- int s = 2*(bitlen(r)-1), sign = n->sign;
- giant tmp, tmp2;
-
- if (isZero(d) || (d->sign < 0))
- {
- exit(SIGN);
- }
-
- tmp = popg();
- tmp2 = popg()
-
- n->sign = abs(n->sign);
- while (1)
- {
- gtog(n, tmp);
- mulg(r, tmp);
- gshiftright(s, tmp);
- mulg(d, tmp);
- subg(tmp, n);
- if (gcompg(n,d) >= 0)
- subg(d,n);
- if (gcompg(n,d) < 0)
- break;
- }
- if (sign >= 0)
- goto done;
- if (isZero(n))
- goto done;
- negg(n);
- addg(d,n);
-done:
- pushg(2);
- return;
-}
-#endif
-
-void
-modg(
- giant d,
- giant n
-)
-/* n becomes n%d. n is arbitrary, but the denominator d must be positive! */
-{
- if (cur_recip == NULL) {
- cur_recip = newgiant(current_max_size);
- cur_den = newgiant(current_max_size);
- gtog(d, cur_den);
- make_recip(d, cur_recip);
- } else if (gcompg(d, cur_den)) {
- gtog(d, cur_den);
- make_recip(d, cur_recip);
- }
- modg_via_recip(d, cur_recip, n);
-}
-
-
-#if 0
-int
-feemulmod (
- giant a,
- giant b,
- int q,
- int k
-)
-/* a becomes (a*b) (mod 2^q-k) where q % 16 == 0 and k is "small" (0 < k < 65535).
- * Returns 0 if unsuccessful, otherwise 1. */
-{
- giant carry, kk, scratch;
- int i, j;
- int asize = abs(a->sign), bsize = abs(b->sign);
- unsigned short *aptr,*bptr,*destptr;
- unsigned int words;
- int kpower, curk;
-
- if ((q % 16) || (k <= 0) || (k >= 65535)) {
- return (0);
- }
-
- carry = popg();
- kk = popg();
- scratch = popg();
-
- for (i=0; i<asize+bsize; i++) scratch->n[i]=0;
-
- words = q >> 4;
-
- bptr = b->n;
- for (i = 0; i < bsize; i++) {
- mult = *bptr++;
- if (mult) {
- kpower = i/words;
-
- if (kpower >= 1) itog (kpower,kk);
- for (j = 1; j < kpower; k++) smulg(kpower,kk);
-
- itog(0,carry);
-
- aptr = a->n;
- for (j = 0; j < bsize; b++) {
- gtog(kk,scratch);
- smulg(*aptr++,scratch);
- smulg(mult,scratch);
- iaddg(*destptr,scratch);
- addg(carry,scratch);
- *destptr++ = scratch->n[0];
- gshiftright(scratch,16);
- gtog(scratch,carry);
- if (destptr - scratch->n >= words) {
- smulg(k, carry);
- smulg(k, kk);
- destptr -= words;
- }
- }
- }
- }
-
- int i,j,m;
- unsigned int prod,carry=0;
- int asize = abs(a->sign), bsize = abs(b->sign);
- unsigned short *aptr,*bptr,*destptr;
- unsigned short mult;
- int words, excess;
- int temp;
- giant scratch = popg(), scratch2 = popg(), scratch3 = popg();
- short *carryptr = scratch->n;
- int kpower,kpowerlimit, curk;
-
- if ((q % 16) || (k <= 0) || (k >= 65535)) {
- return (0);
- }
-
- scratch
-
- for (i=0; i<asize+bsize; i++) scratch->n[i]=0;
-
- words = q >> 4;
-
- bptr = b->n;
- for (i=0; i<bsize; ++i)
- {
- mult = *bptr++;
- if (mult)
- {
- kpower = i/words;
- aptr = a->n;
- destptr = scratch->n + i;
-
- if (kpower == 0) {
- carry = 0;
- } else if (kpower <= kpowerlimit) {
- carry = 0;
- curk = k;
- for (j = 1; j < kpower; j++) curk *= k;
- } else {
- itog (k,scratch);
- for (j = 1; j < kpower; j++) smulg(k,scratch);
- itog(0,scratch2);
- }
-
- for (j = 0; j < asize; j++) {
- if(kpower == 0) {
- prod = *aptr++ * mult + *destptr + carry;
- *destptr++ = (unsigned short)(prod & 0xFFFF);
- carry = prod >> 16;
- } else if (kpower < kpowerlimit) {
- prod = kcur * *aptr++;
- temp = prod >> 16;
- prod &= 0xFFFF;
- temp *= mult;
- prod *= mult;
- temp += prod >> 16;
- prod &= 0xFFFF;
- prod += *destptr + carry;
- carry = prod >> 16 + temp;
- *destptr++ = (unsigned short)(prod & 0xFFFF);
- } else {
- gtog(scratch,scratch3);
- smulg(*aptr++,scratch3);
- smulg(mult,scratch3);
- iaddg(*destptr,scratch3);
- addg(scratch3,scratch2);
- *destptr++ = scratch2->n[0];
- memmove(scratch2->n,scratch2->n+1,2*(scratch2->size-1));
- scratch2->sign--;
- }
- if (destptr - scratch->n > words) {
- if (kpower == 0) {
- curk = k;
- carry *= k;
- } else if (kpower < kpowerlimit) {
- curk *= k;
- carry *= curk;
- } else if (kpower == kpowerlimit) {
- itog (k,scratch);
- for (j = 1; j < kpower; j++) smulg(k,scratch);
- itog(carry,scratch2);
- smulg(k,scratch2);
- } else {
- smulg(k,scratch);
- smulg(k,scratch2);
- }
- kpower++;
- destptr -= words;
- }
- }
-
- /* Next, deal with the carry term. Needs to be improved to
- handle overflow carry cases. */
- if (kpower <= kpowerlimit) {
- iaddg(carry,scratch);
- } else {
- addg(scratch2,scratch);
- }
- while(scratch->sign > q)
- gtog(scratch,scratch2)
- }
- }
- scratch->sign = destptr - scratch->n;
- if (!carry)
- --(scratch->sign);
- scratch->sign *= gsign(a)*gsign(b);
- gtog(scratch,a);
- pushg(3);
- return (1);
-}
-#endif
-
-int
-idivg(
- int divisor,
- giant theg
-)
-{
- /* theg becomes theg/divisor. Returns remainder. */
- int n;
- int base = 1<<(8*sizeof(short));
-
- n = radixdiv(base,divisor,theg);
- return(n);
-}
-
-
-void
-divg(
- giant d,
- giant n
-)
-/* n becomes n/d. n is arbitrary, but the denominator d must be positive! */
-{
- if (cur_recip == NULL) {
- cur_recip = newgiant(current_max_size);
- cur_den = newgiant(current_max_size);
- gtog(d, cur_den);
- make_recip(d, cur_recip);
- } else if (gcompg(d, cur_den)) {
- gtog(d, cur_den);
- make_recip(d, cur_recip);
- }
- divg_via_recip(d, cur_recip, n);
-}
-
-
-void
-powermod(
- giant x,
- int n,
- giant g
-)
-/* x becomes x^n (mod g). */
-{
- giant scratch2 = popg();
- gtog(x, scratch2);
- itog(1, x);
- while (n)
- {
- if (n & 1)
- {
- mulg(scratch2, x);
- modg(g, x);
- }
- n >>= 1;
- if (n)
- {
- squareg(scratch2);
- modg(g, scratch2);
- }
- }
- pushg(1);
-}
-
-
-void
-powermodg(
- giant x,
- giant n,
- giant g
-)
-/* x becomes x^n (mod g). */
-{
- int len, pos;
- giant scratch2 = popg();
-
- gtog(x, scratch2);
- itog(1, x);
- len = bitlen(n);
- pos = 0;
- while (1)
- {
- if (bitval(n, pos++))
- {
- mulg(scratch2, x);
- modg(g, x);
- }
- if (pos>=len)
- break;
- squareg(scratch2);
- modg(g, scratch2);
- }
- pushg(1);
-}
-
-
-void
-fermatpowermod(
- giant x,
- int n,
- int q
-)
-/* x becomes x^n (mod 2^q+1). */
-{
- giant scratch2 = popg();
-
- gtog(x, scratch2);
- itog(1, x);
- while (n)
- {
- if (n & 1)
- {
- mulg(scratch2, x);
- fermatmod(q, x);
- }
- n >>= 1;
- if (n)
- {
- squareg(scratch2);
- fermatmod(q, scratch2);
- }
- }
- pushg(1);
-}
-
-
-void
-fermatpowermodg(
- giant x,
- giant n,
- int q
-)
-/* x becomes x^n (mod 2^q+1). */
-{
- int len, pos;
- giant scratch2 = popg();
-
- gtog(x, scratch2);
- itog(1, x);
- len = bitlen(n);
- pos = 0;
- while (1)
- {
- if (bitval(n, pos++))
- {
- mulg(scratch2, x);
- fermatmod(q, x);
- }
- if (pos>=len)
- break;
- squareg(scratch2);
- fermatmod(q, scratch2);
- }
- pushg(1);
-}
-
-
-void
-mersennepowermod(
- giant x,
- int n,
- int q
-)
-/* x becomes x^n (mod 2^q-1). */
-{
- giant scratch2 = popg();
-
- gtog(x, scratch2);
- itog(1, x);
- while (n)
- {
- if (n & 1)
- {
- mulg(scratch2, x);
- mersennemod(q, x);
- }
- n >>= 1;
- if (n)
- {
- squareg(scratch2);
- mersennemod(q, scratch2);
- }
- }
- pushg(1);
-}
-
-
-void
-mersennepowermodg(
- giant x,
- giant n,
- int q
-)
-/* x becomes x^n (mod 2^q-1). */
-{
- int len, pos;
- giant scratch2 = popg();
-
- gtog(x, scratch2);
- itog(1, x);
- len = bitlen(n);
- pos = 0;
- while (1)
- {
- if (bitval(n, pos++))
- {
- mulg(scratch2, x);
- mersennemod(q, x);
- }
- if (pos>=len)
- break;
- squareg(scratch2);
- mersennemod(q, scratch2);
- }
- pushg(1);
-}
-
-
-void
-gshiftleft(
- int bits,
- giant g
-)
-/* shift g left bits bits. Equivalent to g = g*2^bits. */
-{
- int rem = bits&15, crem = 16-rem, words = bits>>4;
- int size = abs(g->sign), j, k, sign = gsign(g);
- unsigned short carry, dat;
-
- if (!bits)
- return;
- if (!size)
- return;
- if (bits < 0) {
- gshiftright(-bits,g);
- return;
- }
- if (size+words+1 > current_max_size) {
- error = OVFLOW;
- exit(error);
- }
- if (rem == 0) {
- memmove(g->n + words, g->n, size * sizeof(short));
- for (j = 0; j < words; j++) g->n[j] = 0;
- g->sign += (g->sign < 0)?(-words):(words);
- } else {
- k = size+words;
- carry = 0;
- for (j=size-1; j>=0; j--) {
- dat = g->n[j];
- g->n[k--] = (unsigned short)((dat >> crem) | carry);
- carry = (unsigned short)(dat << rem);
- }
- do {
- g->n[k--] = carry;
- carry = 0;
- } while(k>=0);
-
- k = size+words;
- if (g->n[k] == 0)
- --k;
- g->sign = sign*(k+1);
- }
-}
-
-
-void
-gshiftright(
- int bits,
- giant g
-)
-/* shift g right bits bits. Equivalent to g = g/2^bits. */
-{
- register int j,size=abs(g->sign);
- register unsigned int carry;
- int words = bits >> 4;
- int remain = bits & 15, cremain = (16-remain);
-
- if (bits==0)
- return;
- if (isZero(g))
- return;
- if (bits < 0) {
- gshiftleft(-bits,g);
- return;
- }
- if (words >= size) {
- g->sign = 0;
- return;
- }
- if (remain == 0) {
- memmove(g->n,g->n + words,(size - words) * sizeof(short));
- g->sign += (g->sign < 0)?(words):(-words);
- } else {
- size -= words;
-
- if (size)
- {
- for(j=0;j<size-1;++j)
- {
- carry = g->n[j+words+1] << cremain;
- g->n[j] = (unsigned short)((g->n[j+words] >> remain ) | carry);
- }
- g->n[size-1] = (unsigned short)(g->n[size-1+words] >> remain);
- }
-
- if (g->n[size-1] == 0)
- --size;
-
- if (g->sign > 0)
- g->sign = size;
- else
- g->sign = -size;
- }
-}
-
-
-void
-extractbits(
- int n,
- giant src,
- giant dest
-)
-/* dest becomes lowermost n bits of src. Equivalent to dest = src % 2^n. */
-{
- register int words = n >> 4, numbytes = words*sizeof(short);
- register int bits = n & 15;
-
- if (n<=0)
- return;
- if (words >= abs(src->sign))
- gtog(src,dest);
- else
- {
- memcpy((char *)(dest->n), (char *)(src->n), numbytes);
- if (bits)
- {
- dest->n[words] = (unsigned short)(src->n[words] & ((1<<bits)-1));
- ++words;
- }
- while ((dest->n[words-1] == 0) && (words > 0))
- {
- --words;
- }
- if (src->sign<0)
- dest->sign = -words;
- else
- dest->sign = words;
- }
-}
-
-
-int
-allzeros(
- int shorts,
- int bits,
- giant g
-)
-{
- int i=shorts;
-
- while (i>0)
- {
- if (g->n[--i])
- return(0);
- }
- return((int)(!(g->n[shorts] & ((1<<bits)-1))));
-}
-
-
-void
-fermatnegate(
- int n,
- giant g
-)
-/* negate g. g is mod 2^n+1. */
-{
- register int shorts = n>>4,
- bits = n & 15,
- i = shorts,
- mask = 1<<bits;
- register unsigned carry,temp;
-
- for (temp=(unsigned)shorts; (int)temp>g->sign-1; --temp)
- {
- g->n[temp] = 0;
- }
- if (g->n[shorts] & mask)
- { /* if high bit is set, -g = 1. */
- g->sign = 1;
- g->n[0] = 1;
- return;
- }
- if (allzeros(shorts,bits,g))
- return; /* if g=0, -g = 0. */
-
- while (i>0)
- { --i;
- g->n[i] = (unsigned short)(~(g->n[i+1]));
- }
- g->n[shorts] ^= mask-1;
-
- carry = 2;
- i = 0;
- while (carry)
- {
- temp = g->n[i]+carry;
- g->n[i++] = (unsigned short)(temp & 0xffff);
- carry = temp>>16;
- }
- while(!g->n[shorts])
- {
- --shorts;
- }
- g->sign = shorts+1;
-}
-
-
-void
-mersennemod (
- int n,
- giant g
-)
-/* g := g (mod 2^n - 1) */
-{
- int the_sign, s;
- giant scratch3 = popg(), scratch4 = popg();
-
- if ((the_sign = gsign(g)) < 0) absg(g);
- while (bitlen(g) > n) {
- gtog(g,scratch3);
- gshiftright(n,scratch3);
- addg(scratch3,g);
- gshiftleft(n,scratch3);
- subg(scratch3,g);
- }
- if(!isZero(g)) {
- if ((s = gsign(g)) < 0) absg(g);
- itog(1,scratch3);
- gshiftleft(n,scratch3);
- itog(1,scratch4);
- subg(scratch4,scratch3);
- if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
- if (s < 0) {
- g->sign = -g->sign;
- addg(scratch3,g);
- }
- if (the_sign < 0) {
- g->sign = -g->sign;
- addg(scratch3,g);
- }
- }
- pushg(2);
-}
-
-void
-fermatmod (
- int n,
- giant g
-)
-/* g := g (mod 2^n + 1), */
-{
- int the_sign, s;
- giant scratch3 = popg();
-
- if ((the_sign = gsign(g)) < 0) absg(g);
- while (bitlen(g) > n) {
- gtog(g,scratch3);
- gshiftright(n,scratch3);
- subg(scratch3,g);
- gshiftleft(n,scratch3);
- subg(scratch3,g);
- }
- if((bitlen(g) < n) && (the_sign * (g->sign) >= 0)) goto leave;
- if(!isZero(g)) {
- if ((s = gsign(g)) < 0) absg(g);
- itog(1,scratch3);
- gshiftleft(n,scratch3);
- iaddg(1,scratch3);
- if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
- if (s * the_sign < 0) {
- g->sign = -g->sign;
- addg(scratch3,g);
- }
- }
-leave:
- pushg(1);
-
-}
-
-void
-fer_mod (
- int n,
- giant g,
- giant modulus
-)
-/* Same as fermatmod(), except modulus = 2^n should be passed
-if available (i.e. if already allocated and set). */
-{
- int the_sign, s;
- giant scratch3 = popg();
-
- if ((the_sign = gsign(g)) < 0) absg(g);
- while (bitlen(g) > n) {
- gtog(g,scratch3);
- gshiftright(n,scratch3);
- subg(scratch3,g);
- gshiftleft(n,scratch3);
- subg(scratch3,g);
- }
- if((bitlen(g) < n) && (the_sign * (g->sign) >= 0)) goto leave;
- if(!isZero(g)) {
- if ((s = gsign(g)) < 0) absg(g);
- if(gcompg(g,modulus) >= 0) subg(modulus,g);
- if (s * the_sign < 0) {
- g->sign = -g->sign;
- addg(modulus,g);
- }
- }
-leave:
- pushg(1);
-}
-
-
-void
-smulg(
- unsigned short i,
- giant g
-)
-/* g becomes g * i. */
-{
- unsigned short carry = 0;
- int size = abs(g->sign);
- register int j,k,mul = abs(i);
- unsigned short *digit = g->n;
-
- for (j=0; j<size; ++j)
- {
- k = *digit * mul + carry;
- carry = (unsigned short)(k>>16);
- *digit = (unsigned short)(k & 0xffff);
- ++digit;
- }
- if (carry)
- {
- if (++j >= current_max_size)
- {
- error = OVFLOW;
- exit(error);
- }
- *digit = carry;
- }
-
- if ((g->sign>0) ^ (i>0))
- g->sign = -j;
- else
- g->sign = j;
-}
-
-
-void
-squareg(
- giant b
-)
-/* b becomes b^2. */
-{
- auxmulg(b,b);
-}
-
-
-void
-mulg(
- giant a,
- giant b
-)
-/* b becomes a*b. */
-{
- auxmulg(a,b);
-}
-
-
-void
-auxmulg(
- giant a,
- giant b
-)
-/* Optimized general multiply, b becomes a*b. Modes are:
- * AUTO_MUL: switch according to empirical speed criteria.
- * GRAMMAR_MUL: force grammar-school algorithm.
- * KARAT_MUL: force Karatsuba divide-conquer method.
- * FFT_MUL: force floating point FFT method. */
-{
- float grammartime;
- int square = (a==b);
- int sizea, sizeb;
-
- switch (mulmode)
- {
- case GRAMMAR_MUL:
- if (square) grammarsquareg(b);
- else grammarmulg(a,b);
- break;
- case FFT_MUL:
- if (square)
- FFTsquareg(b);
- else
- FFTmulg(a,b);
- break;
- case KARAT_MUL:
- if (square) karatsquareg(b);
- else karatmulg(a,b);
- break;
- case AUTO_MUL:
- sizea = abs(a->sign);
- sizeb = abs(b->sign);
- if((sizea > KARAT_BREAK) && (sizea <= FFT_BREAK) &&
- (sizeb > KARAT_BREAK) && (sizeb <= FFT_BREAK)){
- if(square) karatsquareg(b);
- else karatmulg(a,b);
-
- } else {
- grammartime = (float)sizea;
- grammartime *= (float)sizeb;
- if (grammartime < FFT_BREAK * FFT_BREAK)
- {
- if (square) grammarsquareg(b);
- else grammarmulg(a,b);
- }
- else
- {
- if (square) FFTsquareg(b);
- else FFTmulg(a,b);
- }
- }
- break;
- }
-}
-
-void
-justg(giant x) {
- int s = x->sign, sg = 1;
-
- if(s<0) {
- sg = -1;
- s = -s;
- }
- --s;
- while(x->n[s] == 0) {
- --s;
- if(s < 0) break;
- }
- x->sign = sg*(s+1);
-}
-
-/* Next, improved Karatsuba routines from A. Powell,
- improvements by G. Woltman. */
-
-void
-karatmulg(giant x, giant y)
-/* y becomes x*y. */
-{
- int s = abs(x->sign), t = abs(y->sign), w, bits,
- sg = gsign(x)*gsign(y);
- giant a, b, c, d, e, f;
-
- if((s <= KARAT_BREAK) || (t <= KARAT_BREAK)) {
- grammarmulg(x,y);
- return;
- }
- w = (s + t + 2)/4; bits = 16*w;
- a = popg(); b = popg(); c = popg();
- d = popg(); e = popg(); f = popg();
- gtog(x,a); absg(a); if (w <= s) {a->sign = w; justg(a);}
- gtog(x,b); absg(b);
- gshiftright(bits, b);
- gtog(y,c); absg(c); if (w <= t) {c->sign = w; justg(c);}
- gtog(y,d); absg(d);
- gshiftright(bits,d);
- gtog(a,e); normal_addg(b,e); /* e := (a + b) */
- gtog(c,f); normal_addg(d,f); /* f := (c + d) */
- karatmulg(e,f); /* f := (a + b)(c + d) */
- karatmulg(c,a); /* a := a c */
- karatmulg(d,b); /* b := b d */
- normal_subg(a,f);
- /* f := (a + b)(c + d) - a c */
- normal_subg(b,f);
- /* f := (a + b)(c + d) - a c - b d */
- gshiftleft(bits, b);
- normal_addg(f, b);
- gshiftleft(bits, b);
- normal_addg(a, b);
- gtog(b, y); y->sign *= sg;
- pushg(6);
-
- return;
-}
-
-void
-karatsquareg(giant x)
-/* x becomes x^2. */
-{
- int s = abs(x->sign), w, bits;
- giant a, b, c;
-
- if(s <= KARAT_BREAK) {
- grammarsquareg(x);
- return;
- }
- w = (s+1)/2; bits = 16*w;
- a = popg(); b = popg(); c = popg();
- gtog(x, a); a->sign = w; justg(a);
- gtog(x, b); absg(b);
- gshiftright(bits, b);
- gtog(a,c); normal_addg(b,c);
- karatsquareg(c);
- karatsquareg(a);
- karatsquareg(b);
- normal_subg(b, c);
- normal_subg(a, c);
- gshiftleft(bits, b);
- normal_addg(c,b);
- gshiftleft(bits, b);
- normal_addg(a, b);
- gtog(b, x);
- pushg(3);
-
- return;
-}
-
-void
-grammarmulg(
- giant a,
- giant b
-)
-/* b becomes a*b. */
-{
- int i,j;
- unsigned int prod,carry=0;
- int asize = abs(a->sign), bsize = abs(b->sign);
- unsigned short *aptr,*bptr,*destptr;
- unsigned short mult;
- giant scratch = popg();
-
- for (i=0; i<asize+bsize; ++i)
- {
- scratch->n[i]=0;
- }
-
- bptr = &(b->n[0]);
- for (i=0; i<bsize; ++i)
- {
- mult = *(bptr++);
- if (mult)
- {
- carry = 0;
- aptr = &(a->n[0]);
- destptr = &(scratch->n[i]);
- for (j=0; j<asize; ++j)
- {
- prod = *(aptr++) * mult + *destptr + carry;
- *(destptr++) = (unsigned short)(prod & 0xffff);
- carry = prod >> 16;
- }
- *destptr = (unsigned short)carry;
- }
- }
- bsize+=asize;
- if (!carry)
- --bsize;
- scratch->sign = gsign(a)*gsign(b)*bsize;
- gtog(scratch,b);
- pushg(1);
-}
-
-
-void
-grammarsquareg (
- giant a
-)
-/* a := a^2. */
-{
- unsigned int cur_term;
- unsigned int prod, carry=0, temp;
- int asize = abs(a->sign), max = asize * 2 - 1;
- unsigned short *ptr = a->n, *ptr1, *ptr2;
- giant scratch;
-
- if(asize == 0) {
- itog(0,a);
- return;
- }
-
- scratch = popg();
-
- asize--;
-
- temp = *ptr;
- temp *= temp;
- scratch->n[0] = temp;
- carry = temp >> 16;
-
- for (cur_term = 1; cur_term < max; cur_term++) {
- ptr1 = ptr2 = ptr;
- if (cur_term <= asize) {
- ptr2 += cur_term;
- } else {
- ptr1 += cur_term - asize;
- ptr2 += asize;
- }
- prod = carry & 0xFFFF;
- carry >>= 16;
- while(ptr1 < ptr2) {
- temp = *ptr1++ * *ptr2--;
- prod += (temp << 1) & 0xFFFF;
- carry += (temp >> 15);
- }
- if (ptr1 == ptr2) {
- temp = *ptr1;
- temp *= temp;
- prod += temp & 0xFFFF;
- carry += (temp >> 16);
- }
- carry += prod >> 16;
- scratch->n[cur_term] = (unsigned short) (prod);
- }
- if (carry) {
- scratch->n[cur_term] = carry;
- scratch->sign = cur_term+1;
- } else scratch->sign = cur_term;
-
- gtog(scratch,a);
- pushg(1);
-}
-
-
-/**************************************************************
- *
- * FFT multiply Functions
- *
- **************************************************************/
-
-int initL = 0;
-
-int
-lpt(
- int n,
- int *lambda
-)
-/* Returns least power of two greater than n. */
-{
- register int i = 1;
-
- *lambda = 0;
- while (i<n)
- {
- i<<=1;
- ++(*lambda);
- }
- return(i);
-}
-
-
-void
-addsignal(
- giant x,
- double *z,
- int n
-)
-{
- register int j, k, m, car, last;
- register double f, g,err;
-
- maxFFTerror = 0;
- last = 0;
- for (j=0;j<n;j++)
- {
- f = gfloor(z[j]+0.5);
- if(f != 0.0) last = j;
- if (checkFFTerror)
- {
- err = fabs(f - z[j]);
- if (err > maxFFTerror)
- maxFFTerror = err;
- }
- z[j] =0;
- k = 0;
- do
- {
- g = gfloor(f*TWOM16);
- z[j+k] += f-g*TWO16;
- ++k;
- f=g;
- } while(f != 0.0);
- }
- car = 0;
- for(j=0;j < last + 1;j++)
- {
- m = (int)(z[j]+car);
- x->n[j] = (unsigned short)(m & 0xffff);
- car = (m>>16);
- }
- if (car)
- x->n[j] = (unsigned short)car;
- else
- --j;
-
- while(!(x->n[j])) --j;
-
- x->sign = j+1;
-}
-
-
-void
-FFTsquareg(
- giant x
-)
-{
- int j,size = abs(x->sign);
- register int L;
-
- if (size<4)
- {
- grammarmulg(x,x);
- return;
- }
- L = lpt(size+size, &j);
- if(!z) z = (double *)malloc(MAX_SHORTS * sizeof(double));
- giant_to_double(x, size, z, L);
- fft_real_to_hermitian(z, L);
- square_hermitian(z, L);
- fftinv_hermitian_to_real(z, L);
- addsignal(x,z,L);
- x->sign = abs(x->sign);
-}
-
-
-void
-FFTmulg(
- giant y,
- giant x
-)
-{
- /* x becomes y*x. */
- int lambda, sizex = abs(x->sign), sizey = abs(y->sign);
- int finalsign = gsign(x)*gsign(y);
- register int L;
-
- if ((sizex<=4)||(sizey<=4))
- {
- grammarmulg(y,x);
- return;
- }
- L = lpt(sizex+sizey, &lambda);
- if(!z) z = (double *)malloc(MAX_SHORTS * sizeof(double));
- if(!z2) z2 = (double *)malloc(MAX_SHORTS * sizeof(double));
-
- giant_to_double(x, sizex, z, L);
- giant_to_double(y, sizey, z2, L);
- fft_real_to_hermitian(z, L);
- fft_real_to_hermitian(z2, L);
- mul_hermitian(z2, z, L);
- fftinv_hermitian_to_real(z, L);
- addsignal(x,z,L);
- x->sign = finalsign*abs(x->sign);
-}
-
-
-void
-scramble_real(
- double *x,
- int n
-)
-{
- register int i,j,k;
- register double tmp;
-
- for (i=0,j=0;i<n-1;i++)
- {
- if (i<j)
- {
- tmp = x[j];
- x[j]=x[i];
- x[i]=tmp;
- }
- k = n/2;
- while (k<=j)
- {
- j -= k;
- k>>=1;
- }
- j += k;
- }
-}
-
-
-void
-fft_real_to_hermitian(
- double *z,
- int n
-)
-/* Output is {Re(z^[0]),...,Re(z^[n/2),Im(z^[n/2-1]),...,Im(z^[1]).
- * This is a decimation-in-time, split-radix algorithm.
- */
-{
- register double cc1, ss1, cc3, ss3;
- register int is, id, i0, i1, i2, i3, i4, i5, i6, i7, i8,
- a, a3, b, b3, nminus = n-1, dil, expand;
- register double *x, e;
- int nn = n>>1;
- double t1, t2, t3, t4, t5, t6;
- register int n2, n4, n8, i, j;
-
- init_sinCos(n);
- expand = cur_run/n;
- scramble_real(z, n);
- x = z-1; /* FORTRAN compatibility. */
- is = 1;
- id = 4;
- do
- {
- for (i0=is;i0<=n;i0+=id)
- {
- i1 = i0+1;
- e = x[i0];
- x[i0] = e + x[i1];
- x[i1] = e - x[i1];
- }
- is = (id<<1)-1;
- id <<= 2;
- } while(is<n);
-
- n2 = 2;
- while(nn>>=1)
- {
- n2 <<= 1;
- n4 = n2>>2;
- n8 = n2>>3;
- is = 0;
- id = n2<<1;
- do
- {
- for (i=is;i<n;i+=id)
- {
- i1 = i+1;
- i2 = i1 + n4;
- i3 = i2 + n4;
- i4 = i3 + n4;
- t1 = x[i4]+x[i3];
- x[i4] -= x[i3];
- x[i3] = x[i1] - t1;
- x[i1] += t1;
- if (n4==1)
- continue;
- i1 += n8;
- i2 += n8;
- i3 += n8;
- i4 += n8;
- t1 = (x[i3]+x[i4])*SQRTHALF;
- t2 = (x[i3]-x[i4])*SQRTHALF;
- x[i4] = x[i2] - t1;
- x[i3] = -x[i2] - t1;
- x[i2] = x[i1] - t2;
- x[i1] += t2;
- }
- is = (id<<1) - n2;
- id <<= 2;
- } while(is<n);
- dil = n/n2;
- a = dil;
- for (j=2;j<=n8;j++)
- {
- a3 = (a+(a<<1))&nminus;
- b = a*expand;
- b3 = a3*expand;
- cc1 = s_cos(b);
- ss1 = s_sin(b);
- cc3 = s_cos(b3);
- ss3 = s_sin(b3);
- a = (a+dil)&nminus;
- is = 0;
- id = n2<<1;
- do
- {
- for(i=is;i<n;i+=id)
- {
- i1 = i+j;
- i2 = i1 + n4;
- i3 = i2 + n4;
- i4 = i3 + n4;
- i5 = i + n4 - j + 2;
- i6 = i5 + n4;
- i7 = i6 + n4;
- i8 = i7 + n4;
- t1 = x[i3]*cc1 + x[i7]*ss1;
- t2 = x[i7]*cc1 - x[i3]*ss1;
- t3 = x[i4]*cc3 + x[i8]*ss3;
- t4 = x[i8]*cc3 - x[i4]*ss3;
- t5 = t1 + t3;
- t6 = t2 + t4;
- t3 = t1 - t3;
- t4 = t2 - t4;
- t2 = x[i6] + t6;
- x[i3] = t6 - x[i6];
- x[i8] = t2;
- t2 = x[i2] - t3;
- x[i7] = -x[i2] - t3;
- x[i4] = t2;
- t1 = x[i1] + t5;
- x[i6] = x[i1] - t5;
- x[i1] = t1;
- t1 = x[i5] + t4;
- x[i5] -= t4;
- x[i2] = t1;
- }
- is = (id<<1) - n2;
- id <<= 2;
- } while(is<n);
- }
- }
-}
-
-
-void
-fftinv_hermitian_to_real(
- double *z,
- int n
-)
-/* Input is {Re(z^[0]),...,Re(z^[n/2),Im(z^[n/2-1]),...,Im(z^[1]).
- * This is a decimation-in-frequency, split-radix algorithm.
- */
-{
- register double cc1, ss1, cc3, ss3;
- register int is, id, i0, i1, i2, i3, i4, i5, i6, i7, i8,
- a, a3, b, b3, nminus = n-1, dil, expand;
- register double *x, e;
- int nn = n>>1;
- double t1, t2, t3, t4, t5;
- int n2, n4, n8, i, j;
-
- init_sinCos(n);
- expand = cur_run/n;
- x = z-1;
- n2 = n<<1;
- while(nn >>= 1)
- {
- is = 0;
- id = n2;
- n2 >>= 1;
- n4 = n2>>2;
- n8 = n4>>1;
- do
- {
- for(i=is;i<n;i+=id)
- {
- i1 = i+1;
- i2 = i1 + n4;
- i3 = i2 + n4;
- i4 = i3 + n4;
- t1 = x[i1] - x[i3];
- x[i1] += x[i3];
- x[i2] += x[i2];
- x[i3] = t1 - 2.0*x[i4];
- x[i4] = t1 + 2.0*x[i4];
- if (n4==1)
- continue;
- i1 += n8;
- i2 += n8;
- i3 += n8;
- i4 += n8;
- t1 = (x[i2]-x[i1])*SQRTHALF;
- t2 = (x[i4]+x[i3])*SQRTHALF;
- x[i1] += x[i2];
- x[i2] = x[i4]-x[i3];
- x[i3] = -2.0*(t2+t1);
- x[i4] = 2.0*(t1-t2);
- }
- is = (id<<1) - n2;
- id <<= 2;
- } while (is<n-1);
- dil = n/n2;
- a = dil;
- for (j=2;j<=n8;j++)
- {
- a3 = (a+(a<<1))&nminus;
- b = a*expand;
- b3 = a3*expand;
- cc1 = s_cos(b);
- ss1 = s_sin(b);
- cc3 = s_cos(b3);
- ss3 = s_sin(b3);
- a = (a+dil)&nminus;
- is = 0;
- id = n2<<1;
- do
- {
- for(i=is;i<n;i+=id)
- {
- i1 = i+j;
- i2 = i1+n4;
- i3 = i2+n4;
- i4 = i3+n4;
- i5 = i+n4-j+2;
- i6 = i5+n4;
- i7 = i6+n4;
- i8 = i7+n4;
- t1 = x[i1] - x[i6];
- x[i1] += x[i6];
- t2 = x[i5] - x[i2];
- x[i5] += x[i2];
- t3 = x[i8] + x[i3];
- x[i6] = x[i8] - x[i3];
- t4 = x[i4] + x[i7];
- x[i2] = x[i4] - x[i7];
- t5 = t1 - t4;
- t1 += t4;
- t4 = t2 - t3;
- t2 += t3;
- x[i3] = t5*cc1 + t4*ss1;
- x[i7] = -t4*cc1 + t5*ss1;
- x[i4] = t1*cc3 - t2*ss3;
- x[i8] = t2*cc3 + t1*ss3;
- }
- is = (id<<1) - n2;
- id <<= 2;
- } while(is<n-1);
- }
- }
- is = 1;
- id = 4;
- do
- {
- for (i0=is;i0<=n;i0+=id)
- {
- i1 = i0+1;
- e = x[i0];
- x[i0] = e + x[i1];
- x[i1] = e - x[i1];
- }
- is = (id<<1) - 1;
- id <<= 2;
- } while(is<n);
- scramble_real(z, n);
- e = 1/(double)n;
- for (i=0;i<n;i++)
- {
- z[i] *= e;
- }
-}
-
-
-void
-mul_hermitian(
- double *a,
- double *b,
- int n
-)
-{
- register int k, half = n>>1;
- register double aa, bb, am, bm;
-
- b[0] *= a[0];
- b[half] *= a[half];
- for (k=1;k<half;k++)
- {
- aa = a[k];
- bb = b[k];
- am = a[n-k];
- bm = b[n-k];
- b[k] = aa*bb - am*bm;
- b[n-k] = aa*bm + am*bb;
- }
-}
-
-
-void
-square_hermitian(
- double *b,
- int n
-)
-{
- register int k, half = n>>1;
- register double c, d;
-
- b[0] *= b[0];
- b[half] *= b[half];
- for (k=1;k<half;k++)
- {
- c = b[k];
- d = b[n-k];
- b[n-k] = 2.0*c*d;
- b[k] = (c+d)*(c-d);
- }
-}
-
-
-void
-giant_to_double
-(
- giant x,
- int sizex,
- double *z,
- int L
-)
-{
- register int j;
-
- for (j=sizex;j<L;j++)
- {
- z[j]=0.0;
- }
- for (j=0;j<sizex;j++)
- {
- z[j] = x->n[j];
- }
-}
-
-
-void
-gswap(
- giant *p,
- giant *q
-)
-{
- giant t;
-
- t = *p;
- *p = *q;
- *q = t;
-}
-
-
-void
-onestep(
- giant x,
- giant y,
- gmatrix A
-)
-/* Do one step of the euclidean algorithm and modify
- * the matrix A accordingly. */
-{
- giant s4 = popg();
-
- gtog(x,s4);
- gtog(y,x);
- gtog(s4,y);
- divg(x,s4);
- punch(s4,A);
- mulg(x,s4);
- subg(s4,y);
-
- pushg(1);
-}
-
-
-void
-mulvM(
- gmatrix A,
- giant x,
- giant y
-)
-/* Multiply vector by Matrix; changes x,y. */
-{
- giant s0 = popg(), s1 = popg();
-
- gtog(A->ur,s0);
- gtog( A->lr,s1);
- dotproduct(x,y,A->ul,s0);
- dotproduct(x,y,A->ll,s1);
- gtog(s0,x);
- gtog(s1,y);
-
- pushg(2);
-}
-
-
-void
-mulmM(
- gmatrix A,
- gmatrix B
-)
-/* Multiply matrix by Matrix; changes second matrix. */
-{
- giant s0 = popg();
- giant s1 = popg();
- giant s2 = popg();
- giant s3 = popg();
-
- gtog(B->ul,s0);
- gtog(B->ur,s1);
- gtog(B->ll,s2);
- gtog(B->lr,s3);
- dotproduct(A->ur,A->ul,B->ll,s0);
- dotproduct(A->ur,A->ul,B->lr,s1);
- dotproduct(A->ll,A->lr,B->ul,s2);
- dotproduct(A->ll,A->lr,B->ur,s3);
- gtog(s0,B->ul);
- gtog(s1,B->ur);
- gtog(s2,B->ll);
- gtog(s3,B->lr);
-
- pushg(4);
-}
-
-
-void
-writeM(
- gmatrix A
-)
-{
- printf(" ul:");
- gout(A->ul);
- printf(" ur:");
- gout(A->ur);
- printf(" ll:");
- gout(A->ll);
- printf(" lr:");
- gout(A->lr);
-}
-
-
-void
-punch(
- giant q,
- gmatrix A
-)
-/* Multiply the matrix A on the left by [0,1,1,-q]. */
-{
- giant s0 = popg();
-
- gtog(A->ll,s0);
- mulg(q,A->ll);
- gswap(&A->ul,&A->ll);
- subg(A->ul,A->ll);
- gtog(s0,A->ul);
- gtog(A->lr,s0);
- mulg(q,A->lr);
- gswap(&A->ur,&A->lr);
- subg(A->ur,A->lr);
- gtog(s0,A->ur);
-
- pushg(1);
-}
-
-
-static void
-dotproduct(
- giant a,
- giant b,
- giant c,
- giant d
-)
-/* Replace last argument with the dot product of two 2-vectors. */
-{
- giant s4 = popg();
-
- gtog(c,s4);
- mulg(a, s4);
- mulg(b,d);
- addg(s4,d);
-
- pushg(1);
-}
-
-
-void
-ggcd(
- giant xx,
- giant yy
-)
-/* A giant gcd. Modifies its arguments. */
-{
- giant x = popg(), y = popg();
- gmatrix A = newgmatrix();
-
- gtog(xx,x); gtog(yy,y);
- for(;;)
- {
- fix(&x,&y);
- if (bitlen(y) <= GCDLIMIT )
- break;
- A->ul = popg();
- A->ur = popg();
- A->ll = popg();
- A->lr = popg();
- itog(1,A->ul);
- itog(0,A->ur);
- itog(0,A->ll);
- itog(1,A->lr);
- hgcd(0,x,y,A);
- mulvM(A,x,y);
- pushg(4);
- fix(&x,&y);
- if (bitlen(y) <= GCDLIMIT )
- break;
- modg(y,x);
- gswap(&x,&y);
- }
- bgcdg(x,y);
- gtog(y,yy);
- pushg(2);
- free(A);
-}
-
-
-void
-fix(
- giant *p,
- giant *q
-)
-/* Insure that x > y >= 0. */
-{
- if( gsign(*p) < 0 )
- negg(*p);
- if( gsign(*q) < 0 )
- negg(*q);
- if( gcompg(*p,*q) < 0 )
- gswap(p,q);
-}
-
-
-void
-hgcd(
- int n,
- giant xx,
- giant yy,
- gmatrix A
-)
-/* hgcd(n,x,y,A) chops n bits off x and y and computes th
- * 2 by 2 matrix A such that A[x y] is the pair of terms
- * in the remainder sequence starting with x,y that is
- * half the size of x. Note that the argument A is modified
- * but that the arguments xx and yy are left unchanged.
- */
-{
- giant x, y;
-
- if (isZero(yy))
- return;
-
- x = popg();
- y = popg();
- gtog(xx,x);
- gtog(yy,y);
- gshiftright(n,x);
- gshiftright(n,y);
- if (bitlen(x) <= INTLIMIT )
- {
- shgcd(gtoi(x),gtoi(y),A);
- }
- else
- {
- gmatrix B = newgmatrix();
- int m = bitlen(x)/2;
-
- hgcd(m,x,y,A);
- mulvM(A,x,y);
- if (gsign(x) < 0)
- {
- negg(x); negg(A->ul); negg(A->ur);
- }
- if (gsign(y) < 0)
- {
- negg(y); negg(A->ll); negg(A->lr);
- }
- if (gcompg(x,y) < 0)
- {
- gswap(&x,&y);
- gswap(&A->ul,&A->ll);
- gswap(&A->ur,&A->lr);
- }
- if (!isZero(y))
- {
- onestep(x,y,A);
- m /= 2;
- B->ul = popg();
- B->ur = popg();
- B->ll = popg();
- B->lr = popg();
- itog(1,B->ul);
- itog(0,B->ur);
- itog(0,B->ll);
- itog(1,B->lr);
- hgcd(m,x,y,B);
- mulmM(B,A);
- pushg(4);
- }
- free(B);
- }
- pushg(2);
-}
-
-
-void
-shgcd(
- register int x,
- register int y,
- gmatrix A
-)
-/*
- * Do a half gcd on the integers a and b, putting the result in A
- * It is fairly easy to use the 2 by 2 matrix description of the
- * extended Euclidean algorithm to prove that the quantity q*t
- * never overflows.
- */
-{
- register int q, t, start = x;
- int Aul = 1, Aur = 0, All = 0, Alr = 1;
-
- while (y != 0 && y > start/y)
- {
- q = x/y;
- t = y;
- y = x%y;
- x = t;
- t = All;
- All = Aul-q*t;
- Aul = t;
- t = Alr;
- Alr = Aur-q*t;
- Aur = t;
- }
- itog(Aul,A->ul);
- itog(Aur,A->ur);
- itog(All,A->ll);
- itog(Alr,A->lr);
-}
projects / apple / security.git / blobdiff