znlgis 博客

GIS开发与技术分享

第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位数 ab 各分解为两个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 仍然使用自定义实现,原因:

  1. 跨版本兼容:支持 .NET Framework 和旧版 .NET Core
  2. 特定需求:只需要乘法和比较,不需要完整算术
  3. 可控性:可以针对特定场景优化

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 的高精度运算实现:

  1. 溢出问题:64位乘法可能溢出,需要128位存储
  2. UInt128Struct:自定义128位无符号整数结构
  3. 乘法算法:将64位数分解为32位数,分步计算
  4. 符号处理:分别计算绝对值和符号,最后合并
  5. 性能优化:值类型、内联、避免分支

这些高精度运算是 Clipper2 能够正确处理任意大坐标的关键基础。


上一章:InternalClipper内部工具类 返回目录 下一章:ClipperBase基类详解