PHP FixedPoint64 rough sketch

<?php
declare(strict_types=1);

namespace Kickback\Math;

final class Int64
{
    private const MASK_SIGN = (int)(1 << 63);

    public static function add_unsigned_with_carry(int $a, int $b, int &$carry) : int
    {
        // Optimized version thanks to this page:
        // https://grack.com/blog/2022/12/20/deriving-a-bit-twiddling-hack/
        // (It has been adjusted for 64-bit arithmatic.)
        $c = $a + $b;
        $carry = ((($c ^ $a) & ($c ^ $b)) >> 63);
        // Note that $carry will either be 0 or it will be 0xFFFF...FFFF.
        // But we want it to be 0 or 1.
        // So we can clear the sign extension bits by ANDing with 1:
        $carry &= 1;
        return $c;
    }
    /* TODO: cut?
        // Since 0 is false and anything non-0 is considered "true", then
        // this will return what we want, without having to worry about
        // somehow reducing this to the set of {0,1}.
        */

    public static function unittest_add_unsigned_with_carry() : void
    {
        echo("unittest Int128::unittest_add_unsigned_with_carry");

        // Shorthand to avoid having to write out long 64-bit hex values.
        $x0000 = (int)0;
        $x0001 = (int)1;
        $x0002 = (int)2;
        $x8000 = PHP_INT_MIN;
        $x8001 = PHP_INT_MIN+1;
        $xFFFF = (int)-1;
        $xFFFE = (int)-2;
        $x7FFF = PHP_INT_MAX;
        $x7FFE = PHP_INT_MAX-1;

        $carry = (int)0;

        assert($x0000 === self::add_unsigned_with_carry($x0000,$x0000,$carry)); assert(0 === $carry);
        assert($x0001 === self::add_unsigned_with_carry($x0000,$x0001,$carry)); assert(0 === $carry);
        assert($x0001 === self::add_unsigned_with_carry($x0001,$x0000,$carry)); assert(0 === $carry);
        assert($x0002 === self::add_unsigned_with_carry($x0001,$x0001,$carry)); assert(0 === $carry);

        assert($x0000 === self::add_unsigned_with_carry($x0000,$x0000,$carry)); assert(0 === $carry);
        assert($xFFFF === self::add_unsigned_with_carry($x0000,$xFFFF,$carry)); assert(0 === $carry);
        assert($xFFFF === self::add_unsigned_with_carry($xFFFF,$x0000,$carry)); assert(0 === $carry);
        assert($xFFFE === self::add_unsigned_with_carry($xFFFF,$xFFFF,$carry)); assert(1 === $carry);

        assert($x0000 === self::add_unsigned_with_carry($x0001,$xFFFF,$carry)); assert(1 === $carry);
        assert($x7FFF === self::add_unsigned_with_carry($xFFFF,$x8000,$carry)); assert(1 === $carry);

        assert($x0001 === self::add_unsigned_with_carry($x0002,$xFFFF,$carry)); assert(1 === $carry);
        assert($x7FFE === self::add_unsigned_with_carry($xFFFE,$x8000,$carry)); assert(1 === $carry);

        assert($x8000 === self::add_unsigned_with_carry($x8001,$xFFFF,$carry)); assert(1 === $carry);
        assert($x7FFE === self::add_unsigned_with_carry($x7FFF,$xFFFF,$carry)); assert(1 === $carry);

        echo(" done.\n");
    }

    public static function ensure_ints_are_64bit() : void
    {
        $num_int_bits = PHP_INT_SIZE*8;
        if ($num_int_bits < 64)
        {
            throw new SomethingSomethingException(
                "This system features ${num_int_bits}-bit integers, " .
                "but at least 64-bit wide integers are required for financial calculations.");
        }
    }

    public static function unittest() : void
    {
        ensure_ints_are_64bit();
        unittest_add_unsigned_with_carry();
        echo("\n");
    }

    // Prevent instantiation/construction of the (static/constant) class.
    /** @return never */
    private function __construct() {
        throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
    }
}
?>

<?php
declare(strict_types=1);
namespace Kickback\Math;

use \Kickback\Math\Int64;

final class Int128
{
    private const MASK_LO = (int)((1 << 32) - 1);
    private const MASK_SIGN = (int)(1 << 63);

    // Internal version of `cmp`.
    // This version does not clamp it's outputs to the set {-1, 0, 1}.
    private static function cmp_unsaturated(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
    {
        $res = $a_hi - $b_hi;
        if ( $res ) {
            return $res;
        }

        return ($a_lo - $b_lo);
    }

    /*
    *  Returns -1, 0, or 1 based on a comparison of the given 128-bit integers.
    *
    *  Defined such that:
    *  * a < b implies cmp(a,b) < 0
    *  * a > b implies cmp(a,b) > 0
    *  * a == b implies cmp(a,b) == 0.
    */
    public static function cmp(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
    {
        $res = self::cmp_unsaturated($a_hi,$a_lo,  $b_hi,$b_lo);
        if ( $res === 0x8000_0000_0000_0000 ) {
            return -1;
        }
        $res--;      //    >0    =0    <0
        $res >>= 62; // {0b000,0b111,0b110} = {0,-1,-2}
        $res++;      // {0b001,0b000,0b111} = {1, 0,-1}
        return $res;
/* TODO: cut
        $res--;      // {0b000,0b111,0b110,0b101}
        $res &;      // {0b000,0b000,0b110,0b100}
        want         // {0b011,0b001,0b000,0b011}   {0b11,0b01,0b00,0b11}

        $res += 3    // {0b100,0b011,0b110,0b101}
        $res >>= 1;  // {0b110,0b001,0b111,0b110}
        $res += 1;   // {0b111,0b010,0b000,0b111}
        $res >>= 1;  // {0b111,0b001,0b000,0b111}

        $res >>= 62; // {0b001,0b000,0b111,0b110} ^ {0b00,0b01,0b10,0b11} | {0b11,0b10,0b01,0b10}
        $res += 1    // {0b010,0b001,0b000,0b111} ^ {0b11,0b00,0b01,0b10} | {0b10,0b11,0b10,0b01}
        $res >>= 1;  // {0b001,0b000,0b001,0b111} ^ {0b11,0b10,0b00,0b11}
        $res += 2    // {0b011,0b010,0b011,0b001} ^ {0b10,0b11,0b00,0b01} | {0b01,0b10,0b11,0g10}
        $res += 3    // {0b100,0b011,0b110,0b101} ^ {0b01,0b10,0b11,0b00} | {0}

        $res >>= 62; // {0b001,0b000,0b111,0b110} ^ {0b00,0b01,0b10,0b11} | {0b11,0b10,0b01,0b10}
        $res += 1    // {0b010,0b001,0b000,0b111} ^ {0b11,0b00,0b01,0b10} | {0b10,0b11,0b10,0b01}
        $res += 2    // {0b011,0b010,0b001,0b000} ^ {0b10,0b11,0b00,0b01} | {0b01,0b10,0b11,0g10}
        $res += 3    // {0b100,0b011,0b010,0b001} ^ {0b01,0b10,0b11,0b00} | {0}

        $res >>= 62; // {0b01,0b00,0b11,0b10} - {0b11,0b00,0b01,0b10} | {0b11,0b10,0b01,0b10}
        $res += 1    // {0b10,0b01,0b00,0b11} - {0b10,0b11,0b00,0b01} | {0b10,0b11,0b10,0b01}
        $res += 2    // {0b11,0b10,0b01,0b00} - {0b01,0b10,0b11,0b00} | {0b01,0b10,0b11,0g10}
        $res += 3    // {0b00,0b11,0b10,0b01} - {0b00,0b01,0b10,0b11} | {0}

        $res >>= 62; // {0b01,0b00,0b11,0b10} - {0b11,0b00,0b01,0b10} | {0b11,0b00,0b11,0b10}
        $res += 1    // {0b10,0b01,0b00,0b11} - {0b10,0b11,0b00,0b01} | {0b10,0b11,0b00,0b11}
        $res += 2    // {0b11,0b10,0b01,0b00} - {0b01,0b10,0b11,0b00} | {0b11,0b10,0b11,0b00}
        $res += 3    // {0b00,0b11,0b10,0b01} - {0b00,0b01,0b10,0b11} | {0b00,0b11,0b10,0b11}

        if ( !$res ) {
            return $res;
        }

        $res >>= 63;
        $res <<= 2;
        return 0x
        */

        /*
        if ( $a_hi < $b_hi ) {
            return -1;
        } else
        if ( $a_hi > $b_hi ) {
            return 1;
        }
        assert($a_hi === $b_hi);

        if ( $a_lo < $b_lo ) {
            return -1;
        } else
        if ( $a_lo > $b_lo ) {
            return 1;
        }
        assert($a_lo === $b_lo);

        return 0;
        */
    }

    public static function unittest_cmp()
    {
        echo("unittest Int128::cmp");

        assert(self::cmp( 0, 0,   0, 0) === 0);
        assert(self::cmp( 0, 0,   0, 1)  <  0);
        assert(self::cmp( 0, 1,   0, 0)  >  0);
        assert(self::cmp( 0, 1,   0, 1) === 0);
        assert(self::cmp(-1, 0,  -1, 0) === 0);
        assert(self::cmp(-1, 0,   0, 0)  <  0);
        assert(self::cmp(-1, 0,   1, 0)  <  0);
        assert(self::cmp( 0, 0,  -1, 0)  >  0);
        assert(self::cmp( 0, 0,   0, 0) === 0);
        assert(self::cmp( 0, 0,   1, 0)  <  0);
        assert(self::cmp( 1, 0,  -1, 0)  >  0);
        assert(self::cmp( 1, 0,   0, 0)  >  0);
        assert(self::cmp( 1, 0,   1, 0) === 0);

        // Ensure that it's returning values in the set {-1,0,1}, and not in the set {-N,0,N} or somesuch.
        assert(self::cmp(  0, 0,   16, 0) === -1);
        assert(self::cmp( 16, 0,    0, 0) ===  1);
        assert(self::cmp(-16, 0,    0, 0) === -1);
        assert(self::cmp(  0, 0,  -16, 0) ===  1);

        // Notably, we've carefully avoided negative values in the less-significant parts, so far.
        // That's because the less-significant integers shall be treated as
        // unsigned integers, but PHP only has _signed_ comparison.
        // So we better check for mistakes of that kind!
        $x0000 = (int)0;
        $xFFFF = (int)-1;
        assert(self::cmp($x0000,$x0000,  $x0000,$x0000) === 0); // 0 === 0
        assert(self::cmp($x0000,$x0000,  $x0000,$xFFFF)  <  0); // 0 < (2^32 - 1)
        assert(self::cmp($x0000,$x0000,  $xFFFF,$x0000)  >  0); // 0 > -(2^32)
        assert(self::cmp($x0000,$x0000,  $xFFFF,$xFFFF)  >  0); // 0 > -1
        assert(self::cmp($x0000,$xFFFF,  $x0000,$x0000)  >  0); // (2^32 - 1) > 0
        assert(self::cmp($x0000,$xFFFF,  $x0000,$xFFFF) === 0); // (2^32 - 1) === (2^32 - 1)
        assert(self::cmp($x0000,$xFFFF,  $xFFFF,$x0000)  >  0); // (2^32 - 1) > -(2^32)
        assert(self::cmp($x0000,$xFFFF,  $xFFFF,$xFFFF)  >  0); // (2^32 - 1) > -1
        assert(self::cmp($xFFFF,$x0000,  $x0000,$x0000)  <  0); // -(2^32) < 0
        assert(self::cmp($xFFFF,$x0000,  $x0000,$xFFFF)  <  0); // -(2^32) < (2^32 - 1)
        assert(self::cmp($xFFFF,$x0000,  $xFFFF,$x0000) === 0); // -(2^32) === -(2^32)
        assert(self::cmp($xFFFF,$x0000,  $xFFFF,$xFFFF)  <  0); // -(2^32) < -1
        assert(self::cmp($xFFFF,$xFFFF,  $x0000,$x0000)  <  0); // -1 < 0
        assert(self::cmp($xFFFF,$xFFFF,  $x0000,$xFFFF)  <  0); // -1 < (2^32 - 1)
        assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$x0000)  >  0); // -1 > -(2^32)
        assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$xFFFF) === 0); // -1 === -1

        echo(" done.\n");
    }

    /*
    * Multiplies two 64-bit integers, resulting in a single 128-bit integer.
    *
    * Because PHP does not have a native 128-bit integer type (or the ability
    * to define structs that can be returned from functions), the returned
    * value is placed into two 64-bit integer reference-type parameters.
    */
    public static function multiply_64x64(int $a, int $b, int &$out_hi, int &$out_lo) : void
    {
        // Reorganize the two 64-bit integers into four 32-bit integers.
        // This helps because now the 32-bit integers will have at least
        // 32 bits of "headroom" for doing multiplications.
        $a_lo = $a & self::MASK_LO;
        $a_hi = $a >> 32;
        $b_lo = $b & self::MASK_LO;
        $b_hi = $b >> 32;
        self::multiply_NxN(32, $a_hi,$a_lo,  $b_hi,$b_lo,  $out_hi,$out_lo);
    }

    /*
    * It is probably better to use `multiply_64x64` in most cases.
    *
    * This is the underlying implementation of `multiply_64x64`, and having
    * it separated out as a different function is useful for testing purposes.
    */
    public static function multiply_NxN(int $digit_bit_width,  int $a_hi, int $a_lo, int $b_hi, int $b_lo, int &$out_hi, int &$out_lo) : void
    {
        assert($bit_width <= 32);
        $mask_lo = ((1 << $digit_bit_width) - 1);
        $max_digit = $mask_lo;

        // Multiplies. We need 4 of them for a school-style long multiply.
        // There are faster methods, but this will do for now.
        //
        // Graphically:
        //                a_hi  a_lo
        //             x  b_hi  b_lo
        //             -------------
        //                x_hi  x_lo
        //          y_hi  y_lo
        //          z_hi  z_lo
        // +  w_hi  w_lo
        // -----------------------
        //    ????  ????  ????  ????
        //
        // Note that these multiplications, or at least the ones involving
        // the $*_lo operands, must be _unsigned_ multiplications.
        // However, because these are 32-bit values stored in 64-bit variables,
        // the sign bit will _never_ be set, so it won't matter if we
        // use a 64-bit signed multiplication on them.
        //
        $x = $a_lo * $b_lo;
        $y = $a_lo * $b_hi;
        $z = $a_hi * $b_lo;
        $w = $a_hi * $b_hi;

        // To make the logic more clear, we will also define the
        // variables corresponding to the {x,y,z,w}_{hi,lo} placeholders.
        $x_lo = $x & $mask_lo;
        $y_lo = $y & $mask_lo;
        $z_lo = $z & $mask_lo;
        $w_lo = $w & $mask_lo;

        $x_hi = $x >> $digit_bit_width;
        $y_hi = $y >> $digit_bit_width;
        $z_hi = $z >> $digit_bit_width;
        $w_hi = $w >> $digit_bit_width;

        // PHP doesn't provide an unsigned right-shift, so we must clear
        // any sign-extended bits in things that were right-shifted.
        $x_hi = $x & $mask_lo;
        $y_hi = $y & $mask_lo;
        $z_hi = $z & $mask_lo;
        $w_hi = $w & $mask_lo;

        // Calculate the bottom row of 32-bit integers.
        //
        // Note that some of these additions might "overflow", but
        // because we are storing our 32-bit integers in 64-bit variables,
        // the carry values will be captured.
        //
        // These will get shuffled into the output integers
        // once the accumulation is done.
        //
        // Graphically:
        //                a_hi  a_lo
        //             x  b_hi  b_lo
        //             -------------
        //                x_hi  x_lo
        //          y_hi  y_lo
        //          z_hi  z_lo
        // +  w_hi  w_lo
        // -----------------------
        //    out3  out2  out1  out0
        //
        $out0 = $x_lo;
        $out1 = $x_hi + $y_lo + $z_lo;
        $out2 = $y_hi + $z_hi + $w_lo;
        $out3 = $w_hi;

        // Perform carry operations.
        // (Beware: order-of-operations is important.)
        if ( $out1 > $max_digit )
        {
            $out2 += ($out1 >> $digit_bit_width);
            $out1 &= $mask_lo;
        }

        if ( $out2 > $max_digit )
        {
            $out3 += ($out2 >> $digit_bit_width);
            $out2 &= $mask_lo;
        }
        // Note that $out3 is involved in an addition, but won't generate
        // a carry-out. It won't overflow, for math reasons.

        // Now we have four proper 32-bit integers (two pairs) with no carry bits,
        // so we can concatenate the pairs to get two 64-bit integers and have our 128-bit result.
        $lo = $out0 + ($out1 << $digit_bit_width);
        $hi = $out2 + ($out3 << $digit_bit_width);

        // Return our results from the function.
        $out_lo = $lo;
        $out_hi = $hi;
    }

    // This is not part of the public API because there is no point.
    // It _should_ do exactly what the `*` operator does.
    // It's just designed to use the `multiply_64x64` function, so that
    // we can see if it gives results identical to `*`, at least whenever
    // the two are given operands that don't overflow.
    private static function multiply_NxN_lo(int $bit_width, int $a, int $b) : int
    {
        $out_lo = (int)0;
        $out_hi = (int)0;
        self::multiply_NxN($bit_width,  $a,  $b,  $out_hi,$out_lo);
        return $out_lo;
    }
    /* TODO: cut?
    // Uses the `multiply_32H32Lx32H32L` function to emulate what
    // `multiply_64x64` would do if it had lower bit-width integers available
    // (ex: `multiply_32x32`).
    //
    // Because this is used for testing purposes, and because we've asserted
    // that we'll always use PHP on 64-bit systems and `int` will always have
    // a width of at least 64 bits, then the result of `multiply_32x32` and
    // smaller will always fit in a 64-bit return value. So we don't need
    // to return the output in reference parameters for this one, which is
    // very helpful for testing.
    private static function multiply_NxN_emul(int $bit_width, int $a, int $b) : int
    {
        assert($bit_width <= 32);
        $mask_lo = ((1 << $bit_width) - 1);
        $a_lo = $a & $mask_lo;
        $a_hi = $a >> $bit_width;
        $b_lo = $b & $mask_lo;
        $b_hi = $b >> $bit_width;
        $out_lo = (int)0;
        $out_hi = (int)0;
        self::multiply_32H32Lx32H32L($a_hi,$a_lo,  $b_hi,$b_lo,  $out_hi,$out_lo);

        $int_max = $mask_lo;
        assert($out_hi <= $int_max);
        assert($out_lo <= $int_max);


    }
*/

    public static function unittest_multiply_64x64() : void
    {
        echo("unittest Int128::unittest_multiply_64x64");

        assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x00));
        assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x01));
        assert(0x0000 === self::multiply_NxN_lo(8, 0x01, 0x00));
        assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0xFF));
        assert(0x0000 === self::multiply_NxN_lo(8, 0xFF, 0x00));
        assert(0x0001 === self::multiply_NxN_lo(8, 0x01, 0x01));
        assert(0x000F === self::multiply_NxN_lo(8, 0x01, 0x0F));
        assert(0x000F === self::multiply_NxN_lo(8, 0x0F, 0x01));
        assert(0x00E1 === self::multiply_NxN_lo(8, 0x0F, 0x0F));
        assert(0x0010 === self::multiply_NxN_lo(8, 0x01, 0x10));
        assert(0x0010 === self::multiply_NxN_lo(8, 0x10, 0x01));
        assert(0x0100 === self::multiply_NxN_lo(8, 0x10, 0x10));
        assert(0x4000 === self::multiply_NxN_lo(8, 0x80, 0x80));
        assert(0x3F01 === self::multiply_NxN_lo(8, 0x7F, 0x7F));
        assert(0xFC04 === self::multiply_NxN_lo(8, 0xFE, 0xFE));
        assert(0xFD02 === self::multiply_NxN_lo(8, 0xFE, 0xFF));
        assert(0xFD02 === self::multiply_NxN_lo(8, 0xFF, 0xFE));
        assert(0xFE01 === self::multiply_NxN_lo(8, 0xFF, 0xFF));
        assert(0x7E81 === self::multiply_NxN_lo(8, 0x7F, 0xFF));
        assert(0x7F80 === self::multiply_NxN_lo(8, 0x80, 0xFF));
        assert(0x3F80 === self::multiply_NxN_lo(8, 0x80, 0x7F));

        for ( $i = (int)0; $i < 256; $i++ )
        {
            for ( $j = (int)0; $j < 256; $j++ )
            {
                $expected = (int)($i*$j);
                $got = self::multiply_NxN_lo(8, $i, $j);
                assert($expected === $got, sprintf("Operands: i=%02x; j=%02x; Expected: %04x; Got: %04X", $i, $j, $expected, $got));
            }
        }

        assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0000));
        assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0001));
        assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0001, 0x0000));
        assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0xFFFF));
        assert(0x0000_0000 === self::multiply_NxN_lo(16, 0xFFFF, 0x0000));
        assert(0x0000_0001 === self::multiply_NxN_lo(16, 0x0001, 0x0001));
        assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0x0001, 0xFFFF));
        assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0xFFFF, 0x0001));
        assert(0xFFFE_0001 === self::multiply_NxN_lo(16, 0xFFFF, 0xFFFF));
        assert(0XA518_60E1 === self::multiply_NxN_lo(16, 0xF00F, 0xB00F));

        assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0000));
        assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0001));
        assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0000));
        assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0xFFFF_FFFF));
        assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0000));
        assert(0x0000_0000_0000_0001 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0001));
        assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0x0000_0001, 0xFFFF_FFFF));
        assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0001));
        assert(0xFFFF_FFFE_0000_0001 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0xFFFF_FFFF));

        // Explicit test of 128-bit returning version, and of 64-bit inputs.
        $a64 = (int)0xAceBe11e_CafeBabe;
        $b64 = (int)0xF100fCa7_F1edFa57;
        $out_hi = (int)0;
        $out_lo = (int)0;
        self::multiply_64x64($a64,  $b64,  $out_hi,$out_lo);
        assert(0x8E9C_2945_7ED5_0292 === $out_lo);
        assert(0xA2CA_B997_9FFE_C71C === $out_hi);

        echo(" done.\n");
    }

    /**
    * @param      int  $numerator_hi
    * @param      int  $numerator_lo
    * @param      int  $denominator
    * @param-out  int  $quotient_hi
    * @param-out  int  $quotient_lo
    * @returns    int  The remainder.
    */
    public static function divide_128x64(int $numerator_hi, int $numerator_lo, int $denominator, int &$quotient_hi, int &$quotient_lo) : int
    {
        // TODO: Test what PHP does with integer divide-by-zero.
        // In principle, PHP should throw a DivisionByZeroError whenever this
        // happens. In that case, we can just proceed as normal, because
        // the first time division is used in this function, it'll cause
        // that error to be thrown.
        // See: https://www.php.net/manual/en/class.divisionbyzeroerror.php
        // (But I haven't tested to make _sure_ PHP does this. If it doesn't
        // throw an exception, then we should. Propagating weird values
        // through the calculation could result in hard-to-find bugs.)

        // Special-case the value denominator value of `1`, both to make this
        // identity operation be fast, but also to ensure that we don't have
        // to worry about any weird corner cases for denominators that are
        // less than 2. (It can matter, because binary arithmatic.)
        if ( $denominator === 1 ) {
            $quotient_hi = $numerator_hi;
            $quotient_lo = $numerator_lo;
            return 0;
        }

        // Give these predictable values.
        // This will be helpful in debugging if an error happens:
        // you'll always know what these started as!
        $quotient_hi = (int)0;
        $quotient_lo = (int)0;

        // Store the sign bit that the result will have.
        $sign = ($numerator_hi ^ $denominator) & self::MASK_SIGN;

        // Sign extraction, so it doesn't mess with our division.
        $sign = $hi & self::MASK_SIGN;
        $hi &= ~self::MASK_SIGN;

        // Extract remainders.
        //
        // Right now we are mostly interested in $hi_remainder.
        // This will be necessary later for estimating the low portion of the result.
        //
        // The $lo_remainder will simply be returned by the function as-is,
        // since this may help the caller with rounding logic.
        $lo_remainder = $lo % $denominator;
        $hi_remainder = $hi % $denominator;

        // Component-wise division.
        //
        // (We use `intdiv` because `/` might do a check to see if one or both
        // of its operands are floating point numbers, because it is a floating
        // point operation by default. `intdiv`, on the other hand, is explicitly
        // intended for integer division, so it may perform/behave better for this use.)
        $lo = intdiv($lo, $denominator);
        $hi = intdiv($hi, $denominator);

        // Calculate the carry.
        $min_carry = PHP_INT_MAX; // Basically ((2^64 / 2) - 1). We'd prefer (2^64), but can't represent it.
        $min_carry = intdiv($min_carry, $denominator);
        $min_carry <<= 1; // Undo the earlier divide-by-2 in the ((2^64 / 2) - 1) = PHP_INT_MAX.

        // Perform the carry.
        $lo += ($hi_remainder * $min_carry);

        // For safety reasons, we'll also calculate a "max" carry.
        $max_carry_init = (1 << 62);
        $max_carry = intdiv($max_carry_init, $denominator);
        if ( ($denominator * $max_carry) !== $max_carry_init ) {
            // Always round up.
            $max_carry++;
        }
        $max_carry <<= 2;
        $max_carry += 3; // Always round up.

        // Perform the carry.
        $max_lo += ($hi_remainder * $max_carry);

        // This loop takes our approximation and improves it until we have
        // the correct quotient.
        $test_lo = (int)0;
        $test_hi = (int)0;
        $test_carry = (int)0;
        self::multiply_64x64($denominator,  $lo,  $test_hi,$test_lo);
        $test_hi += ($hi * $denominator);

        while ( true )
        {
            $test_lo = Int64::add_unsigned_with_carry($test_lo, $denominator, $test_carry);
            $test_hi += $test_carry;
            // if ( test > numerator )
            if ( self::cmp($test_hi, $test_lo, $numerator_hi, $numerator_lo) > 0 ) {
                break;
            }

            // The additional increment of $denominator in the $test variable
            // corresponds to incrementing the $lo value by 1.
            $lo++;

            // This handles the aforementioned "safety reasons".
            // It prevents us from infinite-looping, or from trying to
            // iterate most of a 64-bit integer's possible values
            // when it is already futile to get the correct answer.
            // Ideally, this will never be executed.
            // If the surrounding function works as intended, then the loop
            // will ALWAYS converge before this condition becomes true.
            if ($lo > $max_lo) {
                $class_name = __CLASS__;
                $func_name = __FUNCTION__;
                throw new \Exception(sprintf(
                    "Internal error: `$class_name::$func_name` did not converge when it should have. ".
                    "Aborting to prevent excessive looping and CPU drain. ".
                    "This means that the `$func_name` function is broken and may return incorrect results, ".
                    "even when this error isn't thrown. Please report this error. ".
                    "Parameters: {numerator_hi:%016X, numerator_lo:%016X, denominator:%016X}",
                    $numerator_hi, $numerator_lo, $denominator
                    ));
            }
        }

        // Pass results to the caller.
        $quotient_hi = $hi;
        $quotient_lo = $lo;
        return $lo_remainder;

        /*
        TODO: Comment is kinda sorta relevant but also outdated already.

        // This will be multiplied by the $hi_remainder to get the portion
        // of the number that is being carried into the $lo component of the result.
        //
        // Note that this should be `(2^64/$denominator)` EXACTLY.
        // Unfortunately, we can't be exact for two reasons:
        // * Rational numbers (the exact result) can't be represented by finite integers.
        // * The number 2^64 is out of the range of 64-bit integers (one more than the max!)
        //
        // Now we'll get clever, and our function will get slower.
        //
        // We notice that the multiplier doesn't _need_ to be exact.
        // We can make a guess. And the guess can be wrong. As long as it's
        // always wrong in a way that's a little bit low.
        //
        // If it's low, it will result in our end-result being too low.
        // Then we can multiply our guess by the denominator.
        // The result will be less than the numerator.
        // (If we had exact numbers instead of integers, multiplying the quotient
        // by the denominator would result in _exactly_ the numerator.)
        //
        // Then we can just add the $denominator over and over, checking
        // it with a multiplication, until it is too high.
        // Once it is too high, we return the last number that multiplied
        // to a value
        //
        */
    }

    public static function unittest_divide_128x64()
    {
        echo("unittest Int128::unittest_divide_128x64");
        // TODO: Oh god this needs testing BADLY.
        echo(" UNIMPLEMENTED! (Please fix!)\n");
        // echo(" done.\n");
    }

    public static function unittest()
    {
        unittest_cmp();
        unittest_multiply_64x64();
        unittest_divide_128x64();
        echo("\n");
    }

    // Prevent instantiation/construction of the (static/constant) class.
    /** @return never */
    private function __construct() {
        throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
    }
}
?>

<?php
declare(strict_types=1);

namespace Kickback\Math;

// TODO: I forget if it's better to do this, or to use an `enum` type. Whatevs.
final class RoundingMode
{
    // Rounding modes.
    // (I think this is all of the possibilities? TODO: Verify.)
    public const ALL_TOWARDS_NEG_INF  =  0; // Default; mimics x86 integer truncation. (I think? TODO: Double-check this.)
    public const ALL_TOWARDS_POS_INF  =  1;
    public const ALL_TOWARDS_ZERO     =  3;
    public const ALL_AWAY_FROM_ZERO   =  4;
    public const ALL_TOWARDS_EVEN     =  5;
    public const ALL_TOWARDS_ODD      =  6;
    public const HALF_TOWARDS_NEG_INF =  7;
    public const HALF_TOWARDS_POS_INF =  8;
    public const HALF_TOWARDS_ZERO    =  9; // Similar to PHP_ROUND_HALF_DOWN
    public const HALF_AWAY_FROM_ZERO  = 10; // Similar to PHP_ROUND_HALF_UP
    public const HALF_TOWARDS_EVEN    = 11; // Similar to PHP_ROUND_HALF_EVEN; bankers' rounding.
    public const HALF_TOWARDS_ODD     = 12; // Similar to PHP_ROUND_HALF_ODD

    // This shall mimic x86 integer truncation.
    public const DEFAULT = self::ALL_TOWARDS_NEG_INF;

    public static function rounding_mode_to_string(int $rounding_mode) : string
    {
        switch($rounding_mode)
        {
            case self::ALL_TOWARDS_NEG_INF : return "ALL_TOWARDS_NEG_INF";
            case self::ALL_TOWARDS_POS_INF : return "ALL_TOWARDS_POS_INF";
            case self::ALL_TOWARDS_ZERO    : return "ALL_TOWARDS_ZERO";
            case self::ALL_AWAY_FROM_ZERO  : return "ALL_AWAY_FROM_ZERO";
            case self::ALL_TOWARDS_EVEN    : return "ALL_TOWARDS_EVEN";
            case self::ALL_TOWARDS_ODD     : return "ALL_TOWARDS_ODD";
            case self::HALF_TOWARDS_NEG_INF: return "HALF_TOWARDS_NEG_INF";
            case self::HALF_TOWARDS_POS_INF: return "HALF_TOWARDS_POS_INF";
            case self::HALF_TOWARDS_ZERO   : return "HALF_TOWARDS_ZERO";
            case self::HALF_AWAY_FROM_ZERO : return "HALF_AWAY_FROM_ZERO";
            case self::HALF_TOWARDS_EVEN   : return "HALF_TOWARDS_EVEN";
            case self::HALF_TOWARDS_ODD    : return "HALF_TOWARDS_ODD";
        }

        return "Unknown rounding mode ($rounding_mode)";
    }

    // Prevent instantiation/construction of the (static/constant) class.
    /** @return never */
    private function __construct() {
        throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
    }
}
?>

<?php
declare(strict_types=1);

namespace Kickback\Math\Types;

use \Kickback\Math\Int128;
use \Kickback\Math\RoundingMode;

final class FixedPoint64
{
    private int numerator;
    private int divisor;

    private static function banker_round(int $numerator, int $increment) : int
    {
        assert($increment > 1);

        // Even if we don't use $even_increment in most cases,
        // we should still assert that it fits in _all_ cases, just to give
        // this function consistent behavior (e.g. whether it asserts or not
        // should depend on `$increment` and not on `$numerator`.)
        //$even_increment = ($increment<<1);
        assert(($increment<<1) <= PHP_INT_MAX);

        if ( $numerator >= 0 ) {
            return self::banker_round_unsigned_($numerator, $increment);
        }
        else {
            return self::banker_round_signed_($numerator, $increment);
        }
    }

    private static function banker_round_signed_(int $numerator, int $increment) : never
    {
        // TODO!
        assert(false, "Signed rounding is not yet implemented.");
    }

    private static function banker_round_unsigned_(int $numerator, int $increment) : int
    {
        /* TODO: cut?
        // We calculate the $round_down_amount (conventional rounding)
        // from the $round_down_even_amount (bankers' rounding)
        // so that we avoid having to modulus more than once.
        $round_down_even_amount = $numerator % $even_increment;
        $round_down_amount = $round_down_even_amount;
        if ( $round_down_amount > $increment ) {
            $round_down_amount -= $increment;
        }
        */

        // First, attempt conventional rounding (e.g. "round-to-nearest-increment").
        $rounding_amount = (int)0;
        if ( self::get_rounding_amount_($numerator, $increment, $rounding_amount) ) {
            return $rounding_amount;
        } else {
            // If that didn't work, then fall back on the tie breaker.
            return self::banker_round_tie_breaker_($numerator, $increment);
        }
    }

    private static function get_rounding_amount_(int $numerator, int $increment, &$rounding_amount) : bool
    {
        $round_down_amount = $numerator % $increment;
        $round_up_amount = $increment - $round_down_amount;
        if ( $round_down_amount < $round_up_amount ) {
            $rounding_amount = -$round_down_amount;
            return true;
        } else
        if ( $round_down_amount > $round_up_amount ) {
            $rounding_amount = $round_up_amount;
            return true;
        } else {
            // If that didn't work, then there's a problem, and it's probably this tie:
            assert($round_down_amount === $round_up_amount);
            $rounding_amount = 0; // TODO: Is this right?
            return false;
        }
    }

    // This is set out into its own function because the tie breaker
    // is unlikely to be executed very frequently. As such, it might
    // be faster (code-execution-wise) to break this out into it's own function.
    //
    // Reasons why this might be the case:
    // * Interpreter is less likely to need to parse/analyse this code when `banker_round` is called.
    // * This is more cache-friendly: the CPU won't need to load the instructions for the tie breaker every time `banker_round` is called.
    // * This makes `banker_round` smaller and therefore more inline-able.
    //
    // I am not sure how much any of those _actually_ matter, and it's not
    // really worth testing at the moment, but this is _usually_ a good way
    // to optimize code, and it seems to have few or no disadvantages.
    private static function banker_round_tie_breaker_(int $numerator, int $increment) : int
    {
        // Now the bankers' rounding comes into play.
        // We break the tie by rounding to the nearest _even_ increment.
        $even_increment = $increment << 1;

        // This involves just doing the rounding again, but with $increment*2.
        $rounding_amount = (int)0;
        assert(self::get_rounding_amount_($numerator, $even_increment, $rounding_amount));
        return $rounding_amount;
    }
    // TODO: The above bankers' rounding code is pretty sus. I was very sleepy while writing it. o.O

    private static function unittest_banker_round() : void
    {
        echo("unittest FixedPoint64::unittest_banker_round");
        // TODO: This needs testing if you want to trust any financials done with it.
        echo(" UNIMPLEMENTED! (Please fix!)\n");
        // echo(" done.\n");
    }

    /*
    *  Beware: This operation is NOT commutative like normal multiplication.
    *  That's because the choice of `$this` decides which divisor to use
    *  when scaling back the oversized numerator that results from the
    *  initial multiplication. (It will always be the divisor of `$this`,
    *  which is both a source operand and also the destination for the result.)
    *
    *  The `multiply` function has the commutative property, but also
    *  requires the result's divisor to be explicitly stated as an additional
    *  argument.
    */
    public function multiply_by(FixedPoint64 $other, int $rounding_mode = RoundingMode::DEFAULT) : void
    {
        // ASSUMPTION: We wish the output to have the same divisor as $this, not $other (or anything else).
        $this->numerator = self::multiply_raw($this, $other, $this->denominator, $rounding_mode);
    }

    public static function multiply(FixedPoint64 $a, FixedPoint64 $b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : FixedPoint64
    {
        $numerator = self::multiply_raw($a, $b, $divisor, $rounding_mode);
        return new FixedPoint64($numerator, $divisor);
    }

    private static function multiply_raw(FixedPoint64 $fp64a, FixedPoint64 $fp64b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : int
    {
        $a = $fp64a->numerator;
        $b = $fp64b->numerator;

        $lo = (int)0;
        $hi = (int)0;
        Misc::multiply_64x64($a,$b,$hi,$lo);

        // TODO: Implement more rounding modes? (Also `banker_round` is pretty sus at the moment b/c sleep not enough.)
        if ( $rounding_mode == RoundingMode::HALF_TOWARDS_EVEN ) {
            $lo = self::banker_round($lo, $divisor);
        } else
        if ( $rounding_mode != RoundingMode::DEFAULT ) {
            throw new \Exception("Unimplemented rounding mode ".RoundingMode::to_string($rounding_mode));
        }

        $quotient_lo = (int)0;
        $quotient_hi = (int)0;
        $remainder = Int128::divide_128x64($hi,$lo,  $divisor,  $quotient_hi,$quotient_lo);

        if ( $rounding_mode != RoundingMode::DEFAULT ) {
            assert($remainder === 0); // Because rounding should have handled this already.
        }

        if ( 0 < $quotient_hi ) {
            $class_name = __CLASS__;
            $func_name = __FUNCTION__;
            $rounding_mode_str = RoundingMode::to_string($rounding_mode);
            throw new \ArithmeticError(sprintf(
                "Overflow in `$class_name::$func_name`. ".
                "Parameters: {fp64a->numerator:%016X, fp64a->denominator:%016X, fp64b->numerator:%016X, fp64b->denominator:%016x, divisor:%016X, rounding_mode:$rounding_mode_str}",
                $fp64a->numerator, $fp64a->denominator, $fp64b->numerator, $fp64b->denominator, $divisor
                ));
        }

        return $quotient_lo;
    }

    public static function unittest_multiply() : void
    {
        // Left as exercise for reader.
        echo("unittest FixedPoint64->multiply\n");
        echo(" done.\n");
    }

    public static function unittest() : void
    {
        // Put unittests for other methods here...
        self::unittest_banker_round();
        // ... or here ...
        self::unittest_multiply();
        // ... or here.
        echo("\n");
    }
}

?>


<?php
declare(strict_types=1);

use \Kickback\Math\Int64;
use \Kickback\Math\Int128;
use \Kickback\Math\Types\FixedPoint64;

Int64::unittest();
Int128::unittest();
FixedPoint64::unittest();
?>