jdk里的大数模幂运算

荷风君
571 阅读5分钟

专有名词

modpow:模幂运算，即，针对数n的exponent次幂取m模
n：底数；
exponent：模数，缩写为ex；
m：模数；

背景

modpow:模幂运算，即，针对数n的exponent次幂取m模
n：底数；exponent：模数，缩写为ex；m：模数
众所周知，对于超过Long的整型数据一般都会使用大数类（BigInteger），而大数的计算效率是影响性能的因素之一，之前有参与算法竞赛。在竞赛中对于大数的计算中，一般不会采用C++手写，而是选择java的BigInteger和BigDecimal。而经常也需要幂取模的操作，而对比手写的快速幂取模，BigInteger的内部效率相对更高，便怀着好奇写下了这篇文章。

jdk源码

（基于jdk1.8.0_291）

声明定义

public BigInteger modPow(BigInteger exponent, BigInteger m)

源码

    /**
     * Returns a BigInteger whose value is
     * <tt>(this<sup>exponent</sup> mod m)</tt>.  (Unlike {@code pow}, this
     * method permits negative exponents.)
     *
     * @param  exponent the exponent.
     * @param  m the modulus.
     * @return <tt>this<sup>exponent</sup> mod m</tt>
     * @throws ArithmeticException {@code m} &le; 0 or the exponent is
     *         negative and this BigInteger is not <i>relatively
     *         prime</i> to {@code m}.
     * @see    #modInverse
     */
    public BigInteger modPow(BigInteger exponent, BigInteger m) {
        if (m.signum <= 0)
            throw new ArithmeticException("BigInteger: modulus not positive");

        // Trivial cases
        if (exponent.signum == 0)
            return (m.equals(ONE) ? ZERO : ONE);

        if (this.equals(ONE))
            return (m.equals(ONE) ? ZERO : ONE);

        if (this.equals(ZERO) && exponent.signum >= 0)
            return ZERO;

        if (this.equals(negConst[1]) && (!exponent.testBit(0)))
            return (m.equals(ONE) ? ZERO : ONE);

        boolean invertResult;
        if ((invertResult = (exponent.signum < 0)))
            exponent = exponent.negate();

        BigInteger base = (this.signum < 0 || this.compareTo(m) >= 0
                           ? this.mod(m) : this);
        BigInteger result;
        if (m.testBit(0)) { // odd modulus
            result = base.oddModPow(exponent, m);
        } else {
            /*
             * Even modulus.  Tear it into an "odd part" (m1) and power of two
             * (m2), exponentiate mod m1, manually exponentiate mod m2, and
             * use Chinese Remainder Theorem to combine results.
             */

            // Tear m apart into odd part (m1) and power of 2 (m2)
            int p = m.getLowestSetBit();   // Max pow of 2 that divides m

            BigInteger m1 = m.shiftRight(p);  // m/2**p
            BigInteger m2 = ONE.shiftLeft(p); // 2**p

            // Calculate new base from m1
            BigInteger base2 = (this.signum < 0 || this.compareTo(m1) >= 0
                                ? this.mod(m1) : this);

            // Caculate (base ** exponent) mod m1.
            BigInteger a1 = (m1.equals(ONE) ? ZERO :
                             base2.oddModPow(exponent, m1));

            // Calculate (this ** exponent) mod m2
            BigInteger a2 = base.modPow2(exponent, p);

            // Combine results using Chinese Remainder Theorem
            BigInteger y1 = m2.modInverse(m1);
            BigInteger y2 = m1.modInverse(m2);

            if (m.mag.length < MAX_MAG_LENGTH / 2) {
                result = a1.multiply(m2).multiply(y1).add(a2.multiply(m1).multiply(y2)).mod(m);
            } else {
                MutableBigInteger t1 = new MutableBigInteger();
                new MutableBigInteger(a1.multiply(m2)).multiply(new MutableBigInteger(y1), t1);
                MutableBigInteger t2 = new MutableBigInteger();
                new MutableBigInteger(a2.multiply(m1)).multiply(new MutableBigInteger(y2), t2);
                t1.add(t2);
                MutableBigInteger q = new MutableBigInteger();
                result = t1.divide(new MutableBigInteger(m), q).toBigInteger();
            }
        }

        return (invertResult ? result.modInverse(m) : result);
    }

分析

1.校验模数

image.png
模数小于对于0时抛出ArithmeticException异常

2.过特例（Trivial cases）

image.png






3.结果是否需要反转（幂次为负数情况）

image.png
当幂次为负的情况下 幂次转为正数并记录标记invertResult

4. 计算基底

image.png
还是依据这个公式，针对负数或者小于模数的数进行初步的取模（应该是针对不需要取模的情况进行效率优化，不考虑这个可以直接取模不进行判断）

5.对模数奇偶进行分类

image.png
奇数情况进6，偶数情况进7

6.针对奇数模数进行计算

6.1进入oddModPow方法

image.png
oddModPow参考了Colin Plumb的C库，运用了窗口移动来不断计算
针对这一部分，博主找了朋友对其进行算法分析，详见www.cnblogs.com/youchandais…

7.针对偶数模数进行计算

7.1.将模数进行拆分

image.png
通过lowbit获取m从右数第一个1，将模数分为一个奇数m1和2的指数m2的乘积

（易得此时m1和m2互质，满足中国剩余定理）

7.2.计算底数关于m1的基底base2

image.png

7.3.计算base和base2的指数化

image.png
计算
计算

7.4.使用中国剩余定理组合a1，a2

image.png


根据计算长度是否会超过最长长度，判断中间计算过程是否可能超过最长限制分为两种计算方法
本质上都是计算

8.反转返回/直接返回

image.png
根据第三步中计算对标记invertResult对结果进行反转返回或直接返回

涉及算法

朴素模幂运算

针对于最朴素的模幂运算是采用这种

int ans = 1;
for(int i = 1; i <= ex; i++){
 	ans*=n;   
}
ans %= m;

而由于
这种思路也可以在每一步进行取模

int ans = 1;
for(int i = 1; i <= ex; i++){
 	ans*=n;
    ans%=m;
}

快速幂

对于任何一个整数的模幂运算   a^b%c   对于b我们可以拆成二进制的形式   b=b0+b12+b22^2+...+bn2^n   这里我们的b0对应的是b二进制的第一位   那么我们的a^b运算就可以拆解成   a^b0a^b12...a^(bn2^n)   对于b来说，二进制位不是0就是1，那么对于bx为0的项我们的计算结果是1就不用考虑了，我们真正想要的其实是b的非0二进制位     那么假设除去了b的0的二进制位之后我们得到的式子是   a^(bx2^x)...a(bn2^n)   这里我们再应用我们一开始提到的公式，那么我们的a^b%c运算就可以转化为   (a^(bx2^x)%c）...(a^(bn2^n)%c)   这样的话，我们就很接近快速幂的本质了   (a^(bx2^x)%c)...(a^(bn2^n)%c)   我们会发现令   A1=(a^(bx2^x)%c)   ...   An=(a^(bn2^n)%c)   这样的话，An始终是A(n-1)的平方倍（当然加进去了取模匀速那），依次递推  

示例代码

int quick(int a,int b,int c)  
{  
    int ans=1;   //记录结果  
    a=a%c;   //预处理，使得a处于c的数据范围之下  
    while(b!=0)  
    {  
        if(b&1) ans=(ans*a)%c;   //如果b的二进制位不是0，那么我们的结果是要参与运算的  
        b>>=1;    //二进制的移位操作，相当于每次除以2，用二进制看，就是我们不断的遍历b的二进制位  
        a=(a*a)%c;   //不断的加倍  
    }  
    return ans;  
}

快速幂还有衍生的矩阵快速幂等算法，有兴趣可以详细了解

中国剩余定理

中国剩余定理（Chinese remainder theorem(CRT)）又称为孙子定理。
出自《孙子算经》第二十六题“数不知数”
有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？
即，一个整数除以三余二，除以五余三，除以七余二，求这个整数。
《孙子算经》中首次提到了同余方程组问题，以及以上具体问题的解法，因此在中文数学文献中也会将中国剩余定理称为孙子定理。
image.png
（证明可自行查找）

参考资料

Powmod快速幂取模
孙子定理

avatar
文章
阅读
粉丝