【笔记】Big Integer - Lehmer's Euclidean Algorithm

最新推荐文章于 2023-06-07 17:14:52 发布

xueyan0096

最新推荐文章于 2023-06-07 17:14:52 发布

阅读量256

点赞数

本文链接：https://blog.csdn.net/xueyan0096/article/details/84140429

版权

Lehmer's Euclidean Algorithm

给定正整数a, b (a >= b)，重复取二进制前32位进行欧几里得算法，直到a和b的二进制长度小于等于32位，

这时用uint版本的欧几里得算法求得的gcd就是原来的a和b的gcd;

/// <summary>
/// 计算<paramref name="a"/>和<paramref name="b"/>的最大公约数
/// </summary>
public static BigInteger GreatestCommonDivisor(BigInteger a, BigInteger b)
{
    if (a.IsZero) return b;
    if (b.IsZero) return a;
    if (a.Sign == -1) a = -a;
    if (b.Sign == -1) b = -b;            
    if (a < b)
        Swap(ref a, ref b);
            
    var u = new long[3];
    var v = new long[3];
    var t = new long[3];
    var tmp = BigInteger.Zero;
    while (true)
    {
        //D1 init
        if (b.GetByteCount() <= 4)
            return GreatestCommonDivisor((long)a, (long)b);
                
        u[0] = 1; u[1] = 0;
        v[0] = 0; v[1] = 1;
        var k = a.GetByteCount()*8 - 32;
        u[2] = (long)(a >> k);
        v[2] = (long)(b >> k);
                
        while (v[2] + v[0] != 0 && v[2] + v[1] != 0)
        {
            //D2 test quotient
            var q = (u[2] + u[0]) / (v[2] + v[0]);
            if (q != (u[2] + u[1]) / (v[2] + v[1])) break;
            //D3 simulate Euclidean algorithm
            t[0] = u[0] - v[0] * q;
            t[1] = u[1] - v[1] * q;
            t[2] = u[2] - v[2] * q;

            u[0] = v[0]; u[1] = v[1]; u[2] = v[2];

            v[0] = t[0]; v[1] = t[1]; v[2] = t[2];
        }
        //D4 multi-digit calculaton
        if (u[1] == 0)
        {
            tmp = a % b;
            a = b;
            b = tmp;
        }
        else
        {
            tmp = u[0] * a + u[1] * b;
            b = v[0] * a + v[1] * b;
            a = tmp;
        }
    }
}

/// <summary>
/// 计算<paramref name="a"/>和<paramref name="b"/>的最大公约数
/// </summary>
public static long GreatestCommonDivisor(long a, long b)
{
    if (a == 0) return b;
    if (b == 0) return a;

    a = (a ^ (a >> 63)) - (a >> 63);//a = abs(a)
    b = (b ^ (b >> 63)) - (b >> 63);//b = abs(b)

    var aZeros = a.NumberOfTrailingZeros();
    var bZeros = b.NumberOfTrailingZeros();
    a >>= aZeros;
    b >>= bZeros;

    var t = Math.Min(aZeros, bZeros);

    while (a != b)
    {
        if (a > b)
        {
            a -= b;
            a >>= a.NumberOfTrailingZeros();
        }
        else
        {
            b -= a;
            b >>= b.NumberOfTrailingZeros();
        }
    }

    return a << t;
}

/// <summary>
/// <paramref name="n"/>二进制尾数0的个数
/// </summary>
public static int NumberOfTrailingZeros(this long n)
{
    if (n == 0) return 64;

    var c = 63;
    var t = n << 32; if (t != 0L) { c -= 32; n = t; }
    t = n << 16; if (t != 0L) { c -= 16; n = t; }
    t = n << 8; if (t != 0L) { c -= 8; n = t; }
    t = n << 4; if (t != 0L) { c -= 4; n = t; }
    t = n << 2; if (t != 0L) { c -= 2; n = t; }
    t = n << 1; if (t != 0L) { c -= 1; }

    return c;
}