API tools faq
paste
Advertisement
chadjoan

PHP FixedPoint64 rough sketch rev2

chadjoan
Jul 10th, 2024
132
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
PHP 43.68 KB | None | 0 0
raw download clone embed print report
  1. <?php
  2. declare(strict_types=1);
  3.  
  4. namespace Kickback\Math;
  5.  
  6. final class Int64
  7. {
  8.     private const MASK_SIGN = (int)(1 << 63);
  9.  
  10.     public static function add_unsigned_with_carry(int $a, int $b, int &$carry) : int
  11.     {
  12.         // Optimized version thanks to this page:
  13.         // https://grack.com/blog/2022/12/20/deriving-a-bit-twiddling-hack/
  14.         // (It has been adjusted for 64-bit arithmatic.)
  15.         $c = $a + $b;
  16.         $carry = ((($c ^ $a) & ($c ^ $b)) >> 63);
  17.         // Note that $carry will either be 0 or it will be 0xFFFF...FFFF.
  18.         // But we want it to be 0 or 1.
  19.         // So we can clear the sign extension bits by ANDing with 1:
  20.         $carry &= 1;
  21.         return $c;
  22.     }
  23.     /* TODO: cut?
  24.         // Since 0 is false and anything non-0 is considered "true", then
  25.         // this will return what we want, without having to worry about
  26.         // somehow reducing this to the set of {0,1}.
  27.         */
  28.  
  29.     public static function unittest_add_unsigned_with_carry() : void
  30.     {
  31.         echo("unittest Int128::unittest_add_unsigned_with_carry");
  32.  
  33.         // Shorthand to avoid having to write out long 64-bit hex values.
  34.         $x0000 = (int)0;
  35.         $x0001 = (int)1;
  36.         $x0002 = (int)2;
  37.         $x8000 = PHP_INT_MIN;
  38.         $x8001 = PHP_INT_MIN+1;
  39.         $xFFFF = (int)-1;
  40.         $xFFFE = (int)-2;
  41.         $x7FFF = PHP_INT_MAX;
  42.         $x7FFE = PHP_INT_MAX-1;
  43.  
  44.         $carry = (int)0;
  45.  
  46.         assert($x0000 === self::add_unsigned_with_carry($x0000,$x0000,$carry)); assert(0 === $carry);
  47.         assert($x0001 === self::add_unsigned_with_carry($x0000,$x0001,$carry)); assert(0 === $carry);
  48.         assert($x0001 === self::add_unsigned_with_carry($x0001,$x0000,$carry)); assert(0 === $carry);
  49.         assert($x0002 === self::add_unsigned_with_carry($x0001,$x0001,$carry)); assert(0 === $carry);
  50.  
  51.         assert($x0000 === self::add_unsigned_with_carry($x0000,$x0000,$carry)); assert(0 === $carry);
  52.         assert($xFFFF === self::add_unsigned_with_carry($x0000,$xFFFF,$carry)); assert(0 === $carry);
  53.         assert($xFFFF === self::add_unsigned_with_carry($xFFFF,$x0000,$carry)); assert(0 === $carry);
  54.         assert($xFFFE === self::add_unsigned_with_carry($xFFFF,$xFFFF,$carry)); assert(1 === $carry);
  55.  
  56.         assert($x0000 === self::add_unsigned_with_carry($x0001,$xFFFF,$carry)); assert(1 === $carry);
  57.         assert($x7FFF === self::add_unsigned_with_carry($xFFFF,$x8000,$carry)); assert(1 === $carry);
  58.  
  59.         assert($x0001 === self::add_unsigned_with_carry($x0002,$xFFFF,$carry)); assert(1 === $carry);
  60.         assert($x7FFE === self::add_unsigned_with_carry($xFFFE,$x8000,$carry)); assert(1 === $carry);
  61.  
  62.         assert($x8000 === self::add_unsigned_with_carry($x8001,$xFFFF,$carry)); assert(1 === $carry);
  63.         assert($x7FFE === self::add_unsigned_with_carry($x7FFF,$xFFFF,$carry)); assert(1 === $carry);
  64.  
  65.         echo(" done.\n");
  66.     }
  67.  
  68.     public static function ensure_ints_are_64bit() : void
  69.     {
  70.         $num_int_bits = PHP_INT_SIZE*8;
  71.         if ($num_int_bits < 64)
  72.         {
  73.             throw new SomethingSomethingException(
  74.                 "This system features ${num_int_bits}-bit integers, " .
  75.                 "but at least 64-bit wide integers are required for financial calculations.");
  76.         }
  77.     }
  78.  
  79.     public static function unittest() : void
  80.     {
  81.         ensure_ints_are_64bit();
  82.         unittest_add_unsigned_with_carry();
  83.         echo("\n");
  84.     }
  85.  
  86.     // Prevent instantiation/construction of the (static/constant) class.
  87.     /** @return never */
  88.     private function __construct() {
  89.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  90.     }
  91. }
  92. ?>
  93.  
  94. <?php
  95. declare(strict_types=1);
  96. namespace Kickback\Math;
  97.  
  98. use \Kickback\Math\Int64;
  99.  
  100. final class Int128
  101. {
  102.     private const MASK_LO = (int)((1 << 32) - 1);
  103.     private const MASK_SIGN = (int)(1 << 63);
  104.  
  105.     // Internal version of `cmp`.
  106.     // This version does not clamp it's outputs to the set {-1, 0, 1}.
  107.     private static function cmp_unsaturated(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
  108.     {
  109.         $res = $a_hi - $b_hi;
  110.         if ( $res ) {
  111.             return $res;
  112.         }
  113.  
  114.         return ($a_lo - $b_lo);
  115.     }
  116.  
  117.     /*
  118.     *  Returns -1, 0, or 1 based on a comparison of the given 128-bit integers.
  119.     *
  120.     *  Defined such that:
  121.     *  * a < b implies cmp(a,b) < 0
  122.     *  * a > b implies cmp(a,b) > 0
  123.     *  * a == b implies cmp(a,b) == 0.
  124.     */
  125.     public static function cmp(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
  126.     {
  127.         $res = self::cmp_unsaturated($a_hi,$a_lo,  $b_hi,$b_lo);
  128.         if ( $res === 0x8000_0000_0000_0000 ) {
  129.             return -1;
  130.         }
  131.         $res--;      //    >0    =0    <0
  132.         $res >>= 62; // {0b000,0b111,0b110} = {0,-1,-2}
  133.         $res++;      // {0b001,0b000,0b111} = {1, 0,-1}
  134.         return $res;
  135. /* TODO: cut
  136.         $res--;      // {0b000,0b111,0b110,0b101}
  137.         $res &;      // {0b000,0b000,0b110,0b100}
  138.         want         // {0b011,0b001,0b000,0b011}   {0b11,0b01,0b00,0b11}
  139.  
  140.         $res += 3    // {0b100,0b011,0b110,0b101}
  141.         $res >>= 1;  // {0b110,0b001,0b111,0b110}
  142.         $res += 1;   // {0b111,0b010,0b000,0b111}
  143.         $res >>= 1;  // {0b111,0b001,0b000,0b111}
  144.  
  145.         $res >>= 62; // {0b001,0b000,0b111,0b110} ^ {0b00,0b01,0b10,0b11} | {0b11,0b10,0b01,0b10}
  146.         $res += 1    // {0b010,0b001,0b000,0b111} ^ {0b11,0b00,0b01,0b10} | {0b10,0b11,0b10,0b01}
  147.         $res >>= 1;  // {0b001,0b000,0b001,0b111} ^ {0b11,0b10,0b00,0b11}
  148.         $res += 2    // {0b011,0b010,0b011,0b001} ^ {0b10,0b11,0b00,0b01} | {0b01,0b10,0b11,0g10}
  149.         $res += 3    // {0b100,0b011,0b110,0b101} ^ {0b01,0b10,0b11,0b00} | {0}
  150.  
  151.         $res >>= 62; // {0b001,0b000,0b111,0b110} ^ {0b00,0b01,0b10,0b11} | {0b11,0b10,0b01,0b10}
  152.         $res += 1    // {0b010,0b001,0b000,0b111} ^ {0b11,0b00,0b01,0b10} | {0b10,0b11,0b10,0b01}
  153.         $res += 2    // {0b011,0b010,0b001,0b000} ^ {0b10,0b11,0b00,0b01} | {0b01,0b10,0b11,0g10}
  154.         $res += 3    // {0b100,0b011,0b010,0b001} ^ {0b01,0b10,0b11,0b00} | {0}
  155.  
  156.         $res >>= 62; // {0b01,0b00,0b11,0b10} - {0b11,0b00,0b01,0b10} | {0b11,0b10,0b01,0b10}
  157.         $res += 1    // {0b10,0b01,0b00,0b11} - {0b10,0b11,0b00,0b01} | {0b10,0b11,0b10,0b01}
  158.         $res += 2    // {0b11,0b10,0b01,0b00} - {0b01,0b10,0b11,0b00} | {0b01,0b10,0b11,0g10}
  159.         $res += 3    // {0b00,0b11,0b10,0b01} - {0b00,0b01,0b10,0b11} | {0}
  160.  
  161.         $res >>= 62; // {0b01,0b00,0b11,0b10} - {0b11,0b00,0b01,0b10} | {0b11,0b00,0b11,0b10}
  162.         $res += 1    // {0b10,0b01,0b00,0b11} - {0b10,0b11,0b00,0b01} | {0b10,0b11,0b00,0b11}
  163.         $res += 2    // {0b11,0b10,0b01,0b00} - {0b01,0b10,0b11,0b00} | {0b11,0b10,0b11,0b00}
  164.         $res += 3    // {0b00,0b11,0b10,0b01} - {0b00,0b01,0b10,0b11} | {0b00,0b11,0b10,0b11}
  165.  
  166.         if ( !$res ) {
  167.             return $res;
  168.         }
  169.  
  170.         $res >>= 63;
  171.         $res <<= 2;
  172.         return 0x
  173.         */
  174.  
  175.         /*
  176.         if ( $a_hi < $b_hi ) {
  177.             return -1;
  178.         } else
  179.         if ( $a_hi > $b_hi ) {
  180.             return 1;
  181.         }
  182.         assert($a_hi === $b_hi);
  183.  
  184.         if ( $a_lo < $b_lo ) {
  185.             return -1;
  186.         } else
  187.         if ( $a_lo > $b_lo ) {
  188.             return 1;
  189.         }
  190.         assert($a_lo === $b_lo);
  191.  
  192.         return 0;
  193.         */
  194.     }
  195.  
  196.     public static function unittest_cmp()
  197.     {
  198.         echo("unittest Int128::cmp");
  199.  
  200.         assert(self::cmp( 0, 0,   0, 0) === 0);
  201.         assert(self::cmp( 0, 0,   0, 1)  <  0);
  202.         assert(self::cmp( 0, 1,   0, 0)  >  0);
  203.         assert(self::cmp( 0, 1,   0, 1) === 0);
  204.         assert(self::cmp(-1, 0,  -1, 0) === 0);
  205.         assert(self::cmp(-1, 0,   0, 0)  <  0);
  206.         assert(self::cmp(-1, 0,   1, 0)  <  0);
  207.         assert(self::cmp( 0, 0,  -1, 0)  >  0);
  208.         assert(self::cmp( 0, 0,   0, 0) === 0);
  209.         assert(self::cmp( 0, 0,   1, 0)  <  0);
  210.         assert(self::cmp( 1, 0,  -1, 0)  >  0);
  211.         assert(self::cmp( 1, 0,   0, 0)  >  0);
  212.         assert(self::cmp( 1, 0,   1, 0) === 0);
  213.  
  214.         // Ensure that it's returning values in the set {-1,0,1}, and not in the set {-N,0,N} or somesuch.
  215.         assert(self::cmp(  0, 0,   16, 0) === -1);
  216.         assert(self::cmp( 16, 0,    0, 0) ===  1);
  217.         assert(self::cmp(-16, 0,    0, 0) === -1);
  218.         assert(self::cmp(  0, 0,  -16, 0) ===  1);
  219.  
  220.         // Notably, we've carefully avoided negative values in the less-significant parts, so far.
  221.         // That's because the less-significant integers shall be treated as
  222.         // unsigned integers, but PHP only has _signed_ comparison.
  223.         // So we better check for mistakes of that kind!
  224.         $x0000 = (int)0;
  225.         $xFFFF = (int)-1;
  226.         assert(self::cmp($x0000,$x0000,  $x0000,$x0000) === 0); // 0 === 0
  227.         assert(self::cmp($x0000,$x0000,  $x0000,$xFFFF)  <  0); // 0 < (2^32 - 1)
  228.         assert(self::cmp($x0000,$x0000,  $xFFFF,$x0000)  >  0); // 0 > -(2^32)
  229.         assert(self::cmp($x0000,$x0000,  $xFFFF,$xFFFF)  >  0); // 0 > -1
  230.         assert(self::cmp($x0000,$xFFFF,  $x0000,$x0000)  >  0); // (2^32 - 1) > 0
  231.         assert(self::cmp($x0000,$xFFFF,  $x0000,$xFFFF) === 0); // (2^32 - 1) === (2^32 - 1)
  232.         assert(self::cmp($x0000,$xFFFF,  $xFFFF,$x0000)  >  0); // (2^32 - 1) > -(2^32)
  233.         assert(self::cmp($x0000,$xFFFF,  $xFFFF,$xFFFF)  >  0); // (2^32 - 1) > -1
  234.         assert(self::cmp($xFFFF,$x0000,  $x0000,$x0000)  <  0); // -(2^32) < 0
  235.         assert(self::cmp($xFFFF,$x0000,  $x0000,$xFFFF)  <  0); // -(2^32) < (2^32 - 1)
  236.         assert(self::cmp($xFFFF,$x0000,  $xFFFF,$x0000) === 0); // -(2^32) === -(2^32)
  237.         assert(self::cmp($xFFFF,$x0000,  $xFFFF,$xFFFF)  <  0); // -(2^32) < -1
  238.         assert(self::cmp($xFFFF,$xFFFF,  $x0000,$x0000)  <  0); // -1 < 0
  239.         assert(self::cmp($xFFFF,$xFFFF,  $x0000,$xFFFF)  <  0); // -1 < (2^32 - 1)
  240.         assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$x0000)  >  0); // -1 > -(2^32)
  241.         assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$xFFFF) === 0); // -1 === -1
  242.  
  243.         echo(" done.\n");
  244.     }
  245.  
  246.     /*
  247.     * Multiplies two 64-bit integers, resulting in a single 128-bit integer.
  248.     *
  249.     * Because PHP does not have a native 128-bit integer type (or the ability
  250.     * to define structs that can be returned from functions), the returned
  251.     * value is placed into two 64-bit integer reference-type parameters.
  252.     */
  253.     public static function multiply_64x64(int $a, int $b, int &$out_hi, int &$out_lo) : void
  254.     {
  255.         // Reorganize the two 64-bit integers into four 32-bit integers.
  256.         // This helps because now the 32-bit integers will have at least
  257.         // 32 bits of "headroom" for doing multiplications.
  258.         $a_lo = $a & self::MASK_LO;
  259.         $a_hi = $a >> 32;
  260.         $b_lo = $b & self::MASK_LO;
  261.         $b_hi = $b >> 32;
  262.         self::multiply_NxN(32, $a_hi,$a_lo,  $b_hi,$b_lo,  $out_hi,$out_lo);
  263.     }
  264.  
  265.     /*
  266.     * It is probably better to use `multiply_64x64` in most cases.
  267.     *
  268.     * This is the underlying implementation of `multiply_64x64`, and having
  269.     * it separated out as a different function is useful for testing purposes.
  270.     */
  271.     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
  272.     {
  273.         assert($bit_width <= 32);
  274.         $mask_lo = ((1 << $digit_bit_width) - 1);
  275.         $max_digit = $mask_lo;
  276.  
  277.         // Multiplies. We need 4 of them for a school-style long multiply.
  278.         // There are faster methods, but this will do for now.
  279.         //
  280.         // Graphically:
  281.         //                a_hi  a_lo
  282.         //             x  b_hi  b_lo
  283.         //             -------------
  284.         //                x_hi  x_lo
  285.         //          y_hi  y_lo
  286.         //          z_hi  z_lo
  287.         // +  w_hi  w_lo
  288.         // -----------------------
  289.         //    ????  ????  ????  ????
  290.         //
  291.         // Note that these multiplications, or at least the ones involving
  292.         // the $*_lo operands, must be _unsigned_ multiplications.
  293.         // However, because these are 32-bit values stored in 64-bit variables,
  294.         // the sign bit will _never_ be set, so it won't matter if we
  295.         // use a 64-bit signed multiplication on them.
  296.         //
  297.         $x = $a_lo * $b_lo;
  298.         $y = $a_lo * $b_hi;
  299.         $z = $a_hi * $b_lo;
  300.         $w = $a_hi * $b_hi;
  301.  
  302.         // To make the logic more clear, we will also define the
  303.         // variables corresponding to the {x,y,z,w}_{hi,lo} placeholders.
  304.         $x_lo = $x & $mask_lo;
  305.         $y_lo = $y & $mask_lo;
  306.         $z_lo = $z & $mask_lo;
  307.         $w_lo = $w & $mask_lo;
  308.  
  309.         $x_hi = $x >> $digit_bit_width;
  310.         $y_hi = $y >> $digit_bit_width;
  311.         $z_hi = $z >> $digit_bit_width;
  312.         $w_hi = $w >> $digit_bit_width;
  313.  
  314.         // PHP doesn't provide an unsigned right-shift, so we must clear
  315.         // any sign-extended bits in things that were right-shifted.
  316.         $x_hi = $x & $mask_lo;
  317.         $y_hi = $y & $mask_lo;
  318.         $z_hi = $z & $mask_lo;
  319.         $w_hi = $w & $mask_lo;
  320.  
  321.         // Calculate the bottom row of 32-bit integers.
  322.         //
  323.         // Note that some of these additions might "overflow", but
  324.         // because we are storing our 32-bit integers in 64-bit variables,
  325.         // the carry values will be captured.
  326.         //
  327.         // These will get shuffled into the output integers
  328.         // once the accumulation is done.
  329.         //
  330.         // Graphically:
  331.         //                a_hi  a_lo
  332.         //             x  b_hi  b_lo
  333.         //             -------------
  334.         //                x_hi  x_lo
  335.         //          y_hi  y_lo
  336.         //          z_hi  z_lo
  337.         // +  w_hi  w_lo
  338.         // -----------------------
  339.         //    out3  out2  out1  out0
  340.         //
  341.         $out0 = $x_lo;
  342.         $out1 = $x_hi + $y_lo + $z_lo;
  343.         $out2 = $y_hi + $z_hi + $w_lo;
  344.         $out3 = $w_hi;
  345.  
  346.         // Perform carry operations.
  347.         // (Beware: order-of-operations is important.)
  348.         if ( $out1 > $max_digit )
  349.         {
  350.             $out2 += ($out1 >> $digit_bit_width);
  351.             $out1 &= $mask_lo;
  352.         }
  353.  
  354.         if ( $out2 > $max_digit )
  355.         {
  356.             $out3 += ($out2 >> $digit_bit_width);
  357.             $out2 &= $mask_lo;
  358.         }
  359.         // Note that $out3 is involved in an addition, but won't generate
  360.         // a carry-out. It won't overflow, for math reasons.
  361.  
  362.         // Now we have four proper 32-bit integers (two pairs) with no carry bits,
  363.         // so we can concatenate the pairs to get two 64-bit integers and have our 128-bit result.
  364.         $lo = $out0 + ($out1 << $digit_bit_width);
  365.         $hi = $out2 + ($out3 << $digit_bit_width);
  366.  
  367.         // Return our results from the function.
  368.         $out_lo = $lo;
  369.         $out_hi = $hi;
  370.     }
  371.  
  372.     // This is not part of the public API because there is no point.
  373.     // It _should_ do exactly what the `*` operator does.
  374.     // It's just designed to use the `multiply_64x64` function, so that
  375.     // we can see if it gives results identical to `*`, at least whenever
  376.     // the two are given operands that don't overflow.
  377.     private static function multiply_NxN_lo(int $bit_width, int $a, int $b) : int
  378.     {
  379.         $out_lo = (int)0;
  380.         $out_hi = (int)0;
  381.         self::multiply_NxN($bit_width,  $a,  $b,  $out_hi,$out_lo);
  382.         return $out_lo;
  383.     }
  384.     /* TODO: cut?
  385.     // Uses the `multiply_32H32Lx32H32L` function to emulate what
  386.     // `multiply_64x64` would do if it had lower bit-width integers available
  387.     // (ex: `multiply_32x32`).
  388.     //
  389.     // Because this is used for testing purposes, and because we've asserted
  390.     // that we'll always use PHP on 64-bit systems and `int` will always have
  391.     // a width of at least 64 bits, then the result of `multiply_32x32` and
  392.     // smaller will always fit in a 64-bit return value. So we don't need
  393.     // to return the output in reference parameters for this one, which is
  394.     // very helpful for testing.
  395.     private static function multiply_NxN_emul(int $bit_width, int $a, int $b) : int
  396.     {
  397.         assert($bit_width <= 32);
  398.         $mask_lo = ((1 << $bit_width) - 1);
  399.         $a_lo = $a & $mask_lo;
  400.         $a_hi = $a >> $bit_width;
  401.         $b_lo = $b & $mask_lo;
  402.         $b_hi = $b >> $bit_width;
  403.         $out_lo = (int)0;
  404.         $out_hi = (int)0;
  405.         self::multiply_32H32Lx32H32L($a_hi,$a_lo,  $b_hi,$b_lo,  $out_hi,$out_lo);
  406.  
  407.         $int_max = $mask_lo;
  408.         assert($out_hi <= $int_max);
  409.         assert($out_lo <= $int_max);
  410.  
  411.  
  412.     }
  413. */
  414.  
  415.     public static function unittest_multiply_64x64() : void
  416.     {
  417.         echo("unittest Int128::unittest_multiply_64x64");
  418.  
  419.         assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x00));
  420.         assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x01));
  421.         assert(0x0000 === self::multiply_NxN_lo(8, 0x01, 0x00));
  422.         assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0xFF));
  423.         assert(0x0000 === self::multiply_NxN_lo(8, 0xFF, 0x00));
  424.         assert(0x0001 === self::multiply_NxN_lo(8, 0x01, 0x01));
  425.         assert(0x000F === self::multiply_NxN_lo(8, 0x01, 0x0F));
  426.         assert(0x000F === self::multiply_NxN_lo(8, 0x0F, 0x01));
  427.         assert(0x00E1 === self::multiply_NxN_lo(8, 0x0F, 0x0F));
  428.         assert(0x0010 === self::multiply_NxN_lo(8, 0x01, 0x10));
  429.         assert(0x0010 === self::multiply_NxN_lo(8, 0x10, 0x01));
  430.         assert(0x0100 === self::multiply_NxN_lo(8, 0x10, 0x10));
  431.         assert(0x4000 === self::multiply_NxN_lo(8, 0x80, 0x80));
  432.         assert(0x3F01 === self::multiply_NxN_lo(8, 0x7F, 0x7F));
  433.         assert(0xFC04 === self::multiply_NxN_lo(8, 0xFE, 0xFE));
  434.         assert(0xFD02 === self::multiply_NxN_lo(8, 0xFE, 0xFF));
  435.         assert(0xFD02 === self::multiply_NxN_lo(8, 0xFF, 0xFE));
  436.         assert(0xFE01 === self::multiply_NxN_lo(8, 0xFF, 0xFF));
  437.         assert(0x7E81 === self::multiply_NxN_lo(8, 0x7F, 0xFF));
  438.         assert(0x7F80 === self::multiply_NxN_lo(8, 0x80, 0xFF));
  439.         assert(0x3F80 === self::multiply_NxN_lo(8, 0x80, 0x7F));
  440.  
  441.         for ( $i = (int)0; $i < 256; $i++ )
  442.         {
  443.             for ( $j = (int)0; $j < 256; $j++ )
  444.             {
  445.                 $expected = (int)($i*$j);
  446.                 $got = self::multiply_NxN_lo(8, $i, $j);
  447.                 assert($expected === $got, sprintf("Operands: i=%02x; j=%02x; Expected: %04x; Got: %04X", $i, $j, $expected, $got));
  448.             }
  449.         }
  450.  
  451.         assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0000));
  452.         assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0001));
  453.         assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0001, 0x0000));
  454.         assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0xFFFF));
  455.         assert(0x0000_0000 === self::multiply_NxN_lo(16, 0xFFFF, 0x0000));
  456.         assert(0x0000_0001 === self::multiply_NxN_lo(16, 0x0001, 0x0001));
  457.         assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0x0001, 0xFFFF));
  458.         assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0xFFFF, 0x0001));
  459.         assert(0xFFFE_0001 === self::multiply_NxN_lo(16, 0xFFFF, 0xFFFF));
  460.         assert(0XA518_60E1 === self::multiply_NxN_lo(16, 0xF00F, 0xB00F));
  461.  
  462.         assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0000));
  463.         assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0001));
  464.         assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0000));
  465.         assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0xFFFF_FFFF));
  466.         assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0000));
  467.         assert(0x0000_0000_0000_0001 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0001));
  468.         assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0x0000_0001, 0xFFFF_FFFF));
  469.         assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0001));
  470.         assert(0xFFFF_FFFE_0000_0001 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0xFFFF_FFFF));
  471.  
  472.         // Explicit test of 128-bit returning version, and of 64-bit inputs.
  473.         $a64 = (int)0xAceBe11e_CafeBabe;
  474.         $b64 = (int)0xF100fCa7_F1edFa57;
  475.         $out_hi = (int)0;
  476.         $out_lo = (int)0;
  477.         self::multiply_64x64($a64,  $b64,  $out_hi,$out_lo);
  478.         assert(0x8E9C_2945_7ED5_0292 === $out_lo);
  479.         assert(0xA2CA_B997_9FFE_C71C === $out_hi);
  480.  
  481.         echo(" done.\n");
  482.     }
  483.  
  484.     /**
  485.     * @param      int  $numerator_hi
  486.     * @param      int  $numerator_lo
  487.     * @param      int  $denominator
  488.     * @param-out  int  $quotient_hi
  489.     * @param-out  int  $quotient_lo
  490.     * @returns    int  The remainder.
  491.     */
  492.     public static function divide_128x64(int $numerator_hi, int $numerator_lo, int $denominator, int &$quotient_hi, int &$quotient_lo) : int
  493.     {
  494.         // TODO: Test what PHP does with integer divide-by-zero.
  495.         // In principle, PHP should throw a DivisionByZeroError whenever this
  496.         // happens. In that case, we can just proceed as normal, because
  497.         // the first time division is used in this function, it'll cause
  498.         // that error to be thrown.
  499.         // See: https://www.php.net/manual/en/class.divisionbyzeroerror.php
  500.         // (But I haven't tested to make _sure_ PHP does this. If it doesn't
  501.         // throw an exception, then we should. Propagating weird values
  502.         // through the calculation could result in hard-to-find bugs.)
  503.  
  504.         // Special-case the value denominator value of `1`, both to make this
  505.         // identity operation be fast, but also to ensure that we don't have
  506.         // to worry about any weird corner cases for denominators that are
  507.         // less than 2. (It can matter, because binary arithmatic.)
  508.         if ( $denominator === 1 ) {
  509.             $quotient_hi = $numerator_hi;
  510.             $quotient_lo = $numerator_lo;
  511.             return 0;
  512.         }
  513.  
  514.         // Give these predictable values.
  515.         // This will be helpful in debugging if an error happens:
  516.         // you'll always know what these started as!
  517.         $quotient_hi = (int)0;
  518.         $quotient_lo = (int)0;
  519.  
  520.         // Store the sign bit that the result will have.
  521.         $sign = ($numerator_hi ^ $denominator) & self::MASK_SIGN;
  522.  
  523.         // Sign extraction, so it doesn't mess with our division.
  524.         $sign = $hi & self::MASK_SIGN;
  525.         $hi &= ~self::MASK_SIGN;
  526.  
  527.         // Extract remainders.
  528.         //
  529.         // Right now we are mostly interested in $hi_remainder.
  530.         // This will be necessary later for estimating the low portion of the result.
  531.         //
  532.         // The $lo_remainder will simply be returned by the function as-is,
  533.         // since this may help the caller with rounding logic.
  534.         $lo_remainder = $lo % $denominator;
  535.         $hi_remainder = $hi % $denominator;
  536.  
  537.         // Component-wise division.
  538.         //
  539.         // (We use `intdiv` because `/` might do a check to see if one or both
  540.         // of its operands are floating point numbers, because it is a floating
  541.         // point operation by default. `intdiv`, on the other hand, is explicitly
  542.         // intended for integer division, so it may perform/behave better for this use.)
  543.         $lo = intdiv($lo, $denominator);
  544.         $hi = intdiv($hi, $denominator);
  545.  
  546.         // Calculate the carry.
  547.         $min_carry = PHP_INT_MAX; // Basically ((2^64 / 2) - 1). We'd prefer (2^64), but can't represent it.
  548.         $min_carry = intdiv($min_carry, $denominator);
  549.         $min_carry <<= 1; // Undo the earlier divide-by-2 in the ((2^64 / 2) - 1) = PHP_INT_MAX.
  550.  
  551.         // Perform the carry.
  552.         $lo += ($hi_remainder * $min_carry);
  553.  
  554.         // For safety reasons, we'll also calculate a "max" carry.
  555.         $max_carry_init = (1 << 62);
  556.         $max_carry = intdiv($max_carry_init, $denominator);
  557.         if ( ($denominator * $max_carry) !== $max_carry_init ) {
  558.             // Always round up.
  559.             $max_carry++;
  560.         }
  561.         $max_carry <<= 2;
  562.         $max_carry += 3; // Always round up.
  563.  
  564.         // Perform the carry.
  565.         $max_lo += ($hi_remainder * $max_carry);
  566.  
  567.         // This loop takes our approximation and improves it until we have
  568.         // the correct quotient.
  569.         $test_lo = (int)0;
  570.         $test_hi = (int)0;
  571.         $test_carry = (int)0;
  572.         self::multiply_64x64($denominator,  $lo,  $test_hi,$test_lo);
  573.         $test_hi += ($hi * $denominator);
  574.  
  575.         while ( true )
  576.         {
  577.             $test_lo = Int64::add_unsigned_with_carry($test_lo, $denominator, $test_carry);
  578.             $test_hi += $test_carry;
  579.             // if ( test > numerator )
  580.             if ( self::cmp($test_hi, $test_lo, $numerator_hi, $numerator_lo) > 0 ) {
  581.                 break;
  582.             }
  583.  
  584.             // The additional increment of $denominator in the $test variable
  585.             // corresponds to incrementing the $lo value by 1.
  586.             $lo++;
  587.  
  588.             // This handles the aforementioned "safety reasons".
  589.             // It prevents us from infinite-looping, or from trying to
  590.             // iterate most of a 64-bit integer's possible values
  591.             // when it is already futile to get the correct answer.
  592.             // Ideally, this will never be executed.
  593.             // If the surrounding function works as intended, then the loop
  594.             // will ALWAYS converge before this condition becomes true.
  595.             if ($lo > $max_lo) {
  596.                 $class_name = __CLASS__;
  597.                 $func_name = __FUNCTION__;
  598.                 throw new \Exception(sprintf(
  599.                     "Internal error: `$class_name::$func_name` did not converge when it should have. ".
  600.                     "Aborting to prevent excessive looping and CPU drain. ".
  601.                     "This means that the `$func_name` function is broken and may return incorrect results, ".
  602.                     "even when this error isn't thrown. Please report this error. ".
  603.                     "Parameters: {numerator_hi:%016X, numerator_lo:%016X, denominator:%016X}",
  604.                     $numerator_hi, $numerator_lo, $denominator
  605.                     ));
  606.             }
  607.         }
  608.  
  609.         // Pass results to the caller.
  610.         $quotient_hi = $hi;
  611.         $quotient_lo = $lo;
  612.         return $lo_remainder;
  613.  
  614.         /*
  615.         TODO: Comment is kinda sorta relevant but also outdated already.
  616.  
  617.         // This will be multiplied by the $hi_remainder to get the portion
  618.         // of the number that is being carried into the $lo component of the result.
  619.         //
  620.         // Note that this should be `(2^64/$denominator)` EXACTLY.
  621.         // Unfortunately, we can't be exact for two reasons:
  622.         // * Rational numbers (the exact result) can't be represented by finite integers.
  623.         // * The number 2^64 is out of the range of 64-bit integers (one more than the max!)
  624.         //
  625.         // Now we'll get clever, and our function will get slower.
  626.         //
  627.         // We notice that the multiplier doesn't _need_ to be exact.
  628.         // We can make a guess. And the guess can be wrong. As long as it's
  629.         // always wrong in a way that's a little bit low.
  630.         //
  631.         // If it's low, it will result in our end-result being too low.
  632.         // Then we can multiply our guess by the denominator.
  633.         // The result will be less than the numerator.
  634.         // (If we had exact numbers instead of integers, multiplying the quotient
  635.         // by the denominator would result in _exactly_ the numerator.)
  636.         //
  637.         // Then we can just add the $denominator over and over, checking
  638.         // it with a multiplication, until it is too high.
  639.         // Once it is too high, we return the last number that multiplied
  640.         // to a value
  641.         //
  642.         */
  643.     }
  644.  
  645.     public static function unittest_divide_128x64()
  646.     {
  647.         echo("unittest Int128::unittest_divide_128x64");
  648.         // TODO: Oh god this needs testing BADLY.
  649.         echo(" UNIMPLEMENTED! (Please fix!)\n");
  650.         // echo(" done.\n");
  651.     }
  652.  
  653.     public static function unittest()
  654.     {
  655.         unittest_cmp();
  656.         unittest_multiply_64x64();
  657.         unittest_divide_128x64();
  658.         echo("\n");
  659.     }
  660.  
  661.     // Prevent instantiation/construction of the (static/constant) class.
  662.     /** @return never */
  663.     private function __construct() {
  664.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  665.     }
  666. }
  667. ?>
  668.  
  669. <?php
  670. declare(strict_types=1);
  671.  
  672. namespace Kickback\Math;
  673.  
  674. // TODO: I forget if it's better to do this, or to use an `enum` type. Whatevs.
  675. final class RoundingMode
  676. {
  677.     // Rounding modes.
  678.     // (I think this is all of the possibilities? TODO: Verify.)
  679.     public const ALL_TOWARDS_NEG_INF  =  0; // Default; mimics x86 integer truncation. (I think? TODO: Double-check this.)
  680.     public const ALL_TOWARDS_POS_INF  =  1;
  681.     public const ALL_TOWARDS_ZERO     =  3;
  682.     public const ALL_AWAY_FROM_ZERO   =  4;
  683.     public const ALL_TOWARDS_EVEN     =  5;
  684.     public const ALL_TOWARDS_ODD      =  6;
  685.     public const HALF_TOWARDS_NEG_INF =  7;
  686.     public const HALF_TOWARDS_POS_INF =  8;
  687.     public const HALF_TOWARDS_ZERO    =  9; // Similar to PHP_ROUND_HALF_DOWN
  688.     public const HALF_AWAY_FROM_ZERO  = 10; // Similar to PHP_ROUND_HALF_UP
  689.     public const HALF_TOWARDS_EVEN    = 11; // Similar to PHP_ROUND_HALF_EVEN; bankers' rounding.
  690.     public const HALF_TOWARDS_ODD     = 12; // Similar to PHP_ROUND_HALF_ODD
  691.  
  692.     // This shall mimic x86 integer truncation.
  693.     public const DEFAULT = self::ALL_TOWARDS_NEG_INF;
  694.  
  695.     public static function rounding_mode_to_string(int $rounding_mode) : string
  696.     {
  697.         switch($rounding_mode)
  698.         {
  699.             case self::ALL_TOWARDS_NEG_INF : return "ALL_TOWARDS_NEG_INF";
  700.             case self::ALL_TOWARDS_POS_INF : return "ALL_TOWARDS_POS_INF";
  701.             case self::ALL_TOWARDS_ZERO    : return "ALL_TOWARDS_ZERO";
  702.             case self::ALL_AWAY_FROM_ZERO  : return "ALL_AWAY_FROM_ZERO";
  703.             case self::ALL_TOWARDS_EVEN    : return "ALL_TOWARDS_EVEN";
  704.             case self::ALL_TOWARDS_ODD     : return "ALL_TOWARDS_ODD";
  705.             case self::HALF_TOWARDS_NEG_INF: return "HALF_TOWARDS_NEG_INF";
  706.             case self::HALF_TOWARDS_POS_INF: return "HALF_TOWARDS_POS_INF";
  707.             case self::HALF_TOWARDS_ZERO   : return "HALF_TOWARDS_ZERO";
  708.             case self::HALF_AWAY_FROM_ZERO : return "HALF_AWAY_FROM_ZERO";
  709.             case self::HALF_TOWARDS_EVEN   : return "HALF_TOWARDS_EVEN";
  710.             case self::HALF_TOWARDS_ODD    : return "HALF_TOWARDS_ODD";
  711.         }
  712.  
  713.         return "Unknown rounding mode ($rounding_mode)";
  714.     }
  715.  
  716.     // Prevent instantiation/construction of the (static/constant) class.
  717.     /** @return never */
  718.     private function __construct() {
  719.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  720.     }
  721. }
  722. ?>
  723.  
  724. <?php
  725. declare(strict_types=1);
  726.  
  727. namespace Kickback\Math\Types;
  728.  
  729. use \Kickback\Math\Int128;
  730. use \Kickback\Math\RoundingMode;
  731.  
  732. final class FixedPoint64
  733. {
  734.     private int numerator;
  735.     private int divisor;
  736.  
  737.     private static function banker_round(int $numerator, int $increment) : int
  738.     {
  739.         assert($increment > 1);
  740.  
  741.         // Even if we don't use $even_increment in most cases,
  742.         // we should still assert that it fits in _all_ cases, just to give
  743.         // this function consistent behavior (e.g. whether it asserts or not
  744.         // should depend on `$increment` and not on `$numerator`.)
  745.         //$even_increment = ($increment<<1);
  746.         assert(($increment<<1) <= PHP_INT_MAX);
  747.  
  748.         if ( $numerator >= 0 ) {
  749.             return self::banker_round_unsigned_($numerator, $increment);
  750.         }
  751.         else {
  752.             return self::banker_round_signed_($numerator, $increment);
  753.         }
  754.     }
  755.  
  756.     private static function banker_round_signed_(int $numerator, int $increment) : never
  757.     {
  758.         // TODO!
  759.         assert(false, "Signed rounding is not yet implemented.");
  760.     }
  761.  
  762.     private static function banker_round_unsigned_(int $numerator, int $increment) : int
  763.     {
  764.         /* TODO: cut?
  765.         // We calculate the $round_down_amount (conventional rounding)
  766.         // from the $round_down_even_amount (bankers' rounding)
  767.         // so that we avoid having to modulus more than once.
  768.         $round_down_even_amount = $numerator % $even_increment;
  769.         $round_down_amount = $round_down_even_amount;
  770.         if ( $round_down_amount > $increment ) {
  771.             $round_down_amount -= $increment;
  772.         }
  773.         */
  774.  
  775.         // First, attempt conventional rounding (e.g. "round-to-nearest-increment").
  776.         $rounding_amount = (int)0;
  777.         if ( self::get_rounding_amount_($numerator, $increment, $rounding_amount) ) {
  778.             return $rounding_amount;
  779.         } else {
  780.             // If that didn't work, then fall back on the tie breaker.
  781.             return self::banker_round_tie_breaker_($numerator, $increment);
  782.         }
  783.     }
  784.  
  785.     private static function get_rounding_amount_(int $numerator, int $increment, &$rounding_amount) : bool
  786.     {
  787.         $round_down_amount = $numerator % $increment;
  788.         $round_up_amount = $increment - $round_down_amount;
  789.         if ( $round_down_amount < $round_up_amount ) {
  790.             $rounding_amount = -$round_down_amount;
  791.             return true;
  792.         } else
  793.         if ( $round_down_amount > $round_up_amount ) {
  794.             $rounding_amount = $round_up_amount;
  795.             return true;
  796.         } else {
  797.             // If that didn't work, then there's a problem, and it's probably this tie:
  798.             assert($round_down_amount === $round_up_amount);
  799.             $rounding_amount = 0; // TODO: Is this right?
  800.             return false;
  801.         }
  802.     }
  803.  
  804.     // This is set out into its own function because the tie breaker
  805.     // is unlikely to be executed very frequently. As such, it might
  806.     // be faster (code-execution-wise) to break this out into it's own function.
  807.     //
  808.     // Reasons why this might be the case:
  809.     // * Interpreter is less likely to need to parse/analyse this code when `banker_round` is called.
  810.     // * This is more cache-friendly: the CPU won't need to load the instructions for the tie breaker every time `banker_round` is called.
  811.     // * This makes `banker_round` smaller and therefore more inline-able.
  812.     //
  813.     // I am not sure how much any of those _actually_ matter, and it's not
  814.     // really worth testing at the moment, but this is _usually_ a good way
  815.     // to optimize code, and it seems to have few or no disadvantages.
  816.     private static function banker_round_tie_breaker_(int $numerator, int $increment) : int
  817.     {
  818.         // Now the bankers' rounding comes into play.
  819.         // We break the tie by rounding to the nearest _even_ increment.
  820.         $even_increment = $increment << 1;
  821.  
  822.         // This involves just doing the rounding again, but with $increment*2.
  823.         $rounding_amount = (int)0;
  824.         assert(self::get_rounding_amount_($numerator, $even_increment, $rounding_amount));
  825.         return $rounding_amount;
  826.     }
  827.     // TODO: The above bankers' rounding code is pretty sus. I was very sleepy while writing it. o.O
  828.  
  829.     private static function unittest_banker_round() : void
  830.     {
  831.         echo("unittest FixedPoint64::unittest_banker_round");
  832.         // TODO: This needs testing if you want to trust any financials done with it.
  833.         echo(" UNIMPLEMENTED! (Please fix!)\n");
  834.         // echo(" done.\n");
  835.     }
  836.  
  837.  
  838.  
  839.     /*
  840.     * Beware: This operation is NOT commutative like normal multiplication.
  841.     * That's because the choice of destination decides which divisor to use
  842.     * when scaling back the oversized numerator that results from the
  843.     * initial multiplication.
  844.     *
  845.     * The `multiply` and `multiply_into` functions have the commutative property.
  846.     *
  847.     * This method will not perform any explicit memory allocations
  848.     * unless an error has occurred.
  849.     */
  850.     public function multiply_by(FixedPoint64 $other, int $rounding_mode = RoundingMode::DEFAULT) : void
  851.     {
  852.         // ASSUMPTION: We wish the output to have the same divisor as $this, not $other (or anything else).
  853.         $this->numerator = self::multiply_raw($this, $other, $this->denominator, $rounding_mode);
  854.     }
  855.  
  856.     /*
  857.     * Allocates a new FixedPoint64 object, then uses the `multiply_into`
  858.     * function to multiply parameters `$a` and `$b` and store the result
  859.     * into the new FixedPoint64 object, which is then returned.
  860.     *
  861.     * The new FixedPoint64 object will have a $denominator field equal to
  862.     * the `$divisor` parameter.
  863.     *
  864.     * This function could be handy as a quick-and-dirty way to draft code
  865.     * that uses FixedPoint64 arithmetic. It has the disadvantage of doing
  866.     * memory allocation for an operation where it isn't always necessary.
  867.     * For production code, it is recommended to use `multiply_into` or
  868.     * `multiply_by` instead.
  869.     */
  870.     public static function multiply(FixedPoint64 $a, FixedPoint64 $b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : FixedPoint64
  871.     {
  872.         $result = new FixedPoint64(0, $divisor);
  873.         return self::multiply_into($a, $b, $result, $rounding_mode);
  874.     }
  875.  
  876.     /*
  877.     * Multiplies `$fp64a` and `$fp64b` together, then stores the result
  878.     * into the `$destination` object.
  879.     *
  880.     * As a convenience, the `$destination` object is returned.
  881.     *
  882.     * When rounding or truncating the fractional multiplication results,
  883.     * the `$destination` parameter's `$denominator` field is used
  884.     * as a divisor.
  885.     *
  886.     * This operation has commutative property over `$fp64a` and `$fp64b`.
  887.     *
  888.     * This function will not perform any explicit memory allocations
  889.     * unless an error has occurred. (In some cases it is possible for
  890.     * the caller to reuse existing FixedPoint64 objects for the `$destination`
  891.     * parameter, and thus avoid memory allocations entirely.)
  892.     */
  893.     public static function multiply_into(FixedPoint64 $fp64a, FixedPoint64 $fp64b, FixedPoint64 $destination, int $rounding_mode = RoundingMode::DEFAULT) : FixedPoint64
  894.     {
  895.         assert(isset($destination));
  896.         $tmp = FixedPoint64::multiply_raw($fp64a, $fp64b, $destination->denominator, $rounding_mode);
  897.         $destination->numerator = $tmp;
  898.         return $destination;
  899.     }
  900.  
  901.     private static function multiply_raw(FixedPoint64 $fp64a, FixedPoint64 $fp64b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : int
  902.     {
  903.         // TODO: Verify that divisors are being chosen correctly.
  904.         // (Actually the below seems to be pretty well argued, now that it's
  905.         // all been written down. Still, we should make sure it all gets tested.)
  906.         // TODO: Make sure unittests cover all of the branches in this function.
  907.  
  908.         // ---=== Design Notes ===---
  909.         // I am currently reasoning like so:
  910.         //
  911.         // Suppose we multiply two numbers, `a` and `b`.
  912.         // `a` has a denominator of 2^16, giving it 16 bits of precision.
  913.         // `b` has a denominator of 2^8, giving it 8 bits of precision.
  914.         //
  915.         // The product `a * b` will then have 24 bits of precision, with a denominator of 2^24.
  916.         //
  917.         // If we want to store the product into `a'`, which has the same precision as `a`,
  918.         // then we will need to divide the product by 2^8 (the denominator of `b`)
  919.         // to acquire a value with 16 bits of precision (because `16 = 24 - 8`).
  920.         //
  921.         // If we want to store the product into `b'`, which has the same precision as `b`,
  922.         // then we will need to divide the product by 2^16 (the denominator of `a`)
  923.         // to acquire a value with 8 bits of precision (because `8 = 24 - 16`).
  924.         //
  925.         // If we want to store the product into `c`,  which has N bits of precision,
  926.         // then we will need to divide by the number of bits required to reach
  927.         // that: `24 - x = N`, so `x = 24 - N`.
  928.         // We divide by `2^(24 - N)` to reach the precision for `c`.
  929.         //
  930.         // But wait, what if `N>24`? Well, that's a good thing to notice,
  931.         // because the destination may actually have more precision than
  932.         // is given by the product of the multiplicands' denominators.
  933.         // In that case, we will need to multiply instead of divide!
  934.         //
  935.         // Now, what if `a` has `P` bits of precision, and `b` has `Q` bits of precision?
  936.         // We will need to adjust our formula slightly:
  937.         // `x = (P+Q) - N`.
  938.         //
  939.         // Now, what if `a` has an arbitrary denominator `A` and `b`, likewise has `B`?
  940.         // The full (wide) product will have a denominator of `AB = A * B`.
  941.         //
  942.         // To store into `a'`, we'll need to find `A` in terms of `AB` and `B`.
  943.         // We write `A` like so: `A = (A * B) / B`.
  944.         // Thus, `A = AB / B` (by substitution.).
  945.         //
  946.         // To store into `b'`, we'll need to find `B` in terms of `AB` and `A`.
  947.         // We write `B` like so: `B = (A * B) / A`.
  948.         // Thus, `B = AB / A` (by substitution.).
  949.         //
  950.         // To store into `c`, we'll write `C` in terms of `AB` and `X`, then find `X`.
  951.         // `C = AB/X`
  952.         // `C*X = AB`
  953.         // `X = AB/C`
  954.         // In the code, we will already know `C` because it is the denominator of `c`.
  955.         // We will need to find `X` so that we can call `divide_128x64` with it.
  956.         // Now we know how to do that. (I hope.)
  957.         //
  958.         // Another algebraic expression for how `c`'s numerator (`cn`) would be calculated:
  959.         //   cn/cd = (an/ad) * (bn/bd)
  960.         //   cn = (an/ad) * (bn/bd) * cd
  961.         //   cn = ((an*bn) / (ad*bd)) * cd
  962.         //   cn = (an * bn * cd) / (ad * bd)
  963.         //   cn = (an * bn) * (cd / (ad * bd)) <- if( cd < (ad*bd) ) { do this? }
  964.         //   cn = (an * bn) / ((ad * bd) / cd) <- if( (ad*bd) < cd ) { do this? }
  965.  
  966.         assert($divisor > 0);
  967.         assert($divisor <= PHP_INT_MAX);
  968.  
  969.         $a = $fp64a->numerator;
  970.         $b = $fp64b->numerator;
  971.         $a_div = $fp64a->denominator;
  972.         $b_div = $fp64b->denominator;
  973.  
  974.         // ---=== Multiplication ===---
  975.         // Multiply 64bit * 64bit to yield 128bit value.
  976.         $lo = (int)0;
  977.         $hi = (int)0;
  978.         Int128::multiply_64x64($a,  $b,  $hi,$lo);
  979.  
  980.         // Do the same for the denominators/divisor.
  981.         $div_lo = (int)0;
  982.         $div_hi = (int)0;
  983.         Int128::multiply_64x64($a_div,  $b_div,  $div_hi,$div_lo);
  984.  
  985.         // ---=== Division: Denominator ===---
  986.         // We need to figure out what number to divide or new numerator by,
  987.         // so that it ends up with the precision described by `$divisor`.
  988.         // This will be `($a->denominator * $b->denominator) / $divisor`.
  989.         $scale_direction = (int)0;
  990.         $div_quotient_lo = (int)0;
  991.         $div_quotient_hi = (int)0;
  992.         $div_remainder   = (int)0;
  993.         if ( $divisor === 0 || $divisor === 1 ) {
  994.             assert($scale_direction === 0);
  995.             $div_quotient_lo = $div_lo;
  996.             $div_quotient_hi = $div_hi;
  997.             assert($div_remainder === 0);
  998.         } else
  999.         if ( 0 === $div_hi && $divisor > $div_lo ){
  1000.             $scale_direction = 1;
  1001.             // TODO: Division here is not guaranteed to be twos-complement compatible.
  1002.             $div_quotient_lo = intdiv($divisor, $div_lo);
  1003.             assert($div_quotient_hi === 0);
  1004.             $div_remainder = $divisor % $div_lo;
  1005.             // TODO: What to do with $div_remainder?
  1006.         }
  1007.         else {
  1008.             $scale_direction = -1;
  1009.             $div_remainder = Int128::divide_128x64($div_hi,$div_lo,  $divisor,  $div_quotient_hi,$div_quotient_lo);
  1010.  
  1011.             // TODO: This limitation isn't strictly necessary;
  1012.             // this is something that CAN happen with valid inputs,
  1013.             // and there is going to be a way to handle it.
  1014.             // It just seems kinda unlikely, and I don't have the time to write it right now.
  1015.             if ( $div_quotient_hi !== 0 ) {
  1016.                 // TODO: Better error message; put parameters and details about the denominators.
  1017.                 $class_name = __CLASS__;
  1018.                 $func_name = __FUNCTION__;
  1019.                 throw new \Exception(
  1020.                     "$class_name::$func_name: Unimplemented combination of input parameters; ".
  1021.                     "The product of the inputs' denominators divided by the destination's denominator ".
  1022.                     "must be a number that fits within PHP_INT_MAX.");
  1023.             }
  1024.             // TODO: What to do with $div_remainder?
  1025.         }
  1026.         // TODO: $div_remainder is a smell suggesting that we haven't done enough math or something.
  1027.         // In particular, this codepath is making it clear that handling arbitrary
  1028.         // input denominators leads to situations where there is no common factor
  1029.         // to divide by.
  1030.         //
  1031.         // Like, if `A = 2^16` and `B = 2^8`, then either `C = 2^16` or `C = 2^8`
  1032.         // will both work fine because `AB = 2^24` which has even factors of both 2^16 and 2^8.
  1033.         // `C` could be any power-of-2 in that situation, though powers greater than 23
  1034.         // will either cause us to do zero scaling or scale up (multiply the result and zero-fill the LSBs).
  1035.         //
  1036.         // However, if `A = 3` and `B = 5`, then `C` has to be some multiple
  1037.         // of either 3 or 5 in order for the scaling to happen cleanly.
  1038.         // If we plug in `C = 2`, then we get `X = 15/2`, which is clearly not an integer.
  1039.         // (It might be possible to get around this by using those remainders
  1040.         // in the later scaling expressions, but it is clear how in my current sleepy state.)
  1041.  
  1042.         // ---=== Rounding ===---
  1043.         // Round the 128bit value, modulo the `$divisor` parameter.
  1044.         // (We really only need to do this for the lower 64 bits, because $divisor can't be bigger than that.)
  1045.         // TODO: Implement more rounding modes? (Also `banker_round` is pretty sus at the moment b/c sleep not enough.)
  1046.         if ( $scale_direction < 0 )
  1047.         {
  1048.             if ( $rounding_mode == RoundingMode::HALF_TOWARDS_EVEN ) {
  1049.                 $lo = self::banker_round($lo, $divisor);
  1050.             } else
  1051.             if ( $rounding_mode != RoundingMode::DEFAULT ) {
  1052.                 throw new \Exception("Unimplemented rounding mode ".RoundingMode::to_string($rounding_mode));
  1053.             }
  1054.         }
  1055.  
  1056.         // ---=== Division (or multiplication): Numerator ===---
  1057.         // Shrink the 128bit value until it is at the desired precision according to `$divisor`.
  1058.         // This will ready it for overflow checking.
  1059.         $quotient_lo = (int)0;
  1060.         $quotient_hi = (int)0;
  1061.         $remainder   = (int)0;
  1062.         if ( $scale_direction === 0 ) {
  1063.             $quotient_lo = $lo;
  1064.             $quotient_hi = $hi;
  1065.             assert($remainder === 0);
  1066.         } else
  1067.         if ( $scale_direction > 0 ) {
  1068.             // Least likely case of all of them, but also the most complicated.
  1069.             $tmp1_out_lo = (int)0;
  1070.             $tmp1_out_hi = (int)0;
  1071.             $tmp2_out_lo = (int)0;
  1072.             $tmp2_out_hi = (int)0;
  1073.             Int128::multiply_64x64($lo,  $div_quotient_lo,  $tmp1_out_hi,$tmp1_out_lo);
  1074.             Int128::multiply_64x64($hi,  $div_quotient_lo,  $tmp2_out_hi,$tmp2_out_lo);
  1075.             $quotient_lo = $tmp1_out_lo;
  1076.             $quotient_hi = $tmp1_out_hi + $tmp2_out_lo;
  1077.             assert($tmp2_out_hi === 0);
  1078.         } else
  1079.         if ( $scale_direction < 0 ) {
  1080.             $remainder = Int128::divide_128x64($hi,$lo,  $divisor,  $quotient_hi,$quotient_lo);
  1081.         }
  1082.  
  1083.         // ---=== Overflow/Error Handling ===---
  1084.         // Now we check for overflow (and maybe do some validation on rounding logic).
  1085.         // If there is no overflow, then we can safely discard the high part of the quotient.
  1086.         if ( $rounding_mode != RoundingMode::DEFAULT ) {
  1087.             assert($remainder === 0); // Because rounding should have handled this already.
  1088.         }
  1089.  
  1090.         if ( 0 !== $quotient_hi ) {
  1091.             $class_name = __CLASS__;
  1092.             $func_name = __FUNCTION__;
  1093.             $rounding_mode_str = RoundingMode::to_string($rounding_mode);
  1094.             throw new \ArithmeticError(sprintf(
  1095.                 "Overflow in `$class_name::$func_name`. ".
  1096.                 "Parameters: {fp64a->numerator:%016X, fp64a->denominator:%016X, fp64b->numerator:%016X, fp64b->denominator:%016x, divisor:%016X, rounding_mode:$rounding_mode_str}",
  1097.                 $fp64a->numerator, $fp64a->denominator, $fp64b->numerator, $fp64b->denominator, $divisor
  1098.                 ));
  1099.         }
  1100.  
  1101.         // ---=== Return ===---
  1102.         return $quotient_lo;
  1103.     }
  1104.  
  1105.     public static function unittest_multiply() : void
  1106.     {
  1107.         // Left as exercise for reader.
  1108.         echo("unittest FixedPoint64->multiply\n");
  1109.         echo(" done.\n");
  1110.     }
  1111.  
  1112.     public static function unittest() : void
  1113.     {
  1114.         // Put unittests for other methods here...
  1115.         self::unittest_banker_round();
  1116.         // ... or here ...
  1117.         self::unittest_multiply();
  1118.         // ... or here.
  1119.         echo("\n");
  1120.     }
  1121. }
  1122.  
  1123. ?>
  1124.  
  1125.  
  1126. <?php
  1127. declare(strict_types=1);
  1128.  
  1129. // Unittest driver to be run from command-line interface.
  1130. use \Kickback\Math\Int64;
  1131. use \Kickback\Math\Int128;
  1132. use \Kickback\Math\Types\FixedPoint64;
  1133.  
  1134. public function unittest_all() : void
  1135. {
  1136.     Int64::unittest();
  1137.     Int128::unittest();
  1138.     FixedPoint64::unittest();
  1139. }
  1140.  
  1141. unittest_all();
  1142. ?>
  1143.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!