(**************************************************************************)
(*                                                                        *)
(*  PMOS/2 software library                                               *)
(*  Copyright (C) 2017-2025   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 BigNum;

        (********************************************************)
        (*                                                      *)
        (*           Support for multibit integers              *)
        (*                                                      *)
        (*  Programmer:         P. Moylan                       *)
        (*  Started:            12 November 2017                *)
        (*  Last edited:        27 October 2025                 *)
        (*  Status:             OK                              *)
        (*                                                      *)
        (*  To check:                                           *)
        (*      Can we make A MOD B more efficient?             *)
        (*                                                      *)
        (********************************************************)


FROM SYSTEM IMPORT CARD8, CARD32, LOC, ADR;

IMPORT Arith64;

FROM RandCard IMPORT
    (* proc *)  Randomize, RandCardinal;

FROM Storage IMPORT
    (* proc *)  ALLOCATE, DEALLOCATE;

FROM LowLevel IMPORT
    (* proc *)  Copy, IAND, IOR, LS, RS, AddOffset, IXORB;

FROM OS2 IMPORT
    (* const*)  SEM_INDEFINITE_WAIT,
    (* type *)  HMTX,
    (* proc *)  DosCreateMutexSem, DosReleaseMutexSem;

FROM STextIO IMPORT
    (* proc *)  WriteChar, WriteString, WriteLn;

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

<* NEW TRACING - *>

TYPE
    DataSubscript = [0..MAX(CARDINAL) DIV 16];
    DataPtr = POINTER TO ARRAY DataSubscript OF CARD32;

    BNptr = POINTER TO BNrecord;
    BNrecord  = RECORD
                    negative: BOOLEAN;
                    next: BNptr;
                    allocated: CARDINAL;
                    nwords: CARDINAL;
                    val: DataPtr;
                END (*RECORD*);
    BN = BNptr;

    (* A BN is stored as an array of nwords CARD32 words, in        *)
    (* little-endian order.  The array size is allowed to expand and    *)
    (* contract as the calculations proceed.  We use a sign-magnitude   *)
    (* representation, because twos complement arithmetic is harder to  *)
    (* do when the number of bits is not fixed.  When one of the        *)
    (* operations in this module is labelled "unsigned", that means     *)
    (* that that operation either works only for nonnegative arguments, *)
    (* or that it ignores the "negative" field and works on the         *)
    (* magnitude alone.  The "next" field is for memory leak checks.    *)

    (* The "allocated" field holds the number of bytes allocated to     *)
    (* val^.  We have allocated >= bytesperword*nwords, and don't       *)
    (* require strict equality, so that sometimes there is spare        *)
    (* space after the nwords of data to allow for expansion.           *)

CONST
    Nul = CHR(0);
    bitsperword = 32;
    wordsize = SIZE(CARD32);
    bytesperword = wordsize;
    topbit = 80000000H;

<* IF TRACING THEN *>
VAR
    (* Information about memory usage.  We keep a list of all Bignums,  *)
    (* the amount of memory they are using, so we can check for leaks.  *)
    (* But only if TRACING is TRUE.                                     *)

    Usage:  RECORD
                access: HMTX;
                listhead: BNptr;
                BNcount, bytecount: CARDINAL;
            END (*RECORD*);
<* END *>

(************************************************************************)
(*                         TRACKING MEMORY USAGE                        *)
(************************************************************************)

<* IF TRACING THEN *>

PROCEDURE LockList;

    (* Get exclusive access to Usage. *)

    BEGIN
        WaitOnSemaphore (mutex, Usage.access);
    END LockList;

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

PROCEDURE UnlockList;

    (* Relinquish exclusive access to Usage. *)

    BEGIN
        DosReleaseMutexSem (Usage.access);
    END UnlockList;

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

PROCEDURE WriteCard (N: CARDINAL);

    (* Writes N in decimal to standard output. *)

    BEGIN
        IF N > 9 THEN
            WriteCard (N DIV 10);
            N := N MOD 10;
        END (*IF*);
        WriteChar (CHR(N + ORD('0')));
    END WriteCard;

<* END *>

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

PROCEDURE ShowBNUsage;

    (* Reports how much memory is being used by Bignums.  *)

    BEGIN
        <* IF TRACING THEN *>
        LockList;
        WriteString ("BN usage: ");
        WriteCard (Usage.BNcount);
        WriteString (" Bignums, ");
        WriteCard (Usage.bytecount);
        WriteString (" bytes.");
        WriteLn;
        UnlockList;
        <* END *>
    END ShowBNUsage;

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

<* IF TRACING THEN *>

PROCEDURE AddToList (A: BN);

    (* Records that A is using memory. *)

    BEGIN
        <* IF TRACING THEN *>
            LockList;
            A^.next := Usage.listhead;
            Usage.listhead := A;
            INC (Usage.BNcount);
            INC (Usage.bytecount, SIZE(BNrecord) + A^.allocated);
            UnlockList;
        <* END *>;
    END AddToList;

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

PROCEDURE RemoveFromList (A: BN);

    (* Records that A will no longer be using memory once A is discarded. *)

    VAR previous, p: BNptr;

    BEGIN
        <* IF TRACING THEN *>
            LockList;

            (* Find A in the list. *)

            previous := NIL;  p := Usage.listhead;
            WHILE (p <> NIL) AND (p <> A) DO
                previous := p;  p := p^.next;
            END (*WHILE*);
            IF p = NIL THEN
                WriteString ("BIGNUM ERROR: Cannot remove element from list.");
                WriteLn;
            ELSE
                IF previous = NIL THEN
                    Usage.listhead := p^.next;
                ELSE
                    previous^.next := p^.next;
                END (*IF*);
                DEC (Usage.BNcount);
                DEC (Usage.bytecount, SIZE(BNrecord) + A^.allocated);
            END (*IF*);
            UnlockList;
        <* END *>;
    END RemoveFromList;

<* END *>

(************************************************************************)
(*                           OUTPUT FUNCTIONS                           *)
(************************************************************************)

PROCEDURE WriteBignum (A: BN);

    (* Writes N in decimal to standard output. *)

    VAR buffer: ARRAY [0..2047] OF CHAR;
        pos: CARDINAL;

    BEGIN
        pos := 0;
        ToDecimal (A, buffer, pos);
        WriteString (buffer);
    END WriteBignum;

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

PROCEDURE WriteBignumHex (A: BN);

    (* Writes N in Hexadecimal to standard output. *)

    VAR buffer: ARRAY [0..2047] OF CHAR;

    BEGIN
        ToHex (A, buffer, 2048);
        WriteString (buffer);
    END WriteBignumHex;

(************************************************************************)
(*                   CREATING AND DESTROYING BIGNUMS                    *)
(************************************************************************)

PROCEDURE AddShortU (A: BN;  B: CARD32);  FORWARD;          (* UNSIGNED *)
PROCEDURE LSBN (VAR (*INOUT*) val: BN;  N: CARDINAL);  FORWARD;
PROCEDURE Prune (VAR (*INOUT*) A: BN);  FORWARD;

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

PROCEDURE HexToCard8 (ch: CHAR): CARD8;

    (* Binary value of one hexadecimal digit.  Returns 0 on error. *)

    BEGIN
        IF ch <= '0' THEN RETURN 0;
        ELSIF ch <= '9' THEN RETURN ORD(ch) - ORD('0');
        ELSIF ch < 'A' THEN RETURN 0;
        ELSIF ch <= 'Z' THEN RETURN ORD(ch) - ORD('A') + 10;
        ELSIF ch < 'a' THEN RETURN 0;
        ELSIF ch <= 'z' THEN RETURN ORD(ch) - ORD('a') + 10;
        ELSE RETURN 0;
        END (*IF*);
    END HexToCard8;

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

PROCEDURE HexToBN (str: ARRAY OF CHAR): BN;

    (* Converts a hexadecimal string in BigEndian format to a Bignum.   *)
    (* Space characters, if present, are ignored.                       *)

    VAR V: BN;
        pos: CARDINAL;
        ch: CHAR;

    BEGIN
        V := Zero();  pos := 0;
        WHILE pos <= HIGH(str) DO
            ch := str[pos];
            IF ch = Nul THEN
                pos := HIGH(str);
            ELSIF ch <> ' ' THEN
                LSBN (V, 4);
                AddShortU (V, HexToCard8(ch));
            END (*IF*);
            INC(pos);
        END (*WHILE*);
        RETURN V;
    END HexToBN;

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

PROCEDURE SwapIt (VAR (*INOUT*) arg: ARRAY OF LOC;  N: CARDINAL);

    (* Reverses the byte order of its N-byte argument.  Used for        *)
    (* conversion between BigEndian and LittleEndian format.            *)

    VAR j, top: CARDINAL;  temp: LOC;

    BEGIN
        top := N-1;
        FOR j := 0 TO top DIV 2 DO
            temp := arg[j];  arg[j] := arg[top-j];  arg[top-j] := temp;
        END (*FOR*);
    END SwapIt;

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

PROCEDURE BinToBN (VAR (*IN*) bindata: ARRAY OF CARD8;  pos, N: CARDINAL): BN;

    (* Converts a string of N CARD8 values, starting at bindata[pos],   *)
    (* to a Bignum.                                                     *)

    VAR V: BN;
        Np, pad0: CARDINAL;

    BEGIN
        pad0 := N MOD wordsize;
        IF pad0 <> 0 THEN
            pad0 := wordsize - pad0;
        END (*IF*);
        Np := N + pad0;

        NEW (V);
        V^.negative := FALSE;
        V^.nwords := Np DIV wordsize;
        ALLOCATE (V^.val, Np);
        V^.allocated := Np;
        V^.val^[0] := 0;        (* This is the zero padding *)
        Copy (ADR(bindata[pos]), AddOffset(V^.val, pad0), N);
        SwapIt (V^.val^, Np);
        Prune (V);
        <* IF TRACING THEN *> AddToList (V); <* END *>
        RETURN V;
    END BinToBN;

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

PROCEDURE MakeBignum (N: INTEGER): BN;

    (* Converts integer to BN. *)

    VAR result: BN;

    BEGIN
        NEW (result);
        IF N < 0 THEN
            result^.negative := TRUE;
            N := -N;
        ELSE
            result^.negative := FALSE;
        END (*IF*);
        result^.nwords := 1;
        result^.allocated := wordsize;
        ALLOCATE (result^.val, wordsize);
        result^.val^[0] := N;
        <* IF TRACING THEN *> AddToList(result); <* END *>
        RETURN result;
    END MakeBignum;

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

PROCEDURE Zero(): BN;

    (* Returns a BigNum with value zero.  *)

    VAR result: BN;

    BEGIN
        NEW (result);
        result^.negative := FALSE;
        result^.nwords := 1;
        result^.allocated := wordsize;
        ALLOCATE (result^.val, wordsize);
        result^.val^[0] := 0;
        <* IF TRACING THEN *> AddToList(result); <* END *>
        RETURN result;
    END Zero;

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

PROCEDURE LongZero (size: CARDINAL): BN;

    (* Returns a BigNum containing size words of zero.  *)

    VAR i: CARDINAL;  result: BN;

    BEGIN
        NEW (result);
        result^.negative := FALSE;
        result^.nwords := size;
        result^.allocated := size*wordsize;
        ALLOCATE (result^.val, result^.allocated);
        <* IF TRACING THEN *> AddToList(result); <* END *>
        FOR i := 0 TO result^.nwords-1 DO
            result^.val^[i] := 0;
        END (*FOR*);
        RETURN result;
    END LongZero;

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

PROCEDURE Discard (VAR (*INOUT*) num: BN);

    (* Deallocates the storage used by num. *)

    BEGIN
        IF num <> NIL THEN
            <* IF TRACING THEN *> RemoveFromList(num); <* END *>
            IF num^.allocated <> 0 THEN
                DEALLOCATE (num^.val, num^.allocated);
            END (*IF*);
            DISPOSE (num);
        END (*IF*);
    END Discard;

(************************************************************************)
(*                          SIMPLE OPERATIONS                           *)
(************************************************************************)

PROCEDURE Card (num: BN): CARDINAL;

    (* Converts BN to cardinal, saturating result if out of range. *)

    BEGIN
        IF num^.negative THEN RETURN 0;
        ELSIF num^.nwords > 1 THEN RETURN MAX(CARDINAL);
        ELSE RETURN num^.val^[0];
        END (*IF*);
    END Card;

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

PROCEDURE Int (num: BN): INTEGER;

    (* Converts BN to integer, saturating result if out of range. *)

    VAR magnitude: CARDINAL;  overflow: BOOLEAN;

    BEGIN
        overflow := num^.nwords > 1;
        magnitude := num^.val^[0];
        IF num^.negative THEN
            IF overflow THEN RETURN MIN(INTEGER);
            ELSIF magnitude > -MIN(INTEGER) THEN RETURN MIN(INTEGER);
            ELSE RETURN -VAL(INTEGER,magnitude);
            END (*IF*);
        ELSE
            IF overflow THEN RETURN MAX(INTEGER);
            ELSIF magnitude > MAX(INTEGER) THEN RETURN MAX(INTEGER);
            ELSE RETURN magnitude;
            END (*IF*);
        END (*IF*);
    END Int;

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

PROCEDURE SingleWord (A: BN;  VAR (*OUT*) val: CARDINAL): BOOLEAN;

    (* Like Card: sets val to the low-order word of A, and  *)
    (* returns TRUE iff this is the complete value.         *)

    BEGIN
        val := A^.val^[0];
        RETURN (A^.nwords = 1) AND NOT A^.negative;
    END SingleWord;

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

PROCEDURE CopyBN (A: BN): BN;

    (* Returns a copy of A. *)

    VAR tocopy: CARDINAL;
        result: BN;

    BEGIN
        NEW (result);
        IF result = NIL THEN
            WriteLn;
            WriteString ("OUT OF MEMORY");
            WriteLn;
        END (*IF*);
        result^.negative := A^.negative;
        result^.allocated := A^.allocated;
        ALLOCATE (result^.val, result^.allocated);
        tocopy := A^.nwords;
        result^.nwords := tocopy;
        tocopy := wordsize*tocopy;
        Copy (A^.val, result^.val, tocopy);
        <* IF TRACING THEN *> AddToList(result); <* END *>
        RETURN result;
    END CopyBN;

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

PROCEDURE ExpandAllocation (A: BN;  words: CARDINAL);

    (* Increases the allocated memory of A by 'words' words.  This does *)
    (* not change the number held in A, and in particular it does not   *)
    (* change A^.nwords; it simply allocates more space for             *)
    (* future expansion.                                                *)

    VAR oldsize, newsize, tocopy: CARDINAL;
        p: DataPtr;

    BEGIN
        IF words > 0 THEN
            oldsize := A^.allocated;
            newsize := oldsize + bytesperword * words;
            tocopy := bytesperword * A^.nwords;
            ALLOCATE (p, newsize);
            IF tocopy > 0 THEN
                Copy (A^.val, p, tocopy);
            END (*IF*);
            IF A^.allocated > 0 THEN
                DEALLOCATE (A^.val, A^.allocated);
            END (*IF*);
            A^.allocated := newsize;
            A^.val := p;

            <* IF TRACING THEN *>
                (* Update the usage data. *)

                LockList;
                INC(Usage.bytecount, newsize-oldsize);
                UnlockList;
            <* END *>;

        END (*IF*);
    END ExpandAllocation;

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

PROCEDURE CheckAllocation (A: BN;  extra: CARDINAL);

    (* Ensures, by expanding the allocation if necessary, that A has    *)
    (* enough space to hold 'extra' extra words.                        *)

    VAR wanted, available: CARDINAL;

    BEGIN
        available := A^.allocated DIV bytesperword;
        wanted := A^.nwords + extra;
        IF wanted > available THEN
            ExpandAllocation (A, wanted - available);
        END (*IF*);
    END CheckAllocation;

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

PROCEDURE IncreaseSize (A: BN;  newsize: CARDINAL);

    (* Increases the data size of A to newsize words.   *)
    (* Assumption: newsize > A^.nwords.                 *)

    VAR oldsize, j: CARDINAL;

    BEGIN
        oldsize := A^.nwords;
        IF newsize*wordsize > A^.allocated THEN
            ExpandAllocation (A, newsize - oldsize);
        END (*IF*);

        A^.nwords := newsize;
        FOR j := oldsize TO newsize-1 DO
            A^.val^[j] := 0;
        END (*FOR*);

    END IncreaseSize;

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

PROCEDURE Extend (VAR (*INOUT*) A: BN;  newsize: CARDINAL);

    (* Increases the data size of A to newsize words.  Does nothing   *)
    (* if A is already this big or bigger.                            *)

    BEGIN
        IF newsize > A^.nwords THEN
            IncreaseSize (A, newsize);
        END (*IF*);
    END Extend;

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

PROCEDURE Prune (VAR (*INOUT*) A: BN);

    (* Removes redundant higher-order words.  *)

    VAR newsize: CARDINAL;

    BEGIN
        newsize := A^.nwords;

        (* Note: must leave at least one word in the result.  *)

        WHILE (newsize > 1) AND (A^.val^[newsize-1] = 0) DO
            DEC (newsize);
        END (*WHILE*);
        A^.nwords := newsize;
    END Prune;

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

PROCEDURE LowWord (A: BN): CARDINAL;

    (* Returns the least significant word of the value of A. *)

    BEGIN
        RETURN A^.val^[0];
    END LowWord;

(************************************************************************)
(*                               SHIFTS                                 *)
(************************************************************************)

PROCEDURE LSW (VAR (*INOUT*) V: BN;  N: CARDINAL);

    (* Shifts V left by N CARD32 words. *)

    VAR j: CARDINAL;

    BEGIN
        IncreaseSize (V, V^.nwords+N);
        FOR j := V^.nwords - 1 TO N BY -1 DO
            V^.val^[j] := V^.val^[j-N];
        END (*FOR*);
        FOR j := 0 TO N-1 DO
            V^.val^[j] := 0;
        END (*FOR*);
    END LSW;

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

PROCEDURE LSBN (VAR (*INOUT*) V: BN;  N: CARDINAL);

    (* Shifts V left by N binary digits. *)

    CONST Allones = 0FFFFFFFFH;

    VAR maskL, maskR, j, carry, newcarry: CARD32;

    BEGIN
        CheckAllocation (V, N DIV 32 + 1);
        IF N >= 32 THEN
            LSW (V, N DIV 32);
            N := N MOD 32;
        END (*IF*);

        (* Now we may assume that N < 32.  *)

        IF N > 0 THEN
            IncreaseSize (V, V^.nwords+1);
            maskR := RS (Allones, N);
            maskL := Allones - maskR;
            carry := 0;
            FOR j := 0 TO V^.nwords -1 DO
                newcarry := RS (IAND (V^.val^[j], maskL), 32-N);
                V^.val^[j] := LS (IAND(V^.val^[j], maskR), N) + carry;
                carry := newcarry;
            END (*FOR*);
        END (*IF*);
        Prune (V);
    END LSBN;

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

PROCEDURE RSW (VAR (*INOUT*) V: BN;  N: CARDINAL);

    (* Shifts V right by N CARD32 words. *)

    VAR j: CARDINAL;

    BEGIN
        IF N > 0 THEN
            FOR j := 0 TO V^.nwords - N - 1 DO
                V^.val^[j] := V^.val^[j+N];
            END (*FOR*);
            FOR j := V^.nwords-N TO V^.nwords-1 DO
                V^.val^[j] := 0;
            END (*FOR*);
            Prune (V);
        END (*IF*);
    END RSW;

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

PROCEDURE RSBN (VAR (*INOUT*) V: BN;  N: CARDINAL);

    (* Shifts V right by N binary digits. *)

    CONST Allones = 0FFFFFFFFH;

    VAR maskL, maskR, j, carry, newcarry: CARD32;

    BEGIN
        IF N >= 32 THEN
            RSW (V, N DIV 32);
            N := N MOD 32;
        END (*IF*);

        (* Now we may assume that N < 32.  *)

        IF N > 0 THEN
            maskL := LS (Allones, N);
            maskR := Allones - maskL;
            carry := 0;
            FOR j := 0 TO V^.nwords-1 BY -1 DO
                newcarry := LS (IAND (V^.val^[j], maskR), 32-N);
                V^.val^[j] := carry + RS (V^.val^[j], N);
                carry := newcarry;
            END (*FOR*);
            Prune (V);
        END (*IF*);
    END RSBN;

(************************************************************************)
(*                              COMPARISONS                             *)
(************************************************************************)

PROCEDURE CmpU (A, B: BN): INTEGER;              (* UNSIGNED *)

    (* Unsigned comparison.  Looking at magnitudes only, returns -1 if  *)
    (* A < B, 0 if A = B, and +1 if A > B.                              *)

    VAR nA, nB, j: CARDINAL;

    BEGIN
        nA := A^.nwords;  nB := B^.nwords;
        IF nA > nB THEN
            RETURN +1;
        ELSIF nA < nB THEN
            RETURN -1;
        END (*IF*);

        (* That takes care of the obvious cases.  Now we have to deal   *)
        (* with the case where A and B have the same number of bits.    *)

        j := nA - 1;
        LOOP
            IF A^.val^[j] > B^.val^[j] THEN
                RETURN +1;
            ELSIF A^.val^[j] < B^.val^[j] THEN
                RETURN -1;
            ELSIF j = 0 THEN
                RETURN 0;
            ELSE
                DEC (j);
            END (*IF*);
        END (*LOOP*);

    END CmpU;

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

PROCEDURE Cmp (A, B: BN): INTEGER;

    (* Comparison.  Returns -1 if A < B, 0 if A = B, and +1 if A > B.  *)

    BEGIN
        IF A^.negative THEN
            IF B^.negative THEN RETURN CmpU (B, A);
            ELSE RETURN -1;
            END (*IF*);
        ELSIF B^.negative THEN RETURN +1;
        ELSE RETURN CmpU (A, B);
        END (*IF*);
    END Cmp;

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

PROCEDURE IsZero (val: BN): BOOLEAN;

    (* Returns TRUE iff val = 0 *)

    VAR j: CARDINAL;

    BEGIN
        IF (val = NIL) OR (val^.nwords = 0) THEN
            (* Should never happen. *)
            RETURN TRUE;
        END (*IF*);
        FOR j := 0 TO val^.nwords-1 DO
            IF val^.val^[j] <> 0 THEN
                RETURN FALSE;
            END (*IF*);
        END (*FOR*);
        RETURN TRUE;
    END IsZero;

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

PROCEDURE Sign (val: BN): INTEGER;

    (* Returns -1 if val < 0, 0 if val = 0, +1 if val > 0.  *)

    BEGIN
        IF IsZero(val) THEN RETURN 0
        ELSIF val^.negative THEN RETURN -1
        ELSE RETURN +1
        END (*IF*);
    END Sign;

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

PROCEDURE IsOdd (A: BN): BOOLEAN;

    (* Returns TRUE iff A is an odd number.  *)

    BEGIN
        RETURN ODD(A^.val^[0]);
    END IsOdd;

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

PROCEDURE Eq (A: BN;  N: CARD32): BOOLEAN;

    (* Returns TRUE iff A = N *)

    BEGIN
        IF A^.nwords <> 1 THEN RETURN FALSE;
        ELSE RETURN A^.val^[0] = N;
        END (*IF*);
    END Eq;

(************************************************************************)
(*                      ADDITION AND SUBTRACTION                        *)
(************************************************************************)

PROCEDURE Negate (A: BN);

    (* Changes the sign of A.  *)

    BEGIN
        A^.negative := NOT A^.negative;
    END Negate;

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

PROCEDURE SetLowBit (A: BN);

    (* Makes A odd by incrementing if necessary. *)

    BEGIN
        A^.val^[0] := IOR (A^.val^[0], 1);
    END SetLowBit;

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

PROCEDURE SetHighBit (A: BN;  mask: CARDINAL);

    (* Modifies the most significant word in A by retaining the bits    *)
    (* in mask, and then setting the bit corresponding to the           *)
    (* leftmost bit of mask.  We assume that mask defines a contiguous  *)
    (* set of bits including the least significant bit, and that        *)
    (* mask is not zero.                                                *)

    VAR pos, w: CARDINAL;

    BEGIN
        pos := A^.nwords - 1;
        w := A^.val^[pos];
        w := IAND (w, mask);
        mask := mask DIV 2 + 1;
        w := IOR (w, mask);
        A^.val^[pos] := w;
    END SetHighBit;

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

PROCEDURE HighBitSet (A: BN;  mask: CARDINAL): BOOLEAN;

    (* Returns TRUE iff the most significant bit of A is set, after     *)
    (* using mask to clear some high-order bits in the most significant *)
    (* word that should not be set.  The mask defines which bit is most *)
    (* significant.  We assume that mask defines a contiguous set of    *)
    (* bits including the least significant bit, and that mask <> 0.    *)

    VAR pos, w: CARDINAL;

    BEGIN
        pos := A^.nwords - 1;
        w := A^.val^[pos];
        w := IAND (w, mask);
        A^.val^[pos] := w;
        mask := mask DIV 2 + 1;
        RETURN IAND (w, mask) <> 0;
    END HighBitSet;

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

PROCEDURE AddU (A: BN;  B: BN);              (* UNSIGNED *)

    (* Computes  A := A + B.  This differs from Sum (below) by being    *)
    (* an in-place calculation.  The old A is overwritten.  This is an  *)
    (* unsigned calculation; the signs of A and B are ignored.          *)

    VAR nA, nB, carry, j, toadd, newval: CARDINAL;

    BEGIN
        nA := A^.nwords;  nB := B^.nwords;

        (* Make sure that A is big enough to hold the answer. *)

        IF nB > nA THEN
            nA := nB + 1;
        ELSE
            INC (nA);
        END (*IF*);
        Extend (A, nA);

        carry := 0;
        FOR j := 0 TO nA-1 DO
            IF j >= nB THEN
                toadd := 0;
            ELSE
                toadd := B^.val^[j];
            END (*IF*);

            (* Add the three components, taking care of carry. *)

            IF (toadd = MAX(CARDINAL)) AND (carry = 1) THEN
                newval := A^.val^[j];
                carry := 1;
            ELSIF A^.val^[j] > MAX(CARDINAL) - toadd - carry THEN
                newval := A^.val^[j] - (MAX(CARDINAL) - toadd - carry) - 1;
                carry := 1;
            ELSE
                newval := A^.val^[j] + toadd + carry;
                carry := 0;
            END (*IF*);
            A^.val^[j] := newval;
        END (*FOR*);

        (* We are guaranteed that the final carry is zero, because we   *)
        (* made A big enough.  Now see whether it's too big.            *)

        Prune (A);

    END AddU;

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

PROCEDURE AddShortU (A: BN;  B: CARD32);              (* UNSIGNED *)

    (* Computes  A := A + B.  This is a version of AddU (see above)     *)
    (* where the second argument is a CARD32 value.                     *)

    VAR nA, carry, j, newval: CARDINAL;

    BEGIN
        nA := A^.nwords + 1;
        Extend (A, nA);


        (* Calculate the least significant word of the result. *)

        IF A^.val^[0] > MAX(CARD32) - B THEN
            newval := A^.val^[0] - (MAX(CARD32) - B) - 1;
            carry := 1;
        ELSE
            newval := A^.val^[0] + B;
            carry := 0;
        END (*IF*);
        A^.val^[0] := newval;

        (* Propagate carry upwards. *)

        j := 1;
        WHILE carry > 0 DO

            (* Add the three components, taking care of carry. *)

            IF A^.val^[j] = MAX(CARD32) THEN
                A^.val^[j] := 0;
            ELSE
                INC (A^.val^[j]);
                carry := 0;
            END (*IF*);
            INC (j);
        END (*WHILE*);

        (* We are guaranteed that the final carry is zero, because we   *)
        (* made A big enough.  Now see whether it's too big.            *)

        Prune (A);

    END AddShortU;

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

PROCEDURE SubU (A: BN;  B: BN);              (* UNSIGNED *)

    (* Computes  A := A - B.  This differs from Diff (below) by being   *)
    (* an in-place calculation.  The old A is overwritten.  This is an  *)
    (* unsigned calculation; the signs of A and B are ignored.  Also,   *)
    (* we assume that the caller has ensured that A is greater than     *)
    (* or equal to B in magnitude.                                      *)

    VAR nA, nB, borrow, j, tosub, newval: CARDINAL;

    BEGIN
        nA := A^.nwords;  nB := B^.nwords;

        borrow := 0;
        FOR j := 0 TO nA-1 DO
            IF j >= nB THEN
                tosub := 0;
            ELSE
                tosub := B^.val^[j];
            END (*IF*);

            (* Add the three components, taking care of borrow. *)

            IF (tosub = MAX(CARDINAL)) AND (borrow = 1) THEN
                newval := A^.val^[j];
                borrow := 1;
            ELSIF A^.val^[j] < tosub + borrow THEN
                newval := A^.val^[j] + (MAX(CARDINAL) - tosub) + (1 - borrow);
                borrow := 1;
            ELSE
                newval := A^.val^[j] - tosub - borrow;
                borrow := 0;
            END (*IF*);
            A^.val^[j] := newval;
        END (*FOR*);

        Prune (A);

    END SubU;

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

PROCEDURE IncrU (A: BN);              (* UNSIGNED *)

    (* Increases the magnitude of A by 1.  The caller must guarantee    *)
    (* that the magnitude is nonzero.                                   *)

    VAR j: CARDINAL;  carry: BOOLEAN;

    BEGIN
        carry := TRUE;  j := 0;
        WHILE carry DO
            IF j = A^.nwords THEN
                IncreaseSize (A, j+1);
                A^.val^[j] := 1;
                carry := FALSE;
            ELSIF A^.val^[j] = MAX(CARD32) THEN
                A^.val^[j] := 0;
                INC (j);
            ELSE
                INC (A^.val^[j]);
                carry := FALSE;
            END (*IF*);
        END (*WHILE*);
    END IncrU;

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

PROCEDURE DecrU (A: BN);              (* UNSIGNED *)

    (* Reduces the magnitude of A by 1.  The caller must guarantee that *)
    (* the magnitude is nonzero.                                        *)

    VAR j: CARDINAL;  borrow: BOOLEAN;

    BEGIN
        borrow := TRUE;  j := 0;
        WHILE borrow DO
            IF A^.val^[j] = 0 THEN
                A^.val^[j] := MAX(CARD32);
                INC (j);
            ELSE
                DEC (A^.val^[j]);
                borrow := FALSE;
            END (*IF*);
        END (*WHILE*);

        (* The above loop will terminate on the first nonzero word.     *)
        (* It is possible that the top word was reduced to zero, so     *)
        (* we need to check for that case.                              *)

        IF (A^.nwords > 1) AND (A^.val^[A^.nwords-1] = 0) THEN
            Prune (A);
        END (*IF*);

    END DecrU;

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

PROCEDURE Incr (A: BN);

    (* Increases the value of A by 1. *)

    BEGIN
        IF A^.negative THEN
            DecrU (A);
            IF IsZero(A) THEN
                A^.negative := FALSE;
            END (*IF*);
        ELSE
            IncrU (A);
        END (*IF*);
    END Incr;

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

PROCEDURE Decr (A: BN);

    (* Reduces the value of A by 1. *)

    BEGIN
        IF IsZero(A) THEN
            A^.negative := TRUE;
        END (*IF*);
        IF A^.negative THEN
            IncrU (A);
        ELSE
            DecrU (A);
        END (*IF*);
    END Decr;

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

PROCEDURE Sum (A, B: BN): BN;

    (* Returns A+B. *)

    VAR result, other: BN;  AgtB: INTEGER;  negate: BOOLEAN;

    BEGIN
        AgtB := CmpU (A, B);
        IF AgtB >= 0 THEN
            result := CopyBN(A);
            other := B;
        ELSE
            result := CopyBN(B);
            other := A;
        END (*IF*);

        (* The following logic could possibly be simplified, but I've   *)
        (* lost the patience to do it.  It is basically just separating *)
        (* the possibilities into four cases.                           *)

        IF A^.negative = B^.negative THEN
            AddU (result, other);
        ELSE
            SubU (result, other);
        END (*IF*);

        IF A^.negative THEN
            IF B^.negative THEN
                negate := TRUE;
            ELSE
                (* A < 0, B >= 0. *)
                negate := AgtB >= 0;
            END (*IF*);
        ELSE
            IF B^.negative THEN
                (* A >= 0, B < 0. *)
                negate := AgtB < 0;
            ELSE
                negate := FALSE;
            END (*IF*);
        END (*IF*);
        result^.negative := negate;
        RETURN result;

    END Sum;

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

PROCEDURE Diff (A, B: BN): BN;

    (* Returns A-B. *)

    VAR BCopy, result: BN;

    BEGIN
        BCopy := CopyBN (B);
        BCopy^.negative := NOT B^.negative;
        result := Sum (A, BCopy);
        Discard (BCopy);
        RETURN result;
    END Diff;

(************************************************************************)
(*                              MULTIPLICATION                          *)
(************************************************************************)

PROCEDURE Prod (A,B: BN): BN;

    (* Returns A*B.  *)

    VAR result: BN;
        product: Arith64.CARD64LE;

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

    PROCEDURE Addpart (k: CARDINAL);

        (* Adds the partial product into result starting at word k. *)

        VAR carry, j: CARDINAL;
            part: ARRAY [0..1] OF CARD32;

        BEGIN

            part[0] := product.low;  part[1] := product.high;
            carry := 0;
            FOR j := 0 TO 1 DO
                IF (part[j] = MAX(CARD32)) AND (carry = 1) THEN
                    (* Let carry continue to propagate. *)
                ELSIF result^.val^[k] <= MAX(CARD32) - part[j] - carry THEN
                    INC (result^.val^[k], part[j] + carry);
                    carry := 0;
                ELSE
                    DEC (result^.val^[k], (MAX(CARD32) - part[j]) + (1 - carry));
                    carry := 1;
                END (*IF*);
                INC (k);
            END (*FOR*);

            (* We have allocated enough result space to ensure that the *)
            (* following loop won't run off the end of the array.       *)

            WHILE carry > 0 DO
                IF result^.val^[k] = MAX(CARD32) THEN
                    result^.val^[k] := 0;
                ELSE
                    INC (result^.val^[k]);
                    carry := 0;
                END (*IF*);
                INC (k);
            END (*WHILE*);

        END Addpart;

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

    VAR i, j: CARDINAL;

    BEGIN
        (* Create an initially zero result. *)

        result := LongZero (A^.nwords + B^.nwords);
        result^.negative := (A^.negative <> B^.negative);

        (* Collect and add the partial products.  *)

        FOR i := 0 TO A^.nwords-1 DO
            FOR j := 0 TO B^.nwords-1 DO

                product := Arith64.Mul32LE (A^.val^[i], B^.val^[j]);
                Addpart (i+j);

            END (*FOR*);

        END (*FOR*);

        Prune (result);
        RETURN result;

    END Prod;

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

PROCEDURE ScalarProd (A: CARDINAL;  B: BN): BN;

    (* Returns A*B.  *)

    VAR result: BN;
        product: Arith64.CARD64LE;

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

    PROCEDURE Addpart (k: CARDINAL);

        (* Adds the partial product into result starting at word k. *)

        VAR carry, j: CARDINAL;
            part: ARRAY [0..1] OF CARD32;

        BEGIN

            part[0] := product.low;  part[1] := product.high;
            carry := 0;
            FOR j := 0 TO 1 DO
                IF (part[j] = MAX(CARD32)) AND (carry = 1) THEN
                    (* Let carry continue to propagate. *)
                ELSIF result^.val^[k] <= MAX(CARD32) - part[j] - carry THEN
                    INC (result^.val^[k], part[j] + carry);
                    carry := 0;
                ELSE
                    DEC (result^.val^[k], (MAX(CARD32) - part[j]) + (1 - carry));
                    carry := 1;
                END (*IF*);
                INC (k);
            END (*FOR*);

            (* We have allocated enough result space to ensure that the *)
            (* following loop won't run off the end of the array.       *)

            WHILE carry > 0 DO
                IF result^.val^[k] = MAX(CARD32) THEN
                    result^.val^[k] := 0;
                ELSE
                    INC (result^.val^[k]);
                    carry := 0;
                END (*IF*);
                INC (k);
            END (*WHILE*);

        END Addpart;

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

    VAR j: CARDINAL;

    BEGIN
        (* Start with a zero result. *)

        result := LongZero (1 + B^.nwords);

        (* Collect and add the partial products.  *)

        FOR j := 0 TO B^.nwords-1 DO
            product := Arith64.Mul32LE (A, B^.val^[j]);
            Addpart (j);
        END (*FOR*);

        Prune (result);
        RETURN result;
    END ScalarProd;

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

PROCEDURE Double (A: BN);

    (* In-place multiplication by two.  *)

    VAR size, j, carry, nextcarry, x: CARDINAL;

    BEGIN
        size := A^.nwords;
        carry := 0;
        FOR j := 0 TO size-1 DO
            nextcarry := 0;
            x := A^.val^[j];
            IF x >= topbit THEN
                nextcarry := 1;
                DEC (x, topbit);
            END (*IF*);
            x := 2 * x + carry;
            A^.val^[j] := x;
            carry := nextcarry;
        END (*FOR*);

        (* If there was a final carry, we have to extend the size. *)

        IF carry > 0 THEN
            Extend (A, size+1);
            A^.val^[size] := 1;
        END (*IF*);

    END Double;

(************************************************************************)
(*                  MULTIPLICATION - KARATSUBA ALGORITHM                *)
(************************************************************************)

(*
PROCEDURE WriteNum (text: ARRAY OF CHAR;  A: BN);

    (* Debug output. *)

    BEGIN
        WriteString (text);
        WriteBignumHex (A);
        WriteLn;
    END WriteNum;
*)

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

PROCEDURE Split (X: BN;  m: CARDINAL;  VAR (*OUT*) X1, X0: BN);

    (* On return X0 contains the m least significant words of X, and    *)
    (* X1 contains the others, shifted down in the obvious way.         *)

    (* Remark: this will lead to X1 and X0 having unused data in their  *)
    (* val arrays.  The wasted space will be recovered when they are    *)
    (* discarded by the caller.                                         *)

    VAR j, k: CARDINAL;

    BEGIN
        X0 := CopyBN(X);  X0^.negative := FALSE;
        X1 := CopyBN(X0);
        X0^.nwords := m;

        (* Shift the top half of X1 down. *)

        k := X1^.nwords;
        IF k > m THEN
            FOR j := 0 TO k - m - 1 DO
                X1^.val^[j] := X1^.val^[m+j];
            END (*FOR*);
            DEC (X1^.nwords, m);
        ELSE
            X1^.nwords := 1;
            X1^.val^[0] := 0;
        END (*IF*);
    END Split;

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

PROCEDURE ShiftUp (m: CARDINAL;  VAR (*INOUT*) Z: BN);

    (* Inserts m zero words at the least significant end of Z, moving   *)
    (* the existing words upwards.  Note that this involves allocating  *)
    (* extra space for Z.                                               *)

    VAR p: DataPtr;  j, k: CARDINAL;

    BEGIN
        k := Z^.nwords;
        ALLOCATE (p, bytesperword*(k+m));
        FOR j := 0 TO m-1 DO
            p^[j] := 0;
        END (*FOR*);
        Copy (Z^.val, ADR(p^[m]), bytesperword*k);
        IF Z^.allocated > 0 THEN
            DEALLOCATE (Z^.val, Z^.allocated);
        END (*IF*);
        Z^.allocated := bytesperword * (k + m);
        Z^.nwords := k + m;
        Z^.val := p;
    END ShiftUp;

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

PROCEDURE KaraProdU (A,B: BN): BN;

    (* Returns A*B (unsigned) using the Karatsuba algorithm.  *)

    VAR A1, A0, B1, B0: BN;
        z0, z1, z2, z3: BN;
        temp1, temp2: BN;
        j, m, mB: CARDINAL;

    BEGIN
        m := A^.nwords;  mB := B^.nwords;
        IF mB > m THEN
            (* swap the two arguments *)
            temp1 := A;  A := B;  B := temp1;
            j := m;  m := mB;  mB := j;
        END (*IF*);
        IF (m = 0) OR (mB = 0) THEN
            RETURN Zero();
        ELSIF mB = 1 THEN
            RETURN ScalarProd (B^.val^[0], A);
        ELSE
            m := m DIV 2;
            Split (A, m, A1, A0);
            Split (B, m, B1, B0);
            z2 := KaraProdU (A1, B1);
            z0 := KaraProdU (A0, B0);

            (* Compute z3 := (A1 + A0)*(B1 + B0)  *)

            temp1 := CopyBN (A1);
            AddU (temp1, A0);
            temp2 := CopyBN (B1);
            AddU (temp2, B0);
            Discard (A0);  Discard (A1);  Discard (B0);  Discard (B1);
            z3 := KaraProdU (temp1, temp2);
            Discard (temp1);  Discard (temp2);

            (* Compute z1 := z3 - z2 - z0  *)

            z1 := z3;  SubU (z1, z2);  SubU (z1, z0);

            (*  result := ShiftUp (2*m, z2) + ShiftUp (m, z1) + z0;  *)

            ShiftUp (2*m, z2);
            ShiftUp (m, z1);
            AddU (z2, z1);  Discard (z1);
            AddU (z2, z0);  Discard (z0);

            RETURN z2;
        END (*IF*);
    END KaraProdU;

(************************************************************************)
(*                               DIVISION                               *)
(************************************************************************)

PROCEDURE Halve (A: BN);

    (* In-place division by two.  We round towards zero.  *)

    VAR size, j, x: CARDINAL;
        carry, nextcarry: BOOLEAN;

    BEGIN
        size := A^.nwords;
        carry := FALSE;
        FOR j := size-1 TO 0 BY -1 DO
            x := A^.val^[j];
            nextcarry := ODD(x);
            x := x DIV 2;
            IF carry THEN
                INC (x, topbit);
            END (*IF*);
            A^.val^[j] := x;
            carry := nextcarry;
        END (*FOR*);

        IF A^.val^[size-1] = 0 THEN
            Prune (A);
        END (*IF*);

    END Halve;

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

PROCEDURE ModU (A, B: BN): BN;          (* UNSIGNED *)

    (* Returns A MOD B.  This is an unsigned operation, where we        *)
    (* assume that A and B are nonzero and positive.                    *)

    (* The code below is a direct copy of procedure DivideU, except for *)
    (* the elimination of the quotient calculations.  See DivideU for   *)
    (* the justification of how the calculation is done.                *)

    VAR R, BK: BN;
        mask, test: CARD32;  cmp: INTEGER;
        j, shift: CARDINAL;

    BEGIN
        (* Begin by eliminating some trivial cases. *)

        cmp := CmpU (A, B);
        IF cmp = 0 THEN
            RETURN Zero();
        ELSIF cmp < 0 THEN
            RETURN CopyBN(A);
        END (*IF*);

        R := CopyBN(A);

        (* Let K be a sufficiently large number whose bits are all zero *)
        (* except for a single 1 bit, and let BK be B times K.  In this *)
        (* procedure we don't construct K explicitly, but its value is  *)
        (* implied by the variables mask and shift.  In the             *)
        (* main loop below we will keep subtracting BK from R, while    *)
        (* halving BK each time around the loop.                        *)

        j := A^.nwords;
        REPEAT
            DEC(j);
            test := A^.val^[j];
        UNTIL test <> 0;

        (* First approximation to mask.  *)

        IF test < 100H THEN
            mask := 80H;  shift := 7;
        ELSIF test < 10000H THEN
            mask := 8000H;  shift := 15;
        ELSIF test < 1000000H THEN
            mask := 800000H;  shift := 23;
        ELSE
            mask := 80000000H;  shift := 31;
        END (*IF*);

        (* Refinement of mask. *)

        WHILE IAND (mask, test) = 0 DO
            mask := mask DIV 2;  DEC(shift);
        END (*WHILE*);

        (* Recall that A^.val^[j] is the component of A on which    *)
        (* test and mask were based.                                *)

        IF A^.nwords > 1 THEN
            INC (shift, bitsperword*(A^.nwords-1));
        END (*IF*);

        BK := CopyBN (B);
        LSBN (BK, shift);
        BK^.negative := FALSE;     (* to force unsigned calculation *)

        WHILE CmpU (R, B) >= 0 DO
            IF CmpU (R, BK) >= 0 THEN
                SubU (R, BK);
            END (*IF*);
            Halve (BK);
        END (*WHILE*);

        Discard (BK);
        Prune (R);
        RETURN R;

    END ModU;

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

PROCEDURE Mod (A, B: BN): BN;

    (* Returns A MOD B, with signed numbers.  *)

    VAR R: BN;

    BEGIN
        R := ModU (A, B);
        IF IsZero(R) THEN
            R^.negative := FALSE;
        ELSE
            R^.negative := A^.negative;
        END (*IF*);
        RETURN R;
    END Mod;

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

PROCEDURE DivideU (A, B: BN;  VAR (*OUT*) Q, R: BN);      (* UNSIGNED *)

    (* Divides A by B, returns quotient and remainder.  This is the     *)
    (* unsigned case, where we assume that A and B are nonzero and      *)
    (* positive.                                                        *)

    VAR K, BK: BN;
        mask, test: CARD32;
        sign: INTEGER;
        j: CARDINAL;

    BEGIN
        (* Begin by eliminating some trivial cases. *)

        sign := CmpU (A, B);
        IF sign <= 0 THEN
            IF sign = 0 THEN
                Q := MakeBignum(1);  R := Zero();
            ELSE
                Q := Zero();  R := CopyBN(A);
            END (*IF*);
            RETURN;
        END (*IF*);

        Q := Zero();  R := CopyBN(A);

        (* Invariant: A = B*Q + R.  *)

        (* We are going to build up the quotient as the sum of numbers  *)
        (* each of which has a single 1 bit.  In the worst case (B = 1) *)
        (* the quotient can have as many bits as A does.  Let K be a    *)
        (* sufficiently large number whose bits are all zero except for *)
        (* the most significant bit.  My estimate below is a bit        *)
        (* conservative, but we should still get reasonable speed.      *)

        j := A^.nwords;
        REPEAT
            DEC(j);
            test := A^.val^[j];
        UNTIL test <> 0;

        (* First approximation to mask.  *)

        IF test < 100H THEN
            mask := 80H;
        ELSIF test < 10000H THEN
            mask := 8000H;
        ELSIF test < 1000000H THEN
            mask := 800000H;
        ELSE
            mask := 80000000H;
        END (*IF*);

        (* Refinement of mask. *)

        WHILE IAND (mask, test) = 0 DO
            mask := mask DIV 2;
        END (*WHILE*);

        NEW (K);
        K^.negative := FALSE;
        K^.nwords := A^.nwords;
        K^.allocated := K^.nwords * wordsize;
        ALLOCATE (K^.val, K^.allocated);
        IF K^.nwords > 1 THEN
            FOR j := 0 TO K^.nwords - 2 DO
                K^.val^[j] := 0;
            END (*FOR*);
        END (*IF*);
        K^.val^[K^.nwords-1] := mask;
        <* IF TRACING THEN *> AddToList(K); <* END *>

        BK := Prod (B, K);

        WHILE CmpU (R, B) >= 0 DO

            (* For any K, if we set                                 *)
            (*          Qnew = Q + K       Rnew = R - B*K           *)
            (* then the above invariant condition will still be     *)
            (* satisfied.                                           *)

            IF CmpU (R, BK) >= 0 THEN
                AddU (Q, K);  SubU (R, BK);
            END (*IF*);
            Halve (K);  Halve (BK);

        END (*WHILE*);

        Prune (Q);  Prune (R);
        Discard (K);  Discard (BK);

    END DivideU;

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

PROCEDURE Kconst (nwords, mask: CARDINAL): BN;

    (* Returns the constant K used in unsigned division.  This is for   *)
    (* testing, and this procedure will probably be scrapped later.     *)

    VAR K: BN;  j: CARDINAL;

    BEGIN
        NEW (K);
        K^.negative := FALSE;
        K^.nwords := nwords;
        ALLOCATE (K^.val, K^.nwords * wordsize);
        IF K^.nwords > 1 THEN
            FOR j := 0 TO K^.nwords - 2 DO
                K^.val^[j] := 0;
            END (*FOR*);
        END (*IF*);
        K^.val^[K^.nwords-1] := mask;
        RETURN K;
    END Kconst;

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

PROCEDURE Divide (A, B: BN;  VAR (*OUT*) Q, R: BN);

    (* Divides A by B, returns quotient and remainder.  This might      *)
    (* need more work for negative A or B.                              *)

    VAR Aneg, Bneg: BOOLEAN;

    BEGIN
        IF IsZero(B) THEN
            (* Result undefined, pick an arbitrary result. *)
            Q := Zero();  R := CopyBN(A);
        ELSIF IsZero(A) THEN
            Q := Zero();  R := Zero();
        ELSE
            Aneg := A^.negative;  Bneg := B^.negative;
            A^.negative := FALSE;  B^.negative := FALSE;
            DivideU (A, B, Q, R);
            A^.negative := Aneg;  B^.negative := Bneg;

            (* The following rules are the simplest way I know of to    *)
            (* ensure that A=QB+R in the signed case.  There is a       *)
            (* different convention that the remainder should always be *)
            (* nonnegative, but if we used that convention we would     *)
            (* have to adjust the value of Q.                           *)

            IF Aneg <> Bneg THEN
                IF NOT IsZero(Q) THEN
                    Q^.negative := TRUE;
                END (*IF*);
            END (*IF*);
            IF Aneg THEN
                IF NOT IsZero(R) THEN
                    R^.negative := TRUE;
                END (*IF*);
            END (*IF*);
        END (*IF*);
    END Divide;

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

PROCEDURE ModCardU (A: BN;  C: CARDINAL): CARDINAL;

    (* Returns the remainder of A divided by C.  This is a direct copy  *)
    (* of DivideByCardU, below, except that we don't calculate the      *)
    (* quotient.  This is an unsigned operation.                        *)

    VAR k, remainder: CARDINAL;
        current, quotient: Arith64.CARD64LE;

    BEGIN
        remainder := 0;
        FOR k := A^.nwords-1 TO 0 BY -1 DO
            current.high := remainder;
            current.low := A^.val^[k];
            Arith64.LongDiv64LE (current, C, quotient, remainder);
        END (*FOR*);
        RETURN remainder;
    END ModCardU;

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

PROCEDURE DivideByCardU (A: BN;  C: CARDINAL;  VAR (*OUT*) Q: BN;
                              VAR (*OUT*) remainder: CARDINAL);

    (* Divides A by C, giving quotient Q and a remainder.   *)
    (* This is an unsigned operation.                       *)

    (* Justification: using the fact that remainder is always less than *)
    (* C, it is easy to show that the high half of the quotient is 0.   *)

    VAR k: CARDINAL;
        current, quotient: Arith64.CARD64LE;

    BEGIN
        Q := CopyBN (A);
        remainder := 0;
        FOR k := A^.nwords-1 TO 0 BY -1 DO
            current.high := remainder;
            current.low := A^.val^[k];
            Arith64.LongDiv64LE (current, C, quotient, remainder);
            Q^.val^[k] := quotient.low;
        END (*FOR*);
        Prune (Q);
    END DivideByCardU;

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

PROCEDURE DivideBy10U (A: BN;  VAR (*OUT*) Q: BN;
                              VAR (*OUT*) remainder: CARDINAL);

    (* Divides A by 10, giving quotient Q and a remainder.  *)
    (* This is an unsigned operation.                       *)

    (* Justification: our variable "current" is a 36-bit number, with   *)
    (* 32 bits supplied by the current word of A, and the other 4 bits  *)
    (* being the remainder from the previous division by ten.  When we  *)
    (* divide this number by ten, the high half of the quotient = 0.    *)

    VAR k: CARDINAL;
        current, quotient: Arith64.CARD64LE;

    BEGIN
        Q := CopyBN (A);
        remainder := 0;
        FOR k := A^.nwords-1 TO 0 BY -1 DO
            current.high := remainder;
            current.low := A^.val^[k];
            Arith64.Div10LE (current, quotient, remainder);
            Q^.val^[k] := quotient.low;
        END (*FOR*);
        Prune (Q);
    END DivideBy10U;

(************************************************************************)
(*                           MODULAR ARITHMETIC                         *)
(************************************************************************)

PROCEDURE ModularPowerBSB (base: BN;  exponent: CARDINAL;  modulus: BN): BN;

    (* 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.                     *)

    (* This is the version where the second parameter is small. *)

    VAR result, product, R, cbase: BN;

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

        cbase := CopyBN (base);

        result := MakeBignum(1);
        IF cbase^.nwords >= modulus^.nwords THEN
            R := ModU (cbase, modulus);
            Discard (cbase);
            cbase := R;
        END (*IF*);
        WHILE exponent > 0 DO
            IF ODD(exponent) THEN
                product := Prod (result, cbase);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            exponent := exponent DIV 2;
            product := Prod (cbase, cbase);
            Discard (cbase);
            IF product^.nwords >= modulus^.nwords THEN
                cbase := ModU (product, modulus);
                Discard (product);
            ELSE
                cbase := product;
            END (*IF*);
        END (*WHILE*);
        Discard (cbase);
        RETURN result;
    END ModularPowerBSB;

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

PROCEDURE ModularPowerSSB (base, exponent: CARDINAL;  modulus: BN): BN;

    (* 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.                     *)

    (* This is the version where base and exponent are both small.      *)

    VAR result, product, cbase: BN;

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

        result := MakeBignum(1);
        WHILE (base < 65536) AND (exponent <> 0) DO
            IF ODD(exponent) THEN
                product := ScalarProd (base, result);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            exponent := exponent DIV 2;
            base := base*base;
        END (*WHILE*);

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

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

        cbase := MakeBignum (base);
        WHILE exponent <> 0 DO
            IF ODD(exponent) THEN
                product := Prod (result, cbase);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            exponent := exponent DIV 2;
            product := Prod (cbase, cbase);
            Discard (cbase);
            IF product^.nwords >= modulus^.nwords THEN
                cbase := ModU (product, modulus);
                Discard (product);
            ELSE
                cbase := product;
            END (*IF*);
        END (*WHILE*);
        Discard (cbase);
        RETURN result;

    END ModularPowerSSB;

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

PROCEDURE ModularPowerSBB (base: CARDINAL;  exponent, modulus: BN): BN;

    (* 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.                     *)

    (* This is the version where base is small.  Note, however, that    *)
    (* base increases as the calculation proceeds, so at a certain      *)
    (* point we will have to switch to using a BigNum base.             *)

    VAR result, product, cbase, cexp: BN;

    BEGIN
        IF exponent^.nwords = 1 THEN
            RETURN ModularPowerSSB (base, exponent^.val^[0], modulus);
        END(*IF*);

        cexp := CopyBN (exponent);

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

        result := MakeBignum(1);
        WHILE (base < 65536) AND NOT IsZero(cexp) DO
            IF IsOdd (cexp) THEN
                product := ScalarProd (base, result);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            Halve (cexp);
            base := base*base;
        END (*WHILE*);

        IF IsZero(cexp) THEN
            Discard (cexp);
            RETURN result;
        END (*IF*);

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

        (* Make local copies of the arguments we will change, so as not *)
        (* to affect the caller's copies of these variables.            *)

        cbase := MakeBignum (base);
        WHILE NOT IsZero(cexp) DO
            IF IsOdd (cexp) THEN
                product := Prod (result, cbase);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            Halve (cexp);
            product := Prod (cbase, cbase);
            Discard (cbase);
            IF product^.nwords >= modulus^.nwords THEN
                cbase := ModU (product, modulus);
                Discard (product);
            ELSE
                cbase := product;
            END (*IF*);
        END (*WHILE*);
        Discard (cexp);
        Discard (cbase);
        RETURN result;
    END ModularPowerSBB;

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

PROCEDURE ModularPower (base, exponent, modulus: BN): BN;

    (* 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.                     *)

    (* This is the version where all three parameters are Bignums.  *)

    VAR result, product, R, cbase, cexp: BN;

    BEGIN
        (* First check some special cases. *)

        IF Eq(modulus, 1) THEN RETURN MakeBignum(0);
        ELSIF base^.nwords = 1 THEN
            RETURN ModularPowerSBB (base^.val^[0], exponent, modulus);
        ELSIF exponent^.nwords = 1 THEN
            RETURN ModularPowerBSB (base, exponent^.val^[0], modulus);
        END(*IF*);

        (* Make local copies of the arguments we will change, so as not *)
        (* to affect the caller's copies of these variables.            *)

        cbase := CopyBN (base);
        cexp := CopyBN (exponent);

        result := MakeBignum(1);
        IF cbase^.nwords >= modulus^.nwords THEN
            R := ModU (cbase, modulus);
            Discard (cbase);
            cbase := R;
        END (*IF*);
        WHILE NOT IsZero(cexp) DO
            IF IsOdd (cexp) THEN
                product := Prod (result, cbase);
                Discard (result);
                IF product^.nwords >= modulus^.nwords THEN
                    result := ModU (product, modulus);
                    Discard (product);
                ELSE
                    result := product;
                END (*IF*);
            END (*IF*);
            Halve (cexp);
            product := Prod (cbase, cbase);
            Discard (cbase);
            IF product^.nwords >= modulus^.nwords THEN
                cbase := ModU (product, modulus);
                Discard (product);
            ELSE
                cbase := product;
            END (*IF*);
        END (*WHILE*);
        Discard (cexp);
        Discard (cbase);
        RETURN result;
    END ModularPower;

(************************************************************************)
(*                       RANDOM NUMBER GENERATOR                        *)
(************************************************************************)

PROCEDURE RC(): CARDINAL;

    (* RandCardinal() returns results in the range [1..2^31 - 2], which *)
    (* is not the range we want.  Here I'm trying to get around that    *)
    (* by combining two 16-bit halves.                                  *)

    VAR w1, w2: CARDINAL;

    BEGIN
        w1 := IAND (RandCardinal(), 0FFFFH);
        w2 := IAND (RandCardinal(), 0FFFFH);
        RETURN 65536*w1 + w2;
    END RC;

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

PROCEDURE Rand (N: CARDINAL): BN;

    (* Produces a random N-word nonnegative number.  *)

    (* The present implementation probably doesn't have the randomness  *)
    (* properties that we would prefer, but as a proof-of-concept       *)
    (* implementation it should be good enough for now.                 *)

    VAR result: BN;
        j: CARDINAL;

    BEGIN
        NEW (result);
        result^.negative := FALSE;
        result^.nwords := N;
        result^.allocated := wordsize*N;
        ALLOCATE (result^.val, result^.allocated);
        FOR j := 0 TO N-1 DO
            result^.val^[j] := RC();
        END (*FOR*);

        (* The following doubling should be obsolete now that we    *)
        (* are using RC rather than RandCardinal.                   *)

        (*result^.val^[N-1] := 2*result^.val^[N-1];*)
        <* IF TRACING THEN *> AddToList(result); <* END *>
        RETURN result;
    END Rand;

(************************************************************************)
(*                    MISCELLANEOUS OTHER FUNCTIONS                     *)
(************************************************************************)

PROCEDURE NWords (A: BN): CARDINAL;

    (* Returns the number of words used to hold the value of A. *)

    BEGIN
        RETURN A^.nwords;
    END NWords;

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

PROCEDURE WordVal (A: BN;  j: CARDINAL): CARDINAL;

    (* Returns the jth word in the internal representation of A,    *)
    (* where j=0 means the least significant word.                  *)

    BEGIN
        RETURN A^.val^[j];
    END WordVal;

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

PROCEDURE BNtoBytesU (V: BN;  VAR (*OUT*) result: ARRAY OF CARD8): CARDINAL;

    (* Converts an unsigned Bignum to an array of CARD8 values, in      *)
    (* BigEndian order, with a leading zero inserted if the first byte  *)
    (* would otherwise have its most significant bit set.  The function *)
    (* result is the number of bytes in the result.                     *)
    (* This is an unsigned operation.                                   *)

    VAR pos: CARDINAL;
        firstout: BOOLEAN;

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

    PROCEDURE Card32toBytes (A, N: CARDINAL;  suppresszero: BOOLEAN);

        (* Converts A to an N-byte value at result[pos], updates pos.   *)
        (* Special case: if suppresszero = TRUE then we suppress        *)
        (* leading zeroes, so might produce fewer than N bytes.         *)
        (* Another special case: if the first byte to be produced has   *)
        (* its high-order bit set, we insert a leading 0 byte.          *)

        VAR val: CARD8;

        BEGIN
            IF (A > 0) OR NOT suppresszero THEN
                IF N > 1 THEN
                    Card32toBytes (A DIV 256, N-1, suppresszero);
                END (*IF*);
                val := A MOD 256;
                IF firstout THEN
                    IF val >= 80H THEN
                        result[pos] := 0;  INC (pos);
                    END (*IF*);
                    firstout := FALSE;
                END (*IF*);
                result[pos] := val;  INC(pos);
            END (*IF*);
        END Card32toBytes;

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

    VAR j, NW: CARDINAL;

    BEGIN
        pos := 0;  firstout := TRUE;
        NW := V^.nwords;
        IF (NW = 1) AND (V^.val^[0] = 0) THEN
            (* Special case: V is zero. *)
            result[0] := 0;  pos := 1;
        ELSE
            Card32toBytes (V^.val^[NW-1], 4, TRUE);
            IF NW > 1 THEN
                FOR j := NW-2 TO 0 BY -1 DO
                    Card32toBytes (V^.val^[j], 4, FALSE);
                END (*FOR*);
            END (*IF*);
        END (*IF*);
        RETURN pos;
    END BNtoBytesU;

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

PROCEDURE BNtoBytes (V: BN;  VAR (*OUT*) result: ARRAY OF CARD8): CARDINAL;

    (* Converts a signed Bignum to an array of CARD8 values, in         *)
    (* BigEndian twos complement form.  An initial 0 or 0FFH is         *)
    (* inserted if needed to make the leading bit agree with the sign,  *)
    (* and of course a zero value becomes a single zero byte, but       *)
    (* otherwise there are no leading all-zero or al-one bytes.         *)
    (* The function result is the number of bytes in the result.        *)

    VAR j, N: CARDINAL;
        val: CARD8;
        carry: BOOLEAN;

    BEGIN
        N := BNtoBytesU (V, result);
        IF V^.negative THEN
            (* Replace result by its twos complement. *)

            carry := TRUE;
            FOR j := N-1 TO 0 BY -1 DO
                val := IXORB (result[j], 0FFH);
                IF carry THEN
                    IF val = 0FFH THEN
                        val := 0;           (* and leave carry set *)
                    ELSE
                        INC (val);
                        carry := FALSE;
                    END (*IF*);
                END (*IF*);
                result[j] := val;
            END (*FOR*);

        END (*IF*);
        RETURN N;
    END BNtoBytes;

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

PROCEDURE BNtoBytesNLZ (V: BN;  VAR (*OUT*) result: ARRAY OF CARD8): CARDINAL;

    (* Like BNtoBytes, except that this is an unsigned operation, and   *)
    (* there is never a leading zero except in the special case where   *)
    (* V is zero.                                                       *)

    VAR pos: CARDINAL;

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

    PROCEDURE Card32toBytes (A, N: CARDINAL;  suppresszero: BOOLEAN);

        (* Converts A to an N-byte value at result[pos], updates pos.   *)
        (* Special case: if suppresszero = TRUE then we suppress        *)
        (* leading zeroes, so might produce fewer than N bytes.         *)
        (* Another special case: if the first byte to be produced has   *)
        (* its high-order bit set, we insert a leading 0 byte.          *)

        VAR val: CARD8;

        BEGIN
            IF (A > 0) OR NOT suppresszero THEN
                IF N > 1 THEN
                    Card32toBytes (A DIV 256, N-1, suppresszero);
                END (*IF*);
                val := A MOD 256;
                result[pos] := val;  INC(pos);
            END (*IF*);
        END Card32toBytes;

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

    VAR j, NW: CARDINAL;

    BEGIN
        pos := 0;
        NW := V^.nwords;
        IF (NW = 1) AND (V^.val^[0] = 0) THEN
            (* Special case: V is zero. *)
            result[0] := 0;  pos := 1;
        ELSE
            Card32toBytes (V^.val^[NW-1], 4, TRUE);
            IF NW > 1 THEN
                FOR j := NW-2 TO 0 BY -1 DO
                    Card32toBytes (V^.val^[j], 4, FALSE);
                END (*FOR*);
            END (*IF*);
        END (*IF*);
        RETURN pos;
    END BNtoBytesNLZ;

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

PROCEDURE BNtoBytesExact (V: BN;  N: CARDINAL;
                                    VAR (*OUT*) result: ARRAY OF CARD8);

    (* Like BNtoBytes, except that leading zeroes are inserted if       *)
    (* necessary to ensure that the result length is exactly N bytes.   *)

    (* There are no guarantees if the result doesn't fit.   *)

    VAR pos: CARDINAL;

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

    PROCEDURE CardtoBytes (A, M: CARDINAL);

        (* Converts A to an M-byte value at result[pos], updates pos.   *)

        BEGIN
            IF M > 1 THEN
                CardtoBytes (A DIV 256, M-1);
            END (*IF*);
            result[pos] := A MOD 256;  INC(pos);
        END CardtoBytes;

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

    PROCEDURE CardtoBytesSLZ (A: CARDINAL);

        (* Converts A to bytes at result[pos], updates pos.  Leading    *)
        (* zeros are suppressed, but at least one byte is produced.     *)

        BEGIN
            IF A > 255 THEN
                CardtoBytesSLZ (A DIV 256);
            END (*IF*);
            result[pos] := A MOD 256;  INC(pos);
        END CardtoBytesSLZ;

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

    VAR j, k, NW, NB: CARDINAL;

    BEGIN
        pos := 0;
        NB := NbytesNLZ (V);
        NW := V^.nwords;

        (* Insert leading zeros if necessary. *)

        IF N > NB THEN
            k := N - NB;
        ELSE
            k := 0;
        END (*IF*);
        IF k > 0 THEN
            FOR j := 1 TO k DO
                result[pos] := 0;  INC(pos);
            END (*FOR*);
        END (*IF*);

        (* Now the rest of the number. *)

        CardtoBytesSLZ (V^.val^[NW - 1]);
        IF NW >= 2 THEN
            FOR j := NW-2 TO 0 BY -1 DO
                CardtoBytes (V^.val^[j], 4);
            END (*FOR*);
        END (*IF*);

    END BNtoBytesExact;

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

PROCEDURE NbytesNLZ (A: BN): CARDINAL;

    (* Returns the number of bytes needed to express A as a byte string *)
    (* with no leading zeros. (Or 1 for the special case A = 0.)        *)

    VAR N, topword, mask, test: CARDINAL;

    BEGIN
        topword := A^.val^[A^.nwords-1];
        IF topword = 0 THEN
            Prune (A);
            topword := A^.val^[A^.nwords-1];
            IF topword = 0 THEN
                RETURN 1;
            END (*IF*);
        END (*IF*);
        N := 4 * A^.nwords;
        mask := 0FF000000H;
        test := IAND(0FF000000H, topword);
        WHILE test = 0 DO
            mask := mask DIV 256;
            DEC (N);
            test := IAND(mask, topword);
        END (*WHILE*);

        RETURN N;

    END NbytesNLZ;

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

PROCEDURE Nwords (A: BN): CARDINAL;

    (* Returns the number of words in A.  *)

    BEGIN
        RETURN A^.nwords;
    END Nwords;

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

PROCEDURE Nbytes (A: BN): CARDINAL;

    (* Returns the number of bytes that BNtoBytes would return.         *)

    VAR N, msb: CARDINAL;

    BEGIN
        N := NbytesNLZ (A);

        (* Allow for an inserted zero byte if the most significant bit  *)
        (* would otherwise be 1.                                        *)

        CASE N MOD 4 OF
          | 1:  msb := 00000080H;
          | 2:  msb := 00008000H;
          | 3:  msb := 00800000H;
          | 0:  msb := 80000000H;
        END (*CASE*);
        IF IAND (A^.val^[A^.nwords-1], msb) <> 0 THEN
            INC (N);
        END (*IF*);

        RETURN N;

    END Nbytes;

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

PROCEDURE Nbits (A: BN): CARDINAL;

    (* Returns the number of significant bits in A.  Equivalently,      *)
    (* returns an N such that A >= 2^(N-1) but A < 2^N.                 *)
    (* Special case: if A is zero we return 0.                          *)

    VAR N, topword, mask: CARDINAL;

    BEGIN
        topword := A^.val^[A^.nwords-1];
        IF topword = 0 THEN
            Prune (A);
            topword := A^.val^[A^.nwords-1];
            IF topword = 0 THEN
                RETURN 0;
            END (*IF*);
        END (*IF*);
        N := bitsperword * A^.nwords;
        mask := topbit;
        WHILE mask > topword DO
            mask := mask DIV 2;
            DEC (N);
        END (*WHILE*);
        RETURN N;
    END Nbits;

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

PROCEDURE Digits (A: BN): CARDINAL;

    (* An estimate of the number of decimal digits in A.  Since we      *)
    (* don't calculate this precisely, we err if necessary on the high  *)
    (* side.  For negative A, we count the '-' as an extra "digit".     *)

    CONST log2 = 0.30103;

    VAR c: REAL;  d: INTEGER;

    BEGIN
        c := log2 * FLOAT (Nbits(A));
        d := TRUNC (c);
        IF c > FLOAT(d) THEN
            INC (d);
        END (*IF*);
        IF A^.negative THEN
            INC (d);
        END (*IF*);
        RETURN d;
    END Digits;

(************************************************************************)
(*                          STRING CONVERSIONS                          *)
(************************************************************************)

(* CURRENTLY UNUSED *)

(*
PROCEDURE ToHexN (x: CARD32;  N: CARDINAL;  dropLZ: BOOLEAN;
        VAR (*OUT*) buffer: ARRAY OF CHAR;  VAR (*INOUT*) pos: CARDINAL);
                                                            (* UNSIGNED *)

    (* Puts N-digit hexadecimal at buffer[pos], updates pos.  Drops     *)
    (* leading zeros if dropLZ is TRUE.                                 *)

    VAR firstpart: CARD32;

    BEGIN
        IF N > 1 THEN
            firstpart := x DIV 16;
            IF (firstpart > 0) OR NOT dropLZ THEN
                ToHexN (firstpart, N-1, dropLZ, buffer, pos);
            END (*IF*);
            x := x MOD 16;
        END (*IF*);
        IF x > 9 THEN
            buffer[pos] := CHR(ORD('A') - 10 + x);
        ELSE
            buffer[pos] := CHR(ORD('0') + x);
        END (*IF*);
        INC (pos);
    END ToHexN;
*)

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

PROCEDURE Hex1 (val: CARDINAL): CHAR;

    (* Converts a value in range 0..15 to hexadecimal.  *)

    BEGIN
        IF val > 9 THEN
            RETURN CHR(ORD('A') - 10 + val);
        ELSE
            RETURN CHR(ORD('0') + val);
        END (*IF*);
    END Hex1;

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

PROCEDURE ToHexNe (x: CARD32;  N: CARDINAL;  dropLZB: BOOLEAN;
        VAR (*OUT*) buffer: ARRAY OF CHAR;  VAR (*INOUT*) pos: CARDINAL);
                                                            (* UNSIGNED *)

    (* Puts N-digit hexadecimal at buffer[pos], updates pos.  Drops     *)
    (* leading zero bytes if dropLZB is TRUE, but always produces an    *)
    (* even number of characters.                                       *)

    VAR firstpart: CARD32;

    BEGIN
        IF N > 2 THEN
            firstpart := x DIV 256;
            IF (firstpart > 0) OR NOT dropLZB THEN
                ToHexNe (firstpart, N-2, dropLZB, buffer, pos);
            END (*IF*);
            x := x MOD 256;
        END (*IF*);
        buffer[pos] := Hex1(x DIV 16);  INC (pos);
        buffer[pos] := Hex1(x MOD 16);  INC (pos);
    END ToHexNe;

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

(* CURRENTLY UNUSED *)

(*
PROCEDURE ToHex8 (A: CARD32;  dropLZ: BOOLEAN;
                        VAR (*OUT*) buffer: ARRAY OF CHAR;
                        VAR (*INOUT*) pos: CARDINAL);       (* UNSIGNED *)

    (* Puts 8-digit hexadecimal at buffer[pos], updates pos.  Drops     *)
    (* leading zeros if dropLZ is TRUE.                                 *)

    BEGIN
        ToHexN (A, 8, dropLZ, buffer, pos);
    END ToHex8;
*)

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

PROCEDURE ToHex (A: BN;  VAR (*OUT*) buffer: ARRAY OF CHAR;  bufsize: CARDINAL);

    (* Puts A in buffer in hexadecimal. *)

    VAR j, pos: CARDINAL;  dropLZB: BOOLEAN;

    BEGIN
        pos := 0;
        IF A^.negative THEN
            buffer[0] := '-';
            pos := 1;
        END (*IF*);

        (* Check for buffer big enough. *)

        IF 8*A^.nwords + pos > bufsize THEN
            FOR pos := 0 TO bufsize-1 DO
                buffer[pos] := '*';
            END (*FOR*);
            RETURN;
        END (*IF*);

        dropLZB := TRUE;
        FOR j := A^.nwords-1 TO 0 BY -1 DO
            ToHexNe (A^.val^[j], 8, dropLZB, buffer, pos);
            dropLZB := FALSE;
        END (*FOR*);
        IF pos < bufsize THEN
            buffer[pos] := CHR(0);
        END (*IF*);

    END ToHex;

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

PROCEDURE FromHexW (VAR (*IN*) hexdata: ARRAY OF CHAR;
                  VAR (*INOUT*) pos: CARDINAL;  nchars: CARDINAL): CARD32;
                                                            (* UNSIGNED *)

    (* Converts a string of nchars hexadecimal characters, starting at  *)
    (* hexdata[pos], to numeric.  Updates pos.                          *)

    VAR j, result: CARDINAL;  ch: CHAR;

    BEGIN
        result := 0;
        FOR j := 0 TO nchars-1 DO
            ch := hexdata[pos];  INC (pos);
            IF ch > '9' THEN
                ch := CAP(ch);
                result := 16*result + (ORD(ch) - ORD('A') + 10);
            ELSE
                result := 16*result + ORD(ch) - ORD('0');
            END (*IF*);
        END (*FOR*);
        RETURN result;
    END FromHexW;

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

PROCEDURE FromHex (hexdata: ARRAY OF CHAR): BN;

    (* Converts hexadecimal string to BN. *)

    VAR N, nchars, j, pos: CARDINAL;
        negative: BOOLEAN;
        result: BN;

    BEGIN
        pos := 0;  negative := FALSE;
        IF hexdata[0] = '+' THEN
            pos := 1;
        ELSIF hexdata[0] = '-' THEN
            negative := TRUE;
            pos := 1;
        END (*IF*);

        (* There are 8 characters per CARD32 value. *)

        N := HIGH(hexdata) - pos;
        nchars := N MOD 8;
        IF nchars = 0 THEN
            nchars := 8;
        ELSE
            INC (N, 8-nchars);
        END (*IF*);
        N := N DIV 8;
        NEW (result);
        result^.negative := negative;
        result^.nwords := N;
        result^.allocated := N*wordsize;
        ALLOCATE (result^.val, result^.allocated);
        FOR j := N-1 TO 0 BY -1 DO
            result^.val^[j] := FromHexW (hexdata, pos, nchars);
            nchars := 8;
        END (*FOR*);

        <* IF TRACING THEN *> AddToList(result); <* END *>
        RETURN result;

    END FromHex;

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

PROCEDURE ToDecimalU (A, ten: BN;  VAR (*INOUT*) buffer: ARRAY OF CHAR;
                                            VAR (*INOUT*) pos: CARDINAL);

    (* Puts A at buffer[pos] in decimal, updates pos.  The second       *)
    (* parameter is a precomputed 10 decimal.  A must be >= 0.          *)
    (* Parameter ten is not used in this version.                       *)

    VAR digit: CARDINAL;
        Q: BN;

    BEGIN
        IF CmpU (A, ten) >= 0 THEN
            DivideBy10U (A, Q, digit);
            ToDecimalU (Q, ten, buffer, pos);
            Discard (Q);
        ELSE
            digit := A^.val^[0];
        END (*IF*);
        IF pos <= HIGH(buffer) THEN
            buffer[pos] := CHR(ORD('0') + digit);
            INC (pos);
        END (*IF*);
    END ToDecimalU;

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

PROCEDURE ToDecimal (A: BN;  VAR (*OUT*) buffer: ARRAY OF CHAR;
                                            VAR (*INOUT*) pos: CARDINAL);

    (* Puts A at buffer[pos] in decimal, updates pos. *)

    VAR ten: BN;
        Aneg: BOOLEAN;

    BEGIN
        ten := Zero();  ten^.val^[0] := 10;
        Aneg := A^.negative;
        IF Aneg THEN
            buffer[0] := '-';  pos := 1;
            A^.negative := FALSE;
        END (*IF*);
        ToDecimalU (A, ten, buffer, pos);
        IF pos <= HIGH(buffer) THEN
            buffer[pos] := CHR(0);
        END (*IF*);
        A^.negative := Aneg;
        Discard (ten);
    END ToDecimal;

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

PROCEDURE FromDecimal (decdata: ARRAY OF CHAR): BN;

    (* Converts decimal string to BN. *)

    VAR ten, A, result: BN;
        pos: CARDINAL;   negate: BOOLEAN;

    BEGIN
        ten := Zero();  ten^.val^[0] := 10;
        result := Zero();
        negate := FALSE;
        pos := 0;
        IF decdata[0] = '+' THEN
            pos := 1;
        ELSIF decdata[0] = '-' THEN
            negate := TRUE;
            pos := 1;
        END (*IF*);
        WHILE (pos <= HIGH(decdata)) AND (decdata[pos] <> CHR(0)) DO
            A := Prod (result, ten);
            AddShortU (A, ORD(decdata[pos]) - ORD('0'));
            INC (pos);
            Discard (result);
            result := A;
        END (*WHILE*);
        Discard (ten);
        IF negate THEN
            result^.negative := TRUE;
        END (*IF*);
        RETURN result;
    END FromDecimal;

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

BEGIN
    <* IF TRACING THEN *>
    WITH Usage DO
        DosCreateMutexSem (NIL, access, 0, FALSE);
        listhead := NIL;
        BNcount := 0;
        bytecount := 0;
    END (*WITH*);
    <* END *>
    Randomize;
END BigNum.

