Алгоритм Ферма для вычисления простых факторов

Согласно ссылке простые множители нечетного числа могут рассчитываться следующим образом:

а = квадрат (N + b ^ 2)

Ниже я написал программу для этого, но я не получаю простые множители 2345678917. Я знаю, что это простое число, но для других простых чисел программа возвращает 1 и само число, но для этого числа этого не происходит. Почему?

#include <stdio.h>
#include <math.h>

void foo(unsigned long long x)
{
    int i;
    for (i=1;i<x;i++)
        if (fmod(sqrt(x + i*i), 1) == 0) {
            printf("%f %f\n", (sqrt(x + i*i) - i), (sqrt(x + i*i)+i));
            return;
        }
}

int main(void) {
    foo((unsigned long long)2345678917);
    return 0;
}

algorithm math prime-factoring
person noman pouigt    schedule 20.10.2015    source источник
comment
На этот вопрос до сих пор нет ответа, поэтому, пожалуйста, опубликуйте свои решения.   -  person noman pouigt    schedule 21.10.2015

Ответы (1)


arrow_upward
2
arrow_downward

Во-первых, тип i (int) не соответствует x (unsigned long long). Измените тип i на unsigned long long, чтобы начать работу.

Как только вы это сделаете, ваша программа радостно провозгласит, что 14 * 167548494 = 2345678917. Конечно, это неверно, поскольку произведение двух четных чисел не может быть нечетным. Проблема здесь заключается в потере точности, поэтому вам нужно реализовать квадратную функцию проверки целых чисел, а не проверять, является ли квадратный корень с плавающей запятой целым.

#include <stdio.h>
#include <math.h>

unsigned long long find_sqrt(unsigned long long x)
{
    unsigned long long lo = 1;
    while (4 * lo * lo <= x) lo *= 2;
    unsigned long long hi = 2 * lo;
    while (lo + 1 < hi) {
        unsigned long long mid = lo + (hi - lo) / 2;
        if (mid * mid <= x) lo = mid;
        else hi = mid;
    }
    return lo * lo == x ? lo : 0;
}

void foo(unsigned long long x)
{
    unsigned long long i;
    for (i=1;i<x;i++) {
        unsigned long long sqrt_x_ii = find_sqrt(x + i*i);
        if (sqrt_x_ii) {
            printf("%llu = %llu * %llu\n",
                   x, sqrt_x_ii - i, sqrt_x_ii + i);
            return;
        }
    }
}

int main(void) {
    foo((unsigned long long) 2345678917);
    return 0;
}
person Mathias Rav    schedule 20.10.2015
comment
Да, он либо очень медленный, либо застрял в инф-цикле, я не останавливал его в отладчике, чтобы проверить. Возможно, он просто давно не находит точных sqrt. - person Peter Cordes; 21.10.2015
comment
@nomanpouigt: Наименьший i, для которого 2345678917 + i*i является квадратом, равен 1172839458. Вас не должно удивлять, что вычисление более миллиарда целых квадратных корней занимает более 5 секунд. В любом случае этот ответ дает важную информацию: а именно, что вы сталкиваетесь с трудностями с точностью с плавающей запятой. Обратите внимание, что для этого значения i 2345678917 + i*i равно 1375552396587412681, что не может быть точно представлено в двойном бинарном64 IEEE 754 (вероятно, это то, что использует ваша платформа). - person Mark Dickinson; 21.10.2015
comment
@nomanpouigt: Кстати, код в этом ответе отлично компилируется и работает для меня, выдавая результат 2345678917 = 1 * 2345678917. Для получения такого результата требуется пара минут. - person Mark Dickinson; 21.10.2015