(**************************************************************************)
(*                                                                        *)
(*  PMOS/2 software library                                               *)
(*  Copyright (C) 2023   Peter Moylan                                     *)
(*                                                                        *)
(*  This program is free software: you can redistribute it and/or modify  *)
(*  it under the terms of the GNU General Public License as published by  *)
(*  the Free Software Foundation, either version 3 of the License, or     *)
(*  (at your option) any later version.                                   *)
(*                                                                        *)
(*  This program is distributed in the hope that it will be useful,       *)
(*  but WITHOUT ANY WARRANTY; without even the implied warranty of        *)
(*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *)
(*  GNU General Public License for more details.                          *)
(*                                                                        *)
(*  You should have received a copy of the GNU General Public License     *)
(*  along with this program.  If not, see <http://www.gnu.org/licenses/>. *)
(*                                                                        *)
(*  To contact author:   http://www.pmoylan.org   peter@pmoylan.org       *)
(*                                                                        *)
(**************************************************************************)

IMPLEMENTATION MODULE Primes;

        (********************************************************)
        (*                                                      *)
        (*            Generation of prime numbers               *)
        (*                                                      *)
        (*  Programmer:         P. Moylan                       *)
        (*  Started:            15 December 2017                *)
        (*  Last edited:        18 July 2023                    *)
        (*  Status:             OK                              *)
        (*                                                      *)
        (********************************************************)


FROM SYSTEM IMPORT CARD8, CARD32;

IMPORT BigNum;

FROM BigNum IMPORT
    (* type *)  BN,
    (* proc *)  ShowBNUsage, WriteBignum, ModularPower, ModCardU;

FROM Arith64 IMPORT
    (* type *)  CARD64LE,
    (* proc *)  Mul32LE, Mul64LE, ShortMul64LE, ShortModLE;

FROM RandCard IMPORT
    (* proc *)  RandCardinal;

FROM LowLevel IMPORT
    (* proc *)  IAND;

(* FOR TESTING: *)

FROM RandCard IMPORT RandInt;

FROM STextIO IMPORT                 (* for debugging *)
    (* proc *)  WriteString, WriteLn;

FROM SIntIO IMPORT
    (* proc *)  WriteCard;          (* for debugging *)

(************************************************************************)

PROCEDURE IsPerfectSquare (A: BN): BOOLEAN;

    (* Returns TRUE iff A is the square of an integer.  *)

    (* Still needed: an elimination procedure that will let us bypass   *)
    (* these calculations for many cases.                               *)

    VAR xold, xnew, B, Q, R: BN;
        change: INTEGER;
        m, largeMod: CARDINAL;

    BEGIN
        (* Initial check.  Some numbers can be rejected based on simple *)
        (* tests.  I've taken this from                                 *)
        (*      http://mersenneforum.org/showpost.php?p=110896          *)
        (* The multiplications by obscure constants are Bloom filters,  *)
        (* which use hash functions to replace lookup tables.           *)
        (* The point of all these preliminary tests is to avoid using   *)
        (* the slow Newton's method calculation whenever possible.      *)

        (* Start with mod 128 rejection.  82% rejection rate.           *)
        (* VERY fast, can read bits directly.                           *)

        m := IAND (BigNum.LowWord(A), 127);        (* Bottom 7 bits *)

        <* COVERFLOW - *>
        IF IAND(IAND(m * 08BC40D7DH, m * 0A1E2F5D1H), 014020AH) <> 0 THEN
            RETURN FALSE;
        END (*IF*);
        <* COVERFLOW + *>

        (* The next few tests use the remainder computed below.  The    *)
        (* slow division is justified by the high rejection rate.       *)

        xold := BigNum.MakeBignum (63*25*11*17);
        xnew := BigNum.MakeBignum (19*23*31);
        B := BigNum.Prod (xold, xnew);
        BigNum.Discard (xold);
        BigNum.Discard (xnew);
        R := BigNum.ModU (A, B);
        BigNum.Discard (B);
        largeMod := BigNum.LowWord (R);
        BigNum.Discard (R);

        <* COVERFLOW - *>

        (* Residues mod 63. 75% rejection.  *)

        m := largeMod MOD 63;
        IF IAND(IAND (m * 03D491DF7H, m * 0C824A9F9H), 010F14008H) <> 0 THEN
            RETURN FALSE;
        END (*IF*);

        (* Residues mod 25. 56% rejection  *)

        m := largeMod MOD 25;
        IF IAND(IAND(m * 01929FC1BH, m * 04C9EA3B2H), 051001005H) <> 0 THEN RETURN FALSE END(*IF*);

        (* Residues mod 31. 48.4% rejection  *)
        (*  Bloom filter is a little different in this case.  *)

        m := 0D10D829AH * (largeMod MOD 31);
        IF IAND(IAND( m, m + 0672A5354H), 021025115H) <> 0 THEN RETURN FALSE END(*IF*);

        (* Residues mod 23. 47.8% rejection  *)

        m := largeMod MOD 23;
        IF IAND(IAND(m * 07BD28629H, m * 0E7180889H), 0F8300H) <> 0 THEN RETURN FALSE END(*IF*);

        (* Residues mod 19. 47.3% rejection  *)

        m := largeMod MOD 19;
        IF IAND(IAND(m * 01B8BEAD3H, m * 04D75A124H), 04280082BH) <> 0 THEN RETURN FALSE END(*IF*);

        (* Residues mod 17. 47.1% rejection  *)

        m := largeMod MOD 17;
        IF IAND(IAND(m * 06736F323H, m * 09B1D499H), 0C0000300H) <> 0 THEN RETURN FALSE END(*IF*);

        (* Residues mod 11. 45.5% rejection  *)

        m := largeMod MOD 11;
        IF IAND(IAND(m * 0ABF1A3A7H, m * 02612BF93H), 045854000H) <> 0 THEN RETURN FALSE END(*IF*);

        (* Net nonsquare rejection rate: 99.92%  *)

        <* COVERFLOW + *>

        (* End of preliminary tests. *)

        (* Use Newton's method to calculate the square root of A.       *)
        (* With initial guess x = (1+A)/2, we can easily show that:     *)
        (*  (a) if A = N^2-1 for some integer N, the x values decrease  *)
        (*      monotonically to the integer part of sqrt(A), then      *)
        (*      oscillate between two values.  As soon as we see the    *)
        (*      oscillation, we can deduce that A is not a square.      *)
        (*  (b) in all other cases, the x values decrease monotonically *)
        (*      until they remain constant at the integer part of sqrt(A). *)

        xnew := BigNum.CopyBN(A);
        BigNum.Incr (xnew);
        BigNum.Halve (xnew);

        (* The following loop terminates when x stops decreasing.   *)

        REPEAT
            xold := xnew;

            (* xnew := (xold + A/xold)/2 *)

            BigNum.Divide (A, xold, Q, R);
            BigNum.Discard (R);
            xnew := BigNum.Sum (xold, Q);
            BigNum.Discard (Q);
            BigNum.Halve (xnew);
            change := BigNum.Cmp (xnew, xold);
            BigNum.Discard (xold);
        UNTIL change >= 0;

        IF change > 0 THEN
            BigNum.Discard (xnew);
            RETURN FALSE;
        END (*IF*);

        (* Now xnew is the approximate square root of B.  We need   *)
        (* only check whether it is an exact square root.           *)

        Q := BigNum.Prod (xnew, xnew);
        BigNum.Discard (xnew);
        IF BigNum.Cmp (Q, A) = 0 THEN
            BigNum.Discard (Q);
            RETURN TRUE;
        ELSE
            BigNum.Discard (Q);
            RETURN FALSE;
        END (*IF*);

    END IsPerfectSquare;

(************************************************************************)
(*                   TESTS FOR PROBABLE PRIME NUMBERS                   *)
(************************************************************************)

PROCEDURE SimplePrimeTest (B: BN): BOOLEAN;

    (* Returns TRUE iff B does not have a small factor.  Since we are   *)
    (* only interested in large primes, we also rule out 2.             *)

    CONST maxj = 200;

    TYPE PrimeSet = ARRAY [0..maxj] OF CARDINAL;

    CONST
        divisorset = PrimeSet{
                               2,    3,    5,    7,   11,   13,   17,   19,   23,   29,
                              31,   37,   41,   43,   47,   53,   59,   61,   67,   71,
                              73,   79,   83,   89,   97,  101,  103,  107,  109,  113,
                             127,  131,  137,  139,  149,  151,  157,  163,  167,  173,
                             179,  181,  191,  193,  197,  199,  211,  223,  227,  229,
                             233,  239,  241,  251,  257,  263,  269,  271,  277,  281,
                             283,  293,  307,  311,  313,  317,  331,  337,  347,  349,
                             353,  359,  367,  373,  379,  383,  389,  397,  401,  409,
                             419,  421,  431,  433,  439,  443,  449,  457,  461,  463,
                             467,  479,  487,  491,  499,  503,  509,  521,  523,  541,
                             547,  557,  563,  569,  571,  577,  587,  593,  599,  601,
                             607,  613,  617,  619,  631,  641,  643,  647,  653,  659,
                             661,  673,  677,  683,  691,  701,  709,  719,  727,  733,
                             739,  743,  751,  757,  761,  769,  773,  787,  797,  809,
                             811,  821,  823,  827,  829,  839,  853,  857,  859,  863,
                             877,  881,  883,  887,  907,  911,  919,  929,  937,  941,
                             947,  953,  967,  971,  977,  983,  991,  997, 1009, 1013,
                            1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
                            1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
                            1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 0};

    VAR j, k: CARDINAL;

    BEGIN
        IF NOT BigNum.IsOdd(B) THEN
            RETURN FALSE;
        END (*IF*);
        j := 1;
        LOOP
            k := divisorset[j];
            IF k = 0 THEN EXIT(*LOOP*) END(*IF*);
            IF ModCardU (B, k) = 0 THEN
                RETURN FALSE;
            END (*IF*);
            INC (j);
        END (*LOOP*);
        RETURN TRUE;
    END SimplePrimeTest;

(************************************************************************)

PROCEDURE FermatPseudoprime (B: BN): BOOLEAN;

    (* Returns TRUE if there is a high probability that B is prime.     *)
    (* We use the SPRP (strong probable prime) condition, using base 2. *)

    CONST baseval = 2;

    VAR s, j, r: CARDINAL;
        one, base, pwr, Bm1, d, temp, temp2: BN;
        test: INTEGER;
        result: BOOLEAN;

    BEGIN
        one := BigNum.MakeBignum(1);
        IF BigNum.Cmp (B, one) <= 0 THEN
            BigNum.Discard (one);  RETURN FALSE;
        END(*IF*);
        pwr := BigNum.CopyBN (one);

        (* Find an s such that B < 2^s.  That will be the starting  *)
        (* point for the next step.  Also calculate 2^s.            *)

        s := BigNum.Nbits(B) + 1;
        FOR j := 1 TO s DO
            BigNum.Double (pwr);
        END (*FOR*);

        (* Find d and s for which B-1 = d*2^s.   The loop below     *)
        (* will certainly terminate, with s = 0 in the worst case.  *)

        Bm1 := BigNum.CopyBN (B);  BigNum.Decr (Bm1);
        BigNum.Divide (Bm1, pwr, d, temp);
        BigNum.Discard (temp);
        LOOP
            temp := BigNum.Prod (d, pwr);
            test := BigNum.Cmp (temp, Bm1);
            BigNum.Discard (temp);
            IF test = 0 THEN
                EXIT (*LOOP*);
            END (*IF*);
            DEC (s);  BigNum.Halve (pwr);
            BigNum.Discard (d);
            BigNum.Divide (Bm1, pwr, d, temp);
            BigNum.Discard (temp);
        END (*LOOP*);
        BigNum.Discard (pwr);

        (* Check whether base^d = 1 MOD B.  If we pass this test, we    *)
        (* don't need to do the next test.                              *)

        base := BigNum.MakeBignum (baseval);
        pwr := ModularPower (base, d, B);
        test := BigNum.Cmp (pwr, one);
        BigNum.Discard (pwr);
        IF test = 0 THEN
            BigNum.Discard (Bm1);
            BigNum.Discard (one);
            BigNum.Discard (base);
            BigNum.Discard (d);
            RETURN TRUE;
        END(*IF*);

        (* Check whether base^(d*2^r) = B-1 MOD B for some r in [0..s-1].  *)
        (* For this calculation pwr = 2^r.                              *)

        r := 0;  pwr := BigNum.CopyBN (one);
        result := FALSE;
        LOOP
            IF r = s THEN EXIT(*LOOP*) END(*IF*);
            temp2 := BigNum.Prod (d, pwr);
            temp := ModularPower (base, temp2, B);
            BigNum.Discard (temp2);
            test := BigNum.Cmp (temp, Bm1);
            BigNum.Discard (temp);
            IF test = 0 THEN
                result := TRUE;
                EXIT (*LOOP*);
            END(*IF*);
            INC (r);  BigNum.Double (pwr);
        END (*LOOP*);

        BigNum.Discard (pwr);
        BigNum.Discard (Bm1);
        BigNum.Discard (one);
        BigNum.Discard (base);
        BigNum.Discard (d);

        RETURN result;

    END FermatPseudoprime;

(************************************************************************)

PROCEDURE Jacobi (k, n: BN): INTEGER;

    (* Calculates the Jacobi symbol (k/n).  *)

    VAR j: INTEGER;
        a, b, temp, Q, R, one: BN;
        wa, wb: CARDINAL;
        beq1: BOOLEAN;

    BEGIN
        IF (BigNum.Sign(n) <= 0) OR NOT BigNum.IsOdd(n) THEN
            RETURN 0;
        END(*IF*);

        a := BigNum.CopyBN (k);
        b := BigNum.CopyBN (n);
        j := 1;
        IF BigNum.Sign(a) < 0 THEN
            BigNum.Negate (a);
            IF BigNum.LowWord(b) MOD 4 = 3 THEN j := -1 END (*IF*);
        END (*IF*);

        WHILE NOT BigNum.IsZero (a) DO
           WHILE NOT BigNum.IsOdd(a) DO

               (* Process factors of 2: (2,b) = -1 if b = 3 or 5 (mod 8) *)

               BigNum.Halve (a);
               wb := BigNum.LowWord(b) MOD 8;
               IF (wb = 3) OR (wb = 5) THEN j := -j END (*IF*);

           END (*WHILE*);

           (* Quadratic reciprocity: (a,b) = (b,a) if either a or b is  *)
           (* equal to 1 (mod 4), and(a,b) = -(b,a) if both a and b are *)
           (* equal to 3 (mod 4).  The other cases were eliminated      *)
           (* above when we handled the case of even a.                 *)

           temp := a;  a := b;  b := temp;
           wa := BigNum.LowWord (a) MOD 4;
           wb := BigNum.LowWord (b) MOD 4;
           IF (wa = 3) AND (wb = 3) THEN j := -j END (*IF*);
           BigNum.Divide (a, b, Q, R);
           BigNum.Discard (Q);
           BigNum.Discard (a);
           a := R;

        END (*WHILE*);

        one := BigNum.MakeBignum (1);
        beq1 := BigNum.Cmp (b, one) = 0;
        BigNum.Discard (one);
        BigNum.Discard (b);
        BigNum.Discard (a);
        IF beq1 THEN RETURN j
        ELSE RETURN 0
        END (*IF*);

    END Jacobi;

(************************************************************************)

(*
PROCEDURE LucasTest (A, D, P, Q: BN): BOOLEAN;

    (* The Lucas test for probable primality. Returns FALSE if A        *)
    (* cannot be a prime.                                               *)

    (* REMARK: I'm starting to suspect that D, P, and Q can be integers *)
    (* rather than Bignums, but I need to understand the test better.   *)

    BEGIN
        (* Theory: the two relevant Lucas sequences are defined by      *)
        (*                                                              *)
        (*     U[0] = 0         U[1] = 1                                *)
        (*     U[n] = P * U[n-1] - Q * U[n-2]                           *)
        (*                                                              *)
        (*     V[0] = 2         U[1] = P                                *)
        (*     V[n] = P * V[n-1] - Q * V[n-2]                           *)
        (*                                                              *)
        (* The test is, for the choices of D, P, and Q we are using,    *)
        (*      U[A+1] = 0 MOD A                                        *)
        (* There are stronger Lucas conditions that also use the V      *)
        (* values, but I'll ignore those for now.                       *)
        (*                                                              *)
        (* This looks simple, but if A is large then brute force        *)
        (* calculation of U[A+1] would be too slow.  We can speed up    *)
        (* the calculation by using the recurrence relationships        *)
        (*      U[2k] = U[k]*V[k]                                       *)
        (*      V[2k] = V[k]^2 - 2*Q^k                                  *)
        (*      U[2k+1] = (P*U[2k] + V[2k}) / 2                         *)
        (*      V[2k+1] = (D*U[2k] + P*V[2k}) / 2                       *)
        (* If either of those last two numerators is odd, it can be     *)
        (* made even by adding A, because we will carry out all of the  *)
        (* calculations modulo A.  I don't understand this last         *)
        (* comment; is it true that (y MOD A)/2 = y/2 MOD A?            *)

        (* For now, I'm returning a dummy result.  *)

        RETURN TRUE;

    END LucasTest;

(************************************************************************)

PROCEDURE LucasProbablePrime (A: BN): BOOLEAN;

    (* Another probable prime test. Returns FALSE if A cannot be a      *)
    (* prime.  A TRUE result does not guarantee that A is a prime, but  *)
    (* it gives a high probability that A is a prime.                   *)

    VAR D, P, Q: BN;  negD, result: BOOLEAN;

    BEGIN
        IF IsPerfectSquare(A) THEN
            RETURN FALSE;
        END (*IF*);

        (* Find the first D in the sequence 5, -7, 9, -11, 13, -15, ... *)
        (* for which the Jacobi symbol (D/A) is -1. Set P = 1 and       *)
        (* Q = (1 - D) / 4.                                             *)

        negD := FALSE;
        D := BigNum.MakeBignum (5);
        WHILE Jacobi (D, A) <> -1 DO
            IF negD THEN
                BigNum.Negate(D);
                BigNum.Incr(D);
                BigNum.Incr(D);
            ELSE
                BigNum.Incr(D);
                BigNum.Incr(D);
                BigNum.Negate(D);
            END (*IF*);
            negD := NOT negD;
        END (*WHILE*);

        (*
        WriteString ("Lucas: at loop exit, A = ");  BigNum.WriteBignum(A);
        WriteString (", D = ");  BigNum.WriteBignum(D);
        WriteLn;
        *)

        P := BigNum.MakeBignum(1);
        Q := BigNum.Diff (P, D);
        BigNum.Halve(Q);  BigNum.Halve(Q);

        result := LucasTest (A, D, P, Q);
        BigNum.Discard (Q);
        BigNum.Discard (P);
        BigNum.Discard (D);

        (*WriteString ("Leaving LucasTest, ");  ShowBNUsage;*)
        RETURN result;

    END LucasProbablePrime;
*)

(************************************************************************)

(*
PROCEDURE Power (x, y, m: CARDINAL): CARDINAL;

    (* Returns x^y mod m.  *)
    (* function needed only temporarily *)

    VAR result: CARDINAL;

    BEGIN
        result := 1;
        WHILE y > 0 DO
            result := result*x;
            DEC (y);
            result := result MOD m;
        END (*WHILE*);
        RETURN result;
    END Power;
*)

(************************************************************************)

PROCEDURE ModularSquare (base, modulus: BN): BN;

    (* Calculates base^2 MOD modulus.  *)

    VAR R, cbase: BN;

    BEGIN
        (* Make a local copy of the base, so as not *)
        (* to affect the caller's copy.            *)

        cbase := BigNum.CopyBN (base);

        (* First reduce cbase MOD modulus. *)

        R := BigNum.ModU (cbase, modulus);
        BigNum.Discard (cbase);

        (* Compute the square, and reduce again. *)

        cbase := BigNum.Prod (R, R);
        BigNum.Discard (R);
        R := BigNum.ModU (cbase, modulus);
        BigNum.Discard (cbase);

        RETURN R;

    END ModularSquare;

(************************************************************************)

PROCEDURE ModularPowerCCC (base, exponent, modulus: CARDINAL): CARDINAL;

    (* Calculates base^exponent MOD modulus.  We use the traditional    *)
    (* fast squaring method, but with modular reduction to stop the     *)
    (* intermediate results from getting too large.                     *)

    VAR result, cbase: CARDINAL;

    BEGIN
        (* Phase 1: stay in this phase while base is small. *)

        result := 1;
        WHILE (base < 65536) AND (exponent <> 0) DO
            IF ODD(exponent) THEN
                result := (base * result) MOD modulus;
            END (*IF*);
            exponent := exponent DIV 2;
            base := base*base;
        END (*WHILE*);

        IF exponent = 0 THEN
            RETURN result;
        END (*IF*);

        (* Phase 2: base^2 is no longer small enough to fit into a CARDINAL. *)

        cbase := base;
        WHILE exponent <> 0 DO
            IF ODD(exponent) THEN
                result := ShortModLE (Mul32LE (cbase, result), modulus);
            END (*IF*);
            exponent := exponent DIV 2;
            cbase := ShortModLE (Mul32LE (cbase, cbase), modulus);
        END (*WHILE*);
        RETURN result;

    END ModularPowerCCC;

(************************************************************************)

PROCEDURE RabinMillerLong (A: BN;  nrounds: CARDINAL): BOOLEAN;

    (* Returns TRUE if A is a probable prime.  The larger nrounds is,   *)
    (* the lower the risk of returning a false positive.  Assume A odd. *)

    VAR round, ctr2, r, N: CARDINAL;
        One, Am1, d, rnda, x, y: BN;
        endround: BOOLEAN;

    BEGIN
        One := BigNum.MakeBignum (1);

        (* Write A as 2^rúd + 1 with d odd (by factoring out powers of 2 from A - 1)*)

        Am1 := BigNum.CopyBN (A);  BigNum.Decr (Am1);
        d := BigNum.CopyBN (Am1);  r := 0;
        WHILE NOT BigNum.IsOdd(d) DO
            INC (r);  BigNum.Halve(d);
        END (*WHILE*);

        (* Each round uses a different random number rnda.  Using five  *)
        (* rounds seems to give adequately low probability of a false   *)
        (* positive.                                                    *)

        N := BigNum.Nwords(Am1);
        rnda := BigNum.MakeBignum(0);      (* to balance the Discards *)
        FOR round := 1 TO nrounds DO
            endround := FALSE;
            (* Pick a random integer rnda in the range [2, A - 2] *)
            REPEAT
                BigNum.Discard (rnda);
                rnda := BigNum.Rand (N);
            UNTIL (BigNum.Cmp(rnda, One) > 0) AND (BigNum.Cmp(rnda, Am1) < 0);
            x := BigNum.ModularPower (rnda, d, A);
            BigNum.Discard (rnda);

            (* If x = 1 or x = A-1 we have passed this round. *)

            IF (NOT BigNum.Eq (x, 1)) AND (BigNum.Cmp (x, Am1) <> 0) THEN
                ctr2 := 1;
                WHILE NOT endround AND (ctr2 < r) DO
                    y := ModularSquare (x, A);
                    BigNum.Discard (x);
                    x := y;
                    IF BigNum.Cmp (x, One) = 0 THEN
                        BigNum.Discard (d);
                        BigNum.Discard (x);
                        BigNum.Discard (Am1);
                        BigNum.Discard (One);
                        RETURN FALSE;
                    END (*IF*);
                    endround := BigNum.Cmp(x, Am1) = 0;
                    INC (ctr2);
                END (*WHILE*);
                IF NOT endround THEN
                    BigNum.Discard (d);
                    BigNum.Discard (x);
                    BigNum.Discard (Am1);
                    BigNum.Discard (One);
                    RETURN FALSE;
                END (*IF*);
            END (*IF*);
            BigNum.Discard (x);
        END (*FOR*);

        BigNum.Discard (d);
        BigNum.Discard (Am1);
        BigNum.Discard (One);
        RETURN TRUE;

    END RabinMillerLong;

(************************************************************************)

PROCEDURE RabinMillerShort (A, nrounds: CARDINAL): BOOLEAN;

    (* Returns TRUE if A is a probable prime.  The larger nrounds is,   *)
    (* the lower the risk of returning a false positive.  Assume A odd. *)

    (* This is the special case of a one-word number.  *)

    VAR round, ctr2, r, Am1, d, rnda, x: CARDINAL;
        endround: BOOLEAN;

    BEGIN
        (* Write A as 2^rúd + 1 with d odd (by factoring out powers of 2 from A - 1)*)

        Am1 := A-1;
        d := Am1;  r := 0;
        WHILE NOT ODD(d) DO
            INC (r);  d := d DIV 2;
        END (*WHILE*);

        (* Each round uses a different random number rnda.  Using five  *)
        (* rounds seems to give adequately low probability of a false   *)
        (* positive.                                                    *)

        FOR round := 1 TO nrounds DO
            endround := FALSE;
            (* Pick a random integer rnda in the range [2, A - 2] *)
            REPEAT
                rnda := RandCardinal();
            UNTIL (rnda > 1) AND (rnda < Am1);
            x := ModularPowerCCC (rnda, d, A);

            (* If x = 1 or x = A-1 we have passed this round. *)

            IF (x <> 1) AND (x <> Am1) THEN
                ctr2 := 1;
                WHILE NOT endround AND (ctr2 < r) DO
                    x := ShortModLE (Mul32LE(x,x), A);
                    IF x = 1 THEN
                        RETURN FALSE;
                    END (*IF*);
                    endround := (x = Am1);
                    INC (ctr2);
                END (*WHILE*);
                IF NOT endround THEN
                    RETURN FALSE;
                END (*IF*);
            END (*IF*);
        END (*FOR*);

        RETURN TRUE;

    END RabinMillerShort;

(************************************************************************)

PROCEDURE RabinMiller (A: BN;  nrounds: CARDINAL): BOOLEAN;

    (* Returns TRUE if A is a probable prime.  The larger nrounds is,   *)
    (* the lower the risk of returning a false positive.  Assume A odd. *)

    VAR val: CARDINAL;

    BEGIN
        IF BigNum.SingleWord (A, val) THEN
            RETURN RabinMillerShort (val, nrounds);
        ELSE
            RETURN RabinMillerLong (A, nrounds);
        END (*IF*);
    END RabinMiller;

(************************************************************************)

PROCEDURE ProbablePrime (A: BN;  fast: BOOLEAN): BOOLEAN;

    (* Returns TRUE if there is a high probability that A is prime.     *)
    (* We do a Baillie-PSW test, which is actually a combination of     *)
    (* three tests.  First, we eliminate some obvious cases where the   *)
    (* number has a small prime factor.  Next, we use the RabinMiller   *)
    (* test with five rounds.  If "fast" is true we return TRUE at      *)
    (* this point, unless of course the preceding tests have failed.    *)
    (* If "fast" is FALSE we continue with the Fermat SPRP              *)
    (* (strong probable prime) condition, using base 2.  That test      *)
    (* eliminates all but a small proportion of composite numbers.      *)
    (* Finally, if we pass those tests then we apply a Lucas probable   *)
    (* prime test.  A final FALSE result means that the number is       *)
    (* definitely composite.  A TRUE result means that there is a high  *)
    (* probability that the number is prime.                            *)

    (* REMARK: I haven't finished implementing the Lucas test, so it    *)
    (* is being skipped for now.                                        *)

    BEGIN
        IF NOT SimplePrimeTest(A) THEN
            RETURN FALSE;
        (*
        ELSIF NOT FermatPseudoprime (A) THEN    (* temporary arrangement *)
            RETURN FALSE;
        *)
        ELSIF NOT RabinMiller (A, 5) THEN
            RETURN FALSE;

        ELSIF fast THEN
            RETURN TRUE;
        ELSIF NOT FermatPseudoprime (A) THEN
            RETURN FALSE;
        ELSE
            RETURN TRUE;
            (*RETURN LucasProbablePrime(A);*)
        END (*IF*);
    END ProbablePrime;

(************************************************************************)

PROCEDURE RandomPrime (N, mask: CARDINAL;  fast: BOOLEAN): BN;

    (* Returns a random probable prime that is N words long.  This is   *)
    (* just for proof of concept.  I'll probably have to come up with   *)
    (* a better random number generator.                                *)

    (* The mask defines the number of significant bits in the top word, *)
    (* with mask = 0FFFFFFFFH if we use the whole word.                 *)

    (* If fast is FALSE we make the test tougher, but slower.           *)

    VAR p: BN;  highset: BOOLEAN;

    BEGIN
        p := BigNum.MakeBignum(0);
        REPEAT
            BigNum.Discard (p);
            p := BigNum.Rand (N);
            highset := BigNum.HighBitSet(p, mask);
            IF highset THEN
                BigNum.SetLowBit (p);       (* to ensure p is odd *)
            END (*IF*);
        UNTIL highset AND ProbablePrime (p, fast);
        RETURN p;
    END RandomPrime;

(************************************************************************)

END Primes.

