ocaml/otherlibs/num/bignum/c/KerN.c

861 lines
17 KiB
C
Raw Blame History

This file contains invisible Unicode characters!

This file contains invisible Unicode characters that may be processed differently from what appears below. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to reveal hidden characters.

/* Copyright Digital Equipment Corporation & INRIA 1988, 1989 */
/* Last modified_on Thu Feb 20 18:18:12 GMT+1:00 1992 by shand */
/* modified_on Tue Jan 15 19:32:53 GMT+1:00 1991 by herve */
/* KerN.c: the kernel written in C */
/*
* Description of types and constants.
*
* Several conventions are used in the commentary:
* A "BigNum" is the name for an infinite-precision number.
* Capital letters (e.g., "N") are used to refer to the value of BigNums.
* The word "digit" refers to a single BigNum digit.
* The notation "Size(N)" refers to the number of digits in N,
* which is typically passed to the subroutine as "nl".
* The notation "Length(N)" refers to the number of digits in N,
* not including any leading zeros.
* The word "Base" is used for the number 2 ** BN_DIGIT_SIZE, where
* BN_DIGIT_SIZE is the number of bits in a single BigNum digit.
* The expression "BBase(N)" is used for Base ** NumDigits(N).
* The term "leading zeros" refers to any zeros before the most
* significant digit of a number.
*
*
* In the code, we have:
*
* "nn" is a pointer to a big number,
* "nl" is the number of digits from nn,
* "d" is a digit.
*
*/
/* */
#define BNNMACROS_OFF
#include "BigNum.h"
#define NOMEM
/*** copyright ***/
static char copyright[]="@(#)KerN.c: copyright Digital Equipment Corporation & INRIA 1988, 1989\n";
/******* non arithmetic access to digits ********/
#ifndef _NO_PROTO
void BnnSetToZero (BigNum nn, BigNumLength nl)
#else
void BnnSetToZero (nn, nl)
BigNum nn; BigNumLength nl;
#endif
/*
* Sets all the specified digits of the BigNum to 0
*/
{
BigNum nnlim;
if (nl <= 0)
return;
nnlim = nn+nl-1;
do *nn = 0; while(nn++ < nnlim);
}
/***************************************/
#ifndef _NO_PROTO
void BnnAssign (BigNum mm, BigNum nn, BigNumLength nl)
#else /* _NO_PROTO */
void BnnAssign ( mm, nn, nl)
BigNum mm; BigNum nn; BigNumLength nl;
#endif /* _NO_PROTO */
/*
* Copies N => M
*/
{
BigNum nnlim;
if (nl <= 0)
return;
nnlim = nn+nl;
#ifdef MSDOS
if (realaddr(mm) < realaddr(nn) || realaddr(mm) > realaddr(nnlim))
#else
if ((mm < nn) || ( mm > nnlim))
#endif
do *mm++ = *nn++; while(nn < nnlim);
else
#ifdef MSDOS
if (realaddr(mm) > realaddr(nn))
#else
if (mm > nn)
#endif
{
mm += nl;
do *--mm = *--nnlim; while(nn < nnlim);
}
}
/***************************************/
/* */
#ifndef _NO_PROTO
void BnnSetDigit (BigNum nn, BigNumDigit d)
#else /* _NO_PROTO */
void BnnSetDigit ( nn, d)
BigNum nn; BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Sets a single digit of N to the passed value
*/
{
*nn = d;
}
/***************************************/
#ifndef _NO_PROTO
BigNumDigit BnnGetDigit (BigNum nn)
#else /* _NO_PROTO */
BigNumDigit BnnGetDigit ( nn)
BigNum nn;
#endif /* _NO_PROTO */
/*
* Returns the single digit pointed by N
*/
{
return (*nn);
}
/***************************************/
/* */
#ifndef _NO_PROTO
BigNumLength BnnNumDigits (BigNum nn, BigNumLength nl)
#else /* _NO_PROTO */
BigNumLength BnnNumDigits ( nn, nl)
BigNum nn; BigNumLength nl;
#endif /* _NO_PROTO */
/*
* Returns the number of digits of N, not counting leading zeros
*/
{
nn += nl;
while (nl != 0 && *--nn == 0)
nl--;
return (nl == 0 ? 1 : nl);
}
/***************************************/
#ifndef _NO_PROTO
BigNumDigit BnnNumLeadingZeroBitsInDigit (BigNumDigit d)
#else /* _NO_PROTO */
BigNumDigit BnnNumLeadingZeroBitsInDigit ( d)
BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns the number of leading zero bits in a digit
*/
{
register int p = 0;
if (BN_DIGIT_SIZE == 16 || BN_DIGIT_SIZE == 32 || BN_DIGIT_SIZE == 64)
{
register BigNumDigit mask = (~(BigNumDigit)0) << (BN_DIGIT_SIZE/2);
register BigNumLength maskl = BN_DIGIT_SIZE/2;
if (d == 0)
return (BN_DIGIT_SIZE);
while (maskl)
{
if ((d & mask) == 0)
{
p += maskl;
d <<= maskl;
}
maskl >>= 1;
mask <<= maskl;
}
}
else
{
register BigNumDigit mask = ((BigNumDigit)1) << (BN_DIGIT_SIZE-1);
while ((d & mask) == 0)
{
p++;
mask >>= 1;
}
}
return (p);
}
/***************************************/
/* */
/************** Predicates on one digit ***************/
#ifndef _NO_PROTO
Boolean BnnDoesDigitFitInWord (BigNumDigit d)
#else /* _NO_PROTO */
Boolean BnnDoesDigitFitInWord ( d)
BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns TRUE iff the digit can be represented in just BN_WORD_SIZE bits
*/
{
/* The C compiler must evaluate the predicate at compile time */
if (BN_DIGIT_SIZE > BN_WORD_SIZE)
return (d >= ((BigNumDigit)1) << BN_WORD_SIZE ? FALSE : TRUE);
else
return (TRUE);
}
/***************************************/
#ifndef _NO_PROTO
Boolean BnnIsDigitZero (BigNumDigit d)
#else /* _NO_PROTO */
Boolean BnnIsDigitZero ( d)
BigNumDigit d;
#endif /* _NO_PROTO */
/* Returns TRUE iff digit = 0 */
{
return (d == 0);
}
/***************************************/
#ifndef _NO_PROTO
Boolean BnnIsDigitNormalized (BigNumDigit d)
#else /* _NO_PROTO */
Boolean BnnIsDigitNormalized ( d)
BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns TRUE iff Base/2 <= digit < Base
* i.e., if digit's leading bit is 1
*/
{
return (d & (((BigNumDigit)1) << (BN_DIGIT_SIZE - 1)) ? TRUE : FALSE);
}
/***************************************/
#ifndef _NO_PROTO
Boolean BnnIsDigitOdd (BigNumDigit d)
#else /* _NO_PROTO */
Boolean BnnIsDigitOdd ( d)
BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns TRUE iff digit is odd
*/
{
return (d & 1 ? TRUE : FALSE);
}
/***************************************/
#ifndef _NO_PROTO
BigNumCmp BnnCompareDigits (BigNumDigit d1, BigNumDigit d2)
#else /* _NO_PROTO */
BigNumCmp BnnCompareDigits ( d1, d2)
BigNumDigit d1; BigNumDigit d2;
#endif /* _NO_PROTO */
/*
* Returns BN_GREATER if digit1 > digit2
* BN_EQUAL if digit1 = digit2
* BN_LESS if digit1 < digit2
*/
{
return (d1 > d2 ? BN_GT : (d1 == d2 ? BN_EQ : BN_LT));
}
/***************** Logical operations ********************/
#ifndef _NO_PROTO
void BnnComplement (BigNum nn, BigNumLength nl)
#else /* _NO_PROTO */
void BnnComplement ( nn, nl)
BigNum nn; BigNumLength nl;
#endif /* _NO_PROTO */
/*
* Performs the computation BBase(N) - N - 1 => N
*/
{
BigNum nnlim;
if (nl <= 0)
return;
nnlim = nn+nl;
do
{
nn++;
nn[-1] = ~nn[-1];
}
while (nn < nnlim);
}
/***************************************/
/* */
#ifndef _NO_PROTO
void BnnAndDigits (BigNum n, BigNumDigit d)
#else /* _NO_PROTO */
void BnnAndDigits ( n, d)
BigNum n; BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns the logical computation n[0] AND d in n[0]
*/
{
*n &= d;
}
/***************************************/
#ifndef _NO_PROTO
void BnnOrDigits (BigNum n, BigNumDigit d)
#else /* _NO_PROTO */
void BnnOrDigits ( n, d)
BigNum n; BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns the logical computation n[0] OR d2 in n[0].
*/
{
*n |= d;
}
/***************************************/
#ifndef _NO_PROTO
void BnnXorDigits (BigNum n, BigNumDigit d)
#else /* _NO_PROTO */
void BnnXorDigits ( n, d)
BigNum n; BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Returns the logical computation n[0] XOR d in n[0].
*/
{
*n ^= d;
}
/***************************************/
/* */
/****************** Shift operations *******************/
#ifndef _NO_PROTO
BigNumDigit BnnShiftLeft (BigNum mm, BigNumLength ml, int nbits)
#else /* _NO_PROTO */
BigNumDigit BnnShiftLeft ( mm, ml, nbits)
BigNum mm; BigNumLength ml; int nbits;
#endif /* _NO_PROTO */
/*
* Shifts M left by "nbits", filling with 0s.
* Returns the leftmost "nbits" of M in a digit.
* Assumes 0 <= nbits < BN_DIGIT_SIZE.
*/
{
register BigNumDigit res = 0, save;
int rnbits;
if (nbits != 0)
{
rnbits = BN_DIGIT_SIZE - nbits;
while (ml-- > 0)
{
save = *mm;
*mm++ = (save << nbits) | res;
res = save >> rnbits;
}
}
return (res);
}
/***************************************/
/* */
#ifndef _NO_PROTO
BigNumDigit BnnShiftRight (BigNum mm, BigNumLength ml, int nbits)
#else /* _NO_PROTO */
BigNumDigit BnnShiftRight ( mm, ml, nbits)
BigNum mm; BigNumLength ml; int nbits;
#endif /* _NO_PROTO */
/*
* Shifts M right by "nbits", filling with 0s.
* Returns the rightmost "nbits" of M in a digit.
* Assumes 0 <= nbits < BN_DIGIT_SIZE.
*/
{
register BigNumDigit res = 0, save;
int lnbits;
if (nbits != 0)
{
mm += ml;
lnbits = BN_DIGIT_SIZE - nbits;
while (ml-- > 0)
{
save = *(--mm);
*mm = (save >> nbits) | res;
res = save << lnbits;
}
}
return (res);
}
/***************************************/
/* */
/******************* Additions **************************/
#ifndef _NO_PROTO
BigNumCarry BnnAddCarry (BigNum nn, BigNumLength nl, BigNumCarry carryin)
#else /* _NO_PROTO */
BigNumCarry BnnAddCarry ( nn, nl, carryin)
BigNum nn; BigNumLength nl; BigNumCarry carryin;
#endif /* _NO_PROTO */
/*
* Performs the sum N + CarryIn => N.
* Returns the CarryOut.
*/
{
if (carryin == 0)
return (0);
if (nl == 0)
return (1);
while (nl > 0 && !(++(*nn++)))
nl--;
return (nl > 0 ? 0 : 1);
}
/***************************************/
/* */
#ifndef _NO_PROTO
BigNumCarry BnnAdd (BigNum mm, BigNumLength ml, BigNum nn, BigNumLength nl, BigNumCarry carryin)
#else /* _NO_PROTO */
BigNumCarry BnnAdd ( mm, ml, nn, nl, carryin)
BigNum mm; BigNumLength ml; BigNum nn; BigNumLength nl; BigNumCarry carryin;
#endif /* _NO_PROTO */
/*
* Performs the sum M + N + CarryIn => M.
* Returns the CarryOut.
* Assumes Size(M) >= Size(N).
*/
{
register BigNumProduct c = carryin;
register BigNumProduct save;
ml -= nl;
while (nl > 0)
{
save = *mm;
c += save;
if (c < save)
{
*(mm++) = *(nn++);
c = 1;
}
else
{
save = *(nn++);
c += save;
*(mm++) = c;
c = (c < save) ? 1 : 0;
}
nl--;
}
return (BnnAddCarry (mm, ml, (BigNumCarry) c));
}
/***************************************/
/* */
/****************** Subtraction *************************/
#ifndef _NO_PROTO
BigNumCarry BnnSubtractBorrow (BigNum nn, BigNumLength nl, BigNumCarry carryin)
#else /* _NO_PROTO */
BigNumCarry BnnSubtractBorrow ( nn, nl, carryin)
BigNum nn; BigNumLength nl; BigNumCarry carryin;
#endif /* _NO_PROTO */
/*
* Performs the difference N + CarryIn - 1 => N.
* Returns the CarryOut.
*/
{
if (carryin == 1)
return (1);
if (nl == 0)
return (0);
while (nl > 0 && !((*nn++)--))
nl--;
return (nl > 0 ? 1 : 0);
}
/***************************************/
/* */
#ifndef _NO_PROTO
BigNumCarry BnnSubtract (BigNum mm, BigNumLength ml, BigNum nn, BigNumLength nl, BigNumCarry carryin)
#else /* _NO_PROTO */
BigNumCarry BnnSubtract ( mm, ml, nn, nl, carryin)
BigNum mm; BigNumLength ml; BigNum nn; BigNumLength nl; BigNumCarry carryin;
#endif /* _NO_PROTO */
/*
* Performs the difference M - N + CarryIn - 1 => M.
* Returns the CarryOut.
* Assumes Size(M) >= Size(N).
*/
{
register BigNumProduct c = carryin;
register BigNumDigit invn;
register BigNumProduct save;
ml -= nl;
while (nl > 0)
{
save = *mm;
invn = *(nn++) ^ -1;
c += save;
if (c < save)
{
*(mm++) = invn;
c = 1;
}
else
{
c += invn;
*(mm++) = c;
c = (c < invn) ? 1 : 0;
}
nl--;
}
return (BnnSubtractBorrow (mm, ml, (BigNumCarry) c)); }
/***************************************/
/* */
/***************** Multiplication ************************/
#ifndef _NO_PROTO
BigNumCarry BnnMultiplyDigit (BigNum pp, BigNumLength pl, BigNum mm, BigNumLength ml, BigNumDigit d)
#else /* _NO_PROTO */
BigNumCarry BnnMultiplyDigit ( pp, pl, mm, ml, d)
BigNum pp; BigNumLength pl; BigNum mm; BigNumLength ml; BigNumDigit d;
#endif /* _NO_PROTO */
/*
* Performs the product:
* Q = P + M * d
* BB = BBase(P)
* Q mod BB => P
* Q div BB => CarryOut
* Returns the CarryOut.
* Assumes Size(P) >= Size(M) + 1.
*/
{
register BigNumProduct c = 0;
if (d == 0)
return (0);
if (d == 1)
return (BnnAdd (pp, pl, mm, ml, (BigNumCarry) 0));
pl -= ml;
{
/* help for stupid compilers--may actually be counter
productive on pipelined machines with decent register allocation!! */
#define m_digit X0
#define X3 Lm
#define X1 Hm
register BigNumDigit Lm, Hm, Ld, Hd, X0, X2 /*, X1, X3 */;
Ld = d & ((((BigNumDigit)1) << (BN_DIGIT_SIZE / 2)) -1);
Hd = d >> (BN_DIGIT_SIZE / 2);
while (ml != 0)
{
ml--;
m_digit = *mm++;
Lm = m_digit & ((((BigNumDigit)1) << (BN_DIGIT_SIZE / 2)) -1);
Hm = m_digit >> (BN_DIGIT_SIZE / 2);
X0 = Ld * Lm;
X2 = Hd * Lm;
X3 = Hd * Hm;
X1 = Ld * Hm;
if ((c += X0) < X0) X3++;
if ((X1 += X2) < X2) X3 += (((BigNumDigit)1)<<(BN_DIGIT_SIZE / 2));
X3 += (X1 >> (BN_DIGIT_SIZE / 2));
X1 <<= (BN_DIGIT_SIZE / 2);
if ((c += X1) < X1) X3++;
if ((*pp += c) < c) X3++;
pp++;
c = X3;
#undef m_digit
#undef X1
#undef X3
}
X0 = *pp;
c += X0;
*(pp++) = c;
if (c >= X0)
return (0);
pl--;
while (pl != 0 && !(++(*pp++)))
pl--;
return (pl != 0 ? 0 : 1);
}
}
#ifdef mips
#ifndef _NO_PROTO
BigNumCarry BnnMultiply2Digit (BigNum pp, BigNumLength pl, BigNum mm, BigNumLength ml, BigNumDigit d0, BigNumDigit d1)
#else /* _NO_PROTO */
BigNumCarry BnnMultiply2Digit ( pp, pl, mm, ml, d0, d1)
BigNum pp; BigNumLength pl; BigNum mm; BigNumLength ml; BigNumDigit d0; BigNumDigit d1;
#endif /* _NO_PROTO */
/*
* Provided for compatibility with mips assembler implementation.
* Performs the product:
* Q = P + M * d0_d1
* BB = BBase(P)
* Q mod BB => P
* Q div BB => CarryOut
* Returns the CarryOut.
* Assumes Size(P) >= Size(M) + 1.
*/
{
return
BnnMultiplyDigit (pp, pl, mm, ml, d0)
+ BnnMultiplyDigit (pp+1, pl-1, mm, ml, d1);
}
#endif /* mips */
/***************************************/
/* */
/********************** Division *************************/
/* xh:xl -= yh:yl */
#define SUB(xh,xl,yh,yl) if (yl > xl) {xl -= yl; xh -= yh + 1;}\
else {xl -= yl; xh -= yh;}
#define LOW(x) (x & ((((BigNumDigit)1) << (BN_DIGIT_SIZE / 2)) -1))
#define HIGH(x) (x >> (BN_DIGIT_SIZE / 2))
#define L2H(x) (x << (BN_DIGIT_SIZE / 2))
#ifndef _NO_PROTO
BigNumDigit BnnDivideDigit (BigNum qq, BigNum nn, BigNumLength nl, BigNumDigit d)
#else /* _NO_PROTO */
BigNumDigit BnnDivideDigit ( qq, nn, nl, d)
BigNum qq; BigNum nn; BigNumLength nl; BigNumDigit d;
#endif /* _NO_PROTO */
/* Performs the quotient: N div d => Q
* Returns R = N mod d
* Assumes leading digit of N < d, and d > 0.
*/
{
{
int k;
BigNumLength orig_nl;
BigNumDigit rh; /* Two halves of current remainder */
BigNumDigit rl; /* Correspond to quad above */
register BigNumDigit qa; /* Current appr. to quotient */
register BigNumDigit ph, pl; /* product of c and qa */
BigNumDigit ch, cl, prev_qq;
/* Normalize divisor */
k = BnnNumLeadingZeroBitsInDigit (d);
if (k != 0)
{
prev_qq = qq[-1];
orig_nl = nl;
d <<= k;
BnnShiftLeft (nn, nl, k);
}
nn += nl;
nl--;
qq += nl;
ch = HIGH (d);
cl = LOW (d);
rl = *(--nn);
while (nl != 0)
{
nl--;
rh = rl;
rl = *(--nn);
qa = rh / ch; /* appr. quotient */
/* Compute ph, pl */
pl = cl * qa;
ph = ch * qa;
ph += HIGH (pl);
pl = L2H (pl);
/* While ph:pl > rh:rl, decrement qa, adjust qh:ql */
while (ph > rh || ph == rh && pl > rl)
{
qa--;
SUB (ph, pl, ch, L2H (cl));
}
SUB (rh, rl, ph, pl);
/* Top half of quotient is correct; save it */
*(--qq) = L2H (qa);
qa = (L2H (rh) | HIGH (rl)) / ch;
/* Approx low half of q */
/* Compute ph, pl, again */
pl = cl * qa;
ph = ch * qa;
ph += HIGH (pl);
pl = LOW (pl) | L2H (LOW (ph));
ph = HIGH (ph);
/* While ph:pl > rh:rl, decrement qa, adjust qh:ql */
while (ph > rh || ph == rh && pl > rl)
{
qa--;
SUB (ph, pl, 0, d);
}
/* Subtract ph:pl from rh:rl; we know rh will be 0 */
rl -= pl;
*qq |= qa;
}
/* Denormalize dividend */
if (k != 0) {
if((qq > nn) && (qq < &nn[orig_nl])) {
/* Overlap between qq and nn. Care of *qq! */
orig_nl = (qq - nn);
BnnShiftRight (nn, orig_nl, k);
nn[orig_nl - 1] = prev_qq;
} else if(qq == nn) {
BnnShiftRight(&nn[orig_nl - 1], 1, k);
} else {
BnnShiftRight (nn, orig_nl, k);
} }
return (rl >> k);
}
}
/***************************************/