Estou fazendo uma aula BigInt como um exercício de programação. Ele usa um vetor de entradas assinadas de complemento de 2 na base-65536 (para que as multiplicações de 32 bits não excedam. Eu aumentarei a base assim que eu estiver funcionando totalmente).

Todas as operações matemáticas básicas são codificadas, com um problema: a divisão édolorosamente lento com o algoritmo básico que consegui criar. (Funciona como divisão binária para cada dígito do quociente ... Não vou publicá-lo, a menos que alguém queira vê-lo ....)

Em vez do meu algoritmo lento, quero usar Newton-Raphson para encontrar o recíproco (deslocado) e multiplicar (e deslocar). Eu acho que estou pensando no básico: você dá a fórmula(x1 = x0 (2 - x0 * divisor)) um bom palpite inicial e, depois de algumas iterações, x converge para o recíproco. Esta parte parece fácil ... mas estou com alguns problemas ao tentar aplicar esta fórmula a números inteiros grandes:

Problema 1:

Porque estou trabalhando com números inteiros ... bem ... não posso usar frações. Isso parece causar x sempre divergir (x0 * divisor deve ser <2, parece?). Minha intuição me diz que deve haver alguma modificação na equação que permita trabalhar números inteiros (com alguma precisão), mas estou realmente lutando para descobrir o que é. (Minha falta de habilidades em matemática está me superando aqui ....) Acho que preciso encontrar uma equação equivalente onde, em vez ded Há simd * [base ^ somePower]? Pode haver alguma equação como(x1 = x0 (2 - x0 * d)) que funciona com números inteiros?

Problema 2:

Quando uso a fórmula de Newton para encontrar o recíproco de alguns números, o resultado acaba sendo apenas uma pequena facção abaixo da resposta que deveria ser ... ex. ao tentar encontrar recíproco de 4 (em decimal):

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

Se eu estivesse representando números na base 10, desejaria um resultado de 25 (e lembre-se de trocar o produto à direita por 2). Com alguns recíprocos, como 1/3, você pode simplesmente truncar o resultado depois de saber que possui precisão suficiente. Mas como posso extrair a recíproca correta do resultado acima?

Desculpe se isso é muito vago ou se estou pedindo demais. Examinei a Wikipedia e todos os documentos de pesquisa que pude encontrar no Google, mas sinto que estou batendo com a cabeça na parede. Agradeço qualquer ajuda que alguém possa me dar!

...

Edit: O algoritmo está funcionando, embora seja muito mais lento do que eu esperava. Na verdade, perdi muita velocidade em comparação com meu antigo algoritmo, mesmo em números com milhares de dígitos ... Ainda estou perdendo alguma coisa. Não é um problema de multiplicação, que é muito rápido. (De fato, estou usando o algoritmo de Karatsuba).

Para qualquer pessoa interessada, aqui está minha iteração atual do algoritmo 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;
}

E aqui está o meu algoritmo antigo (realmente feio) que é mais rápido:

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;
}