Multiplies two DDValue numbers using double–double arithmetic. It calculates the high product, uses a fused multiply–add (FMA) to capture rounding error, computes the cross products, and then normalizes the result.
(&mut self, a: DDValue, b: DDValue)
| 895 | /// It calculates the high product, uses a fused multiply–add (FMA) to capture rounding error, |
| 896 | /// computes the cross products, and then normalizes the result. |
| 897 | fn dd_mul(&mut self, a: DDValue, b: DDValue) -> DDValue { |
| 898 | // p = a.hi * b.hi (primary product) |
| 899 | let p = self.builder.ins().fmul(a.hi, b.hi); |
| 900 | // err = fma(a.hi, b.hi, -p) recovers the rounding error. |
| 901 | let zero = self.builder.ins().f64const(0.0); |
| 902 | let neg_p = self.builder.ins().fsub(zero, p); |
| 903 | let err = self.builder.ins().fma(a.hi, b.hi, neg_p); |
| 904 | // Compute cross terms: a.hi*b.lo + a.lo*b.hi. |
| 905 | let a_hi_b_lo = self.builder.ins().fmul(a.hi, b.lo); |
| 906 | let a_lo_b_hi = self.builder.ins().fmul(a.lo, b.hi); |
| 907 | let cross = self.builder.ins().fadd(a_hi_b_lo, a_lo_b_hi); |
| 908 | // Sum p and the cross terms. |
| 909 | let s = self.builder.ins().fadd(p, cross); |
| 910 | // Isolate rounding error from the addition. |
| 911 | let t = self.builder.ins().fsub(s, p); |
| 912 | let s_minus_t = self.builder.ins().fsub(s, t); |
| 913 | let part1 = self.builder.ins().fsub(p, s_minus_t); |
| 914 | let part2 = self.builder.ins().fsub(cross, t); |
| 915 | let e = self.builder.ins().fadd(part1, part2); |
| 916 | // Include the error from the low parts multiplication. |
| 917 | let a_lo_b_lo = self.builder.ins().fmul(a.lo, b.lo); |
| 918 | let err_plus_e = self.builder.ins().fadd(err, e); |
| 919 | let lo_sum = self.builder.ins().fadd(err_plus_e, a_lo_b_lo); |
| 920 | // Renormalize the sum. |
| 921 | let hi_new = self.builder.ins().fadd(s, lo_sum); |
| 922 | let hi_new_minus_s = self.builder.ins().fsub(hi_new, s); |
| 923 | let lo_new = self.builder.ins().fsub(lo_sum, hi_new_minus_s); |
| 924 | DDValue { |
| 925 | hi: hi_new, |
| 926 | lo: lo_new, |
| 927 | } |
| 928 | } |
| 929 | |
| 930 | /// Multiplies a DDValue by a regular f64 (Value) using similar techniques as dd_mul. |
| 931 | /// It multiplies both the high and low parts by b, computes the rounding error, |
no test coverage detected