Newton-Raphson分部与大整数

Question

我正在将BigInt课程作为编程练习.它在base-65536中使用2的补码有符号整数的向量(因此32位乘法不会溢出.一旦我完全工作,我将增加基数).

所有基本的数学运算进行编码,其中一个问题:分裂是痛苦的基本算法,我能够创造慢.(对于商的每个数字,它有点像二进制除法......除非有人想看到它,否则我不会发布它.)

而不是我的慢速算法,我想使用Newton-Raphson找到(移位)倒数然后乘以(和移位).我想我已经掌握了基础知识:你给出公式(x1 = x0(2 - x0*除数))一个很好的初始猜测,然后经过一些迭代后,x收敛到倒数.这部分似乎很容易......但是当我尝试将这个公式应用于大整数时,我遇到了一些问题:

问题1:

因为我正在使用整数......好吧......我不能使用分数.这似乎导致x总是发散(x0*除数似乎必须<2)？我的直觉告诉我应该对方程进行一些修改,使其能够整数运算(达到一定的精度),但我真的很难找出它是什么.(我缺乏数学技能在这里打败了我......)我想我需要找到一些等效的等式而不是d有d*[base ^ somePower]？可以有一些方程式(x1 = x0(2 - x0*d))与整数一致吗？

问题2:

当我使用牛顿的公式来找到某些数字的倒数时,结果最终只是一个小部分,低于答案应该是...... ex.当试图找到4的倒数(十进制):

x0 = 0.3
x1 = 0.24
x2 = 0.2496
x3 = 0.24999936
x4 = 0.2499999999983616
x5 = 0.24999999999999999999998926258176

如果我代表基数为10的数字,我希望得到25的结果(并记住将产品右移2).使用一些倒数(例如1/3),您可以在知道足够的准确度后截断结果.但是我怎样才能从上面的结果中得出正确的倒数呢？

对不起,如果这太模糊或者我要求太多了.我查看了维基百科和我在谷歌上可以找到的所有研究论文,但我觉得我正在撞墙.我感谢任何人都能给我的帮助!

...

编辑:算法运行,虽然它比我预期的要慢得多.与我的旧算法相比,我实际上失去了很多速度,即使是数千位的数字......我仍然缺少一些东西.这不是乘法的问题,这是非常快的.(我确实使用Karatsuba的算法).

对于任何感兴趣的人,这是我目前的Newton-Raphson算法的迭代:

bigint operator/(const bigint& lhs, const bigint& rhs) {
    if (rhs == 0) throw overflow_error("Divide by zero exception");
    bigint dividend = lhs;
    bigint divisor = rhs;

    bool negative = 0;
    if (dividend < 0) {
        negative = !negative;
        dividend.invert();
    }
    if (divisor < 0) {
        negative = !negative;
        divisor.invert();
    }

    int k = dividend.numBits() + divisor.numBits();
    bigint pow2 = 1;
    pow2 <<= k + 1;

    bigint x = dividend - divisor;
    bigint lastx = 0;
    bigint lastlastx = 0;
    while (1) {
        x = (x * (pow2 - x * divisor)) >> k;
        if (x == lastx || x == lastlastx) break;
        lastlastx = lastx;
        lastx = x;
    }
    bigint quotient = dividend * x >> k;
    if (dividend - (quotient * divisor) >= divisor) quotient++;
    if (negative)quotient.invert();
    return quotient;
}

这是我的(非常丑陋)旧算法更快:

bigint operator/(const bigint& lhs, const bigint & rhs) {
    if (rhs == 0) throw overflow_error("Divide by zero exception");
    bigint dividend = lhs;
    bigint divisor = rhs;

    bool negative = 0;
    if (dividend < 0) {
        negative = !negative;
        dividend.invert();
    }
    if (divisor < 0) {
        negative = !negative;
        divisor.invert();
    }

    bigint remainder = 0;
    bigint quotient = 0;
    while (dividend.value.size() > 0) {
        remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1));
        remainder.value.push_back(0);
        remainder.unPad();
        dividend.value.pop_back();

        if (divisor > remainder) {
            quotient.value.push_back(0);
        } else {
            int count = 0;
            int i = MSB;
            bigint value = 0;
            while (i > 0) {
                bigint increase = divisor * i;
                bigint next = value + increase;
                if (next <= remainder) {
                    value = next;
                    count += i;
                }
                i >>= 1;
            }
            quotient.value.push_back(count);
            remainder -= value;
        }
    }

    for (int i = 0; i < quotient.value.size() / 2; i++) {
        int swap = quotient.value.at(i);
        quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i);
        quotient.value.at(quotient.value.size() - 1 - i) = swap;
    }

    if (negative)quotient.invert();
    quotient.unPad();
    return quotient;
}

Answer 1

首先,你可以实现时间上的划分O(n^2)和合理的常数,因此它不会比天真的乘法慢得多.但是,如果你使用类似Karatsuba的算法,甚至是基于FFT的乘法算法,那么你确实可以使用Newton-Raphson加速你的除法算法.

甲牛顿-拉夫逊迭代计算的倒数的x是q[n+1]=q[n]*(2-q[n]*x).

假设我们要计算floor(2^k/B)其中B是一个正整数.WLOG , B?2^k; 否则,商是0.x=B/2^k收益率的Newton-Raphson迭代q[n+1]=q[n]*(2-q[n]*B/2^k).我们可以重新安排它

q[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k

这种迭代只需要整数乘法和位移.它会收敛floor(2^k/B)吗？不必要.然而,在最坏的情况下,它最终在floor(2^k/B)和之间交替ceiling(2^k/B)(证明它!).所以你可以使用一些不那么聪明的测试,看看你是否在这种情况下,并提取floor(2^k/B).(这个"不那么聪明的测试"应该比每次迭代中的乘法更快;但是,优化这个东西会很好).

实际上,计算floor(2^k/B)足以计算floor(A/B)任何正整数A,B.就拿k这样A*B?2^k,和验证floor(A/B)=A*ceiling(2^k/B) >> k.

最后,这种方法的一个简单但重要的优化是在Newton-Raphson方法的早期迭代中截断乘法(即仅计算乘积的较高位).这样做的原因是,早期迭代的结果远非商,并且不正确地执行它们并不重要.(优化这个论点,并表明如果你适当地做这件事,你可以及时划分两位?n整数O(M(2n)),假设你可以及时乘以两位?k整数M(k),并且M(x)是一个增加凸函数).