BigCardinalsImpl.mesa
Copyright Ó 1985, 1987, 1988, 1991 by Xerox Corporation. All rights reserved.
Dan Greene, July 2, 1986 2:45:58 pm PDT
Swinehart, July 5, 1985 2:02:26 pm PDT
Russ Atkinson (RRA) September 29, 1988 2:51:09 pm PDT
Michael Plass, December 12, 1988 1:52:27 pm PST
DIRECTORY
Basics USING [CompareCard, Comparison, HighHalf, LongNumber, LowHalf],
BigCardinals USING [BigCARD, BigCARDRec, Zero],
Convert USING [RopeFromCard],
IO USING [PutChar, RopeFromROS, ROS, STREAM],
Basics16 USING [BITAND, BITNOT, BITOR, BITSHIFT, LongDiv, LongDivMod, LongMult],
Random USING [ChooseInt, Create, RandomStream],
Real USING [FScale],
Rope USING [ActionType, Concat, Equal, Fetch, FromChar, FromProc, Map, ROPE, Size, SkipOver, Substr];
BigCardinalsImpl:
CEDAR
PROGRAM
IMPORTS Basics, Convert, IO, Basics16, Random, Real, Rope
EXPORTS BigCardinals = BEGIN OPEN BigCardinals, Basics16;
ROPE: TYPE = Rope.ROPE;
STREAM: TYPE = IO.STREAM;
highBit: Digit = (LAST[Digit]-LAST[Digit]/2);
lastPlus: CARD = CARD32[LAST[Digit]]+1; -- 200000B
bitsPerChar: NAT = BITS[CHAR];
bitsPerDigit: NAT = BITS[Digit];
bytesPerDigit: NAT = BYTES[Digit];
useNewRopeEncoding: BOOL ¬ TRUE;
Most of the algorithms used below can be found in Knuth Volume 2.
BigFail: PUBLIC ERROR [subclass: ATOM ¬ $Unspecified] = CODE;
Basic Arithmetic
BigAdd:
PUBLIC
PROC [in1, in2: BigCARD]
RETURNS [out: BigCARD] ~ {
smaller: BigCARD;
addResults: CARD;
carry: Digit ¬ 0;
Check[$BigAdd, in1, in2];
IF in1.size < in2.size THEN {smaller ¬ in1; out ¬ in2} ELSE {smaller ¬ in2; out ¬ in1};
out ¬ BigCopy[out];
IF (out.size = 0) OR (smaller.size = 0) THEN RETURN;
FOR index:
NAT
IN [0..smaller.size)
DO
addResults ¬ CARD32[out.contents[index]] + CARD32[smaller.contents[index]] + carry;
out.contents[index] ¬ Basics.LowHalf[addResults];
carry ¬ Basics.HighHalf[addResults];
ENDLOOP;
FOR index:
NAT
IN [smaller.size..out.size)
DO
addResults ¬ CARD32[out.contents[index]] + carry;
out.contents[index] ¬ Basics.LowHalf[addResults];
carry ¬ Basics.HighHalf[addResults];
ENDLOOP;
IF carry > 0
THEN {
IF out.size >= out.contents.max THEN out ¬ BigCopy[out]; --Expand `out' if necessary
out.contents[out.size] ¬ carry;
out.size ¬ out.size + 1;
};
CheckResult[$BigAdd, out];
};
BigSubtract:
PUBLIC
PROC [large, small: BigCARD, reuse:
BOOL ¬
FALSE]
RETURNS [out: BigCARD] ~ {
borrow: CARD ¬ 0;
subResults: CARD ¬ 0;
Check[$BigSubtract, large, small];
IF large.size < small.size THEN ERROR BigFail[$NegativeResult];
IF reuse THEN out ¬ large ELSE out ¬ BigCopy[large];
In abuse of type, borrow is conceptually an INT in what follows.
FOR index:
NAT
IN [0..small.size)
DO
subResults ¬ CARD32[out.contents[index]] - CARD32[small.contents[index]] + borrow;
out.contents[index] ¬ Basics.LowHalf[subResults];
borrow ¬
LOOPHOLE[
INT32[
LOOPHOLE[Basics.HighHalf[subResults],
INT16]],
CARD];
The LONG and LOOPHOLEs extend the sign of borrow.
ENDLOOP;
FOR index:
NAT
IN [small.size..out.size)
DO
subResults ¬ CARD32[out.contents[index]] + borrow;
out.contents[index] ¬ Basics.LowHalf[subResults];
borrow ¬ LOOPHOLE[INT32[LOOPHOLE[Basics.HighHalf[subResults], INT16]], CARD];
ENDLOOP;
out ¬ StripLeadingZeros[out];
CheckResult[$BigSubtract, out];
};
StripLeadingZeros:
PROC [in: BigCARD]
RETURNS [out: BigCARD] ~ {
out ¬ in;
FOR index:
NAT
DECREASING
IN [0..in.size)
DO
IF out.contents[index] > 0 THEN RETURN;
out.size ¬ out.size - 1;
ENDLOOP;
out.contents ¬ NIL;
};
BigMultiply:
PUBLIC
PROC [in1, in2: BigCARD]
RETURNS [out: BigCARD ¬ []] ~ {
carry: Digit ¬ 0;
Check[$BigMultiply, in1, in2];
IF in1.size = 0 OR in2.size = 0 THEN RETURN;
out.size ¬ in1.size + in2.size;
out.contents ¬ NEW[BigCARDRec[out.size]];
FOR index1:
NAT
IN [0..in1.size)
DO
FOR index2:
NAT
IN [0..in2.size)
DO
product: CARD ¬ LongMult[in1.contents[index1], in2.contents[index2]]
+ CARD32[out.contents[index1+index2]]+carry;
out.contents[index1+index2] ¬ Basics.LowHalf[product];
carry ¬ Basics.HighHalf[product];
ENDLOOP;
out.contents[index1+in2.size] ¬ carry;
carry ¬ 0;
ENDLOOP;
out ¬ StripLeadingZeros[out];
CheckResult[$BigMultiply, out];
};
BigDivMod:
PUBLIC
PROC [num, den: BigCARD, reuse:
BOOL ¬
FALSE]
RETURNS [quo, rem: BigCARD] ~ {
highDen, nextDen: Digit;
smallNum, readyToShift: CARD;
carry, subtract, borrow, multResults, subResults, addResults: CARD;
normalize, remainder: Digit;
oldNumSize: NAT ¬ num.size;
leftShiftedDigit: Basics.LongNumber ¬ [pair[hi: 0, lo: 0]];
Check[$BigDivMod, num, den];
IF den.size=0 THEN ERROR BigFail[$DivideByZero];
IF den.size=1
THEN {
[quo, remainder] ¬ DivideByDigit[num, den.contents[0], reuse];
rem ¬ BigFromSmall[remainder];
Check[$BigDivModResult, quo, rem];
RETURN;
};
IF NOT reuse THEN num ¬ BigCopy[num];
IF den.size > num.size
THEN {rem ¬ num; Check[$BigDivModResult, quo, rem];
RETURN};
normalize ¬ den.contents[den.size - 1];
CheckLittle[normalize];
SELECT normalize
FROM
0 => ERROR;
>= LAST[Digit]/2 => normalize ¬ 1;
ENDCASE => {
normalize ¬ LongDiv[lastPlus, normalize+1];
CheckLittle[normalize];
num ¬ MultiplyByDigit[num, normalize, TRUE];
};
IF num.size = oldNumSize
THEN {
IF num.size >= num.contents.max THEN num ¬ BigCopy[num];
num.contents[num.size] ¬ 0;
num.size ¬ num.size + 1;
};
IF normalize # 1 THEN den ¬ MultiplyByDigit[den, normalize, reuse];
highDen ¬ den.contents[den.size - 1];
nextDen ¬ den.contents[den.size - 2];
quo.size ¬ num.size - den.size;
quo.contents ¬ NEW[BigCARDRec[quo.size]];
FOR quoIndex:
NAT
DECREASING
IN [0..quo.size)
DO
Try to guess the quotient using the leading two digits of `num' and the first digit of `den'.
guess: Digit ¬ Digit.LAST;
leftShiftedDigit.hi ¬ num.contents[quoIndex + den.size];
smallNum ¬ leftShiftedDigit.card + num.contents[quoIndex + den.size - 1];
IF highDen # leftShiftedDigit.hi THEN guess ¬ LongDiv[smallNum, highDen];
Adjust the guess by using one more digit each of `num' and `den'.
WHILE
TRUE
DO
readyToShift ¬ smallNum - LongMult[guess, highDen];
IF readyToShift >= lastPlus THEN EXIT;
leftShiftedDigit.hi ¬ readyToShift;
IF leftShiftedDigit.card+num.contents[quoIndex + den.size - 2]
>= LongMult[guess, nextDen] THEN EXIT;
guess ¬ guess - 1;
ENDLOOP;
Test the guess with the whole of `den'. This is almost always correct, so den * guess is actually subtracted from `num'.
carry ¬ borrow ¬ 0;
FOR denIndex:
NAT
IN [0..den.size)
DO
multResults ¬ LongMult[den.contents[denIndex], guess] + carry;
carry ¬ Basics.HighHalf[multResults];
subtract ¬ Basics.LowHalf[multResults];
Borrow is treated as an INT, and sign extended with LOOPHOLEs below.
subResults ¬ CARD32[num.contents[quoIndex + denIndex]] - subtract + borrow;
CheckLittle[quoIndex];
CheckLittle[denIndex];
num.contents[quoIndex + denIndex] ¬ Basics.LowHalf[subResults];
CheckLittle[num.contents[quoIndex + denIndex]];
borrow ¬ LOOPHOLE[INT32[LOOPHOLE[Basics.HighHalf[subResults], INT16]], CARD];
ENDLOOP;
In rare cases guess is off by one, so `den' is added back to `num'.
IF carry - borrow > num.contents[quoIndex + den.size]
THEN {
carry ¬ 0;
FOR denIndex:
NAT
IN [0..den.size)
DO
addResults ¬ CARD32[den.contents[denIndex]]
+ num.contents[quoIndex + denIndex] + carry;
carry ¬ Basics.HighHalf[addResults];
CheckLittle[quoIndex];
CheckLittle[denIndex];
num.contents[quoIndex + denIndex] ¬ Basics.LowHalf[addResults];
CheckLittle[num.contents[quoIndex + denIndex]];
ENDLOOP;
guess ¬ guess - 1;
};
CheckLittle[quoIndex];
CheckLittle[guess];
quo.contents[quoIndex] ¬ guess;
ENDLOOP;
quo ¬ StripLeadingZeros[quo];
CheckLittle[num.size];
num.size ¬ den.size;
[rem, remainder] ¬ DivideByDigit[num, normalize, TRUE];
Check[$BigDivModResult, quo, rem];
};
Conversion Routines
BigFromDecimalRope:
PUBLIC
PROC [in: Rope.
ROPE]
RETURNS [out: BigCARD] ~ {
action: Rope.ActionType = {
[c: CHAR] RETURNS [quit: BOOL ← FALSE]
IF pos = 0
THEN
SELECT c
FROM
'0, ' => RETURN;
ENDCASE;
SELECT c
FROM
IN ['0..'9] => {
pos ¬ pos + 1;
chunk ¬ chunk * 10 + CARDINAL[c-'0];
IF (Basics.LowHalf[pos]
MOD 4) = 0
THEN {
out ¬ MultiplyByDigit[out, 10000, TRUE];
out ¬ AddDigit[out, chunk, TRUE];
chunk ¬ 0;
};
};
ENDCASE => BigFail[$InvalidDecimalNumber];
};
pos: CARD ¬ 0;
chunk: Digit ¬ 0;
[] ¬ Rope.Map[base: in, action: action];
SELECT (Basics.LowHalf[pos]
MOD 4)
FROM
0 => RETURN;
1 => out ¬ MultiplyByDigit[out, 10, TRUE];
2 => out ¬ MultiplyByDigit[out, 100, TRUE];
3 => out ¬ MultiplyByDigit[out, 1000, TRUE];
ENDCASE => ERROR;
out ¬ AddDigit[out, chunk, TRUE];
};
StoreBCByte:
PROC [bc: BigCardinals.BigCARD, i:
CARD, byte:
BYTE] ~ {
q: CARDINAL ~ i / bytesPerDigit;
r: CARDINAL ~ i MOD bytesPerDigit;
bc.contents[q] ¬ Basics16.BITOR[Basics16.BITAND[Basics16.BITNOT[Basics16.BITSHIFT[255, (8*r)]], bc.contents[q]], Basics16.BITSHIFT[byte, (8*r)]];
};
BigFromBinaryRope:
PUBLIC PROC [in:
ROPE]
RETURNS [BigCardinals.BigCARD] ~ {
nBytes: CARD ~ Rope.Size[in];
IF nBytes = 0
THEN { RETURN [BigCardinals.Zero] }
ELSE {
nDigits: CARDINAL ~ (nBytes+(bytesPerDigit-1))/bytesPerDigit;
bc: BigCardinals.BigCARD ¬ [size: nDigits, contents: NEW[BigCardinals.BigCARDRec[nDigits]]];
index: CARD ¬ nBytes;
Action:
PROC [c:
CHAR]
RETURNS [quit:
BOOL ¬
FALSE] ~ {
StoreBCByte[bc, index ¬ index - 1, ORD[c]]
};
bc.contents[nDigits-1] ¬ 0;
[] ¬ Rope.Map[base: in, action: Action];
RETURN [bc]
};
};
FetchBCByte:
PROC [in: BigCardinals.BigCARD, i:
CARD]
RETURNS [
BYTE] ~ {
q: NAT ~ i / bytesPerDigit;
r: NAT ~ i MOD bytesPerDigit;
RETURN [VAL[Basics16.BITSHIFT[in.contents[q], -(8*r)] MOD 256]]
};
BigToBinaryRope:
PUBLIC
PROC [in: BigCardinals.BigCARD]
RETURNS [
ROPE] ~ {
length: CARD ¬ CARD32[in.size]*bytesPerDigit;
WHILE length # 0 AND FetchBCByte[in, length-1] = 0 DO length ¬ length - 1 ENDLOOP;
{ index:
CARD ¬ length;
P: PROC RETURNS [CHAR] ~ { RETURN [VAL[FetchBCByte[in, index ¬ index - 1]]] };
RETURN [Rope.FromProc[len: length, p: P]]
};
};
BigFromRope:
PUBLIC
PROC [in: Rope.
ROPE]
RETURNS [out: BigCARD ¬ []] ~ {
IF useNewRopeEncoding
THEN {
len: INT ¬ Rope.Size[in];
IF len # 0
THEN {
pos: INT ¬ Rope.SkipOver[in, 0, "\000"];
IF pos < len
THEN {
bits: INT ¬ (len-pos) * BITS[CHAR];
outlen: NAT ¬ LongDiv[bits+bitsPerDigit-1, bitsPerDigit];
index: NAT ¬ outlen-1;
action: Rope.ActionType = {
[c: CHAR] RETURNS [quit: BOOL ← FALSE]
w: Digit ¬ Basics16.BITSHIFT[ORD[c], shift ¬ shift - BITS[CHAR]];
out.contents[index] ¬ out.contents[index] + w;
IF shift = 0 AND index # 0 THEN {shift ¬ bitsPerDigit; index ¬ index - 1};
};
delta: CARD ¬ len - pos;
shift: NAT ¬ bitsPerDigit - (delta MOD bytesPerDigit)*NAT[BITS[CHAR]];
out.contents ¬ NEW[BigCARDRec[outlen]];
out.size ¬ outlen;
[] ¬ Rope.Map[in, pos, delta, action];
out ¬ Shift[out, pos+1];
};
IF pos # 0
THEN {
There were leading zero chars, so add in the trailing 1 bits
ones: BigCARD ¬ BigSubtract[Shift[BigFromSmall[1], pos], BigFromSmall[1], TRUE];
out ¬ BigAdd[out, ones];
};
};
}
ELSE {
x, y: INT ¬ 0;
trailingOnes: BigCARD ¬ BigFromSmall[1];
j: INT ¬ Rope.Size[in] - 1;
i: NAT;
First count (x) and strip off trailing nulls.
UNTIL j < 0 OR in.Fetch[j] ~= '\000 DO j ¬ j-1; x ¬ x + 1; ENDLOOP;
in ¬ Rope.Substr[in, 0, j+1];
size of bigcardinal is ceiling((1/2)*number of non-trailing-null chars in rope).
out.size ¬ (j + 2)/2;
out.contents ¬ NEW[BigCARDRec[out.size]];
IF ~BigZero[out]
THEN {
i ¬ j ¬ 0;
UNTIL j > Rope.Size[in] - 3
DO
out.contents[i] ¬ (in.Fetch[j] - 0C);
j ¬ j + 1;
out.contents[i] ¬ Basics16.BITOR[out.contents[i], (Rope.Fetch[in, j] - 0C)*256];
j ¬ j + 1;
i ¬ i + 1;
ENDLOOP;
IF j = Rope.Size[in]-2
THEN out.contents[out.size-1] ¬ (Rope.Fetch[in, j+1] - 0C)*256 + (Rope.Fetch[in, j] - 0C)
ELSE out.contents[out.size - 1] ¬ (Rope.Fetch[in, j] - 0C);
};
Now append one 0 and x 1's to the end of this bigcardinal.
UNTIL y = x
DO
trailingOnes ¬ BigMultiply[trailingOnes, BigFromSmall[2]];
y ¬ y + 1;
ENDLOOP;
out ¬ BigMultiply[out, BigMultiply[trailingOnes, BigFromSmall[2]]];
out ¬ BigAdd[out, BigSubtract[trailingOnes, BigFromSmall[1]]];
};
};
MultiplyByDigit:
PUBLIC
PROC [big: BigCARD, small: Digit, reuse:
BOOL ¬
FALSE]
RETURNS [out: BigCARD] ~ {
multResults: CARD;
carry: Digit ¬ 0;
CheckBigLittle[$MultiplyByDigit, big, small];
IF reuse THEN out ¬ big ELSE out ¬ BigCopy[big];
SELECT small
FROM
0 => {
FOR index:
NAT
IN [0..out.size)
DO
out.contents[index] ¬ 0;
ENDLOOP;
};
1 => {
carry ¬ 0; -- for placing counting breakpoints
};
2 => {
cc: Digit ¬ 0;
FOR index:
NAT
IN [0..out.size)
DO
c: CARD ¬ CARD[out.contents[index]]*2+cc;
out.contents[index] ¬ Basics.LowHalf[c];
cc ¬ Basics.HighHalf[c];
ENDLOOP;
carry ¬ cc;
};
4 => {
cc: Digit ¬ 0;
FOR index:
NAT
IN [0..out.size)
DO
c: CARD ¬ CARD[out.contents[index]]*4+cc;
out.contents[index] ¬ Basics.LowHalf[c];
cc ¬ Basics.HighHalf[c];
ENDLOOP;
carry ¬ cc;
};
ENDCASE => {
FOR index:
NAT
IN [0..out.size)
DO
multResults ¬ LongMult[out.contents[index], small] + carry;
out.contents[index] ¬ Basics.LowHalf[multResults];
carry ¬ Basics.HighHalf[multResults];
ENDLOOP;
};
IF carry > 0
THEN {
IF out.size >= out.contents.max THEN out ¬ BigCopy[out];
out.contents[out.size] ¬ carry;
out.size ¬ out.size + 1
};
CheckResult[$MultiplyByDigit, out];
};
AddDigit:
PROC [in: BigCARD, digit: Digit, reuse:
BOOL]
RETURNS [BigCARD] = {
CheckBigLittle[$AddDigit, in, digit];
IF in.size = 0 THEN RETURN [BigFromSmall[digit]];
IF NOT reuse THEN in ¬ BigCopy[in];
IF digit = 0 THEN RETURN [in];
FOR i:
NAT
IN [0..in.size)
DO
addResults: CARD ¬ CARD32[in.contents[i]] + digit;
in.contents[i] ¬ Basics.LowHalf[addResults];
digit ¬ Basics.HighHalf[addResults];
IF digit = 0 THEN RETURN [in];
ENDLOOP;
At this point we need room for carry out
in ¬ BigCopy[in];
in.contents[in.size] ¬ digit;
in.size ¬ in.size + 1;
CheckResult[$AddDigit, in];
RETURN [in];
};
Shift:
PUBLIC
PROC [bc: BigCARD, bits:
INT]
RETURNS [BigCARD] = {
new: BigCARD ¬ [];
Check[$Shift, bc, []];
IF bits = 0 THEN {CheckResult[$Shift, bc]; RETURN [bc]};
IF bc.size # 0
THEN {
bcBits: INT ¬ CARD32[bc.size]*bitsPerDigit + bits;
hiWord: Digit ¬ bc.contents[bc.size-1];
hiLim: Digit ¬ Digit.LAST - Digit.LAST/2;
IF hiWord = 0 THEN ERROR;
WHILE hiWord < hiLim
DO
hiWord ¬ hiWord + hiWord;
bcBits ¬ bcBits - 1;
ENDLOOP;
IF bcBits > 0
THEN {
loWord: Digit ¬ 0;
mod: INTEGER ¬ LOOPHOLE[bits, CARD] MOD bitsPerDigit;
loShift: INTEGER ¬ 0;
hiShift: INTEGER ¬ 0;
srcIndex: NAT ¬ 0;
dstIndex: NAT ¬ 0;
newSize: NAT ¬ (bcBits+bitsPerDigit-1)/bitsPerDigit;
new.contents ¬ NEW[BigCARDRec[newSize]];
new.size ¬ newSize;
SELECT bits
FROM
< 0 => {
Shift Right
div: Digit ¬ LongDiv[-bits, bitsPerDigit];
srcIndex ¬ div;
IF mod = 0
THEN {
The fast case, unshifted word transfers
FOR i:
NAT
IN [0..newSize)
DO
new.contents[dstIndex+i] ¬ bc.contents[srcIndex+i]
ENDLOOP;
}
ELSE {
An offset shift to the right
loShift ¬ mod-bitsPerDigit;
hiShift ¬ mod;
loWord ¬ bc.contents[srcIndex];
THROUGH [1..newSize)
DO
srcIndex ¬ srcIndex + 1;
hiWord ¬ bc.contents[srcIndex];
new.contents[dstIndex] ¬ Basics16.BITSHIFT[hiWord, hiShift]
+ Basics16.BITSHIFT[loWord, loShift];
dstIndex ¬ dstIndex + 1;
loWord ¬ hiWord;
ENDLOOP;
hiWord ¬ 0;
IF srcIndex+1 < bc.size THEN hiWord ¬ bc.contents[srcIndex+1];
new.contents[dstIndex] ¬ Basics16.BITSHIFT[hiWord, hiShift]
+ Basics16.BITSHIFT[loWord, loShift];
};
};
> 0 => {
Shift Left
div: Digit ¬ LongDiv[bits, bitsPerDigit];
nwords: NAT ¬ newSize - div + 1;
dstIndex ¬ div;
IF mod = 0
THEN {
The fast case, unshifted word transfers
FOR i:
NAT
IN [0..bc.size)
DO
new.contents[dstIndex+i] ¬ bc.contents[srcIndex+i]
ENDLOOP;
}
ELSE {
An offset shift to the left
loShift ¬ mod-bitsPerDigit;
hiShift ¬ mod;
loWord ¬ 0;
THROUGH [0..bc.size)
DO
hiWord ¬ bc.contents[srcIndex];
srcIndex ¬ srcIndex + 1;
new.contents[dstIndex] ¬ Basics16.BITSHIFT[hiWord, hiShift]
+ Basics16.BITSHIFT[loWord, loShift];
dstIndex ¬ dstIndex + 1;
loWord ¬ hiWord;
ENDLOOP;
loWord ¬ Basics16.BITSHIFT[loWord, loShift];
IF loWord # 0 THEN new.contents[dstIndex] ¬ loWord;
};
};
ENDCASE => ERROR;
new ¬ StripLeadingZeros[new];
};
};
CheckResult[$Shift, new];
RETURN [new];
};
MulMod:
PROC [x, y, mod: BigCARD]
RETURNS [bc: BigCARD ¬ []] = {
Check[$MulMod, x, y];
IF x.size # 0
AND y.size # 0
AND mod.size # 0
THEN {
SELECT BigCompare[x, mod]
FROM
equal => RETURN;
greater => {x ¬ BigDivMod[x, mod].rem; IF x.size = 0 THEN RETURN};
ENDCASE;
SELECT BigCompare[y, mod]
FROM
equal => RETURN;
greater => {y ¬ BigDivMod[y, mod].rem; IF y.size = 0 THEN RETURN};
ENDCASE;
IF x.size = 1
AND y.size = 1
THEN {
Trivial case can go much faster
prod: CARD ¬ LongMult[x.contents[0], y.contents[0]];
IF mod.size < 3 THEN prod ¬ prod MOD BigToCard[mod];
RETURN [BigFromCard[prod]];
};
bc ¬ BigMultiply[x, y];
SELECT
TRUE
FROM
bc.size < mod.size => RETURN;
bc.size = mod.size AND bc.contents[bc.size-1] < mod.contents[bc.size-1] => RETURN;
ENDCASE;
bc ¬ BigDivMod[bc, BigCopy[mod], TRUE].rem;
};
CheckResult[$MulMod, bc];
};
TwoToTheNth:
PUBLIC
PROC [n:
NAT]
RETURNS [bc: BigCARD] = {
This procedure computes 2 ^ n.
mod: NAT ¬ n MOD bitsPerDigit;
words: NAT ¬ (n+bitsPerDigit) / bitsPerDigit;
bit: Digit ¬ Basics16.BITSHIFT[1, mod];
bc.size ¬ words;
bc.contents ¬ NEW[BigCARDRec[words]]; -- initialized to all zeros by the allocator?
bc.contents[words-1] ¬ bit;
};
TimesTwoToTheNth:
PUBLIC
PROC [in: BigCARD, n:
NAT]
RETURNS [BigCARD] = {
This procedure computes in * 2 ^ n.
RETURN [Shift[in, n]];
};
BigPowerOfBase:
PUBLIC
PROC [base: BigCARD, exponent:
NAT]
RETURNS [BigCARD] ~ {
power: BigCARD ¬ BigFromSmall[1];
IF exponent = 0 THEN RETURN [power];
-- Exponentiation by repeated squaring.
WHILE exponent > 1
DO
IF exponent MOD 2 = 1 THEN power ¬ BigMultiply[power, base];
exponent ¬ exponent / 2;
base ¬ BigMultiply[base, base];
ENDLOOP;
power ¬ BigMultiply[power, base];
RETURN [power];
};
BigPowerOfTen:
PUBLIC
PROC [exponent:
NAT]
RETURNS [BigCARD] ~ {
RETURN [BigPowerOfBase[BigFromSmall[10], exponent]];
};
BigToREAL:
PUBLIC
PROC [in: BigCARD]
RETURNS [
REAL] ~ {
real: REAL ¬ 0.0;
scale: INTEGER ¬ 0;
IF BigZero[in] THEN RETURN [real];
FOR index:
NAT
IN [0..in.size)
DO
real ¬ real + Real.FScale[REAL[in.contents[index]], scale];
scale ¬ scale + BITS[Digit];
ENDLOOP;
RETURN [real]
};
BigOne:
PUBLIC
PROC [in: BigCARD]
RETURNS [
BOOL] ~ {
Test if a BigCARD is the integer one.
IF in.size = 1 THEN RETURN [in.contents[0] = 1] ELSE RETURN [FALSE]
};
FirstOneBit:
PUBLIC PROC [bc: BigCARD]
RETURNS [exp:
INT ¬ -1] = {
IF bc.size > 0
THEN {
temp: Digit ¬ bc.contents[bc.size-1];
exp ¬ CARD[bc.size-1]*CARD[BITS[Digit]];
IF temp = 0 THEN ERROR;
DO
temp ¬ temp / 2;
IF temp = 0 THEN EXIT;
exp ¬ exp + 1;
ENDLOOP;
};
};
LowOrderZeros:
PUBLIC PROC [bc: BigCARD]
RETURNS [
INT] = {
IF bc.size > 0
THEN {
FOR i:
NAT
IN [0..bc.size)
DO
w: Digit ¬ bc.contents[i];
IF w # 0
THEN {
n: INT ¬ INT[i]*BITS[Digit];
WHILE (w MOD 16 = 0) DO w ¬ w / 16; n ¬ n + 4 ENDLOOP;
WHILE (w MOD 2 = 0) DO w ¬ w / 2; n ¬ n + 1 ENDLOOP;
RETURN [n]
};
ENDLOOP;
};
RETURN [-1]
};
BigFactorial:
PUBLIC
PROC [x:
NAT]
RETURNS [out: BigCARD] = {
Returns the factorial of x.
bc: BigCardinals.BigCARD ¬ BigFromSmall[x];
exp: NAT ¬ 0;
WHILE x > 1
DO
digit: Digit ¬ x ¬ x-1;
WHILE (digit MOD 2) = 0 DO digit ¬ digit / 2; exp ¬ exp + 1; ENDLOOP;
IF digit # 1 THEN bc ¬ MultiplyByDigit[bc, digit, TRUE];
ENDLOOP;
IF exp # 0
THEN bc ¬ TimesTwoToTheNth[bc, exp];
Bring in the factors of 2.
RETURN [bc];
};
BigToDecimalRope:
PUBLIC
PROC [in: BigCARD, reuse:
BOOL ¬
FALSE]
RETURNS [out: Rope.
ROPE ¬
NIL] ~ {
SELECT in.size
FROM
0 => RETURN ["0"];
1 => RETURN [Convert.RopeFromCard[in.contents[0]]];
ENDCASE => {
chunk: Digit;
ros: STREAM ¬ IO.ROS[];
pos: INT ¬ 0;
fetch:
PROC
RETURNS [
CHAR] = {
RETURN [Rope.Fetch[out, pos ¬ pos - 1]];
};
IF NOT reuse THEN in ¬ BigCopy[in];
UNTIL in.size = 0
DO
[in, chunk] ¬ DivideByDigit[in, 10, TRUE];
IO.PutChar[ros, '0 + chunk];
ENDLOOP;
out ¬ IO.RopeFromROS[ros];
pos ¬ Rope.Size[out];
RETURN [Rope.FromProc[pos, fetch]];
};
};
BigToRope:
PUBLIC
PROC [in: BigCARD]
RETURNS [out: Rope.
ROPE ¬
NIL] ~ {
IF useNewRopeEncoding
THEN {
The # of low-order ones in the BigCARD determines the number of leading zero bytes in the rope.
IF in.size # 0
THEN {
trailingOnes: INT ¬ 0;
total: INT ¬ trailingOnes+1;
index: NAT ¬ 0;
sample: Digit ¬ 0;
shift: INTEGER ¬ 0;
eachChar:
PROC
RETURNS [c:
CHAR] = {
IF trailingOnes > 0 THEN {trailingOnes ¬ trailingOnes - 1; RETURN [0C]};
c ¬ VAL[Basics16.BITSHIFT[sample, shift] MOD 256];
IF shift = 0
THEN {
IF index = 0 THEN RETURN;
index ¬ index - 1;
sample ¬ in.contents[index];
shift ¬ -bitsPerDigit;
};
shift ¬ shift + BITS[CHAR];
};
FOR i:
NAT
IN [0..in.size)
DO
digit: Digit ¬ in.contents[i];
IF digit # Digit.
LAST
THEN {
WHILE (digit
MOD 2) = 1
DO
trailingOnes ¬ trailingOnes + 1;
digit ¬ digit / 2;
ENDLOOP;
EXIT;
};
trailingOnes ¬ trailingOnes + BITS[Digit];
ENDLOOP;
total ¬ trailingOnes;
in ¬ Shift[in, -1-trailingOnes];
IF in.size # 0
THEN {
chars: INT ¬ in.size * bytesPerDigit;
index ¬ in.size-1;
sample ¬ in.contents[index];
shift ¬ -bitsPerDigit;
WHILE shift # 0
DO
Scan over the leading zero bytes in the digit
shift ¬ shift + BITS[CHAR];
IF Basics.LowHalf[Basics16.BITSHIFT[sample, shift]] # 0 THEN EXIT;
chars ¬ chars - 1;
ENDLOOP;
total ¬ total + chars;
};
out ¬ Rope.FromProc[total, eachChar];
};
}
ELSE {
x, y: INT ¬ 0;
nonTrailingOnes: BigCARD;
i: INT;
trailingOnes: BigCARD ¬ BigFromSmall[1];
trailingNulls, nonTrailingNulls: Rope.ROPE ¬ NIL;
i ¬ 0;
UNTIL i = in.size
DO
mask: Digit ¬ 1B;
thisBit: Digit ¬ Basics16.BITAND[in.contents[i], mask];
The mask will become zero if and when we multiply 2^15 by 2.
WHILE mask > 0B
AND thisBit = mask
DO
x ¬ x + 1;
mask ¬ mask * 2;
thisBit ¬ Basics16.BITAND[in.contents[i], mask];
ENDLOOP;
The IF x = 0 terminates the UNTIL loop if there are no trailing ones.
IF x = 0 OR thisBit = 0 THEN i ¬ in.size ELSE i ¬ i + 1;
ENDLOOP;
UNTIL y = x
DO
trailingOnes ¬ BigMultiply[trailingOnes, BigFromSmall[2]];
trailingNulls ¬ Rope.Concat[trailingNulls, Rope.FromChar['\000]];
y ¬ y + 1;
ENDLOOP;
nonTrailingOnes ¬ BigSubtract[in, BigSubtract[trailingOnes, BigFromSmall[1]]];
nonTrailingOnes ¬ BigDivMod[nonTrailingOnes, BigMultiply[trailingOnes, BigFromSmall[2]]].quo;
IF ~BigZero[nonTrailingOnes]
THEN {
FOR i
IN [0 .. nonTrailingOnes.size - 1)
DO
nonTrailingNulls ¬ Rope.Concat[nonTrailingNulls, Rope.FromChar[nonTrailingOnes.contents[i] MOD 256 + 0C]];
nonTrailingNulls ¬ Rope.Concat[nonTrailingNulls, Rope.FromChar[nonTrailingOnes.contents[i] / 256+ 0C]];
ENDLOOP;
nonTrailingNulls ¬ Rope.Concat[nonTrailingNulls, Rope.FromChar[nonTrailingOnes.contents[nonTrailingOnes.size - 1] MOD 256 + 0C]];
IF nonTrailingOnes.contents[nonTrailingOnes.size - 1] > 257
THEN
nonTrailingNulls ¬ Rope.Concat[nonTrailingNulls, Rope.FromChar[nonTrailingOnes.contents[nonTrailingOnes.size -1] / 256 + 0C]];
};
out ¬ Rope.Concat[nonTrailingNulls, trailingNulls];
};
};
DivideByDigit:
PROC [num: BigCARD, den: Digit, reuse:
BOOL ¬
FALSE]
RETURNS [quo: BigCARD, rem: Digit ¬ 0] ~ {
quotient: Digit;
leftShiftedRem: Basics.LongNumber ¬ [pair[hi: 0, lo: 0]];
CheckBigLittle[$DivideByDigit, num, den];
IF reuse THEN quo ¬ num ELSE quo ¬ BigCopy[num];
SELECT den
FROM
2 => {
high: Digit ¬ 0;
FOR index:
NAT
DECREASING
IN [0..quo.size)
DO
c: Digit ¬ quo.contents[index];
rem ¬ c MOD 2;
quo.contents[index] ¬ c/2 + high;
high ¬ rem * highBit;
ENDLOOP;
};
4 => {
high: Digit ¬ 0;
FOR index:
NAT
DECREASING
IN [0..quo.size)
DO
c: Digit ¬ quo.contents[index];
rem ¬ c MOD 4;
quo.contents[index] ¬ c/4 + high;
high ¬ rem * (highBit/2);
ENDLOOP;
};
ENDCASE => {
FOR index:
NAT
DECREASING
IN [0..quo.size)
DO
leftShiftedRem.hi ¬ rem;
[quotient, rem] ¬ LongDivMod[quo.contents[index]+leftShiftedRem.card, den];
quo.contents[index] ¬ quotient;
ENDLOOP;
};
quo ¬ StripLeadingZeros[quo];
CheckResult[$DivideByDigit, quo];
};
BigFromSmall:
PUBLIC
PROC [in: Digit]
RETURNS [out: BigCARD ¬ []] ~ {
IF in # 0
THEN {
out.contents ¬ NEW[BigCARDRec[1]];
out.size ¬ 1;
out.contents[0] ¬ in;
};
};
BigToSmall:
PUBLIC
PROC [in: BigCARD]
RETURNS [out: Digit] ~ {
SELECT in.size
FROM
0 => RETURN [0];
1 => RETURN [in.contents[0]];
ENDCASE => ERROR BigFail[$TooLarge];
};
BigFromCard:
PUBLIC
PROC [in:
CARD]
RETURNS [out: BigCARD] ~ {
ln: Basics.LongNumber ¬ [card[in]];
SELECT
TRUE
FROM
ln.hi # 0 => {
out.size ¬ 2;
out.contents ¬ NEW[BigCARDRec[2]];
out.contents[0] ¬ ln.lo;
out.contents[1] ¬ ln.hi;
};
ln.lo # 0 => {
out.size ¬ 1;
out.contents ¬ NEW[BigCARDRec[1]];
out.contents[0] ¬ ln.lo;
};
ENDCASE => out.size ¬ 0;
};
BigToCard:
PUBLIC
PROC [in: BigCARD]
RETURNS [out:
CARD] ~ {
ln: Basics.LongNumber ¬ [pair[hi: 0, lo: 0]];
SELECT in.size
FROM
0 => {};
1 => ln.lo ¬ in.contents[0];
2 => {ln.lo ¬ in.contents[0]; ln.hi ¬ in.contents[1]};
ENDCASE => ERROR BigFail[$TooLarge];
RETURN [ln.card];
};
Exotic Functions
BigRandom:
PUBLIC
PROC [length:
NAT ¬ 0, max, min: BigCARD ¬ [0,
NIL] ]
RETURNS [out: BigCARD] ~ {
minMatters, maxMatters: BOOL ¬ TRUE;
lower: Digit ¬ 0;
upper: Digit ¬ 177777B;
IF length > 0
THEN {out.size ¬ length; maxMatters ¬ minMatters ¬ FALSE}
ELSE out.size ¬ max.size;
out.contents ¬ NEW[BigCARDRec[out.size]];
FOR index:
NAT
DECREASING
IN [0..out.size)
DO
IF maxMatters THEN upper ¬ max.contents[index];
IF minMatters AND index < min.size THEN lower ¬ min.contents[index];
out.contents[index] ¬ Random.ChooseInt[rs: markov, min: lower, max: upper];
IF out.contents[index] < upper THEN {maxMatters ¬ FALSE; upper ¬ 177777B};
IF out.contents[index] > lower THEN {minMatters ¬ FALSE; lower ¬ 0};
ENDLOOP;
out ¬ StripLeadingZeros[out];
};
BigExponentiate:
PUBLIC
PROC [base, exponent, modulus: BigCARD, reuse:
BOOL ¬
FALSE]
RETURNS [out: BigCARD] ~ {
Exponentiation by repeated squaring.
out ¬ BigFromSmall[1];
IF exponent.size # 0
THEN {
IF base.size >= modulus.size
THEN
base ¬ BigDivMod[base, modulus, FALSE].rem;
FOR index:
NAT
IN [0..exponent.size)
DO
temporary: Digit ¬ exponent.contents[index];
THROUGH [0..bitsPerDigit)
DO
IF (temporary
MOD 2) = 1
THEN
out ¬ MulMod[out, base, modulus];
temporary ¬ temporary / 2;
IF temporary = 0 AND index+1 = exponent.size THEN RETURN;
base ¬ MulMod[base, base, modulus];
ENDLOOP;
ENDLOOP;
};
};
BigGCD:
PUBLIC
PROC [in1, in2: BigCARD]
RETURNS [out: BigCARD] ~ {
oneLarger: BOOL ¬ (BigCompare[in1, in2] = greater);
DO
IF in1.size = 0 OR in2.size = 0 THEN EXIT;
IF oneLarger THEN in1 ¬ BigDivMod[in1, in2].rem ELSE in2 ¬ BigDivMod[in2, in1].rem;
oneLarger ¬ NOT oneLarger;
ENDLOOP;
out ¬ IF oneLarger THEN in1 ELSE in2;
};
BigInverse:
PUBLIC
PROC [in, modulus: BigCARD]
RETURNS [out: BigCARD] ~ {
This is 1/2 of an extended GCD algorithm. Since the intermediate results can be negative, signs must be maintained: sign1 for temp1 and sign2 for temp2.
sign: INTEGER;
sign1, sign2: INTEGER ¬ 1;
temp1: BigCARD ¬ BigFromSmall[0];
temp2: BigCARD ¬ BigFromSmall[1];
quo, addOrSub: BigCARD;
out ¬ modulus;
IF BigCompare[in, out] = greater THEN in ¬ BigDivMod[in, modulus].rem;
IF in.size # 0
THEN
DO
[quo, out] ¬ BigDivMod[out, in];
IF out.size = 0 THEN {out ¬ temp2; sign ¬ sign2; EXIT};
addOrSub ¬ BigMultiply[temp2, quo];
SELECT
TRUE
FROM
sign1 # sign2 => temp1 ¬ BigAdd[temp1, addOrSub];
BigCompare[temp1, addOrSub] = greater => temp1 ¬ BigSubtract[temp1, addOrSub];
ENDCASE => {sign1 ¬ -sign1; temp1 ¬ BigSubtract[addOrSub, temp1]};
[quo, in] ¬ BigDivMod[in, out];
IF in.size = 0 THEN {out ¬ temp1; sign ¬ sign1; EXIT};
addOrSub ¬ BigMultiply[temp1, quo];
SELECT
TRUE
FROM
sign1 # sign2 => temp2 ¬ BigAdd[temp2, addOrSub];
BigCompare[temp2, addOrSub] = greater => temp2 ¬ BigSubtract[temp2, addOrSub];
ENDCASE => {sign2 ¬ -sign2; temp2 ¬ BigSubtract[addOrSub, temp2]};
ENDLOOP;
IF sign = -1 THEN out ¬ BigSubtract[modulus, out];
};
BigPrimeTest:
PUBLIC
PROC [prime: BigCARD, accuracy:
NAT ¬ 20]
RETURNS [
BOOL] ~ {
one: BigCARD ¬ BigFromSmall[1];
two: BigCARD ¬ BigFromSmall[2];
primeMinusOne: BigCARD ¬ BigSubtract[prime, one];
q: BigCARD ¬ BigCopy[primeMinusOne];
x, y: BigCARD;
j: NAT;
k: NAT ¬ 0;
IF BigOdd[q] THEN RETURN [FALSE];
UNTIL BigOdd[q]
DO
q ¬ DivideByDigit[q, 2, TRUE].quo;
k ¬ k + 1;
ENDLOOP;
THROUGH [0..accuracy/2)
DO
On each iteration through the loop there is a 1/4 chance of falsely believing that `prime' is prime when it is not prime.
x ¬ BigRandom[max: primeMinusOne, min: two];
y ¬ BigExponentiate[x, q, prime];
j ¬ 0;
If `prime' is prime then y should eventually become 1 in the loop below. Furthermore, it should equal -1 MOD prime on the iteration before it becomes 1. If this scenario is not followed then `prime' is definately not a prime.
IF
NOT BigCompare[y, one] = equal
THEN
DO
IF BigCompare[y, primeMinusOne] = equal THEN EXIT;
j ¬ j+1;
y ¬ BigExponentiate[y, two, prime]; -- Now y = x (q 2j)
IF (j = k) OR BigCompare[y, one] = equal THEN RETURN [FALSE];
ENDLOOP;
ENDLOOP;
RETURN [TRUE];
};
Test and Time
Test:
PROC [invTimes:
NAT ¬ 10, primeTest:
BOOL ¬
TRUE]
RETURNS [
NAT] ~ {
one: BigCARD = BigFromSmall[1];
op1: BigCARD ¬ BigFromDecimalRope["222222222222222222"];
op2: BigCARD ¬ BigFromDecimalRope["18446744073709551556"];
op3: BigCARD ¬ BigFromDecimalRope["18446744073709551557"];
res1: BigCARD ¬ BigExponentiate[op1, op2, op3];
res2, res3: BigCARD;
hits, primes: NAT ¬ 0;
IF BigCompare[res1, one] # equal THEN ERROR;
FOR index:
NAT
IN [0..invTimes)
DO
op1 ¬ BigRandom[length: 50];
op2 ¬ BigRandom[length: 40];
res1 ¬ BigSubtract[op1, op2];
res2 ¬ BigAdd[res1, op2];
IF BigCompare[op1, res2] # equal THEN ERROR;
[res1, res2] ¬ BigDivMod[op1, op2];
res3 ¬ BigMultiply[op2, res1];
res3 ¬ BigAdd[res2, res3];
IF BigCompare[op1, res3] # equal THEN ERROR;
IF BigCompare[BigGCD[op1, op2], one] = equal
THEN {
hits ¬ hits + 1;
res1 ¬ BigInverse[op2, op1];
res2 ¬ BigMultiply[res1, op2];
res3 ¬ BigDivMod[res2, op1].rem;
IF BigCompare[one, res3] # equal THEN ERROR;
};
ENDLOOP;
IF primeTest
THEN {
A long arithmetic progression of primes found by Paul Pritchard of Cornell
op1 ¬ BigFromDecimalRope["8297644387"];
op2 ¬ BigFromDecimalRope["4180566390"];
WHILE BigPrimeTest[op1]
DO
primes ¬ primes + 1;
op1 ¬ BigAdd[op1, op2];
ENDLOOP;
IF NOT primes = 19 THEN ERROR;
};
RETURN [hits]
};
TestRopeEncoding:
PROC
[pow2Times:
NAT ¬ 100, randomRopes:
NAT ¬ 100, randomBigs:
NAT ¬ 100] = {
This test runs through some interesting cases for the rope encoding that we use.
First: we test 2**N and 2**N-1 for having invertible encodings as ropes.
Second: we test some random ropes for having invertible encodings as BigCARDs.
Third: we test some random BigCARDs for having invertible encodings as ropes.
one: BigCARD = BigFromSmall[1];
FOR i:
NAT
IN [0..pow2Times)
DO
x: BigCARD ¬ TwoToTheNth[i];
xm: BigCARD ¬ BigSubtract[x, one];
xr: ROPE ¬ BigToRope[x];
xmr: ROPE ¬ BigToRope[xm];
IF BigCompare[x, xm] # greater THEN ERROR;
IF BigCompare[x, BigFromRope[xr]] # equal THEN ERROR;
IF BigCompare[xm, BigFromRope[xmr]] # equal THEN ERROR;
IF BigCompare[x, BigFromDecimalRope[BigToDecimalRope[x]]] # equal THEN ERROR;
IF BigCompare[xm, BigFromDecimalRope[BigToDecimalRope[xm]]] # equal THEN ERROR;
ENDLOOP;
FOR i:
NAT
IN [0..randomRopes)
DO
eachChar:
PROC
RETURNS [
CHAR] = {
RETURN [VAL[BYTE[Random.ChooseInt[rs: markov, min: 0, max: 255]]]];
};
zeroChar:
PROC
RETURNS [
CHAR] = {
RETURN [0C];
};
r: ROPE ¬ Rope.FromProc[len: Random.ChooseInt[rs: markov, min: 0, max: 256], p: eachChar];
z: ROPE ¬ Rope.FromProc[len: Random.ChooseInt[rs: markov, min: 1, max: 256], p: zeroChar];
rz: ROPE ¬ Rope.Concat[z, r];
IF NOT Rope.Equal[r, BigToRope[BigFromRope[r]]] THEN ERROR;
IF NOT Rope.Equal[z, BigToRope[BigFromRope[z]]] THEN ERROR;
IF NOT Rope.Equal[rz, BigToRope[BigFromRope[rz]]] THEN ERROR;
ENDLOOP;
FOR i:
NAT
IN [0..randomBigs)
DO
bc: BigCARD ¬ BigRandom[
length: Random.ChooseInt[rs: markov, min: 0, max: 256],
max: [],
min: TwoToTheNth[Random.ChooseInt[rs: markov, min: 0, max: 1024]]];
bcr: ROPE ¬ BigToRope[bc];
bcbc: BigCARD ¬ BigFromRope[bcr];
IF BigCompare[bc, bcbc] # equal THEN ERROR;
ENDLOOP;
};
Time:
PROC []
RETURNS [Rope.
ROPE] ~ {
bc1: BigCARD ¬ BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234567220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"];
bc2: BigCARD ¬ BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234522220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"];
bc3: BigCARD ¬ BigFromDecimalRope["1234567890994567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"];
RETURN [ BigToDecimalRope[BigExponentiate[bc1, bc2, bc3 ] ] ];
};
markov: Random.RandomStream ¬ Random.Create[range: 0, seed: -1];