DIRECTORY
Basics USING [BITAND, BITOR, BITSHIFT, Comparison, HighByte, HighHalf, LongDiv, LongDivMod, LongMult, LongNumber, LowByte, LowHalf],
BigCardinals,
Convert USING [CardFromRope, RopeFromCard],
IO USING [PutChar, RopeFromROS, ROS, STREAM],
PrincOpsUtils USING [LongCopy],
Random USING [Create, ChooseInt, RandomStream],
Rope USING [Concat, Fetch, FromChar, FromProc, Length, ROPE, Substr];

BigCardinalsImpl: CEDAR PROGRAM
IMPORTS Basics, Convert, IO, PrincOpsUtils, Random, Rope
EXPORTS BigCardinals = {
OPEN BigCardinals;

ROPE: TYPE = Rope.ROPE;
STREAM: TYPE = IO.STREAM;
highBit: CARDINAL = (LAST[CARDINAL]-LAST[CARDINAL]/2);
lastPlus: CARD = LONG[LAST[CARDINAL]]+1; -- 200000B
bitsPerWord: CARDINAL = BITS[WORD];

BigFail: PUBLIC ERROR [subclass: ATOM _ $Unspecified] = CODE;


BigAdd: PUBLIC PROC [in1, in2: BigCARD] RETURNS [out: BigCARD] ~ {
smaller: BigCARD;
addResults: CARD; carry: CARDINAL _ 0;

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: CARDINAL IN [0..smaller.size) DO
addResults _ LONG[out.contents[index]] + LONG[smaller.contents[index]] + carry;
out.contents[index] _ Basics.LowHalf[addResults]; carry _ Basics.HighHalf[addResults];
ENDLOOP;
FOR index: CARDINAL IN [smaller.size..out.size) DO
addResults _ LONG[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;
};
};

BigSubtract: PUBLIC PROC [large, small: BigCARD, reuse: BOOL _ FALSE]
RETURNS [out: BigCARD] ~ {
borrow, subResults: CARD _ 0;

IF large.size < small.size THEN ERROR BigFail[$NegativeResult];
IF reuse THEN out _ large ELSE out _ BigCopy[large];


FOR index: CARDINAL IN [0..small.size) DO
subResults _ LONG[out.contents[index]] - LONG[small.contents[index]] + borrow;
out.contents[index] _ Basics.LowHalf[subResults];
borrow _ LOOPHOLE[LONG[LOOPHOLE[Basics.HighHalf[subResults], INTEGER]], CARD];
ENDLOOP;

FOR index: CARDINAL IN [small.size..out.size) DO
subResults _ LONG[out.contents[index]] + borrow;
out.contents[index] _ Basics.LowHalf[subResults];
borrow _ LOOPHOLE[LONG[LOOPHOLE[Basics.HighHalf[subResults], INTEGER]], CARD];
ENDLOOP;

out _ StripLeadingZeros[out];
};

StripLeadingZeros: PROC [in: BigCARD] RETURNS [out: BigCARD] ~ {
out _ in;
FOR index: CARDINAL 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: CARDINAL _ 0;

IF BigZero[in1] OR BigZero[in2] THEN RETURN;

out.size _ in1.size + in2.size;
out.contents _ NEW[BigCARDRec[out.size]];

FOR index1: CARDINAL IN [0..in1.size) DO
FOR index2: CARDINAL IN [0..in2.size) DO
product: CARD _ Basics.LongMult[in1.contents[index1], in2.contents[index2]]
+ LONG[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];
};

BigDivMod: PUBLIC PROC [num, den: BigCARD, reuse: BOOL _ FALSE]
RETURNS [quo, rem: BigCARD] ~ {
highDen, nextDen: CARDINAL;
smallNum, readyToShift: CARD;
carry, subtract, borrow, multResults, subResults, addResults: CARD;
guess, normalize, remainder: CARDINAL;
oldNumSize: CARDINAL _ num.size;
leftShiftedDigit: Basics.LongNumber _ [lc[0]];

IF den.size=0 THEN ERROR BigFail[$DivideByZero];
IF den.size=1 THEN {
[quo, remainder] _ DivideByDigit[num, den.contents[0], reuse];
rem _ BigFromSmall[remainder];
RETURN;
};

IF NOT reuse THEN num _ BigCopy[num];
IF den.size > num.size THEN {rem _ num; RETURN};

normalize _ den.contents[den.size - 1];
SELECT normalize FROM
0 => ERROR;
>= LAST[CARDINAL]/2 => normalize _ 1;
ENDCASE => {
normalize _ Basics.LongDiv[lastPlus, normalize+1];
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: CARDINAL DECREASING IN [0..quo.size) DO
leftShiftedDigit.hi _ num.contents[quoIndex + den.size];
smallNum _ leftShiftedDigit.lc + num.contents[quoIndex + den.size - 1]; 
guess _ IF highDen = leftShiftedDigit.hi THEN 177777B ELSE Basics.LongDiv[smallNum, highDen];

WHILE TRUE DO
readyToShift _ smallNum - Basics.LongMult[guess, highDen];
IF readyToShift >= 200000B THEN EXIT;
leftShiftedDigit.hi _ readyToShift;
IF leftShiftedDigit.lc+num.contents[quoIndex + den.size - 2]
>= Basics.LongMult[guess, nextDen] THEN EXIT;
guess _ guess - 1;
ENDLOOP;

carry _ borrow _ 0;
FOR denIndex: CARDINAL IN [0..den.size) DO
multResults _ Basics.LongMult[den.contents[denIndex], guess] + carry;
carry _ Basics.HighHalf[multResults]; subtract _ Basics.LowHalf[multResults];
subResults _ LONG[num.contents[quoIndex + denIndex]] - subtract + borrow;
num.contents[quoIndex + denIndex] _ Basics.LowHalf[subResults];
borrow _ LOOPHOLE[LONG[LOOPHOLE[Basics.HighHalf[subResults], INTEGER]], CARD];
ENDLOOP;

IF carry - borrow > num.contents[quoIndex + den.size] THEN {
carry _ 0;
FOR denIndex: CARDINAL IN [0..den.size) DO
addResults _ LONG[den.contents[denIndex]]
+ num.contents[quoIndex + denIndex] + carry;
carry _ Basics.HighHalf[addResults];
num.contents[quoIndex + denIndex] _ Basics.LowHalf[addResults];
ENDLOOP;
guess _ guess - 1;
};

quo.contents[quoIndex] _ guess;
ENDLOOP;

quo _ StripLeadingZeros[quo];

num.size _ den.size;
[rem, remainder] _ DivideByDigit[num, normalize, TRUE];
};


BigFromDecimalRope: PUBLIC PROC [in: Rope.ROPE] RETURNS [out: BigCARD] ~ {
excess: INT _ Rope.Length[in] MOD 4;
chunk: Rope.ROPE _ Rope.Substr[in,0,excess];
temporary: CARD _ IF excess > 0 THEN Convert.CardFromRope[chunk] ELSE 0;

in _ Rope.Substr[in, excess];
out _ BigFromSmall[ temporary ];

UNTIL Rope.Length[in] = 0 DO
chunk _ Rope.Substr[in,0,4];
in _ Rope.Substr[in,4];
out _ MultiplyByDigit[out, 10000, TRUE];
out _ BigAdd[ out, BigFromSmall[Convert.CardFromRope[chunk]]];
ENDLOOP;
};

BigFromHexRope: PUBLIC PROC [in: Rope.ROPE] RETURNS [out: BigCARD] ~ {
excess: INT _ Rope.Length[in] MOD 3;
chunk: Rope.ROPE _ Rope.Substr[in,0,excess];
temporary: CARD _ IF excess > 0 THEN Convert.CardFromRope[chunk, 16] ELSE 0;

in _ Rope.Substr[in, excess];
out _ BigFromSmall[ temporary ];

UNTIL Rope.Length[in] = 0 DO
chunk _ Rope.Substr[in,0,3];
in _ Rope.Substr[in,3];
out _ MultiplyByDigit[out, 1000h, TRUE];
out _ BigAdd[ out, BigFromSmall[Convert.CardFromRope[chunk, 16]]];
ENDLOOP;
};

BigFromRope: PUBLIC PROC [in: Rope.ROPE] RETURNS [out: BigCARD] ~ {
x, y: INT _ 0;
trailingOnes: BigCARD _ BigFromSmall[1];
j: INT _ Rope.Length[in] - 1;
i: CARDINAL;
UNTIL j < 0 OR in.Fetch[j]  ~= '\000 DO j _ j-1; x _ x + 1; ENDLOOP;
in _ Rope.Substr[in, 0, j+1];
out.size _ (j + 2)/2;
out.contents _ NEW[BigCARDRec[out.size]];
IF ~BigZero[out] THEN {
i _ j _ 0;
UNTIL j > Rope.Length[in] - 3 DO
out.contents[i] _ (in.Fetch[j] - 0C);
j _ j + 1;
out.contents[i] _ Basics.BITOR[out.contents[i], (in.Fetch[j] - 0C)*256];
j _ j + 1;
i _ i + 1;
ENDLOOP;
IF j = Rope.Length[in]-2 
THEN out.contents[out.size-1] _ (in.Fetch[j+1] - 0C)*256 + (in.Fetch[j] - 0C)
ELSE out.contents[out.size - 1] _ (in.Fetch[j] - 0C);
};
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: CARDINAL, reuse: BOOL _ FALSE] RETURNS [out: BigCARD] ~ {
multResults: CARD;
carry: CARDINAL _ 0;

IF reuse THEN out _ big ELSE out _ BigCopy[big];

SELECT small FROM
0 => {
FOR index: CARDINAL IN [0..out.size) DO
out.contents[index] _ 0;
ENDLOOP;
};
1 => {
carry _ 0;  -- for placing counting breakpoints
};
2 => {
cc: CARDINAL _ 0;
FOR index: CARDINAL IN [0..out.size) DO
c: CARDINAL _ out.contents[index];
out.contents[index] _ c*2+cc;
cc _ c / highBit;
ENDLOOP;
carry _ cc;
};
4 => {
cc: CARDINAL _ 0;
FOR index: CARDINAL IN [0..out.size) DO
c: CARDINAL _ out.contents[index];
out.contents[index] _ c*4+cc;
cc _ c / (highBit/2);
ENDLOOP;
carry _ cc;
};
ENDCASE => {
FOR index: CARDINAL IN [0..out.size) DO
multResults _ Basics.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
};
};

TwoToTheNth: PUBLIC PROC [n: CARDINAL] RETURNS [bc: BigCARD] = {
mod: CARDINAL _ n MOD bitsPerWord;
words: CARDINAL _ n / bitsPerWord + (IF mod # 0 THEN 1 ELSE 0);
bc.size _ words;
IF words # 0 THEN {
bit: CARDINAL _ 1;
bc.contents _ NEW[BigCARDRec[words]];  -- initialized to all zeros by the allocator
THROUGH [0..mod) DO bit _ bit + bit; ENDLOOP;
bc.contents[words-1] _ bit;
};
};

TimesTwoToTheNth: PUBLIC PROC [in: BigCARD, n: CARDINAL] RETURNS [out: BigCARD] = {
mod: CARDINAL _ n MOD bitsPerWord;
modNotZero: CARDINAL _ IF mod # 0 THEN 1 ELSE 0;
words: CARDINAL _ n / bitsPerWord + modNotZero;
out.size _ 0;
IF words # 0 AND in.size # 0 THEN {
digit: CARDINAL _ 1;
newWords: CARDINAL _ in.size+words;
delta: CARDINAL _ words - modNotZero;
out.contents _ NEW[BigCARDRec[newWords]];  -- initialized to all zeros by the allocator
out.size _ in.size+delta;
TRUSTED {
PrincOpsUtils.LongCopy[
from: @in.contents[0], nwords: in.size, to: @out.contents[delta] ];
};
digit _ Basics.BITSHIFT[digit, mod];
out _ MultiplyByDigit[out, digit, TRUE];
};
};

BigFactorial: PUBLIC PROC [x: CARDINAL] RETURNS [out: BigCARD] = {
bc: BigCardinals.BigCARD _ BigFromSmall[x];
exp: CARDINAL _ 0;
WHILE x > 1 DO
digit: CARDINAL _ 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];
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: CARDINAL;
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 BigZero[in] DO
[in, chunk] _ DivideByDigit[in, 10, TRUE];
IO.PutChar[ros, '0 + chunk];
ENDLOOP;
out _ IO.RopeFromROS[ros];
pos _ Rope.Length[out];
RETURN [Rope.FromProc[pos, fetch]];
};
};

BigToHexRope: 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], 16]];
ENDCASE => {
chunk: CARDINAL;
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 BigZero[in] DO
[in, chunk] _ DivideByDigit[in, 16, TRUE];
IF (chunk < 10) THEN IO.PutChar[ros, '0 + chunk] 
ELSE IO.PutChar[ros, 'a + chunk - 10];
ENDLOOP;
out _ IO.RopeFromROS[ros];
pos _ Rope.Length[out];
RETURN [Rope.FromProc[pos, fetch]];
};
};

BigToRope: PUBLIC PROC [in: BigCARD] RETURNS [out: Rope.ROPE] ~ {
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: CARDINAL _ 1B;
thisBit: CARDINAL _ Basics.BITAND[in.contents[i], mask];
WHILE mask > 0B AND thisBit = mask DO
x _ x + 1;
mask _ mask * 2;
thisBit _ Basics.BITAND[in.contents[i], mask];
ENDLOOP;
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[Basics.LowByte[nonTrailingOnes.contents[i]] + 0C]];
nonTrailingNulls _ Rope.Concat[nonTrailingNulls, Rope.FromChar[Basics.HighByte[nonTrailingOnes.contents[i]]+ 0C]];
ENDLOOP;
nonTrailingNulls _ Rope.Concat[nonTrailingNulls, Rope.FromChar[Basics.LowByte[nonTrailingOnes.contents[nonTrailingOnes.size - 1]] + 0C]];
IF nonTrailingOnes.contents[nonTrailingOnes.size - 1] > 257 THEN nonTrailingNulls _ Rope.Concat[nonTrailingNulls, Rope.FromChar[Basics.HighByte[nonTrailingOnes.contents[nonTrailingOnes.size -1]]+ 0C]];
};
out _ Rope.Concat[nonTrailingNulls, trailingNulls];
};

DivideByDigit: PROC [num: BigCARD, den: CARDINAL, reuse: BOOL _ FALSE]
RETURNS [quo: BigCARD, rem: CARDINAL _ 0] ~ {
quotient: CARDINAL;
leftShiftedRem: Basics.LongNumber _ [lc[0]];

IF reuse THEN quo _ num ELSE quo _ BigCopy[num];

SELECT den FROM
1 => {
quotient _ 0;
};
2 => {
high: CARDINAL _ 0;
FOR index: CARDINAL DECREASING IN [0..quo.size) DO
c: CARDINAL _ quo.contents[index];
rem _ Basics.BITAND[c, 1];
quo.contents[index] _ c/2 + high;
high _ rem * highBit;
ENDLOOP;
};
4 => {
high: CARDINAL _ 0;
FOR index: CARDINAL DECREASING IN [0..quo.size) DO
c: CARDINAL _ quo.contents[index];
rem _ Basics.BITAND[c, 3];
quo.contents[index] _ c/4 + high;
high _ rem * (highBit/2);
ENDLOOP;
};
ENDCASE => {
FOR index: CARDINAL DECREASING IN [0..quo.size) DO
leftShiftedRem.hi _ rem;
[quotient, rem] _ Basics.LongDivMod[quo.contents[index]+leftShiftedRem.lc, den];
quo.contents[index] _ quotient;
ENDLOOP;
};

quo _ StripLeadingZeros[quo];
};

BigFromSmall: PUBLIC PROC [in: CARDINAL] RETURNS [out: BigCARD] ~ {
IF in = 0
THEN out.size _ 0
ELSE {
out.contents _ NEW[BigCARDRec[1]];
out.size _ 1;
out.contents[0] _ in;
};
};

BigToSmall: PUBLIC PROC [in: BigCARD] RETURNS [out: CARDINAL] ~ {
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 _ [lc[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 _ [lc[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.lc];
};



BigCopy: PUBLIC PROC [in: BigCARD] RETURNS [out: BigCARD] ~ {
out.size _ in.size;
out.contents _ NEW[BigCARDRec[in.size + 2]];
IF in.size > 0 THEN TRUSTED {
PrincOpsUtils.LongCopy[from: @in.contents[0], nwords: in.size, to: @out.contents[0] ];
}; 
};


BigCompare: PUBLIC PROC [in1, in2: BigCARD] RETURNS [Basics.Comparison] ~ {

IF in1.size > in2.size THEN RETURN [greater];
IF in1.size < in2.size THEN RETURN [less];

FOR index: CARDINAL DECREASING IN [0..in1.size) DO
IF in1.contents[index] > in2.contents[index] THEN RETURN [greater];
IF in1.contents[index] < in2.contents[index] THEN RETURN [less];
ENDLOOP;

RETURN [equal];
};

BigZero: PUBLIC PROC [big: BigCARD] RETURNS [BOOL] ~ {
RETURN [big.size = 0]
};

BigOdd: PUBLIC PROC [in: BigCARD] RETURNS [BOOL] ~ {
RETURN [in.size # 0 AND in.contents[0] MOD 2 = 1];
};



BigRandom: PUBLIC PROC [length: CARDINAL _ 0, max, min: BigCARD _ [0,NIL] ]
RETURNS [out: BigCARD] ~ {
minMatters, maxMatters: BOOL _ TRUE;
lower: CARDINAL _ 0;
upper: CARDINAL _ 177777B;

IF length > 0
THEN {out.size _ length; maxMatters _ minMatters _ FALSE}
ELSE out.size _ max.size;
out.contents _ NEW[BigCARDRec[out.size]];

FOR index: CARDINAL 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[markov,lower,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] ~ {
temporary: CARDINAL;

out _ BigFromSmall[1];
IF BigZero[exponent] THEN RETURN;
IF NOT reuse THEN base _ BigCopy[base];

FOR index: CARDINAL IN [0..exponent.size) DO
temporary _ exponent.contents[index];
FOR dummy: CARDINAL IN [1..bitsPerWord] DO
IF temporary MOD 2 = 1 THEN
out _ BigDivMod[BigMultiply[out, base], modulus].rem;
temporary _ temporary / 2;
base _ BigDivMod[BigMultiply[base, base], modulus].rem;
ENDLOOP;
ENDLOOP;
};

BigGCD: PUBLIC PROC [in1, in2: BigCARD] RETURNS [out: BigCARD] ~ {
oneLarger: BOOL _ (BigCompare[in1, in2] = greater);
DO
IF BigZero[in1] OR BigZero[in2] 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] ~ {
sign: INTEGER;
sign1, sign2: INTEGER _ 1;
temp1, quo, addOrSub: BigCARD;
temp2: BigCARD _ BigFromSmall[1];

out _ modulus;

IF BigCompare[in, out] = greater THEN in _ BigDivMod[in, modulus].rem;
DO
IF BigZero[in] THEN {out_temp1; sign_sign1; EXIT};
[quo, out] _ BigDivMod[out, in];
addOrSub _ BigMultiply[temp2, quo];

SELECT TRUE FROM
sign1 # sign2 => temp1 _ BigAdd[temp1, addOrSub];
BigCompare[temp1, addOrSub] = greater => temp1 _ BigSubtract[temp1, addOrSub];
ENDCASE => {sign2 _ - sign2; temp1 _ BigSubtract[addOrSub, temp1]};

IF BigZero[out] THEN {out _ temp2; sign _ sign2; EXIT};
[quo, in] _ BigDivMod[in, out];
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: CARDINAL _ 20]
RETURNS [BOOL] ~ {
one: BigCARD _ BigFromSmall[1]; two: BigCARD _ BigFromSmall[2];
primeMinusOne: BigCARD _ BigSubtract[prime, one];
q: BigCARD _ BigCopy[primeMinusOne];
x, y: BigCARD;
j: CARDINAL;
k: CARDINAL _ 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
x _ BigRandom[max: primeMinusOne, min: two];
y _ BigExponentiate[x, q, prime];
j _ 0;

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 2^j)
IF (j = k) OR BigCompare[y, one] = equal THEN RETURN [FALSE];
ENDLOOP;
ENDLOOP;
RETURN [TRUE];
};



Test: PROC [iterations: CARDINAL _ 10] RETURNS [CARDINAL] ~ {
op1, op2, op3, res1, res2, res3, one: BigCARD;
hits, primes: CARDINAL _ 0;

one _ BigFromSmall[1];

op1 _ BigFromDecimalRope["222222222222222222"];
op2 _ BigFromDecimalRope["18446744073709551556"];
op3 _ BigFromDecimalRope["18446744073709551557"];
res1 _ BigExponentiate[op1, op2, op3];
IF NOT(BigCompare[res1, one] = equal) THEN ERROR;

FOR index: CARDINAL IN [0..iterations) DO
op1 _ BigRandom[length: 50]; op2 _ BigRandom[length: 40];

res1 _ BigSubtract[op1, op2]; res2 _ BigAdd[res1, op2];
IF NOT(BigCompare[op1, res2] = equal) THEN ERROR;

[res1, res2] _ BigDivMod[op1, op2];
res3 _ BigMultiply[op2, res1]; res3 _ BigAdd[res2, res3];
IF NOT(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 NOT(BigCompare[one, res3] = equal) THEN ERROR};
ENDLOOP;

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]
};

Time: PROC [] RETURNS [Rope.ROPE] ~ {
RETURN [ BigToDecimalRope[BigExponentiate[BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234567220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"], BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234522220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"], BigFromDecimalRope["1234567890994567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"] ] ] ]
};

markov: Random.RandomStream _ Random.Create[0, -1];

}.
ξBigCardinalsImpl.mesa
Copyright Σ 1985, 1987 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) May 12, 1987 7:18:02 pm PDT
Last Edited by: Ross May 20, 1987 10:43:45 am PDT
Most of the algorithms used below can be found in Knuth Volume 2.

Basic Arithmetic
In abuse of type, borrow is conceptually an INT in what follows.

--The LONG and LOOPHOLEs below extend the sign of borrow.
Try to guess the quotient using the leading two digits of `num' and the first digit of `den'.
Adjust the guess by using one more digit each of `num' and `den'.
Test the guess with the whole of `den'.  This is almost always correct, so  den * guess  is actually subtracted from `num'.

Borrow is treated as an INT, and sign extended with LOOPHOLEs below.
In rare cases guess is off by one, so `den' is added back to `num'.
Conversion Routines

First count (x) and strip off trailing nulls.
size of bigcardinal is ceiling((1/2)*number of non-trailing-null chars in rope).
Now append one 0 and x 1's to the end of this bigcardinal.
This procedure computes 2 ^ n.
This procedure computes in * 2 ^ n.
Returns the factorial of x.
Bring in the factors of 2.
The mask will become zero if and when we multiply 2^15 by 2.
The IF x = 0 terminates the UNTIL loop if there are no trailing ones.
Copy and Comparison
Exotic Functions
Exponentiation by repeated squaring.
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.

On each iteration through the loop there is a 1/4 chance of falsely believing that `prime' is prime when it is not prime.
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.

Test and Time
A long arithmetic progression of primes found by Paul Pritchard of Cornell
Κ˜
codeš
œ™Kš
œB™BKšœ$Οk™'Kšœ#™&Kšœœ™/K™1K˜š
	˜	Kš	œœœœœ_˜„
Kš
œ
˜
Kšœœ˜+Kšœ
œœœ˜-Kšœœ˜Kšœœ#˜/Kšœœ-œ
˜E——K˜šΟnœœ
˜Kšœœ˜8Kšœ˜Kšœ˜K˜Kšœœœ
˜Kšœœœ
œ
˜Kšœ	œœ
œœ
œ˜6Kš	œ
œœ
œ
œ˜3Kšœ
œœ
œ˜#—K˜K™AK™Kš
žœœ
œœœ
˜=K˜™K˜šžœœ
œœ˜BKš
œ˜Kšœœ	œ˜&K˜Kšœœœ˜WK˜Kš
œ˜Kšœœœ
œ
˜4K˜šœœ
œ˜+Kšœ
œœ"˜OK˜VKšœ
˜—šœœ
œ˜2Kšœ
œ˜/K˜VKšœ
˜—K˜šœœ˜KšœœΟc˜UK˜8Kš
œ˜—K˜—K˜šžœœ
œ œœœ˜`Kšœœ˜K˜Kšœœ
œ˜?Kšœœ
œ˜4K˜Kš
Ÿ@™@K˜šœœ
œ˜)Kšœ
œœ!˜NK˜1K™Kš
Ÿ9™9Kšœ	œ
œ
œœœ˜NKšœ
˜—K˜šœœ
œ˜0Kšœ
œ˜0K˜1Kšœ	œ
œ
œœœ˜NKšœ
˜—K˜K˜K˜—K˜šžœœœ˜@K˜	š	œœ

œ
œ˜1Kšœœ
œ
˜'K˜Kšœ
˜—Kšœœ
˜K˜K˜—šžœœ
œœ˜GKšœœ˜K˜Kšœœœ
œ
˜,K˜K˜Kšœœ˜)K˜šœ	œ
œ˜(šœ	œ
œ˜(Kšœ	œAœ$˜vK˜6K˜!Kšœ
˜—K˜1Kšœ
˜—K˜K˜K˜—K˜šž	œœ
œœœœ˜_Kšœœ
˜Kšœœ
˜Kšœ>œ
˜CKšœœ
˜&Kšœœ˜ K˜.K˜Kšœœ
œ˜0šœœ˜K˜>K˜Kšœ
˜Kš
œ˜—K˜Kšœ
œœ˜%šœœ
œ˜0K˜—Kš
œ'˜'šœ˜Kšœœ
˜Kšœœ
œ˜%šœ˜Kš
œ2˜2Kšœ&œ˜,K˜——šœœ˜Kšœœ˜8K˜4Kš
œ˜K˜—Kšœœ.˜CK˜%K˜%K˜K˜Kšœœ˜)K˜š	œœ

œ
œ˜5Kš
Ÿ]™]K˜8K˜HKšœœœ	œ#˜]K˜Kš
ŸA™Ašœ
œ
˜
K˜:Kšœœ
œ
˜%K˜#Kšœ^œ
œ
˜jK˜Kšœ
˜—K˜Kš
Ÿ{™{K˜šœœ
œ˜*K˜EK˜MK™Kš
ŸD™DKšœ
œ8˜IK˜?Kšœ	œ
œ
œœœ˜NKšœ
˜—K˜Kš
ŸC™Cšœ4œ˜<K˜
šœœ
œ˜*Kšœ
œE˜VK˜$K˜?Kšœ
˜—K˜Kš
œ˜—K˜Kš
œ˜Kšœ
˜—K˜K˜K˜Kš
œ˜Kšœ1œ˜7K˜—K˜K˜—™K™š
žœœ
œœœ˜JKšœœœ˜$Kšœœ˜,Kš	œœœœœ˜HK˜K˜Kš
œ ˜ K˜šœ˜K˜K˜Kšœ"œ˜(Kš
œ>˜>Kšœ
˜—K˜K˜—š
žœœ
œœœ˜FKšœœœ˜$Kšœœ˜,Kš	œœœœ!œ˜LK˜K˜Kš
œ ˜ K˜šœ˜K˜K˜Kšœ"œ˜(Kš
œB˜BKšœ
˜—K˜K˜—š
žœœ
œ
œœ˜CKšœœ˜K˜(Kšœœ˜Kšœœ
˜Kš
œ-™-Kšœœœœ
˜DK˜Kš
œP™PK˜Kšœœ˜)šœœ˜K˜
šœ˜ K˜%K˜
Kšœœ*˜HK˜
K˜
Kšœ
˜—šœ˜KšœI˜MKšœ1˜5—Kš
œ˜—Kš
œ:™:šœœ
˜K˜:K˜
Kšœ
˜—K˜CK˜>K˜—K˜šžœœ
œœ	œœœ˜lKšœ
œ
˜Kšœœ˜K˜Kšœœœ˜0K˜šœ˜˜šœœ
œ˜'K˜Kšœ
˜—K˜—˜KšœŸ#˜/Kš
œ˜—˜Kšœœ˜šœœ
œ˜'Kšœœ˜"K˜Kš
œ˜Kšœ
˜—K˜K˜—˜Kšœœ˜šœœ
œ˜'Kšœœ˜"K˜Kš
œ˜Kšœ
˜—K˜K˜—šœ˜šœœ
œ˜'K˜BK˜2K˜%Kšœ
˜—K˜——šœœ˜Kšœœ˜8K˜K˜Kš
œ˜—K˜—K˜š
žœœ
œœœ˜@Kš
œ™Kšœœœ
˜"Kš	œœœ	œœ˜?Kš
œ˜šœœ˜Kšœœ˜KšœœŸ,˜SKšœ
œœ
˜-Kš
œ˜K˜—K˜K˜—š
žœœ
œœœ˜SKš
œ#™#Kšœœœ
˜"Kš	œœœ	œœ˜0Kšœœ ˜/Kš
œ
˜
šœœ
œ˜#Kšœœ˜Kšœ
œ˜#Kšœœ˜%KšœœŸ,˜WKš
œ˜šœ˜	Kš
œ[˜[K˜—Kšœœ
˜$Kšœ"œ˜(K˜—K˜K˜—š
žœœ
œœœ˜BK™Kš
œ+˜+Kšœœ˜šœ˜Kšœœ˜Kšœœœ#œ
˜EKšœœ!œ˜8Kšœ
˜—šœ	œ˜/K™—Kšœ˜K˜K˜—šžœœ
œœœœœœ˜cšœ	˜Kšœœ˜Kšœœ(˜3šœ˜Kšœœ
˜Kšœœœ
œ˜Kšœœ˜
šœœ
œœ˜Kšœ"˜(K˜—Kšœ
œœ˜#šœ
˜Kšœ$œ˜*Kšœ˜Kšœ
˜—Kšœœ˜Kš
œ˜Kšœ˜#K˜——K˜K˜—šžœœ
œœœœœœ˜_šœ	˜Kšœœ˜Kšœœ,˜7šœ˜Kšœœ
˜Kšœœœ
œ˜Kšœœ˜
šœœ
œœ˜Kšœ"˜(K˜—Kšœ
œœ˜#šœ
˜Kšœ$œ˜*šœœ
œ˜1Kšœ
œ˜&—Kšœ
˜—Kšœœ˜Kš
œ˜Kšœ˜#K˜——K˜—K˜š
ž	œœ
œœœ˜AKšœœ˜K˜Kšœœ
˜K˜(Kšœ&œœ
˜1Kš
œ˜šœ
˜Kšœœ˜Kšœ	œ
œ˜8Kš
œ<™<šœœ˜%K˜
K˜Kšœœ˜.Kšœ
˜—Kšœœœ$™EKšœœ
œ
œ˜8Kšœ˜	—šœœ
˜K˜:K˜AK˜
Kšœ
˜—K˜NK˜]šœœ˜#šœœ!˜+K˜rK˜rKšœ
˜—Kš
œ‰
˜‰
Kšœ:œ‰
˜Ι
Kš
œ˜—K˜3K˜—K˜šž
œœœ	œœœœ	˜tKšœ
œ
˜K˜,K˜Kšœœœ˜0K˜šœ˜˜Kš
œ
˜
K˜—˜Kšœœ˜š	œœ

œ
œ˜2Kšœœ˜"Kšœ
œ˜Kš
œ!˜!Kš
œ˜Kšœ
˜—K˜—˜Kšœœ˜š	œœ

œ
œ˜2Kšœœ˜"Kšœ
œ˜Kš
œ!˜!Kš
œ˜Kšœ
˜—K˜—šœ˜š	œœ

œ
œ˜2K˜K˜PK˜Kšœ
˜—K˜——K˜Kš
œ˜K˜—K˜š
žœœ
œœœ˜Cšœ˜	Kšœ
˜šœ˜Kšœœ˜"K˜
K˜K˜——K˜—K˜š
ž
œœ
œœœ˜Ašœ	˜Kšœœ˜Kšœœ˜Kšœœ˜$—K˜K˜—š
žœœ
œœœ˜>Kš
œ!˜!šœ
œ
˜š
œ˜K˜
Kšœœ˜"Kš
œ˜Kš
œ˜K˜—š
œ˜K˜
Kšœœ˜"Kš
œ˜K˜—Kšœ˜—K˜—K˜š
ž	œœ
œœœ˜<K˜ šœ	˜Kš
œ˜Kš
œ˜Kš
œ6˜6Kšœœ˜$—Kšœ	˜K˜—K˜K˜—™K˜šžœœ
œœ˜=K˜Kšœœ˜,šœ
œ
œ˜Kš
œV˜VKš
œ˜—K˜K˜—K˜šž
œœ
œœ˜KK˜Kšœœ
œ˜-Kšœœ
œ˜*K˜š	œœ

œ
œ˜2Kšœ+œ
œ˜CKšœ+œ
œ˜@Kšœ
˜—K˜Kšœ˜K˜—K˜š
žœœ
œœœ˜6Kšœ˜K˜—K˜š
žœœ
œœœ˜4Kšœœœ˜2K˜—K˜K˜—™K˜šž	œœ
œ
œœœ˜fKšœœœ
˜$Kšœœ˜Kšœœ˜K˜šœ˜
Kšœ/œ
˜9Kšœ˜—Kšœœ˜)K˜š	œœ

œ
œ˜2Kšœœ˜/Kšœœœ˜DKš
œ;˜;Kšœœœ˜JKšœœœ
˜DKšœ
˜K˜—K˜K˜—K˜K˜šžœœ
œ+œœœ˜oKšœœ
˜K˜K˜Kšœœ
œ
˜!Kšœ
œœ˜'K˜Kš
Ÿ$™$šœœ
œ˜,K˜%šœœ
œ˜*šœœ˜K˜5—K˜K˜7Kšœ
˜—Kšœ
˜—Kš
œ˜K˜—šžœœ
œœ˜BKšœœ$˜3š
˜Kšœœœ
œ
˜*Kšœœœ˜SKšœœ˜Kšœ
˜—Kšœœœœ˜%K˜—K˜šž
œœ
œœ˜IKš
Ÿ™
™™
—™Kšœœ
˜Kšœœ˜K˜K˜!K˜K˜K˜Kšœœ!˜Fš
˜Kšœ
œœ˜2K˜ K˜#K˜šœ
œ
˜Kš
œ1˜1Kš
œN˜NKšœ<˜C—K˜Kšœœœ˜7K˜K˜#K˜šœ
œ
˜Kš
œ1˜1Kš
œN˜NKšœ<˜C—K˜Kšœ
˜—K˜Kšœœ!˜2K˜K˜—K˜šžœœ
œœœœ˜VK˜?K˜1K˜$K˜Kšœœ
˜Kšœœ˜K˜Kšœœ
œ
œ˜!K˜šœ˜Kšœœ˜"K˜
Kšœ
˜—K˜šœ˜Kš
Ÿy™yK˜,K˜!K˜K˜Kš
Ÿβ
™β
K™šœ
œ˜&š
˜Kšœ&œ
œ
˜2K˜Kšœ&Ÿ˜;Kš
œ	œœ
œ
œ˜=Kšœ
˜——Kšœ
˜—Kšœ
œ˜K˜—K˜K˜—™
K˜š
žœœœœœ˜=K˜.Kšœœ˜K˜K˜K˜K˜/K˜1K˜1K˜&Kšœ
œ œ
œ
˜1K˜šœœ
œ˜)K˜9K˜K˜7Kšœ
œ œ
œ
˜1K˜K˜#K˜9Kšœ
œ œ
œ
˜1K˜šœ+˜1K˜MKš	œ!œ
œ œ
œ˜S—Kšœ
˜K˜—Kš
œJ™JK˜'K˜'šœ˜K˜K˜Kšœ
˜—Kšœ
œ
œ
œ
˜K˜Kšœ˜
K˜K˜—šžœœœœ˜%KšœΖ˜ΝKš
œ˜—K˜K˜3—K˜Kš
œ˜—…—S u(