#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "platform.h"
#include "giantIntegers.h"
#include "feeDebug.h"
#include "ckconfig.h"
#include "ellipticMeasure.h"
#include "falloc.h"
#include "giantPortCommon.h"
#ifdef FEE_DEBUG
#if (GIANT_LOG2_BITS_PER_DIGIT == 4)
#warning Compiling with two-byte giantDigits
#endif
#endif
#if 0
#if FEE_DEBUG
char printbuf1[200];
char printbuf2[200];
char printbuf3[200];
void printGiantBuf(giant x)
{
int i;
sprintf(printbuf2, "sign=%d cap=%d n[]=", x->sign, x->capacity);
for(i=0; i<abs(x->sign); i++) {
sprintf(printbuf3 + 10*i, "%u:", x->n[i]);
}
}
char printbuf4[200];
char printbuf5[200];
void printGiantBuf2(giant x)
{
int i;
sprintf(printbuf4, "sign=%d cap=%d n[]=", x->sign, x->capacity);
for(i=0; i<abs(x->sign); i++) {
sprintf(printbuf5 + 10*i, "%u:", x->n[i]);
}
}
#endif
#endif
#define WARN_UNOPTIMIZE 0
#define LOG_GIANT_STACK 0
#define LOG_GIANT_STACK_OVERFLOW 1
#define WARN_ZERO_GIANT_SIZE FEE_DEBUG
#define GIANT_MAC_DEBUG 0
#if GIANT_MAC_DEBUG
#include <string.h>
#include <TextUtils.h>
static void logCom(unsigned char *str) {
c2pstr((char *)str);
DebugStr(str);
}
void dblog0(const char *str) {
Str255 outStr;
strcpy((char *)outStr, str);
logCom(outStr);
}
#else
#define dblog0(s)
#endif
#ifndef min
#define min(a,b) ((a)<(b)? (a) : (b))
#endif // min
#ifndef max
#define max(a,b) ((a)>(b)? (a) : (b))
#endif // max
#ifndef TRUE
#define TRUE 1
#endif // TRUE
#ifndef FALSE
#define FALSE 0
#endif // FALSE
static void absg(giant g);
#if GIANTS_VIA_STACK
#if LOG_GIANT_STACK
#define gstackDbg(x) printf x
#else // LOG_GIANT_STACK
#define gstackDbg(x)
#endif // LOG_GIANT_STACK
typedef struct {
unsigned numDigits; unsigned numFree; unsigned totalGiants; giant *stack;
} gstack;
static gstack *gstacks = NULL; static unsigned numGstacks = 0; static int gstackInitd = 0;
#define INIT_NUM_GIANTS 16
#define MIN_GIANT_SIZE 4
#define GIANT_SIZE_INCR 2
void initGiantStacks(unsigned maxDigits)
{
unsigned curSize = MIN_GIANT_SIZE;
unsigned sz;
unsigned i;
dblog0("initGiantStacks\n");
if(gstackInitd) {
printf("multiple initGiantStacks calls\n");
return;
}
gstackDbg(("initGiantStacks(%d)\n", maxDigits));
numGstacks = 1;
while(curSize<=maxDigits) {
curSize <<= GIANT_SIZE_INCR;
numGstacks++;
}
sz = sizeof(gstack) * numGstacks;
gstacks = (gstack*) fmalloc(sz);
bzero(gstacks, sz);
curSize = MIN_GIANT_SIZE;
for(i=0; i<numGstacks; i++) {
gstacks[i].numDigits = curSize;
curSize <<= GIANT_SIZE_INCR;
}
gstackInitd = 1;
}
void freeGiantStacks(void)
{
int i;
int j;
gstack *gs;
if(!gstackInitd) {
return;
}
for(i=0; i<numGstacks; i++) {
gs = &gstacks[i];
for(j=0; j<gs->numFree; j++) {
freeGiant(gs->stack[j]);
gs->stack[j] = NULL;
}
if(gs->stack != NULL) {
ffree(gs->stack);
gs->stack = NULL;
}
}
ffree(gstacks);
gstacks = NULL;
gstackInitd = 0;
}
#endif // GIANTS_VIA_STACK
giant borrowGiant(unsigned numDigits)
{
giant result;
#if GIANTS_VIA_STACK
unsigned stackNum;
gstack *gs = gstacks;
#if WARN_ZERO_GIANT_SIZE
if(numDigits == 0) {
printf("borrowGiant(0)\n");
numDigits = gstacks[numGstacks-1].numDigits;
}
#endif
if(numDigits <= MIN_GIANT_SIZE)
stackNum = 0;
else if (numDigits <= (MIN_GIANT_SIZE << GIANT_SIZE_INCR))
stackNum = 1;
else if (numDigits <= (MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR)))
stackNum = 2;
else if (numDigits <= (MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR)))
stackNum = 3;
else if (numDigits <= (MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR)))
stackNum = 4;
else
stackNum = numGstacks;
if(stackNum >= numGstacks) {
#if LOG_GIANT_STACK_OVERFLOW
gstackDbg(("giantFromStack overflow; numDigits %d\n",
numDigits));
#endif return newGiant(numDigits);
}
gs = &gstacks[stackNum];
#if GIANT_MAC_DEBUG
if((gs->numFree != 0) && (gs->stack == NULL)) {
dblog0("borrowGiant: null stack!\n");
}
#endif
if(gs->numFree != 0) {
result = gs->stack[--gs->numFree];
}
else {
result = newGiant(gs->numDigits);
}
#else
result = newGiant(numDigits);
#endif
PROF_INCR(numBorrows);
return result;
}
void returnGiant(giant g)
{
#if GIANTS_VIA_STACK
unsigned stackNum;
gstack *gs;
unsigned cap = g->capacity;
#if FEE_DEBUG
if(!gstackInitd) {
CKRaise("returnGiant before stacks initialized!");
}
#endif
#if GIANT_MAC_DEBUG
if(g == NULL) {
dblog0("returnGiant: null g!\n");
}
#endif
switch(cap) {
case MIN_GIANT_SIZE:
stackNum = 0;
break;
case MIN_GIANT_SIZE << GIANT_SIZE_INCR:
stackNum = 1;
break;
case MIN_GIANT_SIZE << (2 * GIANT_SIZE_INCR):
stackNum = 2;
break;
case MIN_GIANT_SIZE << (3 * GIANT_SIZE_INCR):
stackNum = 3;
break;
case MIN_GIANT_SIZE << (4 * GIANT_SIZE_INCR):
stackNum = 4;
break;
default:
stackNum = numGstacks;
break;
}
if(stackNum >= numGstacks) {
#if LOG_GIANT_STACK_OVERFLOW
gstackDbg(("giantToStack overflow; numDigits %d\n", cap));
#endif freeGiant(g);
return;
}
gs = &gstacks[stackNum];
if(gs->numFree == gs->totalGiants) {
if(gs->totalGiants == 0) {
gstackDbg(("Initial alloc of gstack(%d)\n",
gs->numDigits));
gs->totalGiants = INIT_NUM_GIANTS;
}
else {
gs->totalGiants *= 2;
gstackDbg(("Bumping gstack(%d) to %d\n",
gs->numDigits, gs->totalGiants));
}
gs->stack = (giantstruct**) frealloc(gs->stack, gs->totalGiants*sizeof(giant));
}
g->sign = 0; gs->stack[gs->numFree++] = g;
#if GIANT_MAC_DEBUG
if((gs->numFree != 0) && (gs->stack == NULL)) {
dblog0("borrowGiant: null stack!\n");
}
#endif
#else
freeGiant(g);
#endif
}
void freeGiant(giant x) {
ffree(x);
}
giant newGiant(unsigned numDigits) {
int size;
giant result;
#if WARN_ZERO_GIANT_SIZE
if(numDigits == 0) {
printf("newGiant(0)\n");
#if GIANTS_VIA_STACK
numDigits = gstacks[numGstacks-1].totalGiants;
#else
numDigits = 20;
#endif
}
#endif
size = (numDigits-1) * GIANT_BYTES_PER_DIGIT + sizeof(giantstruct);
result = (giant)fmalloc(size);
if(result == NULL) {
return NULL;
}
result->sign = 0;
result->capacity = numDigits;
return result;
}
giant copyGiant(giant x)
{
int bytes;
giant result = newGiant(x->capacity);
bytes = sizeof(giantstruct) +
((x->capacity - 1) * GIANT_BYTES_PER_DIGIT);
bcopy(x, result, bytes);
return result;
}
unsigned bitlen(giant n) {
unsigned b = GIANT_BITS_PER_DIGIT;
giantDigit c = ((giantDigit)1) << (GIANT_BITS_PER_DIGIT - 1);
giantDigit w;
if (isZero(n)) {
return(0);
}
w = n->n[abs(n->sign) - 1];
if (!w) {
CKRaise("bitlen - no bit set!");
}
while((w&c) == 0) {
b--;
c >>= 1;
}
return(GIANT_BITS_PER_DIGIT * (abs(n->sign)-1) + b);
}
int bitval(giant n, int pos) {
int i = abs(pos) >> GIANT_LOG2_BITS_PER_DIGIT;
giantDigit c = ((giantDigit)1) << (pos & (GIANT_BITS_PER_DIGIT - 1));
return ((0!=((n->n[i]) & c))?1:0);
}
int gsign(giant g)
{
if (isZero(g)) return(0);
if (g->sign > 0) return(1);
return(-1);
}
void gtrimSign(giant g)
{
int numDigits = abs(g->sign);
int i;
for(i=numDigits-1; i>=0; i--) {
if(g->n[i] == 0) {
numDigits--;
}
else {
break;
}
}
if(g->sign < 0) {
g->sign = -numDigits;
}
else {
g->sign = numDigits;
}
}
int isone(giant g) {
return((g->sign==1)&&(g->n[0]==1));
}
int isZero(giant thegiant) {
int count;
int length = abs(thegiant->sign);
giantDigit *numpointer;
if (length) {
numpointer = thegiant->n;
for(count = 0; count<length; ++count,++numpointer)
if (*numpointer != 0 ) return(FALSE);
}
return(TRUE);
}
int gcompg(giant a, giant b)
{
int sa = a->sign;
int j;
int sb = b->sign;
giantDigit va;
giantDigit vb;
int sgn;
if(isZero(a) && isZero(b)) return 0;
if(sa > sb) return(1);
if(sa < sb) return(-1);
if(sa < 0) {
sa = -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);
}
void gtog(giant srcgiant, giant destgiant) {
int numbytes;
CKASSERT(srcgiant != NULL);
numbytes = abs(srcgiant->sign) * GIANT_BYTES_PER_DIGIT;
if (destgiant->capacity < abs(srcgiant->sign))
CKRaise("gtog overflow!!");
memcpy((char *)destgiant->n, (char *)srcgiant->n, numbytes);
destgiant->sign = srcgiant->sign;
}
void int_to_giant(int i, giant g) {
int isneg = (i<0);
unsigned int j = abs(i);
unsigned dex;
g->sign = 0;
if (i==0) {
g->n[0] = 0;
return;
}
if(GIANT_BYTES_PER_DIGIT == sizeof(int)) {
g->n[0] = j;
g->sign = 1;
}
else {
unsigned scnt = GIANT_BITS_PER_DIGIT;
for(dex=0; dex<sizeof(int); dex++) {
g->n[dex] = j & GIANT_DIGIT_MASK;
j >>= scnt;
g->sign++;
if(j == 0) {
break;
}
}
}
if (isneg) {
g->sign = -(g->sign);
}
}
void negg(giant g) {
g->sign = -g->sign;
}
void iaddg(int i, giant g) {
int j;
giantDigit carry;
int size = abs(g->sign);
if (isZero(g)) {
int_to_giant(i,g);
}
else {
carry = i;
for(j=0; ((j<size) && (carry != 0)); j++) {
g->n[j] = giantAddDigits(g->n[j], carry, &carry);
}
if(carry) {
++g->sign;
if (g->sign > (int)g->capacity) CKRaise("iaddg overflow!");
g->n[size] = carry;
}
}
}
void imulg(unsigned n, giant g)
{
giant tmp = borrowGiant(abs(g->sign) + sizeof(int));
int_to_giant(n, tmp);
mulg(tmp, g);
returnGiant(tmp);
}
static void normal_addg(giant a, giant b)
{
giantDigit carry1 = 0;
giantDigit carry2 = 0;
int asize = a->sign, bsize = b->sign;
giantDigit *an = a->n;
giantDigit *bn = b->n;
giantDigit tmp;
int j;
int comSize;
int maxSize;
if(asize < bsize) {
comSize = asize;
maxSize = bsize;
}
else {
comSize = bsize;
maxSize = asize;
}
for(j=0; j<comSize; j++) {
if(carry1 || carry2) {
tmp = giantAddDigits(bn[j], (giantDigit)1, &carry1);
}
else {
carry1 = 0;
tmp = bn[j];
}
bn[j] = giantAddDigits(tmp, an[j], &carry2);
}
if(asize < bsize) {
if(carry2) {
carry1 = 1;
}
if(carry1) {
for(; j<bsize; j++) {
bn[j] = giantAddDigits(bn[j], (giantDigit)1, &carry1);
if(carry1 == 0) {
break;
}
}
}
} else {
if(carry2) {
carry1 = 1;
}
for(; j<asize; j++) {
if(carry1) {
bn[j] = giantAddDigits(an[j], (giantDigit)1, &carry1);
}
else {
bn[j] = an[j];
carry1 = 0;
}
}
}
b->sign = maxSize;
if(carry1) {
bn[j] = 1;
b->sign++;
if (b->sign > (int)b->capacity) CKRaise("iaddg overflow!");
}
}
static void normal_subg(giant a, giant b)
{
int j;
int size = b->sign;
giantDigit tmp;
giantDigit borrow1 = 0;
giantDigit borrow2 = 0;
giantDigit *an = a->n;
giantDigit *bn = b->n;
if(a->sign == 0) {
return;
}
for (j=0; j<a->sign; ++j) {
if(borrow1 || borrow2) {
tmp = giantSubDigits(bn[j], (giantDigit)1, &borrow1);
}
else {
tmp = bn[j];
borrow1 = 0;
}
bn[j] = giantSubDigits(tmp, an[j], &borrow2);
}
if(borrow1 || borrow2) {
borrow1 = 1;
for (j=a->sign; j<size; ++j) {
if(borrow1) {
bn[j] = giantSubDigits(bn[j], (giantDigit)1, &borrow1);
}
else {
break;
}
}
}
while((size-- > 0) && (b->n[size] == 0))
;
b->sign = (b->n[size] == 0)? 0 : size+1;
}
static void reverse_subg(giant a, giant b)
{
int j;
int size = a->sign;
giantDigit tmp;
giantDigit borrow1 = 0;
giantDigit borrow2 = 0;
giantDigit *an = a->n;
giantDigit *bn = b->n;
if(b->sign == 0) {
gtog(a, b);
return;
}
for (j=0; j<b->sign; ++j) {
if(borrow1 || borrow2) {
tmp = giantSubDigits(an[j], (giantDigit)1, &borrow1);
}
else {
tmp = an[j];
borrow1 = 0;
}
bn[j] = giantSubDigits(tmp, bn[j], &borrow2);
}
if(borrow1 || borrow2) {
borrow1 = 1;
}
for (j=b->sign; j<size; ++j) {
if(borrow1) {
bn[j] = giantSubDigits(an[j], (giantDigit)1, &borrow1);
}
else {
bn[j] = an[j];
borrow1 = 0;
}
}
b->sign = size;
while(!b->n[--size]);
b->sign = size+1;
}
void addg(giant a, giant b)
{ 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;
}
negg(a); if(a != b) negg(b); normal_addg(a,b);
negg(a); if(a != b) negg(b); 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)
{
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;
}
static void bdivg(giant v, giant u)
{
int diff = bitlen(u) - bitlen(v);
giant scratch7;
if (diff<0) {
int_to_giant(0,u);
return;
}
scratch7 = borrowGiant(u->capacity);
gtog(v, scratch7);
gshiftleft(diff,scratch7);
if(gcompg(u,scratch7) < 0) diff--;
if(diff<0) {
int_to_giant(0,u);
returnGiant(scratch7);
return;
}
int_to_giant(1,u);
gshiftleft(diff,u);
returnGiant(scratch7);
}
int binvaux(giant p, giant x)
{
giant scratch7;
giant u0;
giant u1;
giant v0;
giant v1;
int result = 1;
int giantSize;
PROF_START;
if(isone(x)) return(result);
giantSize = 4 * abs(p->sign);
scratch7 = borrowGiant(giantSize);
u0 = borrowGiant(giantSize);
u1 = borrowGiant(giantSize);
v0 = borrowGiant(giantSize);
v1 = borrowGiant(giantSize);
int_to_giant(1, v0); gtog(x, v1);
int_to_giant(0,x); gtog(p, u1);
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);
}
if (!isone(u1)) {
gtog(u1,x);
if(x->sign<0) addg(p, x);
result = 0;
goto done;
}
if (x->sign<0) addg(p, x);
done:
returnGiant(scratch7);
returnGiant(u0);
returnGiant(u1);
returnGiant(v0);
returnGiant(v1);
PROF_END(binvauxTime);
return(result);
}
#if 0
int binvg(giant p, giant x)
{
modg(p, x);
return(binvaux(p,x));
}
#endif
static void absg(giant g) {
if (g->sign < 0) g->sign = -g->sign;
}
void gshiftleft(int bits, giant g) {
int rem = bits & (GIANT_BITS_PER_DIGIT - 1);
int crem = GIANT_BITS_PER_DIGIT - rem;
int digits = 1 + (bits >> GIANT_LOG2_BITS_PER_DIGIT);
int size = abs(g->sign);
int j;
int k;
int sign = gsign(g);
giantDigit carry;
giantDigit dat;
#if FEE_DEBUG
if(bits < 0) {
CKRaise("gshiftleft(-bits)\n");
}
#endif
if(!bits) return;
if(!size) return;
if((size+digits) > (int)g->capacity) {
CKRaise("gshiftleft overflow");
}
k = size - 1 + digits; carry = 0;
if(rem == 0) {
g->n[k] = 0; for(j=size-1; j>=0; j--) {
g->n[--k] = g->n[j];
}
do{
g->n[--k] = 0;
} while(k>0);
}
else {
for(j=size-1; j>=0; j--) {
dat = g->n[j];
g->n[k--] = (dat >> crem) | carry;
carry = (dat << rem);
}
do{
g->n[k--] = carry;
carry = 0;
} while(k>=0);
}
k = size - 1 + digits;
if(g->n[k] == 0) --k;
g->sign = sign * (k+1);
if (abs(g->sign) > g->capacity) {
CKRaise("gshiftleft overflow");
}
}
void gshiftright(int bits, giant g) {
int j;
int size=abs(g->sign);
giantDigit carry;
int digits = bits >> GIANT_LOG2_BITS_PER_DIGIT;
int remain = bits & (GIANT_BITS_PER_DIGIT - 1);
int cremain = GIANT_BITS_PER_DIGIT - remain;
#if FEE_DEBUG
if(bits < 0) {
CKRaise("gshiftright(-bits)\n");
}
#endif
if(bits==0) return;
if(isZero(g)) return;
if (digits >= size) {
g->sign = 0;
return;
}
size -= digits;
if(remain == 0) {
if(g->sign > 0) {
g->sign = size;
}
else {
g->sign = -size;
}
for(j=0; j < size; j++) {
g->n[j] = g->n[j+digits];
}
return;
}
for(j=0;j<size;++j) {
if (j==size-1) {
carry = 0;
}
else {
carry = (g->n[j+digits+1]) << cremain;
}
g->n[j] = ((g->n[j+digits]) >> remain ) | carry;
}
if (g->n[size-1] == 0) {
--size;
}
if(g->sign > 0) {
g->sign = size;
}
else {
g->sign = -size;
}
if (abs(g->sign) > g->capacity) {
CKRaise("gshiftright overflow");
}
}
void extractbits(unsigned n, giant src, giant dest) {
int digits = n >> GIANT_LOG2_BITS_PER_DIGIT;
int numbytes = digits * GIANT_BYTES_PER_DIGIT;
int bits = n & (GIANT_BITS_PER_DIGIT - 1);
if (n <= 0) {
return;
}
if (dest->capacity * 8 * GIANT_BYTES_PER_DIGIT < n) {
CKRaise("extractbits - not enough room");
}
if (digits >= abs(src->sign)) {
gtog(src,dest);
}
else {
memcpy((char *)(dest->n), (char *)(src->n), numbytes);
if (bits) {
dest->n[digits] = src->n[digits] & ((1<<bits)-1);
++digits;
}
while((digits > 0) && (dest->n[digits-1] == 0)) {
--digits;
}
if(src->sign < 0) {
dest->sign = -digits;
}
else {
dest->sign = digits;
}
}
if (abs(dest->sign) > dest->capacity) {
CKRaise("extractbits overflow");
}
}
#define NEW_MERSENNE 0
#if NEW_MERSENNE
void
gmersennemod(
int n,
giant g
)
{
int the_sign;
giant scratch3 = borrowGiant(g->capacity);
giant scratch4 = borrowGiant(1);
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)) goto out;
int_to_giant(1,scratch3);
gshiftleft(n,scratch3);
int_to_giant(1,scratch4);
subg(scratch4,scratch3);
if(gcompg(g,scratch3) >= 0) subg(scratch3,g);
if (the_sign < 0) {
g->sign = -g->sign;
addg(scratch3,g);
}
out:
returnGiant(scratch3);
returnGiant(scratch4);
}
#else
void gmersennemod(int n, giant g) {
unsigned bits = n & (GIANT_BITS_PER_DIGIT - 1);
unsigned digits = 1 + ((n-1) >> GIANT_LOG2_BITS_PER_DIGIT);
int isPositive = (g->sign > 0);
int j;
int b;
int size;
int foundzero;
giantDigit mask = (bits == 0) ? GIANT_DIGIT_MASK : (giantDigit)((1<<bits)-1);
giant scratch1;
b = bitlen(g);
if(b < n) {
unsigned numDigits = (n + GIANT_BITS_PER_DIGIT - 1) >>
GIANT_LOG2_BITS_PER_DIGIT;
giantDigit lastWord = 0;
giantDigit bits = 1;
if(g->sign >= 0) return;
scratch1 = borrowGiant(numDigits + 1);
scratch1->sign = numDigits;
for(j=0; j<(int)(numDigits-1); j++) {
scratch1->n[j] = GIANT_DIGIT_MASK;
}
for(j=0; j < (int)(n & (GIANT_BITS_PER_DIGIT-1)); j++) {
lastWord |= bits;
bits <<= 1;
}
scratch1->n[numDigits-1] = lastWord;
addg(g, scratch1);
gtog(scratch1, g);
returnGiant(scratch1);
return;
}
if(b == n) {
for(foundzero=0, j=0; j<b; j++) {
if(bitval(g, j)==0) {
foundzero = 1;
break;
}
}
if (!foundzero) {
int_to_giant(0, g);
return;
}
}
absg(g);
scratch1 = borrowGiant(g->capacity);
while ( ((unsigned)(g->sign) > digits) ||
( ((unsigned)(g->sign)==digits) && (g->n[digits-1] > mask))) {
extractbits(n, g, scratch1);
gshiftright(n, g);
addg(scratch1, g);
}
size = g->sign;
if (!isPositive) {
for(j = digits-1; j >= size; j--) {
g->n[j] = GIANT_DIGIT_MASK;
}
for(j = size-1; j >= 0; j--) {
g->n[j] = ~g->n[j];
}
g->n[digits-1] &= mask;
j = digits-1;
while((g->n[j] == 0) && (j > 0)) {
--j;
}
size = j+1;
}
g->sign = size;
if (abs(g->sign) > g->capacity) {
CKRaise("gmersennemod overflow");
}
if (size < (int)digits) {
goto bye;
}
if (g->n[size-1] != mask) {
goto bye;
}
mask = GIANT_DIGIT_MASK;
for(j=0; j<(size-1); j++) {
if (g->n[j] != mask) {
goto bye;
}
}
g->sign = 0;
bye:
returnGiant(scratch1);
}
#endif
void mulg(giant a, giant b) {
int i;
int asize, bsize;
giantDigit *bptr = b->n;
giantDigit mult;
giant scratch1;
giantDigit carry;
giantDigit *scrPtr;
if (isZero(b)) {
return;
}
if (isZero(a)) {
gtog(a, b);
return;
}
if(a == b) {
grammarSquare(b);
return;
}
bsize = abs(b->sign);
asize = abs(a->sign);
scratch1 = borrowGiant((asize+bsize));
scrPtr = scratch1->n;
for(i=0; i<asize+bsize; ++i) {
scrPtr[i]=0;
}
for(i=0; i<bsize; ++i, scrPtr++) {
mult = bptr[i];
if (mult != 0) {
carry = VectorMultiply(mult,
a->n,
asize,
scrPtr);
scrPtr[asize] += carry;
}
}
bsize+=asize;
if(scratch1->n[bsize - 1] == 0) {
--bsize;
}
scratch1->sign = gsign(a) * gsign(b) * bsize;
if (abs(scratch1->sign) > scratch1->capacity) {
CKRaise("GiantGrammarMul overflow");
}
gtog(scratch1,b);
returnGiant(scratch1);
#if FEE_DEBUG
(void)bitlen(b); #endif
PROF_INCR(numMulg); INCR_MULGS; }
void grammarSquare(giant a) {
giantDigit prodLo;
giantDigit prodHi;
giantDigit carryLo = 0;
giantDigit carryHi = 0;
giantDigit tempLo;
giantDigit tempHi;
unsigned int cur_term;
unsigned asize;
unsigned max;
giantDigit *ptr = a->n;
giantDigit *ptr1;
giantDigit *ptr2;
giant scratch;
if(a->sign == 0) {
goto end;
}
asize = abs(a->sign);
max = asize * 2 - 1;
scratch = borrowGiant(2 * asize);
asize--;
giantMulDigits(*ptr, *ptr, &tempLo, &tempHi);
scratch->n[0] = tempLo;
carryLo = tempHi;
carryHi = 0;
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;
}
prodLo = carryLo;
prodHi = 0;
carryLo = carryHi;
carryHi = 0;
while(ptr1 < ptr2) {
giantMulDigits(*ptr1++, *ptr2--, &tempLo, &tempHi);
giantAddDouble(&prodLo, &prodHi, (tempLo << 1));
giantAddDouble(&carryLo, &carryHi,
(tempLo >> (GIANT_BITS_PER_DIGIT - 1)));
giantAddDouble(&carryLo, &carryHi, (tempHi << 1));
carryHi += (tempHi >> (GIANT_BITS_PER_DIGIT - 1));
}
if (ptr1 == ptr2) {
giantMulDigits(*ptr1, *ptr1, &tempLo, &tempHi);
giantAddDouble(&prodLo, &prodHi, tempLo);
giantAddDouble(&carryLo, &carryHi, tempHi);
}
giantAddDouble(&carryLo, &carryHi, prodHi);
scratch->n[cur_term] = prodLo;
}
if (carryLo) {
scratch->n[cur_term] = carryLo;
scratch->sign = cur_term+1;
} else scratch->sign = cur_term;
gtog(scratch,a);
returnGiant(scratch);
end:
PROF_INCR(numGsquare);
}
void clearGiant(giant g)
{
unsigned i;
for(i=0; i<g->capacity; i++) {
g->n[i] = 0;
}
g->sign = 0;
}
#if ENGINE_127_BITS
int
scompg(int n, giant g) {
if((g->sign == 1) && (g->n[0] == n)) return(1);
return(0);
}
#endif // ENGINE_127_BITS
void
make_recip(giant d, giant r)
{
int b;
int giantSize = 4 * abs(d->sign);
giant tmp = borrowGiant(giantSize);
giant tmp2 = borrowGiant(giantSize);
if (isZero(d) || (d->sign < 0))
{
CKRaise("illegal argument to make_recip");
}
int_to_giant(1, r); subg(r, d); b = bitlen(d); addg(r, d);
gshiftleft(b, r); gtog(r, tmp2);
while(1) {
gtog(r, tmp);
gsquare(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);
}
int_to_giant(1, tmp);
gshiftleft(2*b, tmp);
gtog(r, tmp2); mulg(d, tmp2);
subg(tmp2, tmp);
int_to_giant(1, tmp2);
while(tmp->sign < 0) {
subg(tmp2, r);
addg(d, tmp);
}
returnGiant(tmp);
returnGiant(tmp2);
return;
}
void
divg_via_recip(giant d, giant r, giant n)
{
int s = 2*(bitlen(r)-1), sign = gsign(n);
int giantSize = (4 * abs(d->sign)) + abs(n->sign);
giant tmp = borrowGiant(giantSize);
giant tmp2 = borrowGiant(giantSize);
if (isZero(d) || (d->sign < 0))
{
CKRaise("illegal argument to divg_via_recip");
}
n->sign = abs(n->sign);
int_to_giant(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;
returnGiant(tmp);
returnGiant(tmp2);
return;
}
void
modg_via_recip(
giant d,
giant r,
giant n
)
{
int s = (bitlen(r)-1), sign = n->sign;
int giantSize = (4 * abs(d->sign)) + abs(n->sign);
giant tmp, tmp2;
tmp = borrowGiant(giantSize);
tmp2 = borrowGiant(giantSize);
if (isZero(d) || (d->sign < 0))
{
CKRaise("illegal argument to modg_via_recip");
}
n->sign = abs(n->sign);
while (1)
{
gtog(n, tmp);
if(s == 0) {
gshiftleft(1, tmp);
}
else {
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:
returnGiant(tmp);
returnGiant(tmp2);
return;
}
void
modg(
giant d,
giant n
)
{
giant recip = borrowGiant(2 * d->capacity);
#if WARN_UNOPTIMIZE
dbgLog(("Warning: unoptimized modg!\n"));
#endif
make_recip(d, recip);
modg_via_recip(d, recip, n);
returnGiant(recip);
}
void
divg(
giant d,
giant n
)
{
giant recip = borrowGiant(2 * abs(d->sign));
#if WARN_UNOPTIMIZE
dbgLog(("Warning: unoptimized divg!\n"));
#endif
make_recip(d, recip);
divg_via_recip(d, recip, n);
returnGiant(recip);
}