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
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];
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: 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: BOOLFALSE]
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];
In abuse of type, borrow is conceptually an INT in what follows.
FOR index: CARDINAL IN [0..small.size) DO
subResults ← LONG[out.contents[index]] - LONG[small.contents[index]] + borrow;
out.contents[index] ← Basics.LowHalf[subResults];
--The LONG and LOOPHOLEs below extend the sign of borrow.
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: BOOLFALSE]
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
Try to guess the quotient using the leading two digits of `num' and the first digit of `den'.
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];
Adjust the guess by using one more digit each of `num' and `den'.
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;
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: CARDINAL IN [0..den.size) DO
multResults ← Basics.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 ← LONG[num.contents[quoIndex + denIndex]] - subtract + borrow;
num.contents[quoIndex + denIndex] ← Basics.LowHalf[subResults];
borrow ← LOOPHOLE[LONG[LOOPHOLE[Basics.HighHalf[subResults], INTEGER]], 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: 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];
};
Conversion Routines
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: CARDIF 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: CARDIF 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;
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.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);
};
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: CARDINAL, reuse: BOOLFALSE] 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] = {
This procedure computes 2 ^ n.
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] = {
This procedure computes in * 2 ^ n.
mod: CARDINAL ← n MOD bitsPerWord;
modNotZero: CARDINALIF 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] = {
Returns the factorial of x.
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];
Bring in the factors of 2.
RETURN [bc];
};
BigToDecimalRope: PUBLIC PROC [in: BigCARD, reuse: BOOLFALSE]
RETURNS [out: Rope.ROPENIL] ~ {
SELECT in.size FROM
0 => RETURN ["0"];
1 => RETURN [Convert.RopeFromCard[in.contents[0]]];
ENDCASE => {
chunk: CARDINAL;
ros: STREAMIO.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: BOOLFALSE]
RETURNS [out: Rope.ROPENIL] ~ {
SELECT in.size FROM
0 => RETURN ["0"];
1 => RETURN [Convert.RopeFromCard[in.contents[0], 16]];
ENDCASE => {
chunk: CARDINAL;
ros: STREAMIO.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.ROPENIL;
i ← 0;
UNTIL i = in.size DO
mask: CARDINAL ← 1B;
thisBit: CARDINAL ← Basics.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 ← Basics.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[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: BOOLFALSE]
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];
};
Copy and Comparison
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];
};
Exotic Functions
BigRandom: PUBLIC PROC [length: CARDINAL ← 0, max, min: BigCARD ← [0,NIL] ]
RETURNS [out: BigCARD] ~ {
minMatters, maxMatters: BOOLTRUE;
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: BOOLFALSE] RETURNS [out: BigCARD] ~ {
temporary: CARDINAL;
out ← BigFromSmall[1];
IF BigZero[exponent] THEN RETURN;
IF NOT reuse THEN base ← BigCopy[base];
Exponentiation by repeated squaring.
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] ~ {
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, 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
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 2^j)
IF (j = k) OR BigCompare[y, one] = equal THEN RETURN [FALSE];
ENDLOOP;
ENDLOOP;
RETURN [TRUE];
};
Test and Time
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;
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]
};
Time: PROC [] RETURNS [Rope.ROPE] ~ {
RETURN [ BigToDecimalRope[BigExponentiate[BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234567220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"], BigFromDecimalRope["12345678901234567890123456789012345678901234567890123456789012345678901234522220123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"], BigFromDecimalRope["1234567890994567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890"] ] ] ]
};
markov: Random.RandomStream ← Random.Create[0, -1];
}.