DIRECTORY
Basics USING [HighHalf, LowHalf, LongMult, LongNumber, LongDiv, LongDivMod, BITOR, BITAND, LowByte, HighByte],
BigCardinals,
Convert USING [CardFromRope, RopeFromCard],
PrincOpsUtils USING [LongCopy],
RandomCard USING [Init, Choose],
Rope USING [Concat, Length, ROPE, Substr, Fetch, FromChar];

BigCardinalsImpl: CEDAR PROGRAM
IMPORTS Basics, Convert, PrincOpsUtils, RandomCard, Rope
EXPORTS BigCardinals =
BEGIN OPEN BigCardinals;

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


BigAdd: PUBLIC PROC [in1, in2: BigCARD] RETURNS [out: BigCARD] ~ {
smaller: BigCARD;
addResults: LONG CARDINAL; 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
BEGIN
IF out.size >= out.contents.max THEN out _ BigCopy[out];  --Expand `out' if necessary
out.contents[out.size] _ carry; out.size _ out.size + 1;
END;
};

BigSubtract: PUBLIC PROC [large, small: BigCARD, reuse: BOOLEAN _ FALSE] RETURNS [out: BigCARD] ~ {
borrow, subResults: LONG CARDINAL _ 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]], LONG CARDINAL];
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]], LONG CARDINAL];
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] ~ {
product: LONG CARDINAL;  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 _ 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: BOOLEAN _ FALSE] RETURNS [quo, rem: BigCARD] ~ {
highDen, nextDen: CARDINAL;
smallNum, readyToShift: LONG CARDINAL;
carry, subtract, borrow, multResults, subResults, addResults: LONG CARDINAL;
guess, normalize, remainder: CARDINAL;
oldNumSize: CARDINAL _ num.size;
leftShiftedDigit: Basics.LongNumber _ [num[0,0]];

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

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

normalize _ 200000B/(LONG[den.contents[den.size - 1]] + 1);
num _ MultiplyByDigit[num, normalize, TRUE];
IF num.size = oldNumSize THEN
BEGIN
IF num.size >= num.contents.max THEN num _ BigCopy[num];
num.contents[num.size] _ 0; num.size _ num.size + 1;
END;

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.highbits _ num.contents[quoIndex + den.size];
smallNum _ leftShiftedDigit.lc + num.contents[quoIndex + den.size - 1]; 
guess _ IF highDen = leftShiftedDigit.highbits THEN 177777B ELSE Basics.LongDiv[smallNum, highDen];

WHILE TRUE DO
readyToShift _ smallNum - Basics.LongMult[guess, highDen];
IF readyToShift >= 200000B THEN EXIT;
leftShiftedDigit.highbits _ 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]], LONG CARDINAL];
ENDLOOP;

IF carry - borrow > num.contents[quoIndex + den.size] THEN
BEGIN
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;
END;

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: LONG CARDINAL _ 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;
};

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 
    BEGIN
        i _ 0;
        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);
    END;
    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: PROC [big: BigCARD, small: CARDINAL, reuse: BOOLEAN _ FALSE] RETURNS [out: BigCARD] ~ {
multResults: LONG CARDINAL; carry: LONG CARDINAL _ 0;

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

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 
BEGIN
IF out.size >= out.contents.max THEN out _ BigCopy[out];
out.contents[out.size] _ carry; out.size _ out.size + 1
END;
};

BigToDecimalRope: PUBLIC PROC [in: BigCARD, reuse: BOOLEAN _ FALSE] RETURNS [out: Rope.ROPE _ NIL] ~ {
chunk: CARDINAL;
IF NOT reuse THEN in _ BigCopy[in];
UNTIL BigZero[in] DO
[in, chunk] _ DivideByDigit[in, 10, TRUE];
out _ Rope.Concat[ Convert.RopeFromCard[chunk], out ];
ENDLOOP;
IF out = NIL THEN out _ "0";
};

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 
    BEGIN
        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]];
    END;
    out _ Rope.Concat[nonTrailingNulls, trailingNulls];
};

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

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

FOR index: CARDINAL DECREASING IN [0..quo.size) DO
leftShiftedRem.highbits _ 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];
};



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 [BigRelation] ~ {

IF in1.size > in2.size THEN RETURN [bigGreater];
IF in1.size < in2.size THEN RETURN [bigLess];

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

RETURN[bigEqual];
};

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

BigOdd: PUBLIC PROC [in: BigCARD] RETURNS [BOOLEAN] ~ {
IF in.size=0 THEN RETURN[FALSE] ELSE RETURN[in.contents[0] MOD 2 = 1]
};



BigRandom: PUBLIC PROC [length: CARDINAL _ 0, max, min: BigCARD _ [0,NIL] ] RETURNS [out: BigCARD] ~ {
minMatters, maxMatters: BOOLEAN _ 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] _ RandomCard.Choose[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: BOOLEAN _ 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..16] 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: BOOLEAN _ (BigCompare[in1, in2] = bigGreater);
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] = bigGreater 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];
IF sign1*sign2  = -1 THEN temp1 _ BigAdd[temp1, addOrSub]
ELSE IF BigCompare[temp1,addOrSub]=bigGreater THEN
temp1_BigSubtract[temp1,addOrSub]
ELSE {sign1 _ - sign1; temp1 _ BigSubtract[addOrSub, temp1]};

IF BigZero[out] THEN {out_temp2; sign_sign2; EXIT};
[quo, in] _ BigDivMod[in, out];
addOrSub _ BigMultiply[temp1, quo];
IF sign1*sign2  = -1 THEN temp2 _ BigAdd[temp2, addOrSub]
ELSE IF BigCompare[temp2,addOrSub]=bigGreater THEN
temp2_BigSubtract[temp2,addOrSub]
ELSE {sign2 _ - sign2; temp2 _ BigSubtract[addOrSub, temp2]};

ENDLOOP;

IF sign = -1 THEN out _ BigSubtract[modulus, out];
};


BigPrimeTest: PUBLIC PROC [prime: BigCARD, accuracy: CARDINAL _ 20] RETURNS [BOOLEAN] ~ {
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] = bigEqual THEN
DO
IF BigCompare[y, primeMinusOne] = bigEqual THEN EXIT;
j _ j+1;
y _ BigExponentiate[y, two, prime];   -- Now y = x^ (q 2^j)
IF (j = k) OR BigCompare[y, one] = bigEqual 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] = bigEqual) 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] = bigEqual) THEN ERROR;

[res1, res2] _ BigDivMod[op1, op2];
res3 _ BigMultiply[op2, res1]; res3 _ BigAdd[res2, res3];
IF NOT(BigCompare[op1, res3] = bigEqual) THEN ERROR;

IF BigCompare[BigGCD[op1, op2], one] = bigEqual THEN
{hits _ hits + 1; res1 _ BigInverse[op2, op1]; res2 _ BigMultiply[res1, op2];
res3 _ BigDivMod[res2, op1].rem; IF NOT(BigCompare[one, res3] = bigEqual) 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"] ] ] ] };

[] _ RandomCard.Init[];
END.
*BigCardinalsImpl.mesa
Dan Greene, March 6, 1984 6:29:57 pm

Last Edited by: Feigenbaum, May 31, 1984 3:38:55 pm 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.
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
�Ê2˜Jš
œ™Jšœ$Ïc
™%J™7JšÏk	œžœ@ž
œ,žœ-žœžœžœžœ˜°šœžœ
ž˜Jšžœ1˜8Jšžœ˜—Jšžœ
žœ˜J˜J™AJ™Jš	œ	žœ
žœžœžœ
˜=J˜J™J˜šÏnœžœ
žœžœ˜BJš
œ˜Jšœžœ
žœ	žœ˜/J˜Jšžœžœ˜6Jšžœ˜ J˜Jš
œ˜Jšžœžœžœ
žœ
˜4J˜šžœžœ
žœž˜+Jšœ
žœžœ"˜OJ˜VJšžœ
˜—šžœžœ
žœž˜2Jšœ
žœ˜/J˜VJšžœ
˜—J˜šžœž˜Jš
ž˜Jšžœžœ�˜UJ˜8Jšžœ
˜—J˜—J˜šŸœžœ
žœ žœžœžœ˜cJšœžœ
žœ˜&J˜Jšžœžœ
žœ˜?Jšžœžœ
žœ˜4J˜J™BJ˜šžœžœ
žœž˜)Jšœ
žœžœ!˜NJ˜1J™9Jš
œ	žœ
žœ
žœžœžœ
žœ˜WJšžœ
˜—J˜šžœžœ
žœž˜0Jšœ
žœ˜0J˜1Jš
œ	žœ
žœ
žœžœžœ
žœ˜WJšžœ
˜—J˜J˜J˜—J˜šŸœžœžœ˜@J˜	š	žœžœ
ž
œ
žœž˜1Jšžœžœ
žœ
˜'J˜Jšžœ
˜—Jšœžœ
˜J˜J˜—šŸœžœ
žœžœ˜GJšœ	žœ
žœ
žœ˜-J˜Jšžœžœžœ
žœ
˜,J˜J˜Jšœžœ˜)J˜šžœ	žœ
žœž˜(šžœ	žœ
žœž˜(˜GJšžœ$˜(—J˜XJšžœ
˜—J˜1Jšžœ
˜—J˜J˜J˜—J˜šŸ	œžœ
žœžœžœžœ˜bJšœžœ
˜Jšœžœ
žœ
˜&Jšœ>žœ
žœ
˜LJšœžœ
˜&Jšœžœ˜ J˜1J˜Jšžœžœ
žœ˜0šžœž˜Jš
ž˜J˜>J˜Jšžœ
˜Jšžœ
˜—J˜Jšžœ
žœžœ˜%šžœžœ
žœ˜0J˜—Jšœžœ"˜;Jšœ&žœ˜,šžœž˜Jš
ž˜Jšžœžœ˜8J˜4Jšžœ
˜J˜—J˜-J˜KJ˜J˜Jšœžœ˜)J˜š	žœžœ
ž
œ
žœž˜5J™_J˜>J˜HJšœžœ%žœ	žœ#˜cJ˜J™Cšžœ
žœ
ž˜
J˜:Jšžœžœ
žœ
˜%J˜)Jšžœ^žœ
žœ
˜jJ˜Jšžœ
˜—J˜J™|J˜šžœžœ
žœž˜*J˜EJ˜MJ™FJšœ
žœ8˜IJ˜?Jš
œ	žœ
žœ
žœžœžœ
žœ˜WJšžœ
˜—J˜J™Ešžœ4ž˜:Jš
ž˜J˜
šžœžœ
žœž˜*Jšœ
žœE˜VJ˜dJšžœ
˜—J˜Jšžœ
˜—J˜Jš
œ˜Jšžœ
˜—J˜J˜J˜Jš
œ˜Jšœ1žœ˜7J˜—J˜J˜J™ š
Ÿœžœ
žœžœžœ˜JJšœžœžœ˜$Jšœžœ˜,Jšœžœ
žœžœžœžœ˜QJ˜J˜Jš
œ ˜ J˜šžœž˜J˜J˜Jšœ"žœ˜(Jš
œ>˜>Jšžœ
˜—J˜J˜JšŸœžœ
žœžœ˜AJ˜Jš
œ-™-Jšœ
žœ˜J˜,Jšœžœ˜!Jšœžœ
˜Jš	œžœžœžœžœ
˜HJ˜!Jš
œP™PJ˜Jšœžœ˜-Jšœžœžœ
˜Jšœž˜	J˜J˜Jšœžœž˜(J˜1J˜Jšœ%žœ*˜TJ˜J˜Jšœžœ
˜Jšœžœ˜!JšœžœI˜UJšœžœ1˜=Jšœžœ
˜Jš
œ:™:Jšœžœžœ
˜J˜BJ˜Jšœžœ
˜J˜GJ˜BJ˜J˜—J˜šŸœžœžœ	žœžœžœ˜hJš	œ
žœ
žœ	žœ
žœ˜5J˜Jšžœžœžœ˜0J˜šžœžœ
žœž˜'J˜BJ˜XJšžœ
˜—J˜šžœžœ
˜Jš
ž˜Jšžœžœ˜8J˜7Jšžœ
˜—J˜—J˜šŸœžœ
žœžœžœžœžœžœ˜fJšœžœ
˜Jšžœ
žœžœ˜#šžœ
ž˜Jšœ$žœ˜*J˜6Jšžœ
˜—J˜J˜—J˜JšŸ	œžœžœžœ˜?Jšœ
žœ˜J˜Jšœžœ
˜J˜,J˜5Jš
œ
˜
Jšœžœ
ž˜Jšœžœ˜Jšœžœ
žœ˜@Jš
œ<™<Jšœžœžœž˜-J˜J˜Jšœžœ˜:Jšœžœ
˜Jš
œE™EJš	œžœžœ
žœ
žœ˜@Jšœžœ˜
Jšœžœžœ
˜J˜BJ˜IJ˜Jšœžœ
˜J˜RJ˜aJšœžœžœ
˜&Jšœž˜	Jšœžœžœ!ž˜3J˜~J˜~Jšœžœ
˜Jš
œ�
˜�
Jšœžœ:žœ‰
˜Ð
Jšœžœ
˜J˜7J˜J˜šŸ
œžœžœ	žœžœžœžœ	˜wJšœ
žœ
˜J˜/J˜Jšžœžœžœ˜0J˜š	žœžœ
ž
œ
žœž˜2J˜J˜PJ˜Jšžœ
˜J˜—Jš
œ˜J˜—J˜š
Ÿœžœ
žœžœžœ˜CJšžœžœžœ˜"Jšœžœ˜"J˜%J˜—J˜š
Ÿ
œžœ
žœžœžœ˜Ašžœ	ž˜Jšœžœ˜Jšœžœ˜Jšžœžœ˜$—J˜—J˜J˜J™J˜šŸœžœ
žœžœ˜=J˜Jšœžœ˜,J˜JšžœY˜`J˜J˜—J˜šŸ
œžœ
žœžœ˜EJ˜Jšžœžœ
žœ˜0Jšžœžœ
žœ˜-J˜š	žœžœ
ž
œ
žœž˜2Jšžœ+žœ
žœ˜FJšžœ+žœ
žœ˜CJšžœ
˜—J˜Jšžœ˜J˜—J˜š
Ÿœžœ
žœžœžœ˜9Jšžœ˜J˜—J˜š
Ÿœžœ
žœžœžœ˜7Jšžœžœ
žœ
žœžœ
žœžœ˜EJ˜—J˜J˜J™J˜šŸ	œžœ
žœ
žœžœžœ˜fJšœžœžœ
˜'Jšœžœ˜Jšœžœ˜J˜šžœžœ/žœ
˜GJšžœ˜—Jšœžœ˜)J˜š	žœžœ
ž
œ
žœž˜2Jšžœžœ˜/Jšžœžœžœ˜DJ˜5Jšžœžœžœ˜JJšžœžœžœ
˜DJšžœ
˜J˜—J˜J˜—J˜J˜šŸœžœ
žœ+žœžœžœ˜rJšœžœ
˜J˜J˜Jšžœžœ
žœ
˜!Jšžœ
žœžœ˜'J˜J™&šžœžœ
žœž˜,J˜%šžœžœ
žœ	ž˜!šžœžœž˜J˜5—J˜J˜7Jšžœ
˜—Jšžœ
˜—Jš
œ˜J˜—šŸœžœ
žœžœ˜BJšœžœ'˜9š
ž˜Jšžœžœžœ
žœ
˜*Jšžœžœžœ˜SJšœžœ˜Jšžœ
˜—Jšœžœžœžœ˜%J˜—J˜JšŸ
œžœ
žœžœ˜IJ™š
™Jšœžœžœ˜)J˜J˜!J˜J˜J˜Jšžœ"žœ!˜Iš
ž˜Jšžœ
žœžœ˜2J˜ J˜#Jšžœžœ ˜9šžœ
žœ'ž˜2J˜!—Jšžœ9˜=J˜Jšžœžœžœ˜3J˜J˜#Jšžœžœ ˜9šžœ
žœ'ž˜2J˜!—Jšžœ9˜=J˜Jšžœ
˜—J˜Jšžœžœ!˜2J˜J˜—J˜šŸœžœ
žœžœžœžœ˜YJ˜?J˜1J˜$J˜Jšœžœ
˜Jšœžœ˜J˜Jšžœžœ
žœ
žœ˜ J˜šžœž˜Jšœžœ˜"J˜
Jšžœ
˜—J˜šžœž˜J™{J˜,J˜!J˜J˜J™ä
J™šžœ
žœž˜)š
ž˜Jšžœ)žœ
žœ
˜5J˜Jšœ&�˜;Jš
žœ	žœžœ
žœ
žœ˜?Jšžœ
˜——Jšžœ
˜—Jšžœ
žœ˜
J˜—J˜J˜J™J˜š
Ÿœžœžœžœžœ˜=J˜.Jšœžœ˜J˜J˜J˜J˜/J˜1J˜1J˜&Jšžœ
žœ#žœ
žœ
˜4J˜šžœžœ
žœž˜)J˜9J˜J˜7Jšžœ
žœ#žœ
žœ
˜4J˜J˜#J˜9Jšžœ
žœ#žœ
žœ
˜4J˜šžœ.ž˜4J˜MJš	œ!žœ
žœ#žœ
žœ˜V—Jšžœ
˜J˜—Jš
œJ™JJ˜Ošžœž˜J˜J˜Jšžœ
˜—Jšžœ
žœ
žœ
žœ
˜J˜Jšžœ˜—J˜J˜JšŸœžœžœžœ˜%JšžœÈ˜ÎJ˜J˜Jšžœ
˜—…—BŠ[æ