(**************************************************************************)
(*                                                                        *)
(*  PMOS/2 software library                                               *)
(*  Copyright (C) 2024   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       *)
(*                                                                        *)
(**************************************************************************)

<* M2EXTENSIONS + *>

IMPLEMENTATION MODULE Arith64;

        (********************************************************)
        (*                                                      *)
        (*            Support for 64-bit integers               *)
        (*                                                      *)
        (*  Programmer:         P. Moylan                       *)
        (*  Started:            17 October 2001                 *)
        (*  Last edited:        14 June 2024                    *)
        (*  Status:             Working                         *)
        (*                                                      *)
        (********************************************************)

FROM SYSTEM IMPORT CARD16, CARD32, CAST;

FROM LowLevel IMPORT LS, RS, IAND;

(************************************************************************)
(*                     BIGENDIAN TO/FROM LITTLEENDIAN                   *)
(************************************************************************)

PROCEDURE ToAPI64 (val: CARD64): CARD64API;

    (* Converts a CARD64 value to a CARD64API value. *)

    VAR result: CARD64API;

    BEGIN
        result.low := val.low;
        result.high := val.high;
        RETURN result;
    END ToAPI64;

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

PROCEDURE FromAPI64 (val: CARD64API): CARD64;

    (* Converts a CARD64API value to a CARD64 value. *)

    VAR result: CARD64;

    BEGIN
        result.low := val.low;
        result.high := val.high;
        RETURN result;
    END FromAPI64;

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

PROCEDURE IsZeroLE (val: CARD64LE): BOOLEAN;

    (* Returns TRUE iff val = 0 *)

    BEGIN
        RETURN (val.high = 0) AND (val.low = 0);
    END IsZeroLE;

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

PROCEDURE Compare64LE (A, B: CARD64LE): INTEGER;

    (* Returns >0 if A>B, =0 if A=B, <0 if A<B.  *)

    BEGIN
        IF A.high > B.high THEN
            RETURN +1;
        ELSIF A.high < B.high THEN
            RETURN -1;
        ELSE
            IF A.low > B.low THEN
                RETURN +1;
            ELSIF A.low < B.low THEN
                RETURN -1;
            ELSE
                RETURN 0;
            END (*IF*);
        END (*IF*);
    END Compare64LE;

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

PROCEDURE Negate (VAR (*INOUT*) A: INT64LE);

    (* Negation in place.  *)

    VAR temp: INTEGER;

    BEGIN
        (* Complement and increment. *)

        IF A.low = 0 THEN
            A.high := -A.high;
        ELSE
            A.high := -A.high - 1;
            temp := -CAST(INTEGER, A.low) - 1;
            A.low := CAST(CARDINAL, temp);
        END (*IF*);
    END Negate;

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

PROCEDURE ABS64 (A: INT64LE): CARD64LE;

    (* Absolute value.  *)

    VAR result: CARD64LE;

    BEGIN
        IF A.high < 0 THEN
            Negate (A);
        END (*IF*);
        result.high := A.high;
        result.low := A.low;
        RETURN result;
    END ABS64;

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

PROCEDURE Add64LE (VAR (*INOUT*) A: CARD64LE;  B: CARDINAL);

    (* Computes  A := A + B.  This differs from Sum64 (below) in the    *)
    (* type of B, and in the way the result is returned.                *)

    BEGIN
        IF A.low > MAX(CARDINAL) - B THEN
            A.low := A.low - (MAX(CARDINAL) - B + 1);
            INC (A.high);
        ELSE
            A.low := A.low + B;
        END (*IF*);
    END Add64LE;

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

PROCEDURE Sub64LE (VAR (*INOUT*) A: CARD64LE;  B: CARDINAL);

    (* Computes  A := A - B.  This differs from Diff64 (below) in the   *)
    (* type of B, and in the way the result is returned.                *)

    BEGIN
        IF A.low < B THEN
            A.low := A.low + (MAX(CARDINAL) - B + 1);
            DEC (A.high);
        ELSE
            A.low := A.low - B;
        END (*IF*);
    END Sub64LE;

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

PROCEDURE INC64LE (VAR (*INOUT*) A: CARD64LE);

    (* Increments A by 1.  *)

    BEGIN
        IF A.low = MAX(CARDINAL) THEN
            A.low := 0;
            INC (A.high);
        ELSE
            INC (A.low);
        END (*IF*);
    END INC64LE;

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

PROCEDURE DEC64LE (VAR (*INOUT*) A: CARD64LE);

    (* Decrements A by 1.  *)

    BEGIN
        IF A.low = 0 THEN
            A.low := MAX(CARDINAL);
            DEC (A.high);
        ELSE
            DEC (A.low);
        END (*IF*);
    END DEC64LE;

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

PROCEDURE Sum64LE (A, B: CARD64LE): CARD64LE;

    (* Returns A+B. *)

    VAR result: CARD64LE;

    BEGIN
        result := A;
        IF result.low > MAX(CARDINAL) - B.low THEN
            DEC (result.low, MAX(CARDINAL) - B.low + 1);
            INC (result.high);
        ELSE
            INC (result.low, B.low);
        END (*IF*);
        INC (result.high, B.high);
        RETURN result;
    END Sum64LE;

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

PROCEDURE LongSub64LE (VAR (*INOUT*) A: CARD64LE;  B: CARD64LE);

    (* Computes  A := A - B  *)

    VAR borrow: BOOLEAN;

    BEGIN
        borrow := A.low < B.low;
        IF borrow THEN
            INC (A.low, MAX(CARDINAL) - B.low + 1);
        ELSE
            DEC (A.low, B.low);
        END (*IF*);
        DEC (A.high, B.high);
        IF borrow AND (A.high > 0) THEN
            DEC (A.high);
        END (*IF*);
    END LongSub64LE;

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

PROCEDURE Diff64LE (A, B: CARD64LE): INT64LE;

    (* Returns A-B. *)

    VAR result: INT64LE;  borrow: BOOLEAN;

    BEGIN
        result.high := A.high;
        result.low := A.low;
        borrow := A.low < B.low;
        IF borrow THEN
            INC (result.low, MAX(CARDINAL) - B.low + 1);
        ELSE
            DEC (result.low, B.low);
        END (*IF*);

        DEC (result.high, B.high);
        IF borrow THEN
            DEC (result.high);
        END (*IF*);
        RETURN result;
    END Diff64LE;

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

PROCEDURE ShortSubLE (A, B: CARD64LE): CARDINAL;

    (* Returns A-B as a CARDINAL value, or MAX(CARDINAL) in the case of *)
    (* overflow.                                                        *)

    BEGIN
        LongSub64LE (A, B);
        IF A.high = 0 THEN RETURN A.low
        ELSE RETURN MAX(CARDINAL)
        END (*IF*);
    END ShortSubLE;

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

PROCEDURE Mul32LE (A, B: CARDINAL): CARD64;

    (* Returns A*B. *)

    VAR C: CARD64;

    BEGIN
        (* The compiler firsts generate code to push ebp    *)
        (* and edi, and then makes edi point to where the   *)
        (* function result should go. I believe that these  *)
        (* edi operations will be generated iff the         *)
        (* procedure is a function that returns a result.   *)

        (* Now we get to my code. *)

        ASM
            mov eax, dword ptr [ebp + 0CH]      (* A *)
            mul dword ptr [ebp + 010H]          (* A*B *)
            mov [ebp - 0CH], eax                (* answer to C *)
            mov [ebp - 08H], edx
        END;

        (* Now the compiler-generated code moves the answer *)
        (* from local variable C to where it really should  *)
        (* go.  This double-handling is unfortunate, and if *)
        (* I had my druthers the variable C wouldn't even   *)
        (* exist, but for a function that returns something *)
        (* the compiler won't allow me to avoid having a    *)
        (* high-level RETURN statement.                     *)

        RETURN C;

    END Mul32LE;

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

PROCEDURE OldMul32LE (A, B: CARDINAL): CARD64LE;

    (* Returns A*B. *)

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

    PROCEDURE INC32 (VAR (*INOUT*) X: CARD32; Y: CARD32;
                                        VAR (*OUT*) carry: BOOLEAN);

        (* X := X+Y, but with overflow check. *)

        BEGIN
            IF X > MAX(CARD32) - Y THEN
                X := X - (MAX(CARD32) - Y) - 1;
                carry := TRUE;
            ELSE
                INC (X, Y);  carry := FALSE;
            END (*IF*);
        END INC32;

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

    CONST scale = MAX(CARD16) + 1;

    TYPE twopart =  RECORD
                        CASE :BOOLEAN OF
                            FALSE:  w: CARD32;
                          | TRUE:   h: RECORD
                                           lo, hi: CARD16;
                                       END (*RECORD*);
                        END (*CASE*);
                    END (*RECORD*);

    VAR A2, B2, middle, reslow: twopart;
        result: CARD64LE;  carry: BOOLEAN;

    BEGIN
        A2.w := A;  B2.w := B;
        result.high := A2.h.hi*B2.h.hi;  reslow.w := A2.h.lo*B2.h.lo;

        (* Now add in the "middle" term. *)

        middle.w := reslow.h.hi;
        INC32 (middle.w, A2.h.hi*B2.h.lo, carry);
        IF carry THEN INC (result.high, scale) END (*IF*);
        INC32 (middle.w, A2.h.lo*B2.h.hi, carry);
        IF carry THEN INC (result.high, scale) END (*IF*);

        reslow.h.hi := middle.h.lo;
        result.low := reslow.w;
        INC (result.high, middle.h.hi);

        RETURN result;

    END OldMul32LE;

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

PROCEDURE Mul64LE (A, B: CARD64LE): CARD64LE;

    (* Returns A*B. *)

    CONST scale = 65536;

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

    PROCEDURE ADC (from: CARD16;  VAR (*INOUT*) to: CARD16;
                                  VAR (*INOUT*) carry: BOOLEAN);

        (* Add 'from' to 'to', with carry propagation. *)

        VAR result: CARDINAL;

        BEGIN
            result := from + to;
            IF carry THEN INC(result) END(*IF*);
            carry := result >= scale;
            IF carry THEN DEC(result, scale) END(*IF*);
            to := result;
        END ADC;

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

    VAR ans: ARRAY [0..5] OF CARD16;
        AA, BB: ARRAY [0..3] OF CARD16;
        part: CARDINAL;
        result: CARD64LE;
        j, k, m: [0..5];
        carry: BOOLEAN;

    BEGIN
        (* Clear the ans array. *)

        FOR j := 0 TO 5 DO
           ans[j] := 0;
        END (*FOR*);

        (* Break A and B down into 16-bit sections. *)

        AA[3] := A.high DIV scale;
        AA[2] := A.high MOD scale;
        AA[1] := A.low DIV scale;
        AA[0] := A.low MOD scale;
        BB[3] := B.high DIV scale;
        BB[2] := B.high MOD scale;
        BB[1] := B.low DIV scale;
        BB[0] := B.low MOD scale;

        (* Do the multiplication as a combination of 16-bit     *)
        (* multiplications, since that is what our built-in     *)
        (* multiplication can handle without losing the part    *)
        (* that overflows to the high part.                     *)

        FOR j := 0 TO 3 DO
            FOR k := 0 TO 3-j DO
                part := AA[j]*BB[k];
                carry := FALSE;
                m := j+k;
                ADC (part MOD scale, ans[m], carry);
                ADC (part DIV scale, ans[m+1], carry);
                IF carry THEN INC(ans[m+2]) END(*IF*);
            END (*FOR*);
        END (*FOR*);

        (* Copy the ans array to the final result. *)

        result.low  := scale*ans[1] + ans[0];
        result.high := scale*ans[3] + ans[2];
        RETURN result;

    END Mul64LE;

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

PROCEDURE ShortMul64LE (A: CARD64LE;  B: CARDINAL): CARD64LE;

    (* Returns A*B, for the case where the second operand is a CARDINAL. *)

    VAR BB: CARD64LE;

    BEGIN
        (* It's probably good enough to do this without *)
        (* trying to optimise.                          *)

        BB.low := B;
        BB.high := 0;
        RETURN Mul64LE (A, BB);

    END ShortMul64LE;

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

PROCEDURE LongDiv64LE (A: CARD64LE;  B: CARDINAL;  VAR (*OUT*) quotient: CARD64LE;
                                  VAR (*OUT*) remainder: CARDINAL);

    (* Divides A by B, returns quotient and remainder. *)

    CONST mask = MAX(CARDINAL) DIV 2 + 1;    (* 2^31 *)

    VAR Q, R, BZ: CARD64LE;  Z: CARDINAL;
        carry: BOOLEAN;

    BEGIN
        (* Initial step: get a first approximation to the quotient Q,   *)
        (* and hence an initial value for the remainder R.              *)

        Q.high := A.high DIV B;
        R.high := A.high MOD B;
        Q.low  := 0;
        R.low  := A.low;

        (* Invariant: A = Q*B + R *)
        (* It is not immediately obvious that the above assignments     *)
        (* satisfy the invariant; but I've checked, and they do.        *)

        (* Notation: Z is our current approximation to R DIV B, and     *)
        (* BZ holds the value of B*Z.  If Z is too large we will do     *)
        (* nothing (but reduce the value of Z for the next iteration).  *)
        (* If Z is not too large then we will decrement R by BZ, and    *)
        (* increment Q by Z.                                            *)

        (*  Z := CARD64LE{0,1};  BZ := CARD64LE{0,B};  *)

        (* Observation: at this point in the calculation we             *)
        (* know that R.high < B.  Also, by construction, BZ.high = B.   *)
        (* It follows that R < BZ, so we can pull the first loop        *)
        (* iteration out as a special case.                             *)

        (* Divide Z by 2. *)

        Z := mask;

        (* Divide BZ by 2. *)

        IF ODD(B) THEN
            BZ.low := mask;
        ELSE
            BZ.low := 0;
        END (*IF*);
        BZ.high := B DIV 2;

        (* Now the main reduction loop, where we divide Z by 2 each     *)
        (* time around the loop.                                        *)

        LOOP
            (* Reduction step if R >= BZ.  *)

            IF (R.high > BZ.high) OR
                    ((R.high = BZ.high) AND (R.low >= BZ.low)) THEN
                Add64LE (Q, Z);
                LongSub64LE (R, BZ);
            END (*IF*);

            (* Once R is down to single precision, we can finish the    *)
            (* job outside the loop.                                    *)

            IF R.high = 0 THEN
                EXIT (*LOOP*);
            END (*IF*);

            (* Divide Z by 2. *)

            Z := Z DIV 2;

            (* Divide BZ by 2. *)

            carry   := ODD (BZ.high);
            BZ.high := BZ.high DIV 2;
            BZ.low  := BZ.low DIV 2;
            IF carry THEN
                INC (BZ.low, mask);
            END (*IF*);

        END (*LOOP*);

        (* Final step: once R.high = 0, we can finish the job using     *)
        (* CARD32 arithmetic.  Note that this step also preserves the   *)
        (* invariant.                                                   *)

        Add64LE (Q, R.low DIV B);
        R.low := R.low MOD B;

        quotient := Q;  remainder := R.low;

    END LongDiv64LE;

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

PROCEDURE ShortDivLE (A: CARD64LE;  B: CARDINAL): CARDINAL;

    (* Returns A DIV B as a CARDINAL value, or MAX(CARDINAL) if the     *)
    (* result overflows.                                                *)

    VAR Q: CARD64LE;  R: CARDINAL;

    BEGIN
        LongDiv64LE (A, B, Q, R);
        IF Q.high = 0 THEN RETURN Q.low
        ELSE RETURN MAX(CARDINAL)
        END (*IF*);
    END ShortDivLE;

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

PROCEDURE ShortModLE (A: CARD64LE;  B: CARDINAL): CARDINAL;

    (* Returns A MOD B as a CARDINAL value.  *)

    VAR Q: CARD64LE;  R: CARDINAL;

    BEGIN
        LongDiv64LE (A, B, Q, R);
        RETURN R;
    END ShortModLE;

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

PROCEDURE Div10LE (number: CARD64LE;  VAR (*OUT*) quotient: CARD64LE;
                                  VAR (*OUT*) remainder: CARDINAL);

    (* Divides number by 10, returns quotient and remainder. *)

    BEGIN
        LongDiv64LE (number, 10, quotient, remainder);
    END Div10LE;

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

PROCEDURE RS64LE (VAR (*INOUT*) number: CARD64LE;  N: CARDINAL;
                                         VAR (*OUT*) remainder: CARDINAL);

    (* Right shift by N bits, i.e. division by 2^N, with remainder      *)
    (* holding the lost bits.  N must be in range [0..32].              *)

    VAR maskR: CARDINAL;

    BEGIN
        IF N = 32 THEN
            remainder := number.low;
            number.low := number.high;
            number.high := 0;
        ELSIF N = 0 THEN
            remainder := 0;
        ELSE
            maskR := RS (MAX(CARDINAL), 32-N);
            remainder := IAND (number.low, maskR);
            number.low := RS (number.low, N) + LS (IAND(number.high, maskR), 32-N);
            number.high := RS (number.high, N);
        END (*IF*);
    END RS64LE;

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

PROCEDURE FLOAT64LE (number: CARD64LE): REAL;

    (* Converts CARD64 to REAL. *)

    CONST scale = FLOAT(MAX(CARDINAL)) + 1.0;

    BEGIN
        RETURN scale*FLOAT(number.high) + FLOAT(number.low);
    END FLOAT64LE;

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

PROCEDURE ORDL (number: INT64): CARDINAL;

    (* Converts INT64 to CARDINAL. *)

    BEGIN
        IF number.high <> 0 THEN RETURN MAX(CARDINAL)
        ELSE RETURN number.low;
        END (*IF*)
    END ORDL;

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

END Arith64.

