如何实现高精度运算并掌握128位整数应用？

第7章：高精度运算与128位整数 7.1 概述 在多边形裁剪算法中，精度问题一直是核心挑战。当进行叉积计算时，两个64位整数相乘可能会溢出。Clipper2 通过实现128位整数运算来解决这个问题，确保在任何坐标范围内都能得到精确的结果。 7.2 溢出问题分析 7.2.1 为什么需要128位运算 考虑叉积计算： // 判断三点是否共线或方向关系 // CrossProduct = (x2 - x1) * (y3 - y2) - (y2 - y1) * (x3 - x2) 当坐标值接近 long 的最大值时： long a = 4611686018427387903; // 约 MaxInt64/2 long b = 4611686018427387903; // a * b 将溢出 long 类型！ // 实际结果需要 128 位才能存储 7.2.2 传统解决方案的问题 // 方案1：使用 double - 会丢失精度 double result = (double)a * b; // 精度损失 // 方案2：使用 BigInteger - 性能较差 BigInteger result = new BigInteger(a) * new BigInteger(b); // 堆分配 7.3 UInt128Struct 结构 7.3.1 结构定义 internal struct UInt128Struct { public ulong lo64; // 低64位 public ulong hi64; // 高64位 public UInt128Struct(ulong lo, ulong hi) { lo64 = lo; hi64 = hi; } public static bool operator ==(UInt128Struct a, UInt128Struct b) { return a.lo64 == b.lo64 && a.hi64 == b.hi64; } public static bool operator !=(UInt128Struct a, UInt128Struct b) { return a.lo64 != b.lo64 || a.hi64 != b.hi64; } } 设计原理： 使用两个 ulong（无符号64位整数）组成128位 lo64 存储低64位，hi64 存储高64位 值类型，避免堆分配 7.3.2 128位数值表示 128位整数值 = hi64 * 2^64 + lo64 例如： hi64 = 0x0000000000000001 lo64 = 0x0000000000000000 表示的值 = 1 * 2^64 + 0 = 18446744073709551616 7.4 128位乘法实现 7.4.1 MultiplyUInt64 方法 [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static UInt128Struct MultiplyUInt64(ulong a, ulong b) { // 将64位数分解为高32位和低32位 ulong aLo = a & 0xFFFFFFFF; // a 的低32位 ulong aHi = a >> 32; // a 的高32位 ulong bLo = b & 0xFFFFFFFF; // b 的低32位 ulong bHi = b >> 32; // b 的高32位 // 计算四个部分乘积 ulong axbHi = aHi * bHi; // 高位 × 高位 ulong axbMid = aHi * bLo; // 高位 × 低位 ulong bxaMid = bHi * aLo; // 低位 × 高位 ulong axbLo = aLo * bLo; // 低位 × 低位 // 计算中间结果，处理进位 ulong carry = ((axbMid & 0xFFFFFFFF) + (bxaMid & 0xFFFFFFFF) + (axbLo >> 32)) >> 32; // 组合最终结果 return new UInt128Struct( axbLo + ((axbMid + bxaMid) << 32), // 低64位 axbHi + (axbMid >> 32) + (bxaMid >> 32) + carry // 高64位 ); } 7.4.2 算法原理 将两个64位数 a 和 b 各分解为两个32位数： a = aHi * 2^32 + aLo b = bHi * 2^32 + bLo a * b = (aHi * 2^32 + aLo) * (bHi * 2^32 + bLo) = aHi * bHi * 2^64 + (aHi * bLo + aLo * bHi) * 2^32 + aLo * bLo 图示： aHi (32位) aLo (32位) × bHi (32位) bLo (32位) ───────────────────────────── aLo * bLo → 低64位的低部分 aHi * bLo → 需要左移32位 aLo * bHi → 需要左移32位 aHi * bHi → 需要左移64位 ───────────────────────────── [ hi64 (64位) ][ lo64 (64位) ] 7.4.3 进位处理 关键是正确处理中间乘积的进位： // axbMid 和 bxaMid 相加可能产生65位结果 // 需要将多余的位进位到 hi64 // 计算 lo64： // 1. axbLo 的低32位直接放入 lo64 的低32位 // 2. axbLo 的高32位 + axbMid 的低32位 + bxaMid 的低32位 // = lo64 的高32位 + 可能的进位 ulong carry = ((axbMid & 0xFFFFFFFF) + (bxaMid & 0xFFFFFFFF) + (axbLo >> 32)) >> 32; 7.5 128位比较 7.5.1 等于比较 public static bool operator ==(UInt128Struct a, UInt128Struct b) { return a.lo64 == b.lo64 && a.hi64 == b.hi64; } 7.5.2 大小比较 在 CrossProductSign 中需要比较两个128位数的大小： // 比较 ab 和 cd（两个 UInt128Struct） if (ab.hi64 == cd.hi64) { // 高64位相等，比较低64位 if (ab.lo64 == cd.lo64) return 0; // 完全相等 result = (ab.lo64 > cd.lo64) ? 1 : -1; } else { // 高64位不等，直接根据高64位判断 result = (ab.hi64 > cd.hi64) ? 1 : -1; } 7.6 CrossProductSign 的实现 7.6.1 完整实现 public static int CrossProductSign(Point64 pt1, Point64 pt2, Point64 pt3) { // 计算向量分量 long a = pt2.X - pt1.X; long b = pt3.Y - pt2.Y; long c = pt2.Y - pt1.Y; long d = pt3.X - pt2.X; // 计算 |a| * |b| 和 |c| * |d| 的128位结果 UInt128Struct ab = MultiplyUInt64((ulong)Math.Abs(a), (ulong)Math.Abs(b)); UInt128Struct cd = MultiplyUInt64((ulong)Math.Abs(c), (ulong)Math.Abs(d)); // 计算符号 int signAB = TriSign(a) * TriSign(b); // a*b 的符号 int signCD = TriSign(c) * TriSign(d); // c*d 的符号 // 根据符号和绝对值比较确定叉积符号 if (signAB == signCD) { // 同号时，需要比较绝对值 int result; if (ab.hi64 == cd.hi64) { if (ab.lo64 == cd.lo64) return 0; // |a*b| == |c*d| result = (ab.lo64 > cd.lo64) ? 1 : -1; } else { result = (ab.hi64 > cd.hi64) ? 1 : -1; } // 如果两者都是负数，需要反转结果 return (signAB > 0) ? result : -result; } // 异号时，较大符号决定结果 return (signAB > signCD) ? 1 : -1; } 7.6.2 符号计算逻辑 // 叉积 = a*b - c*d // 我们需要确定 a*b - c*d 的符号 // 情况分析： // 1. signAB > 0, signCD > 0: // 结果符号 = sign(|ab| - |cd|) // 2. signAB < 0, signCD < 0: // 结果符号 = -sign(|ab| - |cd|) // 3. signAB > 0, signCD < 0: // 结果 > 0（正数减负数） // 4. signAB < 0, signCD > 0: // 结果 < 0（负数减正数） // 5. signAB == 0 或 signCD == 0: // 需要特殊处理 7.6.3 TriSign 辅助函数 [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static int TriSign(long x) { return (x < 0) ? -1 : (x > 0) ? 1 : 0; } 返回三种状态：-1（负）、0（零）、1（正）。 7.7 ProductsAreEqual 精确比较 7.7.1 实现代码 internal static bool ProductsAreEqual(long a, long b, long c, long d) { // 将有符号数转换为无符号绝对值 ulong absA = (ulong)Math.Abs(a); ulong absB = (ulong)Math.Abs(b); ulong absC = (ulong)Math.Abs(c); ulong absD = (ulong)Math.Abs(d); // 计算128位乘积 UInt128Struct mul_ab = MultiplyUInt64(absA, absB); UInt128Struct mul_cd = MultiplyUInt64(absC, absD); // 计算符号 int sign_ab = TriSign(a) * TriSign(b); int sign_cd = TriSign(c) * TriSign(d); // 必须绝对值相等且符号相等 return mul_ab.lo64 == mul_cd.lo64 && mul_ab.hi64 == mul_cd.hi64 && sign_ab == sign_cd; } 7.7.2 应用场景 判断三点是否共线： internal static bool IsCollinear(Point64 pt1, Point64 sharedPt, Point64 pt2) { long a = sharedPt.X - pt1.X; long b = pt2.Y - sharedPt.Y; long c = sharedPt.Y - pt1.Y; long d = pt2.X - sharedPt.X; // 当 a*b == c*d 时，叉积为0，三点共线 return ProductsAreEqual(a, b, c, d); } 7.8 性能优化 7.8.1 内联优化 [MethodImpl(MethodImplOptions.AggressiveInlining)] internal static UInt128Struct MultiplyUInt64(ulong a, ulong b) AggressiveInlining 提示JIT编译器内联此方法，减少函数调用开销。 7.8.2 避免分支 // 使用位运算代替条件分支 ulong lo = a & 0xFFFFFFFF; // 而不是 if (a > 0xFFFFFFFF) ... ulong hi = a >> 32; 7.8.3 结构体 vs 元组 // 使用结构体而非元组 internal struct UInt128Struct // 更好 // 而不是 (ulong lo, ulong hi) // 需要装箱 7.9 与 .NET 内置类型的比较 7.9.1 .NET 6+ 的 UInt128 从 .NET 6 开始，.NET 提供了内置的 UInt128 类型。但 Clipper2 仍然使用自定义实现，原因： 跨版本兼容：支持 .NET Framework 和旧版 .NET Core 特定需求：只需要乘法和比较，不需要完整算术 可控性：可以针对特定场景优化 7.9.2 BigInteger 对比 // BigInteger 方式 BigInteger result = new BigInteger(a) * new BigInteger(b); // Clipper2 方式 UInt128Struct result = MultiplyUInt64((ulong)Math.Abs(a), (ulong)Math.Abs(b)); 差异： BigInteger 是引用类型，需要堆分配 BigInteger 功能完整但开销较大 UInt128Struct 是值类型，无GC压力 7.10 实际应用示例 7.10.1 判断点的位置关系 public static int PointRelativeToLine(Point64 p, Point64 lineStart, Point64 lineEnd) { // 使用高精度叉积符号判断 int sign = CrossProductSign(lineStart, lineEnd, p); if (sign > 0) return 1; // 点在线的左侧 if (sign < 0) return -1; // 点在线的右侧 return 0; // 点在线上 } 7.10.2 多边形方向判断 public static bool IsClockwise(Path64 polygon) { // 计算有符号面积 // 使用高精度运算避免大坐标溢出 double area = 0; int cnt = polygon.Count; for (int i = 0; i < cnt; i++) { Point64 pt1 = polygon[i]; Point64 pt2 = polygon[(i + 1) % cnt]; // 注意：这里使用 double 是因为面积累加 // 单次叉积计算使用 CrossProductSign 可以避免溢出 area += (double)(pt1.X + pt2.X) * (pt2.Y - pt1.Y); } return area > 0; // 正面积为顺时针（Y轴向下） } 7.11 边界情况处理 7.11.1 零值处理 // 当 a 或 b 为 0 时 if (a == 0 || b == 0) { // TriSign 返回 0，signAB = 0 // 需要特殊处理 } 7.11.2 最大值测试 // 测试极端情况 long maxCoord = InternalClipper.MaxCoord; Point64 p1 = new Point64(-maxCoord, -maxCoord); Point64 p2 = new Point64(maxCoord, maxCoord); Point64 p3 = new Point64(maxCoord, -maxCoord); // 128位运算确保不会溢出 int sign = CrossProductSign(p1, p2, p3); 7.12 本章小结 本章详细介绍了 Clipper2 的高精度运算实现： 溢出问题：64位乘法可能溢出，需要128位存储 UInt128Struct：自定义128位无符号整数结构 乘法算法：将64位数分解为32位数，分步计算 符号处理：分别计算绝对值和符号，最后合并 性能优化：值类型、内联、避免分支 这些高精度运算是 Clipper2 能够正确处理任意大坐标的关键基础。 上一章：InternalClipper内部工具类 | 返回目录 | 下一章：ClipperBase基类详解