| | | 1 | | //======================================================================= |
| | | 2 | | // EnclosedQ96Coefficient.cs |
| | | 3 | | //======================================================================= |
| | | 4 | | // MIT License, Copyright (c) 2024-present David Oravsky (mrdav30) |
| | | 5 | | // See LICENSE file in the project root for full license information. |
| | | 6 | | //======================================================================= |
| | | 7 | | |
| | | 8 | | using System.Runtime.CompilerServices; |
| | | 9 | | |
| | | 10 | | namespace FixedMathSharp; |
| | | 11 | | |
| | | 12 | | /// <summary> |
| | | 13 | | /// Retains a signed Q96 coefficient and its absolute error bound in integer |
| | | 14 | | /// Q96 units. Geometry preparation normalizes its sign and widths once; |
| | | 15 | | /// multiplying certifies only unique, representable half-even results. |
| | | 16 | | /// </summary> |
| | | 17 | | internal readonly struct EnclosedQ96Coefficient |
| | | 18 | | { |
| | | 19 | | private readonly Signed192 _magnitude; |
| | | 20 | | private readonly uint _errorUnits; |
| | | 21 | | private readonly sbyte _sign; |
| | | 22 | | private readonly bool _fits, _coarse; |
| | | 23 | | |
| | | 24 | | internal EnclosedQ96Coefficient(Signed192 value, uint errorUnits) |
| | | 25 | | { |
| | 1927 | 26 | | WideArithmetic.GetMagnitude(value, out ulong high, out ulong middle, out ulong low); |
| | 1927 | 27 | | _magnitude = new Signed192(high, middle, low); |
| | 1927 | 28 | | _sign = (sbyte)value.Sign; |
| | 1927 | 29 | | _errorUnits = errorUnits; |
| | 1927 | 30 | | _fits = high == 0; |
| | 1927 | 31 | | _coarse = _fits && (middle >> 36) == 0; |
| | 1927 | 32 | | } |
| | | 33 | | |
| | 1921 | 34 | | internal Signed192 Value => _sign < 0 ? WideArithmetic.Negate(_magnitude) : _magnitude; |
| | | 35 | | |
| | | 36 | | internal bool Equals(EnclosedQ96Coefficient other) => |
| | 6 | 37 | | _sign == other._sign && _errorUnits == other._errorUnits && _magnitude.Equals(other._magnitude); |
| | | 38 | | |
| | | 39 | | /// <summary> |
| | | 40 | | /// Returns raw Q32.32 bits on success; returns false/zero for uncertain |
| | | 41 | | /// rounding, a nonzero product whose coefficient magnitude is at least |
| | | 42 | | /// 2^128, or scalar overflow. |
| | | 43 | | /// The default coefficient is zero, including its zero error bound. |
| | | 44 | | /// </summary> |
| | | 45 | | // Keep scalar construction at materialization, outside inlined certificate outputs. |
| | | 46 | | [MethodImpl(MethodImplOptions.AggressiveInlining)] |
| | | 47 | | internal bool TryMultiplyRaw(Fixed64 value, out long result) |
| | | 48 | | { |
| | 1921 | 49 | | result = default; |
| | 1921 | 50 | | if (_sign == 0 || value.m_rawValue == 0) |
| | 165 | 51 | | return true; |
| | 1756 | 52 | | if (!_fits) |
| | 240 | 53 | | return false; |
| | | 54 | | |
| | 1516 | 55 | | ulong magnitude = Fixed64.AbsToUInt64(value.m_rawValue); |
| | 1516 | 56 | | bool negative = (_sign < 0) != (value.m_rawValue < 0); |
| | 1516 | 57 | | if (_coarse && magnitude < (1UL << 58)) |
| | | 58 | | { |
| | | 59 | | // Q60 needs one product. Discarding 36 bits plus any uint error |
| | | 60 | | // loses <2 Q60 units; these widths also prove scalar range safety. |
| | 799 | 61 | | ulong mantissa = (_magnitude.Middle << 28) | (_magnitude.Low >> 36); |
| | 887 | 62 | | if (mantissa == 0) return true; |
| | 711 | 63 | | Fixed64.Multiply64To128(mantissa, magnitude, out ulong upper, out ulong lower); |
| | 711 | 64 | | ulong fraction = lower & ((1UL << 60) - 1); |
| | | 65 | | const ulong half = 1UL << 59; |
| | 711 | 66 | | ulong distance = fraction < half ? half - fraction : fraction - half; |
| | 711 | 67 | | if (distance > 2 * magnitude) |
| | | 68 | | { |
| | | 69 | | // The certificate excludes half ties and proves magnitude <=2^62. |
| | 437 | 70 | | long raw = (long)((upper << 4) | (lower >> 60)) + (fraction > half ? 1L : 0L); |
| | 437 | 71 | | result = negative ? -raw : raw; |
| | 437 | 72 | | return true; |
| | | 73 | | } |
| | | 74 | | } |
| | 991 | 75 | | return TryMultiplyFull(_magnitude.Middle, _magnitude.Low, magnitude, negative, _errorUnits, out result); |
| | | 76 | | } |
| | | 77 | | |
| | | 78 | | private static bool TryMultiplyFull(ulong coefficientMiddle, ulong coefficientLow, |
| | | 79 | | ulong magnitude, bool negative, uint errorUnits, out long result) |
| | | 80 | | { |
| | 991 | 81 | | result = default; |
| | 991 | 82 | | Fixed64.Multiply64To128(coefficientLow, magnitude, out ulong lowHigh, out ulong low); |
| | 991 | 83 | | Fixed64.Multiply64To128(coefficientMiddle, magnitude, out ulong highHigh, out ulong highLow); |
| | 991 | 84 | | ulong middle = unchecked(lowHigh + highLow); |
| | 991 | 85 | | ulong high = highHigh + (middle < lowHigh ? 1UL : 0UL); |
| | 991 | 86 | | if ((high >> 32) != 0) |
| | 56 | 87 | | return false; |
| | | 88 | | |
| | 935 | 89 | | if (errorUnits != 0) |
| | | 90 | | { |
| | | 91 | | // The radius is below 2^95 because |valueRaw| <= 2^63 and |
| | | 92 | | // errorUnits < 2^32. Only the nearest half-integer can change rounding. |
| | 739 | 93 | | Fixed64.Multiply64To128(magnitude, errorUnits, out ulong errorHigh, out ulong errorLow); |
| | 739 | 94 | | ulong fractionHigh = (uint)middle; |
| | | 95 | | const ulong half = 1UL << 31; |
| | 739 | 96 | | bool belowHalf = fractionHigh < half; |
| | 739 | 97 | | ulong distanceHigh = belowHalf |
| | 739 | 98 | | ? half - fractionHigh - (low != 0 ? 1UL : 0UL) |
| | 739 | 99 | | : fractionHigh - half; |
| | 739 | 100 | | ulong distanceLow = belowHalf ? unchecked(0UL - low) : low; |
| | 739 | 101 | | if (distanceHigh < errorHigh || (distanceHigh == errorHigh && distanceLow <= errorLow)) |
| | 190 | 102 | | return false; |
| | | 103 | | } |
| | | 104 | | |
| | 745 | 105 | | ulong quotient = (high << 32) | (middle >> 32); |
| | 745 | 106 | | bool guard = (middle & (1UL << 31)) != 0; |
| | 745 | 107 | | bool sticky = low != 0 || (middle & ((1UL << 31) - 1)) != 0; |
| | 745 | 108 | | Fixed64 rounded = Fixed64.RoundAndApplySign(quotient, guard, sticky, negative, out bool representable); |
| | 745 | 109 | | if (!representable) |
| | 20 | 110 | | return false; |
| | 725 | 111 | | result = rounded.m_rawValue; |
| | 725 | 112 | | return true; |
| | | 113 | | } |
| | | 114 | | } |