API tools faq
paste
Advertisement
chadjoan

PHP FixedPoint64 rough sketch rev3

chadjoan
Jul 13th, 2024
85
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
PHP 77.76 KB | None | 0 0
raw download clone embed print report
  1. <?php
  2. // Notable future direction for the math stuff:
  3. // Int64 turned into IntNN and began supporting both 32-bit and 64-bit operations.
  4. // Int128 doesn't really have an analogous transform, but it should be possible
  5. // to create operations for that too...
  6. // But the thing that REALLY needs to happen is the creation of 3 bins:
  7. // functions that operate on (algebraically) atomic values, ex: cmp($a,$b)
  8. // functions that operate on values represented in two parts, ex: cmp($ah,$al,  $bh,$bl)
  9. // functions that operate on values represented in _four_ parts, ex: cmp($a3,$a2,$a1,$a0,  $b3,$b2,$b1,$b0)
  10. //
  11. // All of those will be "low level", in a sense:
  12. // It's hard to know dynamically which one to call, because the argument counts change.
  13. //
  14. // So there will need to be abstractions that represent numbers with the correct
  15. // magnitude of parts, and then select the appropriate functions to operate
  16. // on those composites. Such an abstraction might be the FixedPoint64 class.
  17. // (I'm starting to be tempted to call it Rational64, because that's more like
  18. // what it _really_ is, though FixedPoint64 still might reflect its _usage_ better,
  19. // hmmmmmmm.)
  20. // The FixedPoint64 class could represent its internals with as many or as few
  21. // integers as it wants, and could switch between operations based on the
  22. // host system's capabilities.
  23. //
  24. // I think that ideally, the file would be templated and preprocessed by
  25. // PHP code into 2 different class files, with each class implementing
  26. // a different magnitude of integer representation (atomic vs 2-part)
  27. // and calling the appropriate functions as needed. (We'll still need the
  28. // 4-part low-level functions, because doing things like multiplying
  29. // two 2-part numbers will require a 4-part multiplication function!)
  30. //
  31. // If that proves too difficult, then we could just always store FixedPoint64's
  32. // internal representation a 2-parts for each integer. In the 64-bit case,
  33. // one of them will always be unused. (Though, oddly enough, this would
  34. // open the possibility of _128-bit_ fixed-point numbers! Albeit, such things
  35. // would then only work on 64-bit systems and not 32-bit systems, and away
  36. // we go again... until we implement 8-part functions to support 128-bit
  37. // numbers on 32-bit systems. We'll be tempted to use the 256-bit math
  38. // that's now available on 64-bit systems, and then need to implement
  39. // 16-part functions as a way to make 256-bit math compatible with 32-bit
  40. // systems. OH NO OWN GOAL (etc)
  41. // Really, we'll probably just say 64-bit FixedPoint/Rational is good enough
  42. // for now, and we'll expand to 128-bit ONLY if we have a reason to.
  43. // But it'll at least be totally doable, with established patterns to guide
  44. // how it'd be done!)
  45. //
  46.  
  47. /*
  48. declare(strict_types=1);
  49.  
  50. namespace Kickback\Math;
  51. */
  52. final class Testing
  53. {
  54.     /*
  55.     * Apparently PHP devs made it _really_ hard for the `assert` statement to
  56.     * actually halt or report errors when triggered:
  57.     * https://stackoverflow.com/questions/40129963/assert-is-not-working-in-php-so-simple-what-am-i-doing-wrong
  58.     *
  59.     * So here's a version that Just Works.
  60.     *
  61.     * We should probably make sure it turns off in production mode.
  62.     * Right now it just always goes.
  63.     */
  64.     public static function _assert(bool $expr, string $msg = "An assertion was triggered.") : void
  65.     {
  66.         if ( !$expr ) {
  67.             throw new \AssertionError($msg);
  68.         }
  69.     }
  70.  
  71.     // Prevent instantiation/construction of the (static/constant) class.
  72.     /** @return never */
  73.     private function __construct() {
  74.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  75.     }
  76. }
  77. ?>
  78.  
  79. <?php
  80. // //*
  81. // /declare(strict_types=1);
  82. // /
  83. // /namespace Kickback\Math\LowLevel;
  84. // /*/
  85. // /
  86. // /We probably don't need this anymore, because I (probably) figured out exactly
  87. // /when PHP converts int64 variables into floats. (It's on signed overflow.)
  88. // /With that knowledge, we can actually use int64 variables normally.
  89. // /Running on a 32-bit system might still present some challenges to make
  90. // /things work correctly, but that's coming along unexpectedly well, to the point
  91. // /where I might just *make it work* instead of trying to error out if we
  92. // /don't have the ints that we want.
  93. // /
  94. // /final class NumberModel
  95. // /{
  96. // /    public const INVALID = 0;
  97. // /    public const INT32 = 1;
  98. // /    public const FLOAT64 = 2;
  99. // /    public const INT64_STORAGE_ONLY = 3;
  100. // /    public const INT64_ALL_OPS_WORK = 4;
  101. // /
  102. // /    public const FIRST_OF_ALL = self::INT32;
  103. // /    public const LAST_OF_ALL = self::INT64_ALL_OPS_WORK;
  104. // /
  105. // /    public const COUNT_OF_ALL = 1+self::LAST_OF_ALL - self::FIRST_OF_ALL;
  106. // /
  107. // /    private function is_available(int $model) : bool
  108. // /    {
  109. // /        switch($model)
  110. // /        {
  111. // /            case self::INVALID: return false;
  112. // /            case self::INT32:   return true;
  113. // /            case self::FLOAT64: return true; TODO: FIGURE OUT WHEN THIS IS TRUE?!
  114. // /            case self::INT64_STORAGE_ONLY:
  115. // /                $num_int_bits = PHP_INT_SIZE*8;
  116. // /                if ($num_int_bits >= 64) {
  117. // /                    return true;
  118. // /                } else {
  119. // /                    return false;
  120. // /                }
  121. // /            case self::INT64_ALL_OPS_WORK: return false; // TODO: How to know?
  122. // /        }
  123. // /        // Anything else isn't even a math model, much less available.
  124. // /        return false;
  125. // /    }
  126. // /
  127. // /    private static function compute_count_available() : int
  128. // /    {
  129. // /        $count = (int)0;
  130. // /        for ($i = self::FIRST_OF_ALL; $i < self::COUNT_OF_ALL; $i++) {
  131. // /            if (self::is_available($i)) {
  132. // /                $count++;
  133. // /            }
  134. // /        }
  135. // /        return $count;
  136. // /    }
  137. // /
  138. // /    private static int $count_available_ = 0;
  139. // /    public static function count_of_available() : int
  140. // /    {
  141. // /        if ( 0 < self::$count_available_ ) {
  142. // /            return self::$count_available_;
  143. // /        } else {
  144. // /            self::$count_available_ = self::compute_count_available();
  145. // /            return self::$count_available_;
  146. // /        }
  147. // /    }
  148. // /
  149. // /    public static function to_string(int $model) : string
  150. // /    {
  151. // /        switch($model)
  152. // /        {
  153. // /            self::INT32             : return __CLASS__."INT32";
  154. // /            self::FLOAT64           : return __CLASS__."FLOAT64";
  155. // /            self::INT64_STORAGE_ONLY: return __CLASS__."INT64_STORAGE_ONLY";
  156. // /            self::INT64_ALL_OPS_WORK: return __CLASS__."INT64_ALL_OPS_WORK";
  157. // /        }
  158. // /    }
  159. // /
  160. // /    // Prevent instantiation/construction of the (static/constant) class.
  161. // /    /** @return never */
  162. // /    private function __construct() {
  163. // /        throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  164. // /    }
  165. // /}
  166. ?>
  167.  
  168. <?php
  169. /*
  170. declare(strict_types=1);
  171.  
  172. namespace Kickback\Math\LowLevel;
  173. */
  174. final class IntNN
  175. {
  176.     // We want to detect overflow during addition or subtraction in some cases.
  177.     //
  178.     // This function also gives us a chance to ensure that (int+int) addition
  179.     // will always return an `int`, because the `+` operator doesn't actually
  180.     // do that by default in PHP! (It usually will, but not if the result overflows.)
  181.     //
  182.     // In cases where there is a signed integer overflow, PHP will convert the
  183.     // result of the addition into a floating point number. This is not what
  184.     // we want when we expect ordinary twos-complement addition. Though we CAN
  185.     // use this fact, along with the `is_int()` function, to detect all
  186.     // signed overflow.
  187.     //
  188.     // In the case of unsigned integer overflow, we will need to detect it
  189.     // separately using bitwise logic exclusively.
  190.     //
  191.     // The bitwise logic for unsigned overflow detection is as follows:
  192.     //
  193.     // In the below table, `a` and `b` are the inputs. `c` is the result of
  194.     // the addition (after correcting any of PHP's signed-overflow-float-conversion).
  195.     // And `r` is the unsigned carry.
  196.     //
  197.     // To derive our truth table, we observe what happens when we add two
  198.     // 2-bit numbers together to yield a 3-bit number. The most-significant bit (msb)
  199.     // of the 3-bit number will tell us if an overflow occurs.
  200.     //
  201.     // To help us analyse a bit further, we also write a 2nd table to the
  202.     // right of the 1st one, and write down the same expression with only
  203.     // the 2nd bit shown for each variable.
  204.     //
  205.     //  a + b =  c  : a+b = c
  206.     //  00+00 = 000 : 0+0 = 0 -> no carry
  207.     //  00+01 = 001 : 0+0 = 0 -> no carry
  208.     //  00+10 = 010 : 0+1 = 1 -> no carry
  209.     //  00+11 = 011 : 0+1 = 1 -> no carry
  210.     //  01+00 = 001 : 0+0 = 0 -> no carry
  211.     //  01+01 = 010 : 0+0 = 1 -> no carry
  212.     //  01+10 = 011 : 0+1 = 1 -> no carry
  213.     //  01+11 = 100 : 0+1 = 0 -> overflow|carry (unsigned)
  214.     //  10+00 = 010 : 1+0 = 1 -> no carry
  215.     //  10+01 = 011 : 1+0 = 1 -> no carry
  216.     //  10+10 = 100 : 1+1 = 0 -> overflow|carry (unsigned)
  217.     //  10+11 = 101 : 1+1 = 0 -> overflow|carry (unsigned)
  218.     //  11+00 = 011 : 1+0 = 1 -> no carry
  219.     //  11+01 = 100 : 1+0 = 0 -> overflow|carry (unsigned)
  220.     //  11+10 = 101 : 1+1 = 0 -> overflow|carry (unsigned)
  221.     //  11+11 = 110 : 1+1 = 1 -> overflow|carry (unsigned)
  222.     //
  223.     // Notably, the less significant bit of these numbers is inconsequental
  224.     // for unsigned carry detection (as long as we have `c`).
  225.     // We only need to concern ourselves with the most significant bits!
  226.     //
  227.     // If we then eliminate the duplicate rows in the msb-only table,
  228.     // we end up with this:
  229.     //
  230.     //  a+b = c
  231.     //  0+0 = 0 -> no carry
  232.     //  0+0 = 1 -> no carry
  233.     //  0+1 = 0 -> overflow|carry (unsigned)
  234.     //  0+1 = 1 -> no carry
  235.     //  1+0 = 0 -> overflow|carry (unsigned)
  236.     //  1+0 = 1 -> no carry
  237.     //  1+1 = 0 -> overflow|carry (unsigned)
  238.     //  1+1 = 1 -> overflow|carry (unsigned)
  239.     //
  240.     // The above generalizes to all integers with larger bit-widths, because
  241.     // we can always use right-bit-shift to discard everything besides the
  242.     // most significant bit.
  243.     //
  244.     // Now we know the relationship between the `a-b-c` msbs, and the
  245.     // carry-or-not status, which we'll call `r` in our truth table.
  246.     //
  247.     // Truth table that we want:
  248.     //
  249.     // : a b c : r :
  250.     // -------------
  251.     // : 0 0 0 : 0 :
  252.     // : 0 0 1 : 0 :
  253.     // : 0 1 0 : 1 :
  254.     // : 0 1 1 : 0 :
  255.     // : 1 0 0 : 1 :
  256.     // : 1 0 1 : 0 :
  257.     // : 1 1 0 : 1 :
  258.     // : 1 1 1 : 1 :
  259.     //
  260.     // Now we'll test expressions in our truth table until we find something
  261.     // that matches `r` by applying bitwise operations to `a`, `b`, and `c`:
  262.     //
  263.     // : a b c : r : ~c : (a|b) : (a&b) : (a|b)&~c : (a&b)|((a|b)&~c) :
  264.     // ----------------------------------------------------------------
  265.     // : 0 0 0 : 0 :  1 :   0   :   0   :     0    :        0         :
  266.     // : 0 0 1 : 0 :  0 :   0   :   0   :     0    :        0         :
  267.     // : 0 1 0 : 1 :  1 :   1   :   0   :     1    :        1         :
  268.     // : 0 1 1 : 0 :  0 :   1   :   0   :     0    :        0         :
  269.     // : 1 0 0 : 1 :  1 :   1   :   0   :     1    :        1         :
  270.     // : 1 0 1 : 0 :  0 :   1   :   0   :     0    :        0         :
  271.     // : 1 1 0 : 1 :  1 :   1   :   1   :     1    :        1         :
  272.     // : 1 1 1 : 1 :  0 :   1   :   1   :     0    :        0         :
  273.     //
  274.     // We have now learned that the expression `(a&b)|((a|b)&~c)` will always
  275.     // predict the correct value for `r`, which is unsigned carry (overflow).
  276.     //
  277.     // This is probably not the most optimal expression, but it may very well
  278.     // be pretty close. We'd need to simplify aggressively to get any better,
  279.     // and it probably won't matter anyways. Most importantly, we can correctly
  280.     // guess the unsigned carry flag after computing an unsigned addition.
  281.     //
  282.     // In the below add and subtract functions, `a`, `b`, and `c` all have
  283.     // the same name, and `r` is represented by the variable `$u_carry`.
  284.     //
  285.     // We can apply the above expression to entire {32|64}-bit integers, and the
  286.     // sign bits in the integers will have the exact same outcomes as in
  287.     // the above truth table.
  288.     //
  289.     // Then we'll need to bit-shift-right by all bits but one
  290.     // (63 bits in the 64-bit case, or 31 bits in the 32-bit case)
  291.     // to get a {0,1} set indicating no carry (0) or carry (1):
  292.     // `$u_carry = (($a&$b)|(($a|$b)&~$c)) >> {31|63};`
  293.     //
  294.     // However, because PHP does not provide an _unsigned_ shift-right,
  295.     // any shift-right on a value with a set sign-bit (sign=1) will result
  296.     // in the variable being filled with 1s.
  297.     //
  298.     // So after the bit-shift, `$u_carry` will either be 0 or it will be 0xFFFF...FFFF.
  299.     // But we want it to be 0 or 1.
  300.     //
  301.     // So we can clear the sign extension bits by ANDing with 1:
  302.     // `$u_carry = (1 & ((($a&$b)|(($a|$b)&~$c)) >> {31|63}));`
  303.  
  304.     private const NUM_BITS_IN_INT = PHP_INT_SIZE*8;
  305.     private const SIGN_SHIFT = (self::NUM_BITS_IN_INT - 1);
  306.     private const HALF_SHIFT = (self::NUM_BITS_IN_INT >> 1);
  307.     private const HALF_MASK  = ((1 << self::HALF_SHIFT)-1);
  308.  
  309.     public static function add(int $a, int $b) : int
  310.     {
  311.         $s_carry = (int)0;
  312.         $u_carry = (int)0;
  313.         return self::subtract_with_borrow($a, $b, $s_carry, $u_carry);
  314.     }
  315.  
  316.     public static function add_with_carry(int $a, int $b, int &$s_carry, &$u_carry) : int
  317.     {
  318.         $c = $a + $b;
  319.  
  320.         if ( is_int($c) ) {
  321.             // Most common branch.
  322.             $u_carry = (1 & ((($a&$b)|(($a|$b)&~$c)) >> self::SIGN_SHIFT));
  323.             $s_carry = 0;
  324.             //echo(sprintf("add_with_carry: a=%016X, b=%016X, c=%016X, (a&b)=%016X, ((a|b)&~c)=%016X, carry=%016X\n", $a, $b, $c, ($a&$b), (($a|$b)&~$c), $u_carry));
  325.             return $c;
  326.         } else {
  327.             // !is_int($c)
  328.             // Uncommon branch: signed overflow.
  329.             $s_carry = 1;
  330.             return self::add_by_parts($a,$b,$u_carry);
  331.         }
  332.     }
  333.  
  334.     // An earlier (and non-working) version of the above function used
  335.     // the below bit-twiddle expressions, from
  336.     //   https://grack.com/blog/2022/12/20/deriving-a-bit-twiddling-hack/
  337.     // to check for overflow. I didn't notice at the time, but these turned
  338.     // out to be for SIGNED overflow detection. Since PHP does that for us
  339.     // (and _to_ us, even), we won't need bit twiddles for the signed case.
  340.     // So we won't be needing these:
  341.     //$carry = ((($c ^ $a) & ($c ^ $b)) >> 63);
  342.     //$carry = ((~($a ^ $b) & ($c ^ $a)) >> 63);
  343.  
  344.     // `add_by_parts(int a, int b, int u_carry):int`
  345.     //
  346.     // Slow-but-accurate addition function that breaks a 64-bit integer
  347.     // into two 32-bit parts (lo and hi), then adds them together.
  348.     // The result of each addition is 33-bits wide, but that fits just
  349.     // fine into a 64-bit integer. After carry propagation, the (lo and hi)
  350.     // parts are reassembled into the final 64-bit integer.
  351.     //
  352.     // This is potentially more complicated and slow than the "just use `+`"
  353.     // function, but it has the tremendous advantage of never converting
  354.     // the final result into a floating point number with incorrect contents.
  355.     // (PHP will automatically convert the result of integer expressions
  356.     // into a floating point value if the integer expression results in
  357.     // a signed overflow.)
  358.     //
  359.     // Since this function already needs to do carry propagation, it can
  360.     // return the unsigned carry status in the `&$u_carry` parameter
  361.     // as a free action.
  362.     //
  363.     // Signed carry is not as trivial, so it is not performed here.
  364.     // (It's still potentially very "easy", but it should probably be
  365.     // computed by the caller only if it's needed.)
  366.     //
  367.     private static function add_by_parts(int $a, int $b, &$u_carry) : int
  368.     {
  369.         // Shorthand:
  370.         $half_mask  = self::HALF_MASK;  // =0xFFFF on 32-bit host, and =0xFFFF_FFFF on 64-bit host.
  371.         $half_shift = self::HALF_SHIFT; // =16 when on 32-bit host, and =32 when on 64-bit host.
  372.  
  373.         // Unpack $a:
  374.         $al = $a & $half_mask;   // $al : Low  {16|32} bits of $a
  375.         $ah = $a >> $half_shift; // $ah : High {16|32} bits of $a
  376.         $ah &= $half_mask;       // Leave space in $ah for carry/overflow detection.
  377.  
  378.         // Unpack $b:
  379.         $bl = $b & $half_mask;   // $bl : Low  {16|32} bits of $b
  380.         $bh = $b >> $half_shift; // $bh : High {16|32} bits of $b
  381.         $bh &= $half_mask;       // Leave space in $bh for carry/overflow detection.
  382.  
  383.         // Part-wise addition:
  384.         $cl = $al + $bl; // $cl : Low  {16|32} bits of result from $a+$b
  385.         $ch = $ah + $bh; // $ch : High {16|32} bits of result from $a+$b, so far unadjusted for carry propagation.
  386.  
  387.         // Propagate low-carry into $ch:
  388.         $carry_lo = ($cl >> $half_shift);
  389.         $ch += $carry_lo;
  390.  
  391.         // Detect unsigned carry from the entire result:
  392.         $carry_hi = ($ch >> $half_shift);
  393.         //echo(sprintf("add_by_parts: ah=%016X, al=%016X, bh=%016X, bl=%016X, ch=%016X, cl=%016X, carry_hi=%016X, carry_lo=%016X\n", $ah,$al,  $bh,$bl,  $ch,$cl,  $carry_hi,$carry_lo));
  394.  
  395.         // Pack the result back into a single integer:
  396.         $c = ($ch << $half_shift) | ($cl & $half_mask);
  397.         //echo(sprintf("add_by_parts: a=%016X, b=%016X, c=%016X, carry=%016X\n", $a, $b, $c, $carry_hi));
  398.  
  399.         // Return results:
  400.         $u_carry = $carry_hi;
  401.         return $c;
  402.     }
  403.  
  404.     public static function unittest_add_with_carry() : void
  405.     {
  406.         echo(__CLASS__."::".__FUNCTION__."...");
  407.  
  408.         // Shorthand to avoid having to write out long 64-bit hex values.
  409.         $x0000 = (int)0;
  410.         $x0001 = (int)1;
  411.         $x0002 = (int)2;
  412.         $x8000 = PHP_INT_MIN;
  413.         $x8001 = PHP_INT_MIN+1;
  414.         $xFFFF = (int)-1;
  415.         $xFFFE = (int)-2;
  416.         $x7FFF = PHP_INT_MAX;
  417.         $x7FFE = PHP_INT_MAX-1;
  418.  
  419.         $s_carry = (int)0;
  420.         $u_carry = (int)0;
  421.  
  422.         Testing::_assert($x0000 === self::add_with_carry($x0000,$x0000,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  423.         Testing::_assert($x0001 === self::add_with_carry($x0000,$x0001,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  424.         Testing::_assert($x0001 === self::add_with_carry($x0001,$x0000,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  425.         Testing::_assert($x0002 === self::add_with_carry($x0001,$x0001,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  426.  
  427.         Testing::_assert($x0000 === self::add_with_carry($x0000,$x0000,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  428.         Testing::_assert($xFFFF === self::add_with_carry($x0000,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  429.         Testing::_assert($xFFFF === self::add_with_carry($xFFFF,$x0000,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(0 === $u_carry);
  430.         Testing::_assert($xFFFE === self::add_with_carry($xFFFF,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(1 === $u_carry);
  431.  
  432.         Testing::_assert($x0000 === self::add_with_carry($x0001,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(1 === $u_carry);
  433.         Testing::_assert($x7FFF === self::add_with_carry($xFFFF,$x8000,$s_carry,$u_carry));
  434.         Testing::_assert(1 === $s_carry);
  435.         Testing::_assert(1 === $u_carry); // PHP Math is b0rked?! It says 0xFFFF...FFFF + 0x8000...0000 === 0x8000...0000, which is WRONG. It should return 0x7FFF...FFFF.
  436.  
  437.         Testing::_assert($x0001 === self::add_with_carry($x0002,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(1 === $u_carry);
  438.         Testing::_assert($x7FFE === self::add_with_carry($xFFFE,$x8000,$s_carry,$u_carry)); Testing::_assert(1 === $s_carry); Testing::_assert(1 === $u_carry); // PHP Math is b0rked?!
  439.  
  440.         Testing::_assert($x8000 === self::add_with_carry($x8001,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(1 === $u_carry);
  441.         Testing::_assert($x7FFE === self::add_with_carry($x7FFF,$xFFFF,$s_carry,$u_carry)); Testing::_assert(0 === $s_carry); Testing::_assert(1 === $u_carry);
  442.  
  443.         echo(" done.\n");
  444.     }
  445.  
  446.     public static function subtract(int $a, int $b) : int
  447.     {
  448.         $s_borrow = (int)0;
  449.         $u_borrow = (int)0;
  450.         return self::subtract_with_borrow($a, $b, $s_borrow, $u_borrow);
  451.     }
  452.  
  453.     public static function subtract_with_borrow(int $a, int $b, int &$s_borrow, &$u_borrow) : int
  454.     {
  455.         $c = $a - $b;
  456.  
  457.         //echo(sprintf("a=%016X, b=%016X, c=%016X, (a&b)=%016X, ((a|b)&~c)=%016X, carry=%016X\n", $a, $b, $c, ($a&$b), (($a|$b)&~$c), $carry));
  458.  
  459.         if ( is_int($c) ) {
  460.             // Most common branch.
  461.             $u_borrow = (1 & ((($a&$b)|(($a|$b)&~$c)) >> self::SIGN_SHIFT)); // TODO: CHECK THIS LOGIC. It might be different for subtraction?
  462.             $s_borrow = 0;
  463.             return $c;
  464.         } else {
  465.             // !is_int($c)
  466.             // Uncommon branch: signed overflow.
  467.             $s_borrow = 1;
  468.             return self::add_by_parts($a,$b,$u_borrow);
  469.         }
  470.     }
  471.  
  472.     // Pretty much the same thing as `add_by_parts`, except for subtraction.
  473.     private static function subtract_by_parts(int $a, int $b, &$u_borrow) : int
  474.     {
  475.         // Shorthand:
  476.         $half_mask  = self::HALF_MASK;  // =0xFFFF on 32-bit host, and =0xFFFF_FFFF on 64-bit host.
  477.         $half_shift = self::HALF_SHIFT; // =16 when on 32-bit host, and =32 when on 64-bit host.
  478.  
  479.         // Unpack $a:
  480.         $al = $a & $half_mask;   // $al : Low  {16|32} bits of $a
  481.         $ah = $a >> $half_shift; // $ah : High {16|32} bits of $a
  482.         $ah &= $half_mask;       // Leave space in $ah for borrow/overflow detection.
  483.  
  484.         // Unpack $b:
  485.         $bl = $b & $half_mask;   // $bl : Low  {16|32} bits of $b
  486.         $bh = $b >> $half_shift; // $bh : High {16|32} bits of $b
  487.         $bh &= $half_mask;       // Leave space in $bh for borrow/overflow detection.
  488.  
  489.         // Part-wise subtraction:
  490.         $cl = $al - $bl; // $cl : Low  {16|32} bits of result from $a+$b
  491.         $ch = $ah - $bh; // $ch : High {16|32} bits of result from $a+$b, so far unadjusted for borrow propagation.
  492.  
  493.         // Propagate low-borrow into $ch:
  494.         $borrow_lo = (($cl >> $half_shift) & 1);
  495.         $ch -= $borrow_lo;
  496.  
  497.         // Detect unsigned borrow from the entire result:
  498.         $borrow_hi = (($ch >> $half_shift) & 1);
  499.  
  500.         // Pack the result back into a single integer:
  501.         $c = ($ch << $half_shift) | ($cl & $half_mask);
  502.  
  503.         // Return results:
  504.         $u_borrow = $borrow_hi;
  505.         return $c;
  506.     }
  507.  
  508.  
  509.     public static function unittest() : void
  510.     {
  511.         echo("Unittests for ".__CLASS__."\n");
  512.         self::unittest_add_with_carry();
  513.         echo("\n");
  514.     }
  515. }
  516. ?>
  517.  
  518. <?php
  519. /*
  520. declare(strict_types=1);
  521.  
  522. namespace Kickback\Math\LowLevel;
  523. */
  524. final class NumberModelsIterator implements \Iterator, \Countable
  525. {
  526.     private int $num_passed = 0;
  527.     private int $which = NumberModel::INVALID;
  528.     private bool $include_unavailable = false;
  529.  
  530.     public function count() : int {
  531.         if ( $this->include_unavailable ) {
  532.             return NumberModel::COUNT_OF_ALL - $this->num_passed;
  533.         } else {
  534.             return NumberModel::count_of_available() - $this->num_passed;
  535.         }
  536.     }
  537.  
  538.     public function current() : int {
  539.         return $this->which;
  540.     }
  541.  
  542.     public function key() : int {
  543.         return $this->num_passed;
  544.     }
  545.  
  546.     private function scan_until_current_elem_is_available() : void
  547.     {
  548.         if ( $this->include_unavailable ) {
  549.             return;
  550.         }
  551.  
  552.         while (
  553.             ($this->which < NumberModel::count_of_available()) &&
  554.             (!NumberModel::is_available($this->which)) )
  555.         {
  556.             $this->which++;
  557.         }
  558.         return;
  559.     }
  560.  
  561.     public function next() : void
  562.     {
  563.         if ($this->which < $this->count())
  564.         {
  565.             $this->num_passed++;
  566.             $this->which++;
  567.             $this->scan_until_current_elem_is_available();
  568.         }
  569.     }
  570.  
  571.     public function rewind() : void
  572.     {
  573.         $this->num_passed = 0;
  574.         $this->which = NumberModel::FIRST_OF_ALL;
  575.         $this->scan_until_current_elem_is_available();
  576.     }
  577.  
  578.     public function valid() : bool {
  579.         return (0 < $this->count());
  580.     }
  581. }
  582. ?>
  583.  
  584. <?php
  585. // /
  586. // / This class is almost certainly obsolete now.
  587. // / Will probably be removed in later revisions.
  588. // / But hey, it was very useful: its failing unittests taught us about
  589. // / how PHP will convert integer expression into float values _sometimes_
  590. // / instead of just using ordinary twos-complement arithmetic.
  591. // / Yeah... that, uh, matters!
  592. // /
  593. // /declare(strict_types=1);
  594. // /
  595. // /namespace Kickback\Math\LowLevel;
  596. // /
  597. // /final class Logic64
  598. // /{
  599. // /    // IDEA: The int+int->float problem is probably PHP detecting overflow and
  600. // /    // converting the value to floating point when an overflow occurs.
  601. // /    // If true, then we can do 63-bit math safely. We'll even have efficient
  602. // /    // carry extraction for the 63-bit ints, because PHP may very well
  603. // /    // respect the int+int->int rule when doing such things.
  604. // /    // However, the extra bit will be costly!
  605. // /
  606. // /    private static \SplFixedArray $math_models_ = new \SplFixedArray();
  607. // /    public static function math_models() : \SplFixedArray
  608. // /    {
  609. // /        if ( 0 < $math_models->count() ) {
  610. // /            return self::$math_models_;
  611. // /        }
  612. // /
  613. // /        self::$math_models_->append(MATH_MODEL_FLOAT64);
  614. // /        self::$math_models_->append(MATH_MODEL_INT32);
  615. // /
  616. // /        $num_int_bits = PHP_INT_SIZE*8;
  617. // /        if ($num_int_bits >= 64) {
  618. // /            self::$math_models_->append(MATH_MODEL_INT64_STORAGE_ONLY);
  619. // /        }
  620. // /        return self::$math_models_;
  621. // /    }
  622. // /
  623. // /    private const MASK_SIGN = (1 << 63);
  624. // /
  625. // /    public static function add(int $number_model,   int $ah, int $al,   int $bh, int $bl) : int
  626. // /    {
  627. // /        $carry = (int)0;
  628. // /        return self::add_with_carry($number_model,  $ah,$al,  $bh,$bl,  $carry);
  629. // /    }
  630. // /
  631. // /    public static function add_with_carry(int $a, int $b, int &$carry) : int
  632. // /    {
  633. // /        $al = $a & 0xFFFF_FFFF;
  634. // /        $ah = $a >> 32;
  635. // /        $ah &= 0xFFFF_FFFF;
  636. // /
  637. // /        $bl = $b & 0xFFFF_FFFF;
  638. // /        $bh = $b >> 32;
  639. // /        $bh &= 0xFFFF_FFFF;
  640. // /
  641. // /        $cl = $al + $bl;
  642. // /        $ch = $ah + $bh;
  643. // /
  644. // /        $carry_lo = (($cl >> 32) & 1);
  645. // /        $ch += $carry_lo;
  646. // /
  647. // /        $carry_hi = (($ch >> 32) & 1);
  648. // /        $c = ($ch << 32) | ($cl & 0xFFFF_FFFF);
  649. // /
  650. // /        $carry = $carry_hi;
  651. // /        return $c;
  652. // /    }
  653. // /
  654. // /    public static function subtract(int $a, int $b) : int
  655. // /    {
  656. // /        $carry = (int)0;
  657. // /        return self::subtract_with_carry($a, $b, $carry);
  658. // /    }
  659. // /
  660. // /    public static function subtract_with_borrow(int $a, int $b, int &$borrow) : int
  661. // /    {
  662. // /        $al = $a & 0xFFFF_FFFF;
  663. // /        $ah = $a >> 32;
  664. // /        $ah &= 0xFFFF_FFFF;
  665. // /
  666. // /        $bl = $b & 0xFFFF_FFFF;
  667. // /        $bh = $b >> 32;
  668. // /        $bh &= 0xFFFF_FFFF;
  669. // /
  670. // /        $cl = $al - $bl;
  671. // /        $ch = $ah - $bh;
  672. // /
  673. // /        $borrow_lo = (($cl >> 32) & 1);
  674. // /        $ch -= $borrow_lo;
  675. // /
  676. // /        $borrow_hi = (($ch >> 32) & 1);
  677. // /        $c = ($ch << 32) | ($cl & 0xFFFF_FFFF);
  678. // /
  679. // /        $borrow = $borrow_hi;
  680. // /        return $c;
  681. // /    }
  682. // /
  683. // /    public static function unittest_add_with_carry() : void
  684. // /    {
  685. // /        echo(__CLASS__."::".__FUNCTION__."...");
  686. // /
  687. // /        // Shorthand to avoid having to write out long 64-bit hex values.
  688. // /        $x0000 = (int)0;
  689. // /        $x0001 = (int)1;
  690. // /        $x0002 = (int)2;
  691. // /        $x8000 = PHP_INT_MIN;
  692. // /        $x8001 = PHP_INT_MIN+1;
  693. // /        $xFFFF = (int)-1;
  694. // /        $xFFFE = (int)-2;
  695. // /        $x7FFF = PHP_INT_MAX;
  696. // /        $x7FFE = PHP_INT_MAX-1;
  697. // /
  698. // /        $carry = (int)0;
  699. // /
  700. // /        Testing::_assert($x0000 === self::add_with_carry($x0000,$x0000,$carry)); Testing::_assert(0 === $carry);
  701. // /        Testing::_assert($x0001 === self::add_with_carry($x0000,$x0001,$carry)); Testing::_assert(0 === $carry);
  702. // /        Testing::_assert($x0001 === self::add_with_carry($x0001,$x0000,$carry)); Testing::_assert(0 === $carry);
  703. // /        Testing::_assert($x0002 === self::add_with_carry($x0001,$x0001,$carry)); Testing::_assert(0 === $carry);
  704. // /
  705. // /        Testing::_assert($x0000 === self::add_with_carry($x0000,$x0000,$carry)); Testing::_assert(0 === $carry);
  706. // /        Testing::_assert($xFFFF === self::add_with_carry($x0000,$xFFFF,$carry)); Testing::_assert(0 === $carry);
  707. // /        Testing::_assert($xFFFF === self::add_with_carry($xFFFF,$x0000,$carry)); Testing::_assert(0 === $carry);
  708. // /        Testing::_assert($xFFFE === self::add_with_carry($xFFFF,$xFFFF,$carry)); Testing::_assert(1 === $carry);
  709. // /
  710. // /        Testing::_assert($x0000 === self::add_with_carry($x0001,$xFFFF,$carry)); Testing::_assert(1 === $carry);
  711. // /        // Testing::_assert($x7FFF === self::add_with_carry($xFFFF,$x8000,$carry)); Testing::_assert(1 === $carry); PHP Math is b0rked?! It says 0xFFFF...FFFF + 0x8000...0000 === 0x8000...0000, which is WRONG. It should return 0x7FFF...FFFF.
  712. // /
  713. // /        Testing::_assert($x0001 === self::add_with_carry($x0002,$xFFFF,$carry)); Testing::_assert(1 === $carry);
  714. // /        // Testing::_assert($x7FFE === self::add_with_carry($xFFFE,$x8000,$carry)); Testing::_assert(1 === $carry); PHP Math is b0rked?!
  715. // /
  716. // /        Testing::_assert($x8000 === self::add_with_carry($x8001,$xFFFF,$carry)); Testing::_assert(1 === $carry);
  717. // /        Testing::_assert($x7FFE === self::add_with_carry($x7FFF,$xFFFF,$carry)); Testing::_assert(1 === $carry);
  718. // /
  719. // /        echo(" done.\n");
  720. // /    }
  721. // /
  722. // /    public static function ensure_ints_are_64bit() : void
  723. // /    {
  724. // /        $num_int_bits = PHP_INT_SIZE*8;
  725. // /        if ($num_int_bits < 64)
  726. // /        {
  727. // /            throw new SomethingSomethingException(
  728. // /                "This system features ${num_int_bits}-bit integers, " .
  729. // /                "but at least 64-bit wide integers are required for financial calculations.");
  730. // /        }
  731. // /    }
  732. // /
  733. // /    public static function unittest() : void
  734. // /    {
  735. // /        echo("Unittests for ".__CLASS__."\n");
  736. // /        self::ensure_ints_are_64bit();
  737. // /        self::unittest_add_with_carry();
  738. // /        echo("\n");
  739. // /    }
  740. // /
  741. // /    // Prevent instantiation/construction of the (static/constant) class.
  742. // /    /// @return never
  743. // /    private function __construct() {
  744. // /        throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  745. // /    }
  746. // /}
  747. // /*/
  748. /*
  749. ?>
  750.  
  751. <?php
  752. //declare(strict_types=1);
  753. namespace Kickback\Math;
  754.  
  755. use \Kickback\Math\Int64;
  756. */
  757. final class Int128
  758. {
  759.     private const MASK_LO = ((1 << 32) - 1);
  760.     private const MASK_SIGN = (1 << 63);
  761. /*
  762.     // Internal version of `cmp`.
  763.     // This version does not clamp it's outputs to the set {-1, 0, 1}.
  764.     private static function cmp_unsaturated(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
  765.     {
  766.         $res = $a_hi - $b_hi;
  767.         if ( $res ) {
  768.             return $res;
  769.         }
  770.  
  771.         return ($a_lo - $b_lo);
  772.     }
  773. */
  774.     // Below is a branchless implementation of a comparison algorithm for
  775.     // a comparison operation on 128-bit integers that area each specified
  776.     // as `hi` and `lo` 64-bit integers.
  777.     //
  778.     // This probably doesn't need to be branchless, and I doubt that PHP code
  779.     // will benefit from that very much.
  780.     //
  781.     // I'm mostly just practicing my skill at writing branchless code, because
  782.     // it is a useful skill to have in other places.
  783.     //
  784.     // I'm still going to use the branchless version because (1) it works
  785.     // and (2) it _might_ help somewhere. It's probably the best way to do things.
  786.     // The only disadvantage would be that it just took longer to write.
  787.     //
  788.     // -- Lily Joan  2024-07-11
  789.     //
  790.     //
  791.     // With the personal statement out of the way, here's out it actually works:
  792.     //
  793.     // The code looks like this:
  794.     //
  795.     // $d1 = $a_hi - $b_hi;
  796.     // $d0 = $a_lo - $b_lo;
  797.     //
  798.     // // Set all bits of `$mask` to 1 if p is nonzero.
  799.     // // Clear all bits if p is zero.
  800.     // $mask = ($d1 | -$d1) >> 63;
  801.     //
  802.     // $r = $d1 | ($d0 & ~$mask);
  803.     //
  804.     // $p = $r >> 63;
  805.     // $q = ~($r-1) >> 63
  806.     //
  807.     // return ($p-$q)|$p;
  808.     //
  809.     // =========================================================================
  810.     // ================== $mask = ($d1 | -$d1) >> 63; ==========================
  811.     // The statement `$mask = ($d1 | -$d1) >> 63;` calculates a mask that is all 1s
  812.     // whenever `$d1` is non-zero, and all 0s whenever `$d1` is zero.
  813.     //
  814.     // The first step to understanding this is to realize that the subexpression
  815.     // `($d1 | -$d1)` will, for non-zero inputs, ALWAYS have its sign-bit set.
  816.     // For the input of zero, it will have its sign-bit clear.
  817.     //
  818.     // It does this by exploiting the fact that the number 0 is neither positive nor
  819.     // negative. Here are the possibilities:
  820.     // :    $d1    : sign($d1) :   -$d1     : sign(-$d1) : sign($d1) | sign(-$d1) :
  821.     // ----------------------------------------------------------------------------
  822.     // : $d1  >  0 :    0      : -$d1  <  0 :     1      :           1            :
  823.     // : $d1 === 0 :    0      : -$d1 === 0 :     0      :           0            :
  824.     // : $d1  <  0 :    1      : -$d1  >  0 :     0      :           1            :
  825.     //
  826.     // Once we have a sign bit that is 1 for all (and only) non-zero inputs,
  827.     // we can shift the value to the right by at least 63 bits to sign-extend the
  828.     // the sign bit to fill the entire variable.
  829.     //
  830.     // Here's another truth table showing what happens for 2-bit integers:
  831.     //
  832.     // :             $d1             : 00 : 01 : 10 : 11 :
  833.     // ---------------------------------------------------
  834.     // :            -$d1             : 00 : 10 : 10 : 01 :
  835.     // :          $d1 | -$d1         : 00 : 11 : 10 : 11 :
  836.     // : ($d1 | -$d1) >> (#bits - 1) : 00 : 11 : 11 : 11 :
  837.     //
  838.     // =========================================================================
  839.     // ================== $r = $d1 | ($d0 & ~$mask); ===========================
  840.     // The statement `$r = $d1 | ($d0 & ~$mask);` selects which "digit" of
  841.     // the 128-bit operands is being used for comparison. It uses information
  842.     // in `$mask` to ensure it picks the high parts if they differ, and
  843.     // picks the low parts if the high parts are equal.
  844.     //
  845.     // The full expression would look more like this:
  846.     //   `($d1 & $mask) | ($d0 & ~$mask);`
  847.     //
  848.     // We can simplify it to
  849.     //   `$d1 | ($d0 & ~$mask)`
  850.     // by noting that $mask is derived from $d1 and will always be 0 when
  851.     // $d1 is 0. As such, there is no reason to AND the two together,
  852.     // because it will never alter the resulting value.
  853.     //
  854.     // (thought discontinued)
  855.  
  856.  
  857.  
  858.     // Other notes, most about how to derive p and q logic:
  859.     // :    p    : 00 : 01 : 10 : 11 : 00 : 01 : 10 : 11 : 00 : 01 : 10 : 11 : 00 : 01 : 10 : 11 :
  860.     // :    q    : 00 : 00 : 00 : 00 : 01 : 01 : 01 : 01 : 10 : 10 : 10 : 10 : 11 : 11 : 11 : 11 :
  861.     // :    m    : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 : 00 :
  862.     // -------------------------------------------------------------------------------------------
  863.     // :  p ^ q  : 00 : 01 : 10 : 11 : 01 : 00 : 11 : 10 : 10 : 11 : 00 : 01 : 11 : 10 : 01 : 00 :
  864.     //
  865.     //
  866.     // : p | -p  : 00 : 11 : 10 : 11 :
  867.     // : >> bits : 00 : 11 : 11 : 11 :
  868.     //
  869.     // $p = $r >> 63
  870.     //
  871.     // >0 =0 <0 min
  872.     // 00 00 11 11 r >> 63 == p
  873.     // //00 11 11 00 r-1 >> 63
  874.     // 11 00 00 11 ~(r-1) >> 63 == q
  875.     // //11 00 11 00 p^q
  876.     // //00 00 00 11 p&q
  877.     // //11 00 11 11 p|q
  878.     // 01 00 11 00 p-q
  879.     // 01 00 11 11 (p-q)|p
  880.  
  881.     /*
  882.     *  Returns -1, 0, or 1 based on a comparison of the given 128-bit integers.
  883.     *
  884.     *  Defined such that:
  885.     *  * a < b implies cmp(a,b) < 0
  886.     *  * a > b implies cmp(a,b) > 0
  887.     *  * a == b implies cmp(a,b) == 0.
  888.     */
  889.     public static function cmp(int $a_hi, int $a_lo,   int $b_hi, int $b_lo) : int
  890.     {
  891.         // TODO: BUG: Needs work to support integers with parts >PHP_INT_MAX
  892.         // (Thus it currently does not pass unittests.)
  893.         // (We'll probably end up looking at u_carry flags. This is a common
  894.         // technique for assembly/CPUs to detect if one number is larger
  895.         // than another: just subtract them and see if a borrow flag is set!)
  896.         // TODO: Also, we should have signed+unsigned versions of this!
  897.         $d1 = IntNN::subtract($a_hi, $b_hi);
  898.         $d0 = IntNN::subtract($a_lo, $b_lo);
  899.  
  900.         // Set all bits of `$mask` to 1 if p is nonzero.
  901.         // Clear all bits if p is zero.
  902.         $mask = ($d1 | -$d1) >> 63;
  903.  
  904.         $r = $d1 | ($d0 & ~$mask);
  905.  
  906.         $p = $r >> 63;
  907.         $q = ~($r-1) >> 63;
  908.         $r2 = (($p-$q)|$p);
  909.         echo(sprintf("\nd1=%016X, d0=%016X, mask=%016X, r=%016X, p=%016X, q=%016X, r2=%016X", $d1,$d0, $mask, $r, $p, $q, $r2));
  910.  
  911.         return $r2;
  912.     }
  913. /*
  914.     public static function cmp(int $a_hi, int $a_lo, int $b_hi, int $b_lo) : int
  915.     {
  916.         // First, we do the comparison with clamping the results to {-1, 0, 1}
  917.         $r = self::cmp_unsaturated($a_hi,$a_lo,  $b_hi,$b_lo);
  918.  
  919.         // Then we clamp the results to {-1, 0, 1} so that we fulfill the function's contract.
  920.         // We'll do this by noting that if we shift the result _right_ by 62 bits
  921.         // then we have 2 bits left with 4 possible states, and 3 of those states
  922.         // are the results we are trying to acquire.
  923.         //
  924.         // The possible values of $r at the start look like this:
  925.         //
  926.         //  case # : target values :   $r, as a bit pattern    :   $r >> 62    : $r >> 62, low 2 bits.
  927.         // ----------------------------------------------------------------------------------------------------
  928.         //    1    : all positives : 0b0???_????_..._????_???? : 0b0000_..._000? :   0b0?
  929.         //    2    : zero, 0       : 0b0000_0000_..._0000_0000 : 0b0000_..._0000 :   0b00
  930.         //    3    : all negatives : 0b1???_????_..._????_???? : 0b1111_..._111? :   0b1?
  931.         //    4    : corner case   : 0b1000_0000_..._0000_0000 : 0b1111_..._1111 :   0b11
  932.         //
  933.         // But we already have a problem, because cases (1) and (2) aren't always
  934.         // distinct (currently, 0 is _actually_ in case (1)).
  935.         //
  936.         // So we'll subtract by 1 before doing our bitshift, and then this happens:
  937.         //
  938.         //  case # : target values :  $r-1, as a bit pattern   :   $r-1 >> 62  : $r-1 >> 62, low 2 bits.
  939.         // ----------------------------------------------------------------------------------
  940.         //    1    : all positives : 0b0???_????_..._????_???? : 0b0000...000? :   0b0?
  941.         //    2    : zero, 0       : 0b1111_1111_..._1111_1111 : 0b1111...1111 :   0b11
  942.         //    3    : all negatives : 0b1???_????_..._????_???? : 0b1111...111? :   0b1?
  943.         //    4    : corner case   : 0b0111_1111_..._1111_1111 : 0b0000...0000 :   0b00
  944.         //
  945.  
  946.         // OLD--------------------
  947.         //
  948.         //  case # : target values :   $r, as a bit pattern    :   $r >> 62    : $r >> 62, low 2 bits.
  949.         // ----------------------------------------------------------------------------------------------------
  950.         //    1a   : all pos > 1   : 0b0???_????_..._????_???0 : 0b0000_..._0000 :   0b00
  951.         //    1b   : all pos > 1   : 0b0???_????_..._????_???1 : 0b0000_..._0001 :   0b01
  952.         //    1.5  : one, 1        : 0b0000_0000_..._0000_0001 : 0b0000_..._0001 :   0b01
  953.         //    2    : zero, 0       : 0b0000_0000_..._0000_0000 : 0b0000_..._0000 :   0b00
  954.         //    3a   : all neg > MIN : 0b1???_????_..._????_???0 : 0b1111_..._1110 :   0b10
  955.         //    3b   : all neg > MIN : 0b1???_????_..._????_???1 : 0b1111_..._1111 :   0b11
  956.         //    4    : corner case   : 0b1000_0000_..._0000_0000 : 0b1111_..._1111 :   0b11
  957.         //
  958.         // But we already have a problem, because cases (1) and (2) aren't always
  959.         // distinct (currently, 0 is _actually_ in case (1)). In the same way,
  960.         // cases (3) and (4) also aren't distinct.
  961.         //
  962.         // So we'll subtract by 1 before doing our bitshift, and then this happens:
  963.         //
  964.         //  case # : target values :  $r-1, as a bit pattern   :   $r-1 >> 62  : $r-1 >> 62, low 2 bits.
  965.         // ----------------------------------------------------------------------------------
  966.         //    1a   : all pos > 1   : 0b0???_????_..._????_???1 : 0b0000_..._0001 :   0b01
  967.         //    1b   : all pos > 1   : 0b0???_????_..._????_???0 : 0b0000_..._0000 :   0b00
  968.         //    1.5  : one, 1        : 0b0000_0000_..._0000_0000 : 0b0000_..._0000 :   0b00
  969.         //    2    : zero, 0       : 0b1111_1111_..._1111_1111 : 0b1111_..._1111 :   0b11
  970.         //    3a   : all neg > MIN : 0b1???_????_..._????_???1 : 0b1111_..._1111 :   0b11
  971.         //    3b   : all neg > MIN : 0b1???_????_..._????_???0 : 0b1111_..._1110 :   0b10
  972.         //    4    : corner case   : 0b1111_1111_..._1111_1111 : 0b1111_..._1111 :   0b11
  973.         // TODO: Did my earlier self fnck up?
  974.         //
  975.         if ( $r === 0x8000_0000_0000_0000 ) {
  976.             return -1;
  977.         }
  978.         $r--;      //    >0    =0    <0
  979.         $r >>= 62; // {0b000,0b111,0b110} = {0,-1,-2}
  980.         $r++;      // {0b001,0b000,0b111} = {1, 0,-1}
  981.         return $r;
  982.     }
  983.     */
  984.  
  985.     public static function unittest_cmp()
  986.     {
  987.         echo(__CLASS__."::".__FUNCTION__."...");
  988.  
  989.         Testing::_assert(self::cmp( 0, 0,   0, 0) === 0);
  990.         Testing::_assert(self::cmp( 0, 0,   0, 1)  <  0);
  991.         Testing::_assert(self::cmp( 0, 1,   0, 0)  >  0);
  992.         Testing::_assert(self::cmp( 0, 1,   0, 1) === 0);
  993.         Testing::_assert(self::cmp(-1, 0,  -1, 0) === 0);
  994.         Testing::_assert(self::cmp(-1, 0,   0, 0)  <  0);
  995.         Testing::_assert(self::cmp(-1, 0,   1, 0)  <  0);
  996.         Testing::_assert(self::cmp( 0, 0,  -1, 0)  >  0);
  997.         Testing::_assert(self::cmp( 0, 0,   0, 0) === 0);
  998.         Testing::_assert(self::cmp( 0, 0,   1, 0)  <  0);
  999.         Testing::_assert(self::cmp( 1, 0,  -1, 0)  >  0);
  1000.         Testing::_assert(self::cmp( 1, 0,   0, 0)  >  0);
  1001.         Testing::_assert(self::cmp( 1, 0,   1, 0) === 0);
  1002.  
  1003.         // Ensure that it's returning values in the set {-1,0,1}, and not in the set {-N,0,N} or somesuch.
  1004.         Testing::_assert(self::cmp(  0, 0,   16, 0) === -1);
  1005.         Testing::_assert(self::cmp( 16, 0,    0, 0) ===  1);
  1006.         Testing::_assert(self::cmp(-16, 0,    0, 0) === -1);
  1007.         Testing::_assert(self::cmp(  0, 0,  -16, 0) ===  1);
  1008.  
  1009.         // Also test another even value (looking for edge cases).
  1010.         Testing::_assert(self::cmp( 0, 0,   0, 0) === 0);
  1011.         Testing::_assert(self::cmp( 0, 0,   0, 2)  <  0);
  1012.         Testing::_assert(self::cmp( 0, 2,   0, 0)  >  0);
  1013.         Testing::_assert(self::cmp( 0, 2,   0, 2) === 0);
  1014.         Testing::_assert(self::cmp(-2, 0,  -2, 0) === 0);
  1015.         Testing::_assert(self::cmp(-2, 0,   0, 0)  <  0);
  1016.         Testing::_assert(self::cmp(-2, 0,   2, 0)  <  0);
  1017.         Testing::_assert(self::cmp( 0, 0,  -2, 0)  >  0);
  1018.         Testing::_assert(self::cmp( 0, 0,   0, 0) === 0);
  1019.         Testing::_assert(self::cmp( 0, 0,   2, 0)  <  0);
  1020.         Testing::_assert(self::cmp( 2, 0,  -2, 0)  >  0);
  1021.         Testing::_assert(self::cmp( 2, 0,   0, 0)  >  0);
  1022.         Testing::_assert(self::cmp( 2, 0,   2, 0) === 0);
  1023.  
  1024.         // Notably, we've carefully avoided negative values in the less-significant parts, so far.
  1025.         // That's because the less-significant integers shall be treated as
  1026.         // unsigned integers, but PHP only has _signed_ comparison.
  1027.         // So we better check for mistakes of that kind!
  1028.         $x0000 = (int)0;
  1029.         $xFFFF = (int)-1;
  1030.         /*
  1031.         Testing::_assert(self::cmp($x0000,$x0000,  $x0000,$x0000) === 0); // 0 === 0
  1032.         Testing::_assert(self::cmp($x0000,$x0000,  $x0000,$xFFFF)  <  0); // 0 < (2^32 - 1)
  1033.         Testing::_assert(self::cmp($x0000,$x0000,  $xFFFF,$x0000)  >  0); // 0 > -(2^32)
  1034.         Testing::_assert(self::cmp($x0000,$x0000,  $xFFFF,$xFFFF)  >  0); // 0 > -1
  1035.         Testing::_assert(self::cmp($x0000,$xFFFF,  $x0000,$x0000)  >  0); // (2^32 - 1) > 0
  1036.         Testing::_assert(self::cmp($x0000,$xFFFF,  $x0000,$xFFFF) === 0); // (2^32 - 1) === (2^32 - 1)
  1037.         Testing::_assert(self::cmp($x0000,$xFFFF,  $xFFFF,$x0000)  >  0); // (2^32 - 1) > -(2^32)
  1038.         Testing::_assert(self::cmp($x0000,$xFFFF,  $xFFFF,$xFFFF)  >  0); // (2^32 - 1) > -1
  1039.         Testing::_assert(self::cmp($xFFFF,$x0000,  $x0000,$x0000)  <  0); // -(2^32) < 0
  1040.         Testing::_assert(self::cmp($xFFFF,$x0000,  $x0000,$xFFFF)  <  0); // -(2^32) < (2^32 - 1)
  1041.         Testing::_assert(self::cmp($xFFFF,$x0000,  $xFFFF,$x0000) === 0); // -(2^32) === -(2^32)
  1042.         Testing::_assert(self::cmp($xFFFF,$x0000,  $xFFFF,$xFFFF)  <  0); // -(2^32) < -1
  1043.         Testing::_assert(self::cmp($xFFFF,$xFFFF,  $x0000,$x0000)  <  0); // -1 < 0
  1044.         Testing::_assert(self::cmp($xFFFF,$xFFFF,  $x0000,$xFFFF)  <  0); // -1 < (2^32 - 1)
  1045.         Testing::_assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$x0000)  >  0); // -1 > -(2^32)
  1046.         Testing::_assert(self::cmp($xFFFF,$xFFFF,  $xFFFF,$xFFFF) === 0); // -1 === -1
  1047. */
  1048.         echo(" done.\n");
  1049.     }
  1050.  
  1051.     /*
  1052.     * Multiplies two 64-bit integers, resulting in a single 128-bit integer.
  1053.     *
  1054.     * Because PHP does not have a native 128-bit integer type (or the ability
  1055.     * to define structs that can be returned from functions), the returned
  1056.     * value is placed into two 64-bit integer reference-type parameters.
  1057.     */
  1058.     public static function multiply_64x64(int $a, int $b, int &$out_hi, int &$out_lo) : void
  1059.     {
  1060.         self::multiply_NxN(32, $a, $b, $out_hi,$out_lo);
  1061.     }
  1062.  
  1063.     /*
  1064.     * It is probably better to use `multiply_64x64` in most cases.
  1065.     *
  1066.     * This is the underlying implementation of `multiply_64x64`, and having
  1067.     * it separated out as a different function is useful for testing purposes.
  1068.     */
  1069.     public static function multiply_NxN(int $digit_bit_width,  int $a, int $b,  int &$out_hi, int &$out_lo) : void
  1070.     {
  1071.         Testing::_assert($digit_bit_width <= 32);
  1072.         $mask_lo = ((1 << $digit_bit_width) - 1);
  1073.         $max_digit = $mask_lo;
  1074.  
  1075.         // Reorganize the two 64-bit integers into four 32-bit integers.
  1076.         // This helps because now the 32-bit integers will have at least
  1077.         // 32 bits of "headroom" for doing multiplications.
  1078.         $a_lo = $a & $mask_lo;
  1079.         $a_hi = $a >> $digit_bit_width;
  1080.         $b_lo = $b & $mask_lo;
  1081.         $b_hi = $b >> $digit_bit_width;
  1082.  
  1083.         // Multiplies. We need 4 of them for a school-style long multiply.
  1084.         // There are faster methods, but this will do for now.
  1085.         //
  1086.         // Graphically:
  1087.         //                a_hi  a_lo
  1088.         //             x  b_hi  b_lo
  1089.         //             -------------
  1090.         //                x_hi  x_lo
  1091.         //          y_hi  y_lo
  1092.         //          z_hi  z_lo
  1093.         // +  w_hi  w_lo
  1094.         // -----------------------
  1095.         //    ????  ????  ????  ????
  1096.         //
  1097.         // Note that these multiplications, or at least the ones involving
  1098.         // the $*_lo operands, must be _unsigned_ multiplications.
  1099.         // However, because these are 32-bit values stored in 64-bit variables,
  1100.         // the sign bit will _never_ be set, so it won't matter if we
  1101.         // use a 64-bit signed multiplication on them.
  1102.         //
  1103.         $x = $a_lo * $b_lo;
  1104.         $y = $a_lo * $b_hi;
  1105.         $z = $a_hi * $b_lo;
  1106.         $w = $a_hi * $b_hi;
  1107.         // TODO: BUG? If the above multiplies cause a signed overflow, PHP will probably convert the result to a `float`!
  1108.  
  1109.         // To make the logic more clear, we will also define the
  1110.         // variables corresponding to the {x,y,z,w}_{hi,lo} placeholders.
  1111.         $x_lo = $x & $mask_lo;
  1112.         $y_lo = $y & $mask_lo;
  1113.         $z_lo = $z & $mask_lo;
  1114.         $w_lo = $w & $mask_lo;
  1115.  
  1116.         $x_hi = $x >> $digit_bit_width;
  1117.         $y_hi = $y >> $digit_bit_width;
  1118.         $z_hi = $z >> $digit_bit_width;
  1119.         $w_hi = $w >> $digit_bit_width;
  1120.  
  1121.         // PHP doesn't provide an unsigned right-shift, so we must clear
  1122.         // any sign-extended bits in things that were right-shifted.
  1123.         $x_hi = $x & $mask_lo;
  1124.         $y_hi = $y & $mask_lo;
  1125.         $z_hi = $z & $mask_lo;
  1126.         $w_hi = $w & $mask_lo;
  1127.  
  1128.         // Calculate the bottom row of 32-bit integers.
  1129.         //
  1130.         // Note that some of these additions might "overflow", but
  1131.         // because we are storing our 32-bit integers in 64-bit variables,
  1132.         // the carry values will be captured.
  1133.         //
  1134.         // These will get shuffled into the output integers
  1135.         // once the accumulation is done.
  1136.         //
  1137.         // Graphically:
  1138.         //                a_hi  a_lo
  1139.         //             x  b_hi  b_lo
  1140.         //             -------------
  1141.         //                x_hi  x_lo
  1142.         //          y_hi  y_lo
  1143.         //          z_hi  z_lo
  1144.         // +  w_hi  w_lo
  1145.         // -----------------------
  1146.         //    out3  out2  out1  out0
  1147.         //
  1148.         $out0 = $x_lo;
  1149.         $out1 = $x_hi + $y_lo + $z_lo;
  1150.         $out2 = $y_hi + $z_hi + $w_lo;
  1151.         $out3 = $w_hi;
  1152.  
  1153.         // Perform carry operations.
  1154.         // (Beware: order-of-operations is important.)
  1155.         if ( $out1 > $max_digit )
  1156.         {
  1157.             $out2 += ($out1 >> $digit_bit_width);
  1158.             $out1 &= $mask_lo;
  1159.         }
  1160.  
  1161.         if ( $out2 > $max_digit )
  1162.         {
  1163.             $out3 += ($out2 >> $digit_bit_width);
  1164.             $out2 &= $mask_lo;
  1165.         }
  1166.         // Note that $out3 is involved in an addition, but won't generate
  1167.         // a carry-out. It won't overflow, for math reasons.
  1168.  
  1169.         // Now we have four proper 32-bit integers (two pairs) with no carry bits,
  1170.         // so we can concatenate the pairs to get two 64-bit integers and have our 128-bit result.
  1171.         $lo = $out0 | ($out1 << $digit_bit_width);
  1172.         $hi = $out2 | ($out3 << $digit_bit_width);
  1173.  
  1174.         // Return our results from the function.
  1175.         $out_lo = $lo;
  1176.         $out_hi = $hi;
  1177.     }
  1178.  
  1179.     // This is not part of the public API because there is no point.
  1180.     // It _should_ do exactly what the `*` operator does.
  1181.     // It's just designed to use the `multiply_64x64` function, so that
  1182.     // we can see if it gives results identical to `*`, at least whenever
  1183.     // the two are given operands that don't overflow.
  1184.     private static function multiply_NxN_lo(int $bit_width, int $a, int $b) : int
  1185.     {
  1186.         $out_lo = (int)0;
  1187.         $out_hi = (int)0;
  1188.         self::multiply_NxN($bit_width, $a, $b,  $out_hi,$out_lo);
  1189.         return $out_lo;
  1190.     }
  1191.     /* TODO: cut?
  1192.     // Uses the `multiply_32H32Lx32H32L` function to emulate what
  1193.     // `multiply_64x64` would do if it had lower bit-width integers available
  1194.     // (ex: `multiply_32x32`).
  1195.     //
  1196.     // Because this is used for testing purposes, and because we've asserted
  1197.     // that we'll always use PHP on 64-bit systems and `int` will always have
  1198.     // a width of at least 64 bits, then the result of `multiply_32x32` and
  1199.     // smaller will always fit in a 64-bit return value. So we don't need
  1200.     // to return the output in reference parameters for this one, which is
  1201.     // very helpful for testing.
  1202.     private static function multiply_NxN_emul(int $bit_width, int $a, int $b) : int
  1203.     {
  1204.         Testing::_assert($bit_width <= 32);
  1205.         $mask_lo = ((1 << $bit_width) - 1);
  1206.         $a_lo = $a & $mask_lo;
  1207.         $a_hi = $a >> $bit_width;
  1208.         $b_lo = $b & $mask_lo;
  1209.         $b_hi = $b >> $bit_width;
  1210.         $out_lo = (int)0;
  1211.         $out_hi = (int)0;
  1212.         self::multiply_32H32Lx32H32L($a_hi,$a_lo,  $b_hi,$b_lo,  $out_hi,$out_lo);
  1213.  
  1214.         $int_max = $mask_lo;
  1215.         Testing::_assert($out_hi <= $int_max);
  1216.         Testing::_assert($out_lo <= $int_max);
  1217.  
  1218.  
  1219.     }
  1220. */
  1221.  
  1222.     public static function unittest_multiply_64x64() : void
  1223.     {
  1224.         echo(__CLASS__."::".__FUNCTION__."...");
  1225.  
  1226.         Testing::_assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x00));
  1227.         Testing::_assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0x01));
  1228.         Testing::_assert(0x0000 === self::multiply_NxN_lo(8, 0x01, 0x00));
  1229.         Testing::_assert(0x0000 === self::multiply_NxN_lo(8, 0x00, 0xFF));
  1230.         Testing::_assert(0x0000 === self::multiply_NxN_lo(8, 0xFF, 0x00));
  1231.         Testing::_assert(0x0001 === self::multiply_NxN_lo(8, 0x01, 0x01));
  1232.         Testing::_assert(0x000F === self::multiply_NxN_lo(8, 0x01, 0x0F));
  1233.         Testing::_assert(0x000F === self::multiply_NxN_lo(8, 0x0F, 0x01));
  1234.         Testing::_assert(0x00E1 === self::multiply_NxN_lo(8, 0x0F, 0x0F));
  1235.         Testing::_assert(0x0010 === self::multiply_NxN_lo(8, 0x01, 0x10));
  1236.         Testing::_assert(0x0010 === self::multiply_NxN_lo(8, 0x10, 0x01));
  1237.         Testing::_assert(0x0100 === self::multiply_NxN_lo(8, 0x10, 0x10));
  1238.         Testing::_assert(0x4000 === self::multiply_NxN_lo(8, 0x80, 0x80));
  1239.         Testing::_assert(0x3F01 === self::multiply_NxN_lo(8, 0x7F, 0x7F));
  1240.         Testing::_assert(0xFC04 === self::multiply_NxN_lo(8, 0xFE, 0xFE));
  1241.         Testing::_assert(0xFD02 === self::multiply_NxN_lo(8, 0xFE, 0xFF));
  1242.         Testing::_assert(0xFD02 === self::multiply_NxN_lo(8, 0xFF, 0xFE));
  1243.         Testing::_assert(0xFE01 === self::multiply_NxN_lo(8, 0xFF, 0xFF));
  1244.         Testing::_assert(0x7E81 === self::multiply_NxN_lo(8, 0x7F, 0xFF));
  1245.         Testing::_assert(0x7F80 === self::multiply_NxN_lo(8, 0x80, 0xFF));
  1246.         Testing::_assert(0x3F80 === self::multiply_NxN_lo(8, 0x80, 0x7F));
  1247.  
  1248.         for ( $i = (int)0; $i < 256; $i++ )
  1249.         {
  1250.             for ( $j = (int)0; $j < 256; $j++ )
  1251.             {
  1252.                 $expected = (int)($i*$j);
  1253.                 $got = self::multiply_NxN_lo(8, $i, $j);
  1254.                 Testing::_assert($expected === $got, sprintf("Operands: i=%02x; j=%02x; Expected: %04x; Got: %04X", $i, $j, $expected, $got));
  1255.             }
  1256.         }
  1257.  
  1258.         Testing::_assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0000));
  1259.         Testing::_assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0x0001));
  1260.         Testing::_assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0001, 0x0000));
  1261.         Testing::_assert(0x0000_0000 === self::multiply_NxN_lo(16, 0x0000, 0xFFFF));
  1262.         Testing::_assert(0x0000_0000 === self::multiply_NxN_lo(16, 0xFFFF, 0x0000));
  1263.         Testing::_assert(0x0000_0001 === self::multiply_NxN_lo(16, 0x0001, 0x0001));
  1264.         Testing::_assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0x0001, 0xFFFF));
  1265.         Testing::_assert(0x0000_FFFF === self::multiply_NxN_lo(16, 0xFFFF, 0x0001));
  1266.         Testing::_assert(0xFFFE_0001 === self::multiply_NxN_lo(16, 0xFFFF, 0xFFFF));
  1267.         Testing::_assert(0XA518_60E1 === self::multiply_NxN_lo(16, 0xF00F, 0xB00F));
  1268.  
  1269.         Testing::_assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0000));
  1270.         Testing::_assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0x0000_0001));
  1271.         Testing::_assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0000));
  1272.         Testing::_assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0x0000_0000, 0xFFFF_FFFF));
  1273.         Testing::_assert(0x0000_0000_0000_0000 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0000));
  1274.         Testing::_assert(0x0000_0000_0000_0001 === self::multiply_NxN_lo(32, 0x0000_0001, 0x0000_0001));
  1275.         Testing::_assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0x0000_0001, 0xFFFF_FFFF));
  1276.         Testing::_assert(0x0000_0000_FFFF_FFFF === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0x0000_0001));
  1277.         Testing::_assert(0xFFFF_FFFE_0000_0001 === self::multiply_NxN_lo(32, 0xFFFF_FFFF, 0xFFFF_FFFF));
  1278.  
  1279.         // Explicit test of 128-bit returning version, and of 64-bit inputs.
  1280.         $a64 = (int)0xAceBe11e_CafeBabe;
  1281.         $b64 = (int)0xF100fCa7_F1edFa57;
  1282.         $out_hi = (int)0;
  1283.         $out_lo = (int)0;
  1284.         self::multiply_64x64($a64,  $b64,  $out_hi,$out_lo);
  1285.         Testing::_assert(0x8E9C_2945_7ED5_0292 === $out_lo);
  1286.         Testing::_assert(0xA2CA_B997_9FFE_C71C === $out_hi);
  1287.  
  1288.         echo(" done.\n");
  1289.     }
  1290.  
  1291.     /**
  1292.     * @param      int  $numerator_hi
  1293.     * @param      int  $numerator_lo
  1294.     * @param      int  $denominator
  1295.     * @param-out  int  $quotient_hi
  1296.     * @param-out  int  $quotient_lo
  1297.     * @returns    int  The remainder.
  1298.     */
  1299.     public static function divide_128x64(int $numerator_hi, int $numerator_lo, int $denominator, int &$quotient_hi, int &$quotient_lo) : int
  1300.     {
  1301.         // TODO: Test what PHP does with integer divide-by-zero.
  1302.         // In principle, PHP should throw a DivisionByZeroError whenever this
  1303.         // happens. In that case, we can just proceed as normal, because
  1304.         // the first time division is used in this function, it'll cause
  1305.         // that error to be thrown.
  1306.         // See: https://www.php.net/manual/en/class.divisionbyzeroerror.php
  1307.         // (But I haven't tested to make _sure_ PHP does this. If it doesn't
  1308.         // throw an exception, then we should. Propagating weird values
  1309.         // through the calculation could result in hard-to-find bugs.)
  1310.  
  1311.         // Special-case the value denominator value of `1`, both to make this
  1312.         // identity operation be fast, but also to ensure that we don't have
  1313.         // to worry about any weird corner cases for denominators that are
  1314.         // less than 2. (It can matter, because binary arithmatic.)
  1315.         if ( $denominator === 1 ) {
  1316.             $quotient_hi = $numerator_hi;
  1317.             $quotient_lo = $numerator_lo;
  1318.             return 0;
  1319.         }
  1320.  
  1321.         // Give these predictable values.
  1322.         // This will be helpful in debugging if an error happens:
  1323.         // you'll always know what these started as!
  1324.         $quotient_hi = (int)0;
  1325.         $quotient_lo = (int)0;
  1326.  
  1327.         // Store the sign bit that the result will have.
  1328.         $sign = ($numerator_hi ^ $denominator) & self::MASK_SIGN;
  1329.  
  1330.         // Sign extraction, so it doesn't mess with our division.
  1331.         $sign = $hi & self::MASK_SIGN;
  1332.         $hi &= ~self::MASK_SIGN;
  1333.  
  1334.         // Extract remainders.
  1335.         //
  1336.         // Right now we are mostly interested in $hi_remainder.
  1337.         // This will be necessary later for estimating the low portion of the result.
  1338.         //
  1339.         // The $lo_remainder will simply be returned by the function as-is,
  1340.         // since this may help the caller with rounding logic.
  1341.         $lo_remainder = $lo % $denominator;
  1342.         $hi_remainder = $hi % $denominator;
  1343.  
  1344.         // Component-wise division.
  1345.         //
  1346.         // (We use `intdiv` because `/` might do a check to see if one or both
  1347.         // of its operands are floating point numbers, because it is a floating
  1348.         // point operation by default. `intdiv`, on the other hand, is explicitly
  1349.         // intended for integer division, so it may perform/behave better for this use.)
  1350.         $lo = intdiv($lo, $denominator);
  1351.         $hi = intdiv($hi, $denominator);
  1352.  
  1353.         // Calculate the carry.
  1354.         $min_carry = PHP_INT_MAX; // Basically ((2^64 / 2) - 1). We'd prefer (2^64), but can't represent it.
  1355.         $min_carry = intdiv($min_carry, $denominator);
  1356.         $min_carry <<= 1; // Undo the earlier divide-by-2 in the ((2^64 / 2) - 1) = PHP_INT_MAX.
  1357.  
  1358.         // Perform the carry.
  1359.         $lo += ($hi_remainder * $min_carry);
  1360.  
  1361.         // For safety reasons, we'll also calculate a "max" carry.
  1362.         $max_carry_init = (1 << 62);
  1363.         $max_carry = intdiv($max_carry_init, $denominator);
  1364.         if ( ($denominator * $max_carry) !== $max_carry_init ) {
  1365.             // Always round up.
  1366.             $max_carry++;
  1367.         }
  1368.         $max_carry <<= 2;
  1369.         $max_carry += 3; // Always round up.
  1370.  
  1371.         // Perform the carry.
  1372.         $max_lo += ($hi_remainder * $max_carry);
  1373.  
  1374.         // This loop takes our approximation and improves it until we have
  1375.         // the correct quotient.
  1376.         $test_lo = (int)0;
  1377.         $test_hi = (int)0;
  1378.         $test_carry = (int)0;
  1379.  
  1380.         // TODO: BUG? Inspect the below `+=` operators and the IntNN:add_with_carry
  1381.         // to ensure that this function will be handling overflow conditions
  1382.         // correctly. Right now, it probably _isn't_ (though we could get lucky?).
  1383.  
  1384.         self::multiply_64x64($denominator,  $lo,  $test_hi,$test_lo);
  1385.         $test_hi += ($hi * $denominator);
  1386.  
  1387.         while ( true )
  1388.         {
  1389.             $test_lo = IntNN::add_with_carry($test_lo, $denominator, $test_carry);
  1390.             $test_hi += $test_carry;
  1391.             // if ( test > numerator )
  1392.             if ( self::cmp($test_hi, $test_lo, $numerator_hi, $numerator_lo) > 0 ) {
  1393.                 break;
  1394.             }
  1395.  
  1396.             // The additional increment of $denominator in the $test variable
  1397.             // corresponds to incrementing the $lo value by 1.
  1398.             $lo++;
  1399.  
  1400.             // This handles the aforementioned "safety reasons".
  1401.             // It prevents us from infinite-looping, or from trying to
  1402.             // iterate most of a 64-bit integer's possible values
  1403.             // when it is already futile to get the correct answer.
  1404.             // Ideally, this will never be executed.
  1405.             // If the surrounding function works as intended, then the loop
  1406.             // will ALWAYS converge before this condition becomes true.
  1407.             if ($lo > $max_lo) {
  1408.                 $class_name = __CLASS__;
  1409.                 $func_name = __FUNCTION__;
  1410.                 throw new \Exception(sprintf(
  1411.                     "Internal error: `$class_name::$func_name` did not converge when it should have. ".
  1412.                     "Aborting to prevent excessive looping and CPU drain. ".
  1413.                     "This means that the `$func_name` function is broken and may return incorrect results, ".
  1414.                     "even when this error isn't thrown. Please report this error. ".
  1415.                     "Parameters: {numerator_hi:%016X, numerator_lo:%016X, denominator:%016X}",
  1416.                     $numerator_hi, $numerator_lo, $denominator
  1417.                     ));
  1418.             }
  1419.         }
  1420.  
  1421.         // Pass results to the caller.
  1422.         $quotient_hi = $hi;
  1423.         $quotient_lo = $lo;
  1424.         return $lo_remainder;
  1425.  
  1426.         /*
  1427.         TODO: Comment is kinda sorta relevant but also outdated already.
  1428.  
  1429.         // This will be multiplied by the $hi_remainder to get the portion
  1430.         // of the number that is being carried into the $lo component of the result.
  1431.         //
  1432.         // Note that this should be `(2^64/$denominator)` EXACTLY.
  1433.         // Unfortunately, we can't be exact for two reasons:
  1434.         // * Rational numbers (the exact result) can't be represented by finite integers.
  1435.         // * The number 2^64 is out of the range of 64-bit integers (one more than the max!)
  1436.         //
  1437.         // Now we'll get clever, and our function will get slower.
  1438.         //
  1439.         // We notice that the multiplier doesn't _need_ to be exact.
  1440.         // We can make a guess. And the guess can be wrong. As long as it's
  1441.         // always wrong in a way that's a little bit low.
  1442.         //
  1443.         // If it's low, it will result in our end-result being too low.
  1444.         // Then we can multiply our guess by the denominator.
  1445.         // The result will be less than the numerator.
  1446.         // (If we had exact numbers instead of integers, multiplying the quotient
  1447.         // by the denominator would result in _exactly_ the numerator.)
  1448.         //
  1449.         // Then we can just add the $denominator over and over, checking
  1450.         // it with a multiplication, until it is too high.
  1451.         // Once it is too high, we return the last number that multiplied
  1452.         // to a value
  1453.         //
  1454.         */
  1455.     }
  1456.  
  1457.     public static function unittest_divide_128x64()
  1458.     {
  1459.         echo(__CLASS__."::".__FUNCTION__."...");
  1460.         // TODO: Oh god this needs testing BADLY.
  1461.         echo(" UNIMPLEMENTED! (Please fix!)\n");
  1462.         // echo(" done.\n");
  1463.     }
  1464.  
  1465.     public static function unittest()
  1466.     {
  1467.         echo("Unittests for ".__CLASS__."\n");
  1468.         self::unittest_cmp();
  1469.         self::unittest_multiply_64x64();
  1470.         self::unittest_divide_128x64();
  1471.         echo("\n");
  1472.     }
  1473.  
  1474.     // Prevent instantiation/construction of the (static/constant) class.
  1475.     /** @return never */
  1476.     private function __construct() {
  1477.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  1478.     }
  1479. }
  1480.  
  1481. /*
  1482. ?>
  1483.  
  1484. <?php
  1485. declare(strict_types=1);
  1486.  
  1487. namespace Kickback\Math;
  1488. */
  1489. // TODO: I forget if it's better to do this, or to use an `enum` type. Whatevs.
  1490. final class RoundingMode
  1491. {
  1492.     // Rounding modes.
  1493.     // (I think this is all of the possibilities? TODO: Verify.)
  1494.     public const ALL_TOWARDS_NEG_INF  =  0; // Default; mimics x86 integer truncation. (I think? TODO: Double-check this.)
  1495.     public const ALL_TOWARDS_POS_INF  =  1;
  1496.     public const ALL_TOWARDS_ZERO     =  3;
  1497.     public const ALL_AWAY_FROM_ZERO   =  4;
  1498.     public const ALL_TOWARDS_EVEN     =  5;
  1499.     public const ALL_TOWARDS_ODD      =  6;
  1500.     public const HALF_TOWARDS_NEG_INF =  7;
  1501.     public const HALF_TOWARDS_POS_INF =  8;
  1502.     public const HALF_TOWARDS_ZERO    =  9; // Similar to PHP_ROUND_HALF_DOWN
  1503.     public const HALF_AWAY_FROM_ZERO  = 10; // Similar to PHP_ROUND_HALF_UP
  1504.     public const HALF_TOWARDS_EVEN    = 11; // Similar to PHP_ROUND_HALF_EVEN; bankers' rounding.
  1505.     public const HALF_TOWARDS_ODD     = 12; // Similar to PHP_ROUND_HALF_ODD
  1506.  
  1507.     // This shall mimic x86 integer truncation.
  1508.     public const DEFAULT = self::ALL_TOWARDS_NEG_INF;
  1509.  
  1510.     public static function rounding_mode_to_string(int $rounding_mode) : string
  1511.     {
  1512.         switch($rounding_mode)
  1513.         {
  1514.             case self::ALL_TOWARDS_NEG_INF : return "ALL_TOWARDS_NEG_INF";
  1515.             case self::ALL_TOWARDS_POS_INF : return "ALL_TOWARDS_POS_INF";
  1516.             case self::ALL_TOWARDS_ZERO    : return "ALL_TOWARDS_ZERO";
  1517.             case self::ALL_AWAY_FROM_ZERO  : return "ALL_AWAY_FROM_ZERO";
  1518.             case self::ALL_TOWARDS_EVEN    : return "ALL_TOWARDS_EVEN";
  1519.             case self::ALL_TOWARDS_ODD     : return "ALL_TOWARDS_ODD";
  1520.             case self::HALF_TOWARDS_NEG_INF: return "HALF_TOWARDS_NEG_INF";
  1521.             case self::HALF_TOWARDS_POS_INF: return "HALF_TOWARDS_POS_INF";
  1522.             case self::HALF_TOWARDS_ZERO   : return "HALF_TOWARDS_ZERO";
  1523.             case self::HALF_AWAY_FROM_ZERO : return "HALF_AWAY_FROM_ZERO";
  1524.             case self::HALF_TOWARDS_EVEN   : return "HALF_TOWARDS_EVEN";
  1525.             case self::HALF_TOWARDS_ODD    : return "HALF_TOWARDS_ODD";
  1526.         }
  1527.  
  1528.         return "Unknown rounding mode ($rounding_mode)";
  1529.     }
  1530.  
  1531.     // Prevent instantiation/construction of the (static/constant) class.
  1532.     /** @return never */
  1533.     private function __construct() {
  1534.         throw new \Exception("Instantiation of static class " . get_class($this) . "is not allowed.");
  1535.     }
  1536. }/*
  1537. ?>
  1538.  
  1539. <?php
  1540. declare(strict_types=1);
  1541.  
  1542. namespace Kickback\Math\Types;
  1543.  
  1544. use \Kickback\Math\Int128;
  1545. use \Kickback\Math\RoundingMode;
  1546. */
  1547. final class FixedPoint64
  1548. {
  1549.     private int $numerator;
  1550.     private int $divisor;
  1551.  
  1552.     private static function banker_round(int $numerator, int $increment) : int
  1553.     {
  1554.         Testing::_assert($increment > 1);
  1555.  
  1556.         // Even if we don't use $even_increment in most cases,
  1557.         // we should still assert that it fits in _all_ cases, just to give
  1558.         // this function consistent behavior (e.g. whether it asserts or not
  1559.         // should depend on `$increment` and not on `$numerator`.)
  1560.         //$even_increment = ($increment<<1);
  1561.         Testing::_assert(($increment<<1) <= PHP_INT_MAX);
  1562.  
  1563.         if ( $numerator >= 0 ) {
  1564.             return self::banker_round_unsigned_($numerator, $increment);
  1565.         }
  1566.         else {
  1567.             return self::banker_round_signed_($numerator, $increment);
  1568.         }
  1569.     }
  1570.  
  1571.     private static function banker_round_signed_(int $numerator, int $increment) : never
  1572.     {
  1573.         // TODO!
  1574.         Testing::_assert(false, "Signed rounding is not yet implemented.");
  1575.     }
  1576.  
  1577.     private static function banker_round_unsigned_(int $numerator, int $increment) : int
  1578.     {
  1579.         /* TODO: cut?
  1580.         // We calculate the $round_down_amount (conventional rounding)
  1581.         // from the $round_down_even_amount (bankers' rounding)
  1582.         // so that we avoid having to modulus more than once.
  1583.         $round_down_even_amount = $numerator % $even_increment;
  1584.         $round_down_amount = $round_down_even_amount;
  1585.         if ( $round_down_amount > $increment ) {
  1586.             $round_down_amount -= $increment;
  1587.         }
  1588.         */
  1589.  
  1590.         // First, attempt conventional rounding (e.g. "round-to-nearest-increment").
  1591.         $rounding_amount = (int)0;
  1592.         if ( self::get_rounding_amount_($numerator, $increment, $rounding_amount) ) {
  1593.             return $rounding_amount;
  1594.         } else {
  1595.             // If that didn't work, then fall back on the tie breaker.
  1596.             return self::banker_round_tie_breaker_($numerator, $increment);
  1597.         }
  1598.     }
  1599.  
  1600.     private static function get_rounding_amount_(int $numerator, int $increment, &$rounding_amount) : bool
  1601.     {
  1602.         $round_down_amount = $numerator % $increment;
  1603.         $round_up_amount = $increment - $round_down_amount;
  1604.         if ( $round_down_amount < $round_up_amount ) {
  1605.             $rounding_amount = -$round_down_amount;
  1606.             return true;
  1607.         } else
  1608.         if ( $round_down_amount > $round_up_amount ) {
  1609.             $rounding_amount = $round_up_amount;
  1610.             return true;
  1611.         } else {
  1612.             // If that didn't work, then there's a problem, and it's probably this tie:
  1613.             Testing::_assert($round_down_amount === $round_up_amount);
  1614.             $rounding_amount = 0; // TODO: Is this right?
  1615.             return false;
  1616.         }
  1617.     }
  1618.  
  1619.     // This is set out into its own function because the tie breaker
  1620.     // is unlikely to be executed very frequently. As such, it might
  1621.     // be faster (code-execution-wise) to break this out into it's own function.
  1622.     //
  1623.     // Reasons why this might be the case:
  1624.     // * Interpreter is less likely to need to parse/analyse this code when `banker_round` is called.
  1625.     // * This is more cache-friendly: the CPU won't need to load the instructions for the tie breaker every time `banker_round` is called.
  1626.     // * This makes `banker_round` smaller and therefore more inline-able.
  1627.     //
  1628.     // I am not sure how much any of those _actually_ matter, and it's not
  1629.     // really worth testing at the moment, but this is _usually_ a good way
  1630.     // to optimize code, and it seems to have few or no disadvantages.
  1631.     private static function banker_round_tie_breaker_(int $numerator, int $increment) : int
  1632.     {
  1633.         // Now the bankers' rounding comes into play.
  1634.         // We break the tie by rounding to the nearest _even_ increment.
  1635.         $even_increment = $increment << 1;
  1636.  
  1637.         // This involves just doing the rounding again, but with $increment*2.
  1638.         $rounding_amount = (int)0;
  1639.         Testing::_assert(self::get_rounding_amount_($numerator, $even_increment, $rounding_amount));
  1640.         return $rounding_amount;
  1641.     }
  1642.     // TODO: The above bankers' rounding code is pretty sus. I was very sleepy while writing it. o.O
  1643.  
  1644.     private static function unittest_banker_round() : void
  1645.     {
  1646.         echo(__CLASS__."::".__FUNCTION__."...");
  1647.         // TODO: This needs testing if you want to trust any financials done with it.
  1648.         echo(" UNIMPLEMENTED! (Please fix!)\n");
  1649.         // echo(" done.\n");
  1650.     }
  1651.  
  1652.  
  1653.  
  1654.     /*
  1655.     * Beware: This operation is NOT commutative like normal multiplication.
  1656.     * That's because the choice of destination decides which divisor to use
  1657.     * when scaling back the oversized numerator that results from the
  1658.     * initial multiplication.
  1659.     *
  1660.     * The `multiply` and `multiply_into` functions have the commutative property.
  1661.     *
  1662.     * This method will not perform any explicit memory allocations
  1663.     * unless an error has occurred.
  1664.     */
  1665.     public function multiply_by(FixedPoint64 $other, int $rounding_mode = RoundingMode::DEFAULT) : void
  1666.     {
  1667.         // ASSUMPTION: We wish the output to have the same divisor as $this, not $other (or anything else).
  1668.         $this->numerator = self::multiply_raw($this, $other, $this->denominator, $rounding_mode);
  1669.     }
  1670.  
  1671.     /*
  1672.     * Allocates a new FixedPoint64 object, then uses the `multiply_into`
  1673.     * function to multiply parameters `$a` and `$b` and store the result
  1674.     * into the new FixedPoint64 object, which is then returned.
  1675.     *
  1676.     * The new FixedPoint64 object will have a $denominator field equal to
  1677.     * the `$divisor` parameter.
  1678.     *
  1679.     * This function could be handy as a quick-and-dirty way to draft code
  1680.     * that uses FixedPoint64 arithmetic. It has the disadvantage of doing
  1681.     * memory allocation for an operation where it isn't always necessary.
  1682.     * For production code, it is recommended to use `multiply_into` or
  1683.     * `multiply_by` instead.
  1684.     */
  1685.     public static function multiply(FixedPoint64 $a, FixedPoint64 $b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : FixedPoint64
  1686.     {
  1687.         $result = new FixedPoint64(0, $divisor);
  1688.         return self::multiply_into($a, $b, $result, $rounding_mode);
  1689.     }
  1690.  
  1691.     /*
  1692.     * Multiplies `$fp64a` and `$fp64b` together, then stores the result
  1693.     * into the `$destination` object.
  1694.     *
  1695.     * As a convenience, the `$destination` object is returned.
  1696.     *
  1697.     * When rounding or truncating the fractional multiplication results,
  1698.     * the `$destination` parameter's `$denominator` field is used
  1699.     * as a divisor.
  1700.     *
  1701.     * This operation has commutative property over `$fp64a` and `$fp64b`.
  1702.     *
  1703.     * This function will not perform any explicit memory allocations
  1704.     * unless an error has occurred. (In some cases it is possible for
  1705.     * the caller to reuse existing FixedPoint64 objects for the `$destination`
  1706.     * parameter, and thus avoid memory allocations entirely.)
  1707.     */
  1708.     public static function multiply_into(FixedPoint64 $fp64a, FixedPoint64 $fp64b, FixedPoint64 $destination, int $rounding_mode = RoundingMode::DEFAULT) : FixedPoint64
  1709.     {
  1710.         Testing::_assert(isset($destination));
  1711.         $tmp = FixedPoint64::multiply_raw($fp64a, $fp64b, $destination->denominator, $rounding_mode);
  1712.         $destination->numerator = $tmp;
  1713.         return $destination;
  1714.     }
  1715.  
  1716.     private static function multiply_raw(FixedPoint64 $fp64a, FixedPoint64 $fp64b, int $divisor, int $rounding_mode = RoundingMode::DEFAULT) : int
  1717.     {
  1718.         // TODO: Verify that divisors are being chosen correctly.
  1719.         // (Actually the below seems to be pretty well argued, now that it's
  1720.         // all been written down. Still, we should make sure it all gets tested.)
  1721.         // TODO: Make sure unittests cover all of the branches in this function.
  1722.  
  1723.         // ---=== Design Notes ===---
  1724.         // I am currently reasoning like so:
  1725.         //
  1726.         // Suppose we multiply two numbers, `a` and `b`.
  1727.         // `a` has a denominator of 2^16, giving it 16 bits of precision.
  1728.         // `b` has a denominator of 2^8, giving it 8 bits of precision.
  1729.         //
  1730.         // The product `a * b` will then have 24 bits of precision, with a denominator of 2^24.
  1731.         //
  1732.         // If we want to store the product into `a'`, which has the same precision as `a`,
  1733.         // then we will need to divide the product by 2^8 (the denominator of `b`)
  1734.         // to acquire a value with 16 bits of precision (because `16 = 24 - 8`).
  1735.         //
  1736.         // If we want to store the product into `b'`, which has the same precision as `b`,
  1737.         // then we will need to divide the product by 2^16 (the denominator of `a`)
  1738.         // to acquire a value with 8 bits of precision (because `8 = 24 - 16`).
  1739.         //
  1740.         // If we want to store the product into `c`,  which has N bits of precision,
  1741.         // then we will need to divide by the number of bits required to reach
  1742.         // that: `24 - x = N`, so `x = 24 - N`.
  1743.         // We divide by `2^(24 - N)` to reach the precision for `c`.
  1744.         //
  1745.         // But wait, what if `N>24`? Well, that's a good thing to notice,
  1746.         // because the destination may actually have more precision than
  1747.         // is given by the product of the multiplicands' denominators.
  1748.         // In that case, we will need to multiply instead of divide!
  1749.         //
  1750.         // Now, what if `a` has `P` bits of precision, and `b` has `Q` bits of precision?
  1751.         // We will need to adjust our formula slightly:
  1752.         // `x = (P+Q) - N`.
  1753.         //
  1754.         // Now, what if `a` has an arbitrary denominator `A` and `b`, likewise has `B`?
  1755.         // The full (wide) product will have a denominator of `AB = A * B`.
  1756.         //
  1757.         // To store into `a'`, we'll need to find `A` in terms of `AB` and `B`.
  1758.         // We write `A` like so: `A = (A * B) / B`.
  1759.         // Thus, `A = AB / B` (by substitution.).
  1760.         //
  1761.         // To store into `b'`, we'll need to find `B` in terms of `AB` and `A`.
  1762.         // We write `B` like so: `B = (A * B) / A`.
  1763.         // Thus, `B = AB / A` (by substitution.).
  1764.         //
  1765.         // To store into `c`, we'll write `C` in terms of `AB` and `X`, then find `X`.
  1766.         // `C = AB/X`
  1767.         // `C*X = AB`
  1768.         // `X = AB/C`
  1769.         // In the code, we will already know `C` because it is the denominator of `c`.
  1770.         // We will need to find `X` so that we can call `divide_128x64` with it.
  1771.         // Now we know how to do that. (I hope.)
  1772.         //
  1773.         // Another algebraic expression for how `c`'s numerator (`cn`) would be calculated:
  1774.         //   cn/cd = (an/ad) * (bn/bd)
  1775.         //   cn = (an/ad) * (bn/bd) * cd
  1776.         //   cn = ((an*bn) / (ad*bd)) * cd
  1777.         //   cn = (an * bn * cd) / (ad * bd)
  1778.         //   cn = (an * bn) * (cd / (ad * bd)) <- if( cd < (ad*bd) ) { do this? }
  1779.         //   cn = (an * bn) / ((ad * bd) / cd) <- if( (ad*bd) < cd ) { do this? }
  1780.  
  1781.         Testing::_assert($divisor > 0);
  1782.         Testing::_assert($divisor <= PHP_INT_MAX);
  1783.  
  1784.         $a = $fp64a->numerator;
  1785.         $b = $fp64b->numerator;
  1786.         $a_div = $fp64a->denominator;
  1787.         $b_div = $fp64b->denominator;
  1788.  
  1789.         // ---=== Multiplication ===---
  1790.         // Multiply 64bit * 64bit to yield 128bit value.
  1791.         $lo = (int)0;
  1792.         $hi = (int)0;
  1793.         Int128::multiply_64x64($a,  $b,  $hi,$lo);
  1794.  
  1795.         // Do the same for the denominators/divisor.
  1796.         $div_lo = (int)0;
  1797.         $div_hi = (int)0;
  1798.         Int128::multiply_64x64($a_div,  $b_div,  $div_hi,$div_lo);
  1799.  
  1800.         // ---=== Division: Denominator ===---
  1801.         // We need to figure out what number to divide or new numerator by,
  1802.         // so that it ends up with the precision described by `$divisor`.
  1803.         // This will be `($a->denominator * $b->denominator) / $divisor`.
  1804.         $scale_direction = (int)0;
  1805.         $div_quotient_lo = (int)0;
  1806.         $div_quotient_hi = (int)0;
  1807.         $div_remainder   = (int)0;
  1808.         if ( $divisor === 0 || $divisor === 1 ) {
  1809.             Testing::_assert($scale_direction === 0);
  1810.             $div_quotient_lo = $div_lo;
  1811.             $div_quotient_hi = $div_hi;
  1812.             Testing::_assert($div_remainder === 0);
  1813.         } else
  1814.         if ( 0 === $div_hi && $divisor > $div_lo ){
  1815.             $scale_direction = 1;
  1816.             // TODO: Division here is not guaranteed to be twos-complement compatible.
  1817.             $div_quotient_lo = intdiv($divisor, $div_lo);
  1818.             Testing::_assert($div_quotient_hi === 0);
  1819.             $div_remainder = $divisor % $div_lo;
  1820.             // TODO: What to do with $div_remainder?
  1821.         }
  1822.         else {
  1823.             $scale_direction = -1;
  1824.             $div_remainder = Int128::divide_128x64($div_hi,$div_lo,  $divisor,  $div_quotient_hi,$div_quotient_lo);
  1825.  
  1826.             // TODO: This limitation isn't strictly necessary;
  1827.             // this is something that CAN happen with valid inputs,
  1828.             // and there is going to be a way to handle it.
  1829.             // It just seems kinda unlikely, and I don't have the time to write it right now.
  1830.             if ( $div_quotient_hi !== 0 ) {
  1831.                 // TODO: Better error message; put parameters and details about the denominators.
  1832.                 $class_name = __CLASS__;
  1833.                 $func_name = __FUNCTION__;
  1834.                 throw new \Exception(
  1835.                     "$class_name::$func_name: Unimplemented combination of input parameters; ".
  1836.                     "The product of the inputs' denominators divided by the destination's denominator ".
  1837.                     "must be a number that fits within PHP_INT_MAX.");
  1838.             }
  1839.             // TODO: What to do with $div_remainder?
  1840.         }
  1841.         // TODO: $div_remainder is a smell suggesting that we haven't done enough math or something.
  1842.         // In particular, this codepath is making it clear that handling arbitrary
  1843.         // input denominators leads to situations where there is no common factor
  1844.         // to divide by.
  1845.         //
  1846.         // Like, if `A = 2^16` and `B = 2^8`, then either `C = 2^16` or `C = 2^8`
  1847.         // will both work fine because `AB = 2^24` which has even factors of both 2^16 and 2^8.
  1848.         // `C` could be any power-of-2 in that situation, though powers greater than 23
  1849.         // will either cause us to do zero scaling or scale up (multiply the result and zero-fill the LSBs).
  1850.         //
  1851.         // However, if `A = 3` and `B = 5`, then `C` has to be some multiple
  1852.         // of either 3 or 5 in order for the scaling to happen cleanly.
  1853.         // If we plug in `C = 2`, then we get `X = 15/2`, which is clearly not an integer.
  1854.         // (It might be possible to get around this by using those remainders
  1855.         // in the later scaling expressions, but it is clear how in my current sleepy state.)
  1856.  
  1857.         // ---=== Rounding ===---
  1858.         // Round the 128bit value, modulo the `$divisor` parameter.
  1859.         // (We really only need to do this for the lower 64 bits, because $divisor can't be bigger than that.)
  1860.         // TODO: Implement more rounding modes? (Also `banker_round` is pretty sus at the moment b/c sleep not enough.)
  1861.         if ( $scale_direction < 0 )
  1862.         {
  1863.             if ( $rounding_mode == RoundingMode::HALF_TOWARDS_EVEN ) {
  1864.                 $lo = self::banker_round($lo, $divisor);
  1865.             } else
  1866.             if ( $rounding_mode != RoundingMode::DEFAULT ) {
  1867.                 throw new \Exception("Unimplemented rounding mode ".RoundingMode::to_string($rounding_mode));
  1868.             }
  1869.         }
  1870.  
  1871.         // ---=== Division (or multiplication): Numerator ===---
  1872.         // Shrink the 128bit value until it is at the desired precision according to `$divisor`.
  1873.         // This will ready it for overflow checking.
  1874.         $quotient_lo = (int)0;
  1875.         $quotient_hi = (int)0;
  1876.         $remainder   = (int)0;
  1877.         if ( $scale_direction === 0 ) {
  1878.             $quotient_lo = $lo;
  1879.             $quotient_hi = $hi;
  1880.             Testing::_assert($remainder === 0);
  1881.         } else
  1882.         if ( $scale_direction > 0 ) {
  1883.             // Least likely case of all of them, but also the most complicated.
  1884.             $tmp1_out_lo = (int)0;
  1885.             $tmp1_out_hi = (int)0;
  1886.             $tmp2_out_lo = (int)0;
  1887.             $tmp2_out_hi = (int)0;
  1888.             Int128::multiply_64x64($lo,  $div_quotient_lo,  $tmp1_out_hi,$tmp1_out_lo);
  1889.             Int128::multiply_64x64($hi,  $div_quotient_lo,  $tmp2_out_hi,$tmp2_out_lo);
  1890.             $quotient_lo = $tmp1_out_lo;
  1891.             $quotient_hi = $tmp1_out_hi + $tmp2_out_lo;
  1892.             Testing::_assert($tmp2_out_hi === 0);
  1893.         } else
  1894.         if ( $scale_direction < 0 ) {
  1895.             $remainder = Int128::divide_128x64($hi,$lo,  $divisor,  $quotient_hi,$quotient_lo);
  1896.         }
  1897.  
  1898.         // ---=== Overflow/Error Handling ===---
  1899.         // Now we check for overflow (and maybe do some validation on rounding logic).
  1900.         // If there is no overflow, then we can safely discard the high part of the quotient.
  1901.         if ( $rounding_mode != RoundingMode::DEFAULT ) {
  1902.             Testing::_assert($remainder === 0); // Because rounding should have handled this already.
  1903.         }
  1904.  
  1905.         if ( 0 !== $quotient_hi ) {
  1906.             $class_name = __CLASS__;
  1907.             $func_name = __FUNCTION__;
  1908.             $rounding_mode_str = RoundingMode::to_string($rounding_mode);
  1909.             throw new \ArithmeticError(sprintf(
  1910.                 "Overflow in `$class_name::$func_name`. ".
  1911.                 "Parameters: {fp64a->numerator:%016X, fp64a->denominator:%016X, fp64b->numerator:%016X, fp64b->denominator:%016x, divisor:%016X, rounding_mode:$rounding_mode_str}",
  1912.                 $fp64a->numerator, $fp64a->denominator, $fp64b->numerator, $fp64b->denominator, $divisor
  1913.                 ));
  1914.         }
  1915.  
  1916.         // ---=== Return ===---
  1917.         return $quotient_lo;
  1918.     }
  1919.  
  1920.     public static function unittest_multiply() : void
  1921.     {
  1922.         // Left as exercise for reader.
  1923.         echo(__CLASS__."::".__FUNCTION__."...");
  1924.         echo(" done.\n");
  1925.     }
  1926.  
  1927.     public static function unittest() : void
  1928.     {
  1929.         echo("Unittests for ".__CLASS__."\n");
  1930.         // Put unittests for other methods here...
  1931.         self::unittest_banker_round();
  1932.         // ... or here ...
  1933.         self::unittest_multiply();
  1934.         // ... or here.
  1935.         echo("\n");
  1936.     }
  1937. }
  1938.  
  1939. ?>
  1940.  
  1941.  
  1942. <?php
  1943. //declare(strict_types=1);
  1944.  
  1945. // Unittest driver to be run from command-line interface.
  1946. /*use \Kickback\Math\IntNN;
  1947. use \Kickback\Math\Int128;
  1948. use \Kickback\Math\Types\FixedPoint64;
  1949. */
  1950.  
  1951. function unittest_all() : void
  1952. {
  1953.     IntNN::unittest();
  1954.     Int128::unittest();
  1955.     FixedPoint64::unittest();
  1956. }
  1957.  
  1958. unittest_all();
  1959. ?>
  1960.  
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!