Опровержение предложения Эйлера грубой силой в C

Я написал немного кода для грубой силы, опровергая предложение Эйлера , в котором говорится

  

$$ a ^ 4 + b ^ 4 + c ^ 4 = d ^ 4 $$

     

не имеет решения, когда \ $ a \ $, \ $ b \ $, \ $ c \ $, \ $ d \ $ - целые положительные числа.

Я не математик, но, читая это, по крайней мере одно решение было найдено Noam Elkies в 1987 году (a = 95800; b = 217519; c = 414560; d = 422481). Я хотел получить представление о том, сколько огневой мощи потребуется решить с помощью грубой силы, поэтому я написал следующее в C:

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

int prop(long int A, long int B, long int C, long int D) {
    return (pow(A, 4) + pow(B, 4) + pow(C, 4) == pow(D, 4));
}

int main() {
    long int a, b, c, d;
    clock_t t;
    t = clock();

    for (a = 1; a < 100000; a++) {
        for (b = 1; b < 300000; b++) {
            for (c = 1; c < 500000; c++) {
                for (d = 1; d < 500000; d++) {
                    if (prop(a, b, c, d))
                        printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
                }
                if (!(c%1000))
                    printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, ((double)(clock() - t))/CLOCKS_PER_SEC);
            }
            printf("a = %ld, b = %ld, time = %fs\n", a, b, ((double)(clock() - t))/CLOCKS_PER_SEC);
        }
        printf("a = %ld, time = %fs\n", a, ((double)(clock() - t))/CLOCKS_PER_SEC);
    }

}

Я запустил его некоторое время и решил, что для этого потребуется примерно $ 85 \ times 10 ^ 6 \ $ лет, чтобы получить ответ выше на моем текущем компьютере.

Я имею в виду, может быть, если бы я ждал несколько тысяч лет, у меня была бы немного лучшая машина, но мой вопрос (и я новичок в C и компьютерной науке в целом - так что будьте осторожны); какие стратегии я мог бы предпринять, чтобы сделать вышеприведенный код быстрее? Я думал о потоке и использовал вместо смещение бит вместо вызовов pow().

Будет ли какой-либо способ (на моей нынешней машине) заставить его работать в моей жизни?

123 голоса | спросил Aidenhjj 25 +03002016-10-25T20:20:08+03:00312016bEurope/MoscowTue, 25 Oct 2016 20:20:08 +0300 2016, 20:20:08

17 ответов


128

Примечание: В какой-то момент этот обзор перешел в сферу ассемблера и GMP. Фактический обзор находится в конце этого сообщения, тогда как в первом разделе обсуждаются проблемы времени выполнения pow, неправильные типы данных и произвольные большие целые числа.

Нет времени жизни для времени выполнения

  

Будет ли какой-либо способ (на моей нынешней машине) заставить его работать в моей жизни?

В столярной мастерской есть великая поговорка: измерьте дважды, вырезайте один раз. Это касается резки дерева или другого материала, где вам нужно выбросить свои ресурсы, если вы случайно порежете не в том месте.

Аналогичное высказывание существует для инженеров-программистов: вы не можете оптимизировать то, что вы не можете измерить. Существует несколько способов измерения вашего кода, например. бенчмаркинг, профилирование или просмотр сгенерированного ассемблера, чтобы узнать, сколько команд потребуется определенной части вашего кода.

Здесь мы рассмотрим последний маршрут, начнем с ассемблера, шаг за шагом рассмотрим, где мы закончим.

Исследование в сборке

Давайте посмотрим на ваш код. Ну, а не ваш, но ассемблер генерирует компилятор. Вы можете использовать gcc -S -O3. На моей платформе это приводит к следующему «горячему» разделу в main:

.L6:
    movapd  %xmm6, %xmm1
    movapd  %xmm10, %xmm0
    call    pow
    movapd  %xmm6, %xmm1
    movapd  %xmm0, %xmm13
    movapd  %xmm11, %xmm0
    call    pow
    movapd  %xmm6, %xmm1
    movapd  %xmm0, %xmm8
    movapd  %xmm9, %xmm0
    call    pow
    movapd  %xmm0, %xmm7
    pxor    %xmm0, %xmm0
    movapd  %xmm6, %xmm1
    cvtsi2sdq   %rbx, %xmm0
    call    pow
    addsd   %xmm13, %xmm8
    addsd   %xmm8, %xmm7
    ucomisd %xmm0, %xmm7

Даже если вы можете не знать ассемблера, но вы можете увидеть эти четыре вызова pow. Первое, что вам нужно знать, это то, что call медленнее по сравнению с другими операциями. Эти четыре вызова call происходят в самом внутреннем цикле. Компилятор удалил вызов call в prop и вместо этого заменил его своим кодом (это быстрее).

mov* назначает значения для регистров, add* добавляет значения и т. д. Регистры с xmm* - это регистры с двойной точностью, предназначенные для переменных double. Поэтому мы в основном вызываем pow с правильными значениями, а затем добавляем, вычитаем и изменяем наши маленькие маленькие двойные значения.

Двойная проблема

Но подождите секунду. Мы пытаемся решить полностью интегральную задачу! Почему наша сгенерированная программа использует эти регистры вообще?

Это должно поднять красный флаг. И действительно, если мы помним подпись pow , это должно ясно, что это не правильный инструмент. Он принимает двойную базу и показатель степени, что указывает на то, что он подходит для таких терминов, как \ $ 15.12151 ^ {3.1415926} \ $. Это полный избыток вашей проблемы.

Используя собственные функции

Итак, давайте вместо этого используем другую версию pow:

long int pow4(long int x){
    return x * x * x * x;
}

Обратите внимание, что ваш компилятор должен создать что-то вроде этого:

movq    %rdi, %rax
imulq   %rdi, %rax
imulq   %rax, %rax
ret

, но если ваш компилятор не распознает эту потенциальную (микро) оптимизацию, вы можете использовать

long int pow4(long int x){
    const long int t = x * x;
    return t * t;
}

вместо.

Нам также нужно изменить prop:

int prop(long int A, long int B, long int C, long int D) {
  return (pow4(A) + pow4(B) + pow4(C) == pow4(D));
}

Allright. Теперь, прежде чем я покажу время новой программы, давайте проверим вывод вашего старого:

a = 1, b = 1, c = 1000, время = 114.156000s

Вот тогда я ударил ^ C. Как поддерживается pow4?

a = 1, b = 1, c = 1000, время = 0.296000s
a = 1, b = 1, c = 2000, время = 0,578000 с
a = 1, b = 1, c = 3000, время = 0,859000s
a = 1, b = 1, c = 4000, время = 1,140000s
a = 1, b = 1, c = 5000, время = 1.421000s
a = 1, b = 1, c = 6000, время = 1,703000s
a = 1, b = 1, c = 7000, время = 1,984000s
a = 1, b = 1, c = 8000, время = 2,265000s
a = 1, b = 1, c = 9000, время = 2,546000s
a = 1, b = 1, c = 10000, время = 2,828000s
a = 1, b = 1, c = 11000, время = 3.109000s
a = 1, b = 1, c = 12000, время = 3,390000s
a = 1, b = 1, c = 13000, время = 3,687000s
a = 1, b = 1, c = 14000, время = 3,968000s
a = 1, b = 1, c = 15000, время = 4.250000s
a = 1, b = 1, c = 16000, время = 4,531000 с

Это 0,2% от вашего первоначального времени или 500-кратное ускорение.

Однако , это происходит за счет: pow4(500000) слишком велико для int64_t, так как \ $ \ log_2 (500000 ^ 4) \ approx 76 \ $. Наибольшее число, которое вы можете проверить с помощью uint64_t, - это 65535, \ $ 2 ^ {16} -1 \ $,что не должно быть удивительно. Поскольку стандарт не предоставляет int128_t или аналогичный, вы должны убедиться, что ваши номера не превышают эти границы.

Вы можете либо написать свою собственную большую целую логику для этого, либо использовать GMP.

Собственные оценки и оценки параметров

Далее вы можете увеличить нижние границы b и c, чтобы \ $ a \ le b \ le c \ $. И для d, ну, если у нас есть a, b, c, то есть только одно решение для d. Мы можем напрямую искать это решение с бинарным поиском.

Двоичный поиск делает алгоритм \ $ \ mathcal O (n ^ 3 \ log n) \ $ из вашего текущего \ $ \ mathcal O (n ^ 4) \ $ one, который должен обеспечивать намного большую скорость, чем предыдущее ускорение.

Еще лучше, если вы использовали соответствующие ограничения для a, b и c, мы можем связать d by

$$ d ^ 4 = a ^ 4 + b ^ 4 + c ^ 4 \ le c ^ 4 + c ^ 4 + c ^ 4 = 3c ^ 4 $$

и, следовательно, получим

$$ c \ le d \ le \ sqrt [4] {3} \, c. $$

С помощью правильного бинарного алгоритма вы можете быстро завершить первый a=1, b=1:

â € |
a = 1, b = 1, c = 463000, время = 0,031000s
a = 1, b = 1, c = 464000, время = 0,031000s
a = 1, b = 1, c = 465000, время = 0,031000s
a = 1, b = 1, c = 466000, время = 0,031000s
a = 1, b = 1, c = 467000, время = 0,031000s
a = 1, b = 1, c = 468000, время = 0,031000s
a = 1, b = 1, c = 469000, время = 0,031000s
a = 1, b = 1, c = 470000, время = 0,031000s
a = 1, b = 1, c = 471000, время = 0,031000s
a = 1, b = 1, c = 472000, время = 0,031000s
a = 1, b = 1, c = 473000, время = 0,031000s
a = 1, b = 1, c = 474000, время = 0,031000s
a = 1, b = 1, c = 475000, время = 0,031000s
a = 1, b = 1, c = 476000, время = 0,031000s
a = 1, b = 1, c = 477000, время = 0,031000s
a = 1, b = 1, c = 478000, время = 0,031000s
a = 1, b = 1, c = 479000, время = 0,031000s
a = 1, b = 1, c = 480000, время = 0,031000s
a = 1, b = 1, c = 481000, время = 0,031000s
a = 1, b = 1, c = 482000, время = 0,031000s
a = 1, b = 1, c = 483000, время = 0,031000s
a = 1, b = 1, c = 484000, время = 0,031000s
a = 1, b = 1, c = 485000, время = 0,031000s
a = 1, b = 1, c = 486000, время = 0,031000s
a = 1, b = 1, c = 487000, время = 0,031000s
a = 1, b = 1, c = 488000, время = 0,031000s
a = 1, b = 1, c = 489000, время = 0,031000s
a = 1, b = 1, c = 490000, время = 0,031000s
a = 1, b = 1, c = 491000, время = 0,031000s
a = 1, b = 1, c = 492000, время = 0,031000s
a = 1, b = 1, c = 493000, время = 0,031000s
a = 1, b = 1, c = 494000, время = 0,031000s
a = 1, b = 1, c = 495000, время = 0,031000s
a = 1, b = 1, c = 496000, время = 0,031000s
a = 1, b = 1, c = 497000, время = 0,031000s
a = 1, b = 1, c = 498000, время = 0,031000s
a = 1, b = 1, c = 499000, время = 0,031000s
a = 1, b = 1, время = 0,031000 с

Это возвращает нас в сферу вашей жизни.

Упражнение

Напишите функцию, которая задана a, b и c проверяет, существует ли d, такой что ваша собственность держится. Он должен возвращать -1, если не существует такого d и d в противном случае.

Используйте эту функцию в своем коде. Убедитесь, что вам нужна грубая \ $ \ log d _ {\ text {max}} \ $ итерация в этой функции.

Важное замечание о размерах целых

Имейте в виду, что long int обычно представляет собой просто 64-битное целое число, что означает, что наибольшее целое число, которое вы можете сохранить, - \ $ 2 ^ {63} -1 \ $. Целочисленные типы с большим количеством битов имеют большие границы, но специфичны для платформы. Кроме того, умножение может быть медленнее, так как умножение 128-битных чисел не так просто, как умножение 64-битных чисел.

См. следующий раздел, как получить умножения вниз

Фактический обзор

Наш pow4 теперь представляет собой по существу два умножения. Однако мы по-прежнему часто используем pow4. В конце концов, нам не нужно пересчитывать \ $ a ^ 4 \ $ на каждой итерации. Компилятор с радостью делает, так как он недостаточно оптимизирован.

Это подводит нас к фактическому обзору: ваш код написан просто, легко читается и понимается. К сожалению, хорошо написанный модульный код часто не выжимает последний бит (хе) из вашего оборудования, если ваш компилятор /время выполнения не очень умное (и, следовательно, часто дорого).

Итак, давайте вернемся к чертежной доске для окончательного обзора вашего кода:

Включает

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

Я бы отсортировал их по имени, но все в порядке. Вы не включаете ничего, что вам не нужно, и вы ничего не забыли (и были спасены нестандартным компилятором).

Декларация

int main() {
    long int a, b, c, d;
    clock_t t;
    t = clock();

В зависимости от того, будете ли вы писать ANSI-C или C99, я отложил бы объявление переменных как можно дольше. Например, на данный момент легкослучайно измените c на какое-то фиктивное значение или забудьте { и случайно проверьте prop после циклов или подобных:

for (a = 1; a < 100000; a++) 
    for (b = 1; b < 300000; b++) 
        for (c = 1; c < 500000; c++) 
            for (d = 1; d < 500000; d++) 
                printf("inner loop");
                if (prop(a,b,c,d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);

Упс. if не проверяется, и вы не получите предупреждение (в старых версиях компилятора, новые предупреждают о возможных проблемах с пробелами). Если вы позже объявите свои переменные (например, C99-стиль), таких ошибок не может быть (хотя он вводит возможное затенение):

for (long int a = 1; a < 100000; a++) 
    for (long int b = 1; b < 300000; b++) 
        for (long int c = 1; c < 500000; c++) 
            for (long int d = 1; d < 500000; d++) 
                printf("inner loop");
                if (prop(a,b,c,d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd %ld\n", a, b, c, d);

Теперь это приведет к ошибке компилятора, поскольку a, b и т. д. выходят за рамки. В любом случае, это зависит от языка, который вы хотите использовать. Некоторые люди предпочитают один путь, другие - другой. Выберите свой.

Результаты кэша (вручную)

Одна вещь, которая больше всего меня поражает, заключается в том, что вы каждый раз пересчитываете \ $ a ^ 4 \ $ и так далее. Мы можем легко обработать это с помощью большего количества переменных (используя стиль объявления):

int main() {
    long int a, b, c, d;
    long int a4, b4, c4, d4; // new variables
    clock_t t;
    t = clock();

    for (a = 1; a < 100000; a++) {
        a4 = pow4(a);                             // remember
        for (b = a; b < 300000; b++) {
            b4 = pow4(b);                         // remember
            for (c = b; c < 500000; c++) {
                c4 = pow4(c);                     // the fourth power
                for (d = c; d < 500000; d++) {
                    d4 = pow4(d);                 // of this member
                    if (a4 + b4 + c4 == d4)
                        printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
…

Помните, как я сказал, что красиво написанный модульный код не всегда оптимален? Это один из тех неудачных примеров, когда вам нужно помочь компилятору (если вы точно не знаете, какие флаги оптимизации вы должны использовать или ваш компилятор слишком агрессивен). prop отсутствует, вызовы pow4 теперь находятся в вашем цикле.

Но компилятор здесь больше не может совершить ошибку: теперь very ясно, что a4 не нужно пересчитывать 300000 * 500000 * 500000 раз.

Пришло время

Хотя наш код теперь более подробный, есть один маленький фрагмент кода, который три раза повторяется в вашем main:

((double)(clock() - t))/CLOCKS_PER_SEC)

Это довольно трудно читать, не так ли? Это идеальный кандидат на выполнение функции:

static inline seconds_since(clock_t t){
    return ((double)(clock() - t))/CLOCKS_PER_SEC;
}

Это изменяет ваш printf на

printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, ((double)(clock() - t))/CLOCKS_PER_SEC);

к

printf("a = %ld, b = %ld, c = %ld, time = %fs\n", a, b, c, seconds_since(t));

Ах. Намного легче читать. Вот для чего предназначены функции inline. Обратите внимание, что любой сложный компилятор должен встроить эту функцию в любом случае, поэтому вы можете также отказаться от inline, если вы не хотите использовать C99.

ответил Zeta 25 +03002016-10-25T21:20:34+03:00312016bEurope/MoscowTue, 25 Oct 2016 21:20:34 +0300 2016, 21:20:34
73

Хорошо, давайте начнем с этого:

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {

Давайте теперь будем игнорировать d. Что ты здесь делаешь? Вы проверяете 1, 1, 1, затем 1, 1, 2, затем 1, 1, 3, ... до 1, 1, 499999. Затем вы начинаете с 1, 2, 1. Но вы уже проверено 1, 1, 2 , так почему вы проверяете 1, 2, 1? Вы можете перейти прямо к 1, 2, 2. Это не спасает вас от этих низких чисел, но поверьте мне, когда вы дойдете до больших чисел, это добавит.

Короче: a, b, c должны быть неубывающие . Мы можем достичь этого, начав b при a и начиная c с b, так что b никогда не меньше a, а c никогда не меньше b. Поэтому сразу вы можете избавиться от примерно половины работы, которую вы делаете с помощью

for (a = 1; a < 100000; a++) {
    for (b = a; b < 300000; b++) {
        for (c = b; c < 500000; c++) {

Затем рассмотрим d.

Конечно, мы можем сразу сделать ту же оптимизацию, которую мы только что сделали для a, b и c, так как d никогда не будет меньше , чем любой из них, а c всегда самый большой.

Кроме того, как только d 4 будет больше , чем сумма, мы можем остановить приращение d, потому что это только увеличится.

Так что это сэкономит много времени прямо там. Но мы можем сделать way лучше этого.

Вопрос, который вы задаете, заключается в следующем: «У этих четырех чисел есть свойство суммы?» но вопрос, который вы должны задавать, - это «a 4 + b 4 + c 4 равный любой четвертый уровень ?» Если да, то вы можете легко вычислить d намного быстрее, чем пытаться использовать все возможные четвертые полномочия. Итак, можете ли вы написать быстрый метод, который говорит вам, является ли конкретная сумма четвертой властью или нет?

Если вы знаете, как взять квадратный корень, что вы можете сделать, это взять сумму, дважды взять квадратный корень и затем дважды выровнять результат. Если вы вернете первоначальную сумму, то это будет четвертая власть, а если нет, то это не так.

Если вы не знаете, как взять квадратный корень, вы можете сделать следующую логику: у нас есть сумма; составляет 128 4 , равное сумме? Нет. Это больше? Да. Затем попробуйте 64 4 . Он равен сумме? Нет. Это меньше? Да. Затем попробуйте 96 4 и т. Д. Двоичный поиск результата до тех пор, пока вы его не найдете, или пока в вашем диапазоне не будет больше номеров для проверки. Это делает лишь крошечную часть вычислений, по сравнению с попытками тысяч и тысяч возможных четвертых сил.

ответил Eric Lippert 25 +03002016-10-25T21:13:35+03:00312016bEurope/MoscowTue, 25 Oct 2016 21:13:35 +0300 2016, 21:13:35
31

Все ответы (когда этот «ответ» был впервые написан) игнорируют важное ограничение: родные типы. В C код long int - это 32-разрядный (или более) подписанный тип, что означает, что наибольшее положительное значение, которое может быть (считается) представлено: \ $ 2 ^ {31} -1 \ $. Наибольший возможный вход для функции, которая вычисляет четвертую мощность и [[может быть засчитана]], возвращает правильный длинный int, это \ $ \ left \ lfloor \ sqrt [4] {2 ^ {31} -1} \ право \ rfloor = 215 \ $. Даже если вы перейдете к unsigned long long int, 64-разрядному (или более) неподписанному типу, \ $ \ left \ lfloor \ sqrt [4] {2 ^ {64} -1} \ right \ rfloor = 65 \, 535 \ $. Тип с плавающей точкой не приведет к дальнейшему: только биты мантиссы ценны из-за необходимой точности. (Если у вас есть доступ к 128-битовому типу, этого было бы достаточно, как \ $ \ left \ lfloor \ sqrt [4] {2 ^ {128} -1} \ right \ rfloor = 4 \, 294 \, 967 \ , 295 \ $).

Если родные типы не будут полностью достаточными, вам понадобится «большая библиотека номеров», которая обрабатывает такие значения. Предупреждение: необходимость вызова в библиотеку, скорее всего, приведет к хиту производительности.

Все, что было сказано, все еще возможно написать код в основном с родными типами (до 128 бит) таким образом, чтобы добавлять и сравнивать только с младшими битами, но вы бы должны быть очень осторожными, и у вас есть возможность ложных срабатываний. В этот момент вы можете позвонить в «большую библиотеку номеров», чтобы проверить нахождение.

ответил whydoubt 26 +03002016-10-26T00:22:31+03:00312016bEurope/MoscowWed, 26 Oct 2016 00:22:31 +0300 2016, 00:22:31
30

Некоторые интересные вещи никто не упомянул об этом, но это улучшение при поиске наименьшего решения:

Если \ $ \ gcd (a, b, c) \ neq 1 \ $, то \ $ a ^ 4 \ $, \ $ b ^ 4 \ $, \ $ c ^ 4 \ $ и \ $ d ^ 4 \ $ делятся на gcd на 4-ю степень, давая меньшее решение.

Следовательно, хотя бы один из \ $ a \ $, \ $ b \ $ и \ $ c \ $ должен быть нечетным, так как если они все четные, то они имеют общий делитель 2, что означает \ $ \ gcd (a, b, c) \ geq 2 \ $. Проверяя, что они не все четные, перед выполнением цикла для d вы сократите свой возможный поиск на большое количество.

В то же время, если \ $ d \ $ четно, то 2 должно делить \ $ a ^ 4 + b ^ 4 + c ^ 4 \ $, но мы уже говорили, что они не могут быть даже четными быть наименьшим решением. Если два (или ноль) из трех переменных четны, тогда сумма их 4-го уровня будет нечетной (четная + четная + нечетная = нечетная, нечетная + нечетная + нечетная = нечетная), поэтому, если d четно, тогда 1 число должно быть четным, а другое 2 нечетным.

Если \ $ d \ $ нечетно, то 2 не может делить \ $ a ^ 4 + b ^ 4 + c ^ 4 \ $. Если ровно 1 из чисел, то сумма их 4-й степени будет четной (четной + нечетной + нечетной = четной). При этом либо они все нечетные, либо только 1 нечетно.

Используя это, вы можете решить, является ли D нечетным или даже основанным на наборе \ $ a, b, c \ $, сокращающем ваше решение, заданное пополам, соответствующим образом адаптируя самый внутренний цикл. В сочетании с возрастающими значениями d> c> b> a это должно улучшить обработку с большим отрывом.

ответил mascoj 26 +03002016-10-26T03:17:44+03:00312016bEurope/MoscowWed, 26 Oct 2016 03:17:44 +0300 2016, 03:17:44
29

Я не могу комментировать решение vnp, но vnp был прав в первый раз: вы можете использовать его в \ $ O (n ^ 2 \ log (n)) \ $ time и \ $ O (n ) \ $. Вам не нужно \ $ O (n ^ 2) \ $ пространство, потому что вам не нужно хранить весь список \ $ a ^ 4 + b ^ 4 \ $ или \ $ d ^ 4-c ^ 4 \ $ upfront. Вместо этого вам нужно только иметь возможность перечислять значения \ $ a ^ 4 + b ^ 4 \ $ и \ $ d ^ 4-c ^ 4 \ $ в порядке возрастания, а затем совпадать с двумя значениями.

Чтобы перечислить значения \ $ a ^ 4 + b ^ 4 \ $, вы можете сохранить для каждого \ $ a \ $ наименьшее значение \ $ b \ $, назовите его \ $ b_a \ $, такой, что \ $ a ^ 4 + b_a ^ 4 \ ge k \ $, где \ $ k \ $ - предыдущее значение \ $ a ^ 4 + b ^ 4 \ $, которое вы читаете. Чтобы прочитать следующее значение, вам нужно найти значение \ $ a \ $, так что \ $ a ^ 4 + b_a ^ 4 \ $ является наименьшим, что вы можете сделать в \ $ O (\ log (n) ) \ $ время, используя кучу или что-то в этом роде. Затем вы увеличиваете \ $ b_a \ $, и вы готовы прочитать следующее значение \ $ a ^ 4 + b ^ 4 \ $.

Это пример программы C ++ 11 (с использованием расширения GNU для получения 128-битных целых чисел), который находит решение примерно через 7 часов на моем ПК. (Это просто для иллюстрации метода - это не очень эффективно другими способами.)

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <queue>
#include <vector>
using std::priority_queue;
using std::vector;

//typedef long long int bigint;
typedef __int128 bigint;

// Look for a^4+b^4+c^4 = d^4 by comparing
// ascending list of a^4+b^4 with ascending list of d^4-c^4

vector<bigint> list;
vector<int> diffptr,sumptr;

bool sumcmp(int a0,int a1){
  int b0=sumptr[a0];
  int b1=sumptr[a1];
  return list[a0]+list[b0]>list[a1]+list[b1];
}

bool diffcmp(int c0,int c1){
  int d0=diffptr[c0];
  int d1=diffptr[c1];
  return list[d0]-list[c0]>list[d1]-list[c1];
}

int main(int ac,char**av){
  int a,c,i,n;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint sumval,diffval;

  list.resize(n);
  diffptr.resize(n);
  sumptr.resize(n);

  priority_queue<int,vector<int>,decltype(&sumcmp)> sumnext(sumcmp);
  priority_queue<int,vector<int>,decltype(&diffcmp)> diffnext(diffcmp);

  for(i=0;i<n;i++)list[i]=bigint(i)*i*i*i;
  for(a=1;a<n;a++){sumptr[a]=a;sumnext.push(a);}
  for(c=0;c<n-1;c++){diffptr[c]=c+1;diffnext.push(c);}

  while(!sumnext.empty()&&!diffnext.empty()){
    a=sumnext.top();sumval=list[sumptr[a]]+list[a];
    c=diffnext.top();diffval=list[diffptr[c]]-list[c];
    if(sumval==diffval){printf("BINGO %d^4+%d^4+%d^4=%d^4\n",a,sumptr[a],c,diffptr[c]);fflush(stdout);}
    if(sumval<=diffval){
      sumnext.pop();
      sumptr[a]++;if(sumptr[a]<n)sumnext.push(a);
    }else{
      diffnext.pop();
      diffptr[c]++;if(diffptr[c]<n)diffnext.push(c);
    }
  }

}
ответил Alex Selby 26 +03002016-10-26T20:56:46+03:00312016bEurope/MoscowWed, 26 Oct 2016 20:56:46 +0300 2016, 20:56:46
21

Считать источник

Поскольку вы уже определили автора, почему бы не прочитать оригинальная бумага ?

Некоторые заблуждения

Фактически, числа, найденные Elkies, не были указаны в вопросе. Они были найдены Роджером Фрием с помощью компьютера. Как он это сделал еще в 1988 году? Это упоминается в документе:

  

Postscript.

     

В то время как наш первый контрпример   $$ (A, B, C; D) = (2682440,15365639,18796760; 20615673) $$   к гипотезе Эйлера по-прежнему выходит за рамки разумного исчерпывающего компьютерного поиска, осталась возможность того, что такие поиски могут найти более мелкие решения. Вскоре после того, как услышал о первом решении, Роджер Фрай из корпорации «Мышление Машины» спросил, было ли это минимально; Я не знал, но предложил, как можно исчерпывающе искать меньшие решения: устраняя общие факторы и переставляя \ $ A \ $, \ $ B \ $, \ $ C \ $, если необходимо, мы можем взять \ $ D \ $ нечетные и не делится на 5, а \ $ C <D \ $, что \ $ D ^ 4 - C ^ 4 \ $ делится на 625 и удовлетворяет нескольким другим свойствам конгруэнции и делимости, и для каждого такого \ $ D \ $ и \ $ C \ $ искать представление \ $ D ^ 4 - C ^ 4 \ $ as \ $ A ^ 4 + B ^ 4 \ $ с \ $ A \ $, \ $ B \ $ делимым на 5. Фрай перевел это в компьютерную программу и запустил ее на различных машинах соединения в течение примерно 100 часов, чтобы найти минимальный контрпример к гипотезе Эйлера:   $$ 95800 ^ 4 + 217519 ^ 4 + 414560 ^ 4 = 422481 ^ 4. $$   Он продолжил поиск и обнаружил, что это решение уникально в диапазоне \ $ D <10 ^ 6 \ $. Это решение появляется при параметризации (6) с \ $ (m, n) = (20, -9) \ $. Мы включаем результат Фрая с его разрешения.

Источник: Elkies, Noam D. «On \ $ A ^ 4 + B ^ 4 + C ^ 4 = D ^ 4 \ $." Математика вычислений (1988): 825-835.

Для тех, кто не знаком с поздним большим компанией Thinking Machines Corporation , это была компания в 1980-х и 1990-х годах, которая производила массово параллельных компьютеров.

Взятые вместе, эти намеки показывают, что многопоточный подход, использующий исчерпывающий поиск с этими теоретическими теоретическими числами, определенно может быть выполнен на современных машинах в течение нескольких часов, а не лет.

ответил Alex Selby 26 +03002016-10-26T20:56:46+03:00312016bEurope/MoscowWed, 26 Oct 2016 20:56:46 +0300 2016, 20:56:46
18

Я рекомендую переписать проблему как нахождение \ $ a, b, c, d \ $ такой, что \ $ a ^ 4 + b ^ 4 = d ^ 4 - c ^ 4 \ $. Теперь вы можете действовать только парами сил. Построение таблицы сумм принимает \ $ O (n ^ 2) \ $. Построение таблицы различий также принимает \ $ O (n ^ 2) \ $. Сортировка их принимает \ $ O (n ^ 2 \ log {n}) \ $. Соответствующие отсортированные таблицы принимают \ $ O (n ^ 2) \ $ в линейном слиянии.

Общая сложность, очевидно, \ $ O (n ^ 2 \ log {n}) \ $.

ответил vnp 25 +03002016-10-25T22:45:08+03:00312016bEurope/MoscowTue, 25 Oct 2016 22:45:08 +0300 2016, 22:45:08
12

Это другой метод из другого ответа (приоритетной очереди), который я представил ранее. Это проще и достигает \ $ O (n ^ 2 \ log (n)) \ $ time и \ $ O (n) \ $ памяти с лучшим постоянным коэффициентом. Вы выбираете некоторые интервалы \ $ [A_i, A_ {i + 1}) \ $, покрывающие \ $ [0, n ^ 4) \ $, а затем отдельно (для каждого \ $ i \ $) находим \ $ a, b \ $ такие, что \ $ a ^ 4 + b ^ 4 \ в [A_i, A_ {i + 1}) \ $ и \ $ c, d \ $ такие, что \ $ d ^ 4-c ^ 4 \ в [A_i, A_ {я + 1}) \ $. Затем вы сортируете два списка из \ $ a ^ 4 + b ^ 4 \ $ и \ $ d ^ 4-c ^ 4 \ $ и проверяете их для пересечения.

Чтобы убедиться, что он работает в правильном временном порядке, вам нужно убедиться, что число \ $ a, b \ $ с \ $ a ^ 4 + b ^ 4 \ в [A_i, A_ {i + 1}) \ $ пропорциональна \ $ n \ $ и грубо не зависит от \ $ i \ $. Вы можете сделать это, используя \ $ A_i = C \ cdot n ^ 2 \ cdot i ^ 2 \ $ для некоторой константы \ $ C \ $. Оптимальное значение \ $ C \ $ будет зависеть от ваших конкретных характеристик оборудования.

Если вы это сделаете (первая программа ниже), то время работы для поиска вышеуказанного решения \ $ 95800 ^ 4 + 414560 ^ 4 + 217519 ^ 4 = 422481 ^ 4 \ $ кажется немного более часа (на моем настольный компьютер). Это использует __int128, расширение GCC, дающее 128-битные целые числа.

Для удовольствия я также модифицировал его, чтобы использовать простейшие (и наиболее полезные) модульные отношения: \ $ a, b \ $ - кратны \ $ 5 \ $ и \ $ d ^ 4-c ^ 4 \ $ кратно \ $ 625 \ $. Это вторая программа ниже, и она находит решение через четыре минуты.

Я также изменил его на многопоточность. Это просто, потому что куски \ $ [A_i, A_ {i + 1}) \ $ можно рассматривать отдельно. Это третья программа ниже, и она находит решение (на моей шестиядерной машине) примерно через 32 секунды.

Еще одна техника, которая ускоряет работу, заключается в том, чтобы свести элементы двух списков по модулю \ $ 2 ^ {64} \ $ до их сортировки и сравнения. Если существует совпадение между \ $ a ^ 4 + b ^ 4 \ $ и \ $ d ^ 4-c ^ 4 \ $ mod \ $ 2 ^ {64} \ $, тогда вам нужна дополнительная проверка, что она подлинная, что означает он возникает из-за нередуцированного матча. (Вероятно, он будет подлинным, поэтому дополнительная проверка занимает незначительное дополнительное время.) Поскольку большую часть времени тратится на сортировку, лучше сортировать 64-битные величины, чем 128 бит. Это уменьшает потребность в памяти в два раза (приблизительно) 2 и время работы примерно на 30%, что составляет около 23 секунд для первого решения на шестиугольной машине. Это четвертая программа ниже.

(Он делает предположения, такие как кастинг 128-битных чисел до 64 бит, совпадает с принятием нижних 64 бит, а long long int - 64 бита, и я оставил некоторый специальный код синхронизации Unix, так что это не все это красиво или портативно.)

Дополнительные примечания:

  • Программы 2 и 3 отредактированы для улучшения Эдварда.

  • Если вы запустите программу 4 с помощью n=3000000, она найдет новое решение \ $ 673865 ^ 4 + 1390400 ^ 4 + 2767624 ^ 4 = 2813001 ^ 4 \ $ примерно через 5 минут (на шестигранной машине). Если вы запустите его с помощью n=9000000 и MAXMEM=12e9, он найдет \ $ 1705575 ^ 4 + 5507880 ^ 4 + 8332208 ^ 4 = 8707481 ^ 4 \ $ после примерно 80 мин.

  • Может быть возможно усилить этот метод интервалов для запуска в \ $ O (n ^ 2) \ $ time и \ $ O (n ^ {2/3}) \ $ памяти с использованием непрерывных дробей. Это довольно сложно (что означает раздражающие накладные расходы с постоянным коэффициентом) и, вероятно, будет полезно только для \ $ n \ ge10 ^ 8 \ $ или так, предполагая память порядка 32 ГБ.

  • Современное состояние этой проблемы, разработанное Робертом Гербичем, записано здесь и здесь (на русском). Их программа использует гораздо больше «модульных фильтров» и использует китайскую теорему о остатках вместо «посередине», что означает, что она использует гораздо меньше памяти, а также быстрее. Они нашли решения до \ $ d = 1871713857 \ $.

Программа 1 - базовое разделение

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4).

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

vector<bigint> pow4;
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,s0,s1;
  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=bigint(min(5.*n,MAXMEM/26.));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){
    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)

    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=1;a<n&&pow4[a]<ub;a++){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=min(int(sqrt(sqrt(ub-pow4[a]))),n);
      for(b=bmin;b<bmax;b++)list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    for(c=0;c<n-1;c++){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(d=dmin;d<dmax;d++)list1.push_back(pow4[d]-pow4[c]);
    }
    std::sort(list1.begin(),list1.end());

    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }
  }
}

Программа 2 - в качестве программы 1, но также используйте сокращение по модулю 5

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

#define MAXMEM 4e9 // Allow the use of up to (approx) this memory in bytes

vector<bigint> pow4;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(int n,bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

int main(int ac,char**av){
  int a,b,c,d,i,j,n,r,s0,s1;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  bigint chunk,chunksize;
  vector<bigint> list0,list1;

  pow4.resize(n);

  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint((MAXMEM-2*n*sizeof(int))*2.8));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));
  for(chunk=0;chunk<bigint(n)*n;chunk+=chunksize){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=chunk*chunk;
    bigint ub=min(chunk+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(chunk/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(n,list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

Программа 3 - как программа 2, но также использует многопоточность

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.
// Parallelise over [A_i,B_i).

// Compile with -pthread flag

#include <cstdio>
#include <cmath>
#include <cassert>
#include <algorithm>
#include <vector>
#include <thread>
#include <mutex>
using std::vector;
using std::min;
using std::max;

//typedef long long int bigint;
typedef __int128 bigint;

int n;// Maximum value of d
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
const int nthreads=12;
bigint chunk,chunksize;
vector<bigint> pow4;
std::mutex chunkserver,printer;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}

void printdecode(bigint cv){// Find a,b,c,d given the common value, cv=a^4+b^4=d^4-c^4
  int a,b,c,d;
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  assert(a<n);
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(pow4[d]-pow4[c]==cv)break;
  }
  assert(c<n-1);
  std::lock_guard<std::mutex> lock(printer);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4\n",a,b,c,d);fflush(stdout);
}

bigint nextchunk(){// Returns next chunk, or -1 if none left
  std::lock_guard<std::mutex> lock(chunkserver);
  if(chunk>=bigint(n)*n)return -1;
  bigint chunkcopy=chunk;
  chunk+=chunksize;
  return chunkcopy;
}

void searcher(){
  int a,b,c,d,i,j,r,s0,s1;
  bigint ch;
  vector<bigint> list0,list1;
  while((ch=nextchunk())>=0){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=ch*ch;
    bigint ub=min(ch+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=min(int(sqrt(sqrt(lb+pow4[c]))),n);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(ch/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){
        printdecode(list0[i]);
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

int main(int ac,char**av){
  int i,t;

  n=500000;
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}

  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(300)*n,bigint(MAXMEM*2.8/nthreads));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  vector<std::thread> thr(nthreads);
  for(t=0;t<nthreads;t++)thr[t]=std::thread(searcher);
  for(t=0;t<nthreads;t++)thr[t].join();
}

Программа 4 - как программа 3, но уменьшает 128-битные целые числа до 64 бит

// Look for a^4+b^4+c^4 = d^4 by comparing the sorted list of a^4+b^4 in the range [A,B)
// with the sorted list of d^4-c^4 in the range [A,B). Use a suitable collection of
// [A_i,B_i) covering the range [0,n^4). Use a=b=0 mod 5, d^4==c^4 mod 625.
// Parallelise over [A_i,B_i).

// http://codereview.stackexchange.com/questions/145221/disproving-euler-proposition-by-brute-force-in-c

// Compile with -pthread flag

#include <cstdio>
#include <cmath>
#include <ctime>
#include <sys/time.h>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <vector>
#include <thread>
#include <mutex>
using std::vector;
using std::min;
using std::max;

typedef __int128 bigint;// GCC extension: 128 bit integer
typedef unsigned long long int medint;// Need this to be 64 bits

int n;// Maximum value of d
double tim0,tim1,dtim;
const double MAXMEM=4e9; // Allow the use of up to (approx) this memory in bytes
const int nthreads=12;
bigint chunk,chunksize;
vector<bigint> pow4;
std::mutex chunkserver,printer;

int gcd(int x,int y){if(y==0)return x; return gcd(y,x%y);}
double cpu(){return clock()/double(CLOCKS_PER_SEC);}
double realtime(){
  struct timeval tv;
  gettimeofday(&tv,0);
  return tv.tv_sec+tv.tv_usec/1e6;
}

void printdecode(bigint cv){// Find a,b,c,d given a possible common value, cv=a^4+b^4=d^4-c^4
  int a,b=0,c,d=0;// Initialise to prevent unnecessary compiler warnings
  for(a=1;a<n;a++){
    bigint b4=cv-pow4[a];
    if(pow4[a]<=b4){
      b=int(round(sqrt(sqrt(b4))));
      if(pow4[b]==b4)break;
    }
  }
  if(a==n)return;
  for(c=0;c<n-1;c++){
    d=int(round(sqrt(sqrt(cv+pow4[c]))));
    if(d>=n)return;
    if(pow4[d]-pow4[c]==cv)break;
  }
  if(c==n-1)return;
  std::lock_guard<std::mutex> lock(printer);
  if(gcd(a,gcd(b,gcd(c,d)))==1)printf("BINGO"); else printf("MULTIPLE");
  printf(" %d^4+%d^4+%d^4=%d^4 @ CPU=%.2fs Real=%.2fs\n",a,b,c,d,cpu(),realtime()-tim0);fflush(stdout);
}

bigint nextchunk(){// Returns next chunk, or -1 if none left
  std::lock_guard<std::mutex> lock(chunkserver);
  if(chunk>=bigint(n)*n)return -1;
  bigint chunkcopy=chunk;
  chunk+=chunksize;
  {
    double tim2=realtime()-tim0;
    if(tim2>=tim1){
      tim1=tim2+dtim;dtim=min(1.2*dtim,100.);
      double adjchunk=chunk-(nthreads+1)/2.*chunksize;// Estimated chunkage completed
      double p=adjchunk/double(bigint(n)*n);// Approx proportion complete (not strictly correct)
      printf("%11.2f   %8.3f%%   %12g   %11.2f  %11.2f\n",tim2,p*100,sqrt(max(adjchunk,0.)),tim2/p,tim2/p*(1-p));fflush(stdout);
    }
  }
  return chunkcopy;
}

void searcher(){
  int a,b,c,d,i,j,r,s0,s1;
  bigint ch;
  vector<medint> list0,list1;
  while((ch=nextchunk())>=0){

    // We will try a^4+b^4 in the range [chunk^2,(chunk+chunksize)^2)
    bigint lb=ch*ch;
    bigint ub=min(ch+chunksize,bigint(n)*n);ub*=ub;

    list0.resize(0);
    for(a=5;a<n&&pow4[a]<ub;a+=5){
      int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
      int bmax=int(sqrt(sqrt(ub-pow4[a])));
      r=(-bmin)%5;if(r<0)r+=5;
      for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1))list0.push_back(pow4[a]+pow4[b]);
    }
    std::sort(list0.begin(),list0.end());

    list1.resize(0);
    int rou[4]={1,182,624,443};// Fourth roots of unity mod 625
    for(c=1;c<n-1;c++)if(c%5){
      int dmin=max(int(sqrt(sqrt(lb+pow4[c]))),c+1);
      int dmax=min(int(sqrt(sqrt(ub+pow4[c]))),n);
      for(i=0;i<4;i++){
        // Need to find d in [dmin,dmax) congruent to rou[i]*c mod 625
        r=(rou[i]*(c%625)-dmin)%625;
        if(r<0)r+=625;
        for(d=dmin+r;d<dmax;d+=625)list1.push_back(pow4[d]-pow4[c]);
      }
    }
    std::sort(list1.begin(),list1.end());

    //printf("%10d   %10lu   %10lu\n",int(ch/chunksize),list0.size(),list1.size());
    i=j=0;
    s0=list0.size();s1=list1.size();
    while(i<s0&&j<s1){
      if(list0[i]==list1[j]){// Match modulo 2^64: find full 128-bit value to see if this is a genuine match
        for(a=5;a<n&&pow4[a]<ub;a+=5){// Ugly code repeat
          int bmin=max(int(sqrt(sqrt(max(lb-pow4[a],bigint(0))))),a);
          int bmax=int(sqrt(sqrt(ub-pow4[a])));
          r=(-bmin)%5;if(r<0)r+=5;
          for(b=bmin+r;b<bmax;b+=5)if(!(a&b&1)&&medint(pow4[a]+pow4[b])==list0[i])printdecode(pow4[a]+pow4[b]);
        }
        j++;continue;
      }
      if(list0[i]<list1[j])i++; else j++;
    }

  }
}

int main(int ac,char**av){
  int i,t;

  n=(ac>1?atoi(av[1]):500000);
  printf("Using n = %d\n",n);
  if(4*log(n)/log(2)+2>sizeof(bigint)*8){fprintf(stderr,"Error: %lu-bit type not large enough\n",sizeof(bigint)*8);exit(1);}
  assert(sizeof(medint)>=8);

  tim0=realtime();
  tim1=0;dtim=1;
  pow4.resize(n);
  for(i=0;i<n;i++)pow4[i]=bigint(i)*i*i*i;

  chunk=0;chunksize=min(bigint(600)*n,bigint(MAXMEM*5.8/nthreads));
  printf("Using chunksize = %lld\n",(long long int)(chunksize));

  vector<std::thread> thr(nthreads);
  for(t=0;t<nthreads;t++)thr[t]=std::thread(searcher);
  for(t=0;t<nthreads;t++)thr[t].join();
  printf("Total CPU %.2fs\n",cpu());
  printf("Total real %.2fs\n",realtime()-tim0);
}
ответил Alex Selby 2 32016vEurope/Moscow11bEurope/MoscowWed, 02 Nov 2016 04:10:05 +0300 2016, 04:10:05
8

Умножение не требуется.

Рассмотрим, как вычислить серию квадратов 1,4,9,16

unsigned sq = 1;
for (unsigned x=1;; x++) {
  sq += x + x + 1;
  printf("%u\n", sq);
}

Рассмотрим, как вычислить серию кубов 1,8,27,64

unsigned sq = 1;
unsigned cube = 1;
for (unsigned x=1;; x++) {
  cube += (sq + sq + sq) + (x + x + x) + 1;
  sq += x + x + 1;
  printf("%u\n", cube);
}

Аналогичная операция для 4-й мощности.


Более широкая, но не произвольно более широкая математическая необходимость

Чтобы достичь a = 95800; b = 217519; c = 414560; d = 422481, требуется 80-битная математика. Предложите использовать компилятор, который поддерживает uint128_t или эквивалент.


Просто integer требуется дополнение и сравнение

Я бы постепенно вычислял a**4 + b**4 + c**4 и d**4, увеличивая меньшую из двух. Конечно, когда они равны - вуаля!


Пример кода. Используйте код factor = 1 для окончательного кода. Используйте factor >= 7 для тестирования.

Использование коэффициента factor позволяет коду запускать только часть значений, необходимых для тестирования. Фактор factor из 50 занял несколько минут. Фактор factor из 20 занял 4 часа. Ожидайте, что коэффициент factor из 1 займет годы (5?). Когда коэффициент factor < 7, вам понадобится тип uint128t, который на самом деле составляет 128 бит, в отличие от используемого здесь владельца места uintmax_t.

#include <stdint.h>
#include <stdlib.h>
#include <stdio.h>

// Adjust to some 128 bit type
// The below is a place-holder for a true 128-bit type
typedef uintmax_t uint_fast128t;

typedef struct {
  uint_fast128t x4;  // 4th power of the value
  uint_fast128t x3;  // cube of the value
  uint_fast64_t x2;  // square of the value
  uint_fast32_t x1;  // location of OP's value a,b,c or d
                     // Assume a,b,c,d are 20 bits or less.
}quad;

#define quad_next(x) { \
  x.x4 += (uint_fast128t) x.x3*4 + x.x2*6 + 4*x.x1 + 1; \
  x.x3 += (uint_fast128t) x.x2*3 + x.x1*3 + 1; \
  x.x2 += (uint_fast64_t) x.x1*2 + 1; \
  x.x1++; \
  }

#define factor 20
uint_fast32_t a_max = 95800 / factor;
uint_fast32_t b_max = 217519 / factor;
uint_fast32_t c_max = 414560 / factor;
uint_fast32_t d_max = 422481 / factor;

void Elkies_B(quad a, quad b) {
  quad c = b;
  quad d = c;
  c.x4 += b.x4;

  // Could use a more efficient step here
  while (d.x4 < c.x4) {
    quad_next(d);
  }
  while (c.x1 <= c_max) {
    while (d.x4 < c.x4) {
      quad_next(d);
    }
    if (d.x4 == c.x4) {
      // Success! Print a,b,c,d
      printf("+ a:%lu b:%lu c:%lu d:%lu\n", //
              (unsigned long) a.x1, (unsigned long) b.x1, //
              (unsigned long) c.x1, (unsigned long) d.x1);
      exit(0);
    }
    quad_next(c);
  }
}

void Elkies_A(quad a) {
  quad b = a;
  b.x4 += a.x4;
  while (b.x1 <= b_max) {
    Elkies_B(a, b);
    quad_next(b);
  }
}

void Elkies(void) {
  quad a = { 1, 1, 1, 1 };
  while (a.x1 <= a_max) {
    printf("- a:%lu %llu %ju %ju\n", //
            (unsigned long) a.x1, (unsigned long long) a.x2, //
            a.x3, a.x4);
    fflush( stdout);
    Elkies_A(a);
    quad_next(a);
  }
  printf("Failed\n");
}

int main() {
  Elkies();
}
ответил chux 26 +03002016-10-26T02:04:55+03:00312016bEurope/MoscowWed, 26 Oct 2016 02:04:55 +0300 2016, 02:04:55
7

Я не уверен, сколько практического использования это, но так же как d и точно один из \ $ a \ $, \ $ b \ $ и \ $ c \ $, который должен быть нечетным, мы можем использовать modulo 5 и 10 свойств.

  1. Modulo 5, четвертые степени равны либо 0, либо 1. Единственный способ объединить сумму из трех из них, чтобы сделать четвертый, состоит в том, что два из трех на LHS равны 0 и с обоими \ $ d ^ 4 \ $ и один из \ $ a ^ 4 \ $, \ $ b ^ 4 \ $ и \ $ c ^ 4 \ $, являющийся 1 по модулю 5 (т. Е. D заканчивается только в 1, 3, 7 или 9).

  2. В базе 10 четвертые полномочия заканчиваются либо в 0000, E1, U6 или 625, где E - четная цифра, а U - нечетная цифра. Рассматривая способы их комбинирования и исключая \ $ 0 + 0 + 5 = 5 \ $ и \ $ 0 + 0 + 6 = 6 \ $, как не в их младших членах, имеем, что конечные разряды четвертых степеней удовлетворяют либо некоторым перестановка \ $ 0 + 0 + 1 = 1 \ $ (как в решении Ноама Элкиса) или некоторая перестановка \ $ 0 + 5 + 6 = 1 \ $.

Таким образом, на LHS должно быть два кратных 5, и по крайней мере один из них должен быть кратным 10 (в решении Elkies они оба). Это тоже должно сократить поиск.

Я не вижу преимущества при просмотре любых других модулей.

Также может быть небольшой выигрыш в тестировании, только если RHS - это просто идеальный квадрат (необходимое условие четвертой мощности), и только если он удовлетворяет этому тесту, тогда тестирование, если оно также является квадратом квадрата , Любое преимущество может зависеть от того, будет ли ваш подход к тестированию с \ $ d \ $ и \ $ d ^ 4 \ $ сверху вниз (от root до root) или снизу вверх (root до его мощности).

ответил Staycator 28 +03002016-10-28T05:36:15+03:00312016bEurope/MoscowFri, 28 Oct 2016 05:36:15 +0300 2016, 05:36:15
7

Решение, о котором вы говорите, не было найдено Elkies, но Фрай. Это описано в [Elkies] :

  

Эйлер предположил в 1769 году, что диофантово уравнение   \ $ A ^ 4 + B ^ 4 + C ^ 4 = D ^ 4 \ $, (...) не имеет решения в натуральных числах.

  

(...), устраняя общие факторы и, если необходимо, переставляя \ $ A, B, C \ $, мы можем взять \ $ D \ $ нечетные и не делимые на \ $ 5 \ $ и \ $ C \ lt D \ $ такое, что \ $ D ^ 4-C ^ 4 \ $ делится на \ $ 625 \ $ и удовлетворяет нескольким другим свойствам конгруэнции и делимости, и для каждого такого D и C ищут представление D ^ 4-C ^ 4 как \ $ A ^ 4 + B ^ 4 \ $, причем \ A, B \ $ делится на \ $ 5 \ $. Фрай перевел это в компьютерную программу и запустил (...) около 100 часов [в 1988 году], чтобы найти минимальный контрпример к гипотезе Эйлера

Роджер Фрай опубликовал подробное описание своего алгоритма здесь: [Frye] .

Первое исследование (для значений менее 10000) было выполнено [Уорд] . Он не использовал компьютер. Статья не доступна бесплатно. Но первую страницу, содержащую важную теорему 1 и две, можно найти здесь и можно найти резюме, на которое ссылается здесь , здесь .

Здесь было опубликовано другое расследование (теперь компьютер) для значений менее 220 000: [Lander, Parkin, Selfridge]

Более поздние исследования: [Ekl] и [Bernstein] .

Самый продвинутый поиск (значения меньше 2 000 000), который я нашел, опубликован на странице страницы Роберта Гербича : это C-программа euler413.c .

В настоящий момент (2016-11-04) файл ошибочно содержит программу дважды. Удалите все строки из строки 1478 до конца. Суб>

Примечание от @AlexSelby: Чтобы заставить код работать в 64-битной системе, вам нужно изменить пять строк multipliers17=(unsigned int**) (malloc) (17*17*sizeof(unsigned int)); и так далее , заменив sizeof(unsigned int) на sizeof(unsigned int*) (и аналогично для следующих четырех строк)

[Bernstein] , Bernstein, DJ; ЭНЕРГЕТИЧЕСКИЕ РЕШЕНИЯ К p (a) + q (b) = r (c) + s (d); МАТЕМАТИКА ВЫЧИСЛЕНИЯ, том 70, N0 233, стр. 389-394

[Ekl] , Ekl, RL; НОВЫЕ РЕЗУЛЬТАТЫ В РАВНОМЕРНЫХ СУММАХ КАК ПОЛНОМОЧИЙ; МАТЕМАТИКА ВЫЧИСЛЕНИЯ, том 67, № 223, июль 1998, стр. 1309-1315

[Elkies] Elkies, Noam D .; On \ $ A ^ 4 + B ^ 4 + C ^ 4 = D ^ 4 \ $; МАТЕМАТИКА ВЫЧИСЛЕНИЯ, 51/184, ОКТ. 1988, 825-835

[Фрай] Фрай, Роджер; Поиск $ 95800 ^ 4 + 217519 ^ 4 + 414560 ^ 4 = 422481 ^ 4 $ на машине подключения

[Lander, Parkin, Selfridge] Lander, L. J .; Parkin T.R .; Selfridge, J.L .; Обзор равных сумм подобных держав; Математика Вычисления, том 21, 1967

[Уорд] Уорд, Морган. Эйлера по сумме трех четвертых степеней. Герцог Матх. J. 15 (1948), вып. 3, 827--837
здесь - сводка и здесь является первой страницей статьи Wards


[Robert Gerbicz 'euler413.c] https://sites.google.com/сайт /robertgerbicz /euler413.c

[Robert Gerbicz 'euler413.exe] https://sites.google.com/site/robertgerbicz/euler413.exe



Просто для удовольствия!

Скомпилированную версию для Windows можно также найти на странице Robert Gerbicz : euler413.exe . Я попробовал. Я ничего не знаю о программе, ее параметрах и ее выходе. Я получил следующее:

C: & gt; euler413.exe
Я не нашел незавершенной работы!
Дайте параметр R, Range будет R * 10240000: 10
Вы хотите начать исчерпывающий поиск в этом диапазоне или
только для поиска специальных решений?
(y = полный поиск, n = специальный поиск) n
Дайте начальное значение a0. Он должен быть 0 & lt = start_a0 <lt = 16384: 0
Дайте конечное значение a0. Это должно быть start_a0 & lt = end_a0 <lt = 16384: 16384
Создание таблиц
Готово
Тестирование: a0 = 1
Закончено: a0 = 1, Диапазон = 102400000, тип = 0, Время: 0h0m0s, Дата: Пт ноя 04 18:20:33 2016
Тестирование: a0 = 9
Закончено: a0 = 9, Диапазон = 102400000, тип = 0, Время: 0h0m1s, Дата: Пт ноя 04 18:20:34 2016
....

Программа создает много вывода на экран и следующие файлы:

euler413work.txt
stat_euler (4,3,1) .txt
results_euler (4,1,3) .txt

euler413work.txt - это файл состояния. Если вы остановите программу и запустите ее позже в том же каталоге, она продолжит работу в тот момент, когда вы ее остановили. stat_euler(4,3,1).txt записывает вывод времени, который вы также можете увидеть на экране. В results_euler(4,1,3).txt найденное решение сохраняется.

Примерно через 9 минут файл результатов содержал следующие строки:

Решение найдено! 98438073 ^ 4 = 50681927 ^ 4 + 96592480 ^ 4 + 22321400 ^ 4
Решение найдено! 17321721 ^ 4 = 8918279 ^ 4 + 16996960 ^ 4 + 3927800 ^ 4
Решение найдено! 40980657 ^ 4 = 21099343 ^ 4 + 40212320 ^ 4 + 9292600 ^ 4
Решение найдено! 20615673 ^ 4 = 15365639 ^ 4 + 18796760 ^ 4 + 2682440 ^ 4 (примитивный)
Решение найдено! 20615673 ^ 4 = 15365639 ^ 4 + 2682440 ^ 4 + 18796760 ^ 4 (примитивный)
Решение найдено! 7182177 ^ 4 = 3697823 ^ 4 + 7047520 + ^ 4 ^ 4 1628600
Решение найдено! 30841113 ^ 4 = 15878887 ^ 4 + 30262880 ^ 4 + 6993400 ^ 4
Решение найдено! 101817921 ^ 4 = 52422079 ^ 4 + 99908960 ^ 4 + 23087800 ^ 4
Решение найдено! 12197457 ^ 4 = 8282543 ^ 4 + 11289040 ^ 4 + 5870000 ^ 4 (примитивный)
Решение найдено! 34220961 ^ 4 = 17619039 ^ 4 + 33579360 ^ 4 + 7759800 ^ 4
Решение найдено! 57879897 ^ 4 = 29800103 ^ 4 + 56794720 ^ 4 + 13124600 ^ 4
Решение найдено! 81538833 ^ 4 = 41981167 ^ 4 + 80010080 ^ 4 + 18489400 ^ 4
Решение найдено! 16003017 ^ 4 = 4479031 ^ 4 + 14173720 ^ 4 + 12552200 ^ 4 (примитивный)
Решение найдено! 422481 ^ 4 = 217519 ^ 4 + 414560 ^ 4 + 95800 ^ 4 (примитивный)
Решение найдено! 16430513 ^ 4 = 16281009 ^ 4 + 3642840 ^ 4 + 7028600 ^ 4 (примитивный)
Решение найдено! 47740353 ^ 4 = 24579647 ^ 4 + 46845280 ^ 4 + 10825400 ^ 4
Решение найдено! 37600809 ^ ​​4 = 19359191 ^ 4 + 36895840 ^ 4 + 8526200 ^ 4
ответил miracle173 30 +03002016-10-30T22:48:36+03:00312016bEurope/MoscowSun, 30 Oct 2016 22:48:36 +0300 2016, 22:48:36
6

Я думаю, что у меня есть решение \ $ O (n ^ 2) \ $, по существу нуждающееся в \ $ O (n ^ 2) \ $ пространстве, но получающее очень большой постоянный коэффициент (~ 100) усадки. Это должно сделать его реалистичным для первого \ $ d = 422481 \ $.

Использование нечетных и четных пределов Дан Узнанского

$$ a ^ 4 + b ^ 4 + c ^ 4 = d ^ 4 $$

(\ $ a \ $ и \ $ b \ $ четные, \ $ c \ $ и \ $ d \ $ нечетны)

Использование материалов стиля Staycator

\ $ 0 + 0 + 1 = 1 \ $, \ $ a ^ 4 \ $ и \ $ b ^ 4 \ $ end в 0000 и \ $ c \ $ и \ $ d \ $ end в 1,3, 7 или 9.

Так как \ $ a ^ 4 + b ^ 4 = d ^ 4 - c ^ 4 \ $, \ $ c \ $ и \ $ d \ $ заканчиваются теми же четырьмя цифрами.

\ $ 0 + 6 + 5 = 1 \ $, \ $ a ^ 4 \ $ заканчивается в 0000 и \ $ c ^ 4 \ $ в 0625. \ $ b \ $ заканчивается на 2, 4, 6 или 8 и d в 1, 3, 7 или 9.

Так как \ $ a ^ 4 + c ^ 4 = d ^ 4 - b ^ 4 \ $, то последние четыре цифры \ $ d ^ 4 - b ^ 4 \ $ равны 0625.

Решение

(Я еще не помещал это в код, только некоторые заметки в Блокноте)

  • Предварительно расчитайте \ $ n ^ 4 \ $ для \ $ n \ $ 1 до 500k
  • Бин по последней цифре \ $ n \ $
  • Bin \ $ n ^ 4 \ $, где n заканчивается на 1, 3, 7 или 9 четырьмя четырьмя цифрами \ $ n ^ 4 \ $
  • Bin \ $ n ^ 4 \ $, где \ $ n \ $ заканчивается на 2, 4, 6 или 8 четырьмя четырьмя цифрами \ $ n ^ 4 \ $
  • Сохранить хэш-таблицу \ $ a ^ 4 + b ^ 4 \ $, где оба \ $ a \ $ и \ $ b \ $ заканчиваются в нуле
  • Зациклируйте последние четыре бина и сравните возможные хэши \ $ d ^ 4 - c ^ 4 \ $ с сохраненными \ $ a ^ 4 + b ^ 4 \ $
  • Сохранить хэш-таблицу \ $ a ^ 4 + c ^ 4 \ $, где \ $ a ^ 4 \ $ заканчивается в 0000 и \ $ c ^ 4 \ $ в 0625
  • Найти бины, где нечетные \ $ n \ $ минус даже \ $ n \ $ bins = 0625 (обертывание). сравнить возможно \ $ d ^ 4 - b ^ 4 \ $ хэши для хранения \ $ a ^ 4 + c ^ 4 \ $.
ответил JollyJoker 28 +03002016-10-28T20:13:04+03:00312016bEurope/MoscowFri, 28 Oct 2016 20:13:04 +0300 2016, 20:13:04
4

Возьмем идею вычисления \ $ a ^ 4 + b ^ 4 \ $ и \ $ d ^ 4-c ^ 4 \ $ и нахождения значений, равных. Предположим, что \ $ a \ ge b \ $ и \ $ d \ gt c \ $.

Создайте две очереди приоритетов: одну для значений \ $ a ^ 4 + b ^ 4 \ $ и одну для значений \ $ d ^ 4-c ^ 4 \ $, но также отслеживание, которое \ $ b \ $ и \ $ c \ $. Заполните очереди приоритетов значениями, где \ $ b = 0 \ $, \ $ c = d-1 \ $ и для каждого \ $ a, d \ $ отслеживать, какое число вы поместили в очереди (изначально \ $ b = 0 \ $, \ $ c = d-1 \ $). Если наименьшее число в каждой очереди одинаково, у вас есть решение. Если нет, вы удаляете, например, (\ $ a, b \ $) из очереди, и вставляете (\ $ a, b + 1 \ $), если \ $ b + 1 \ gt a \ $; или вы удаляете (\ $ d, c \ $) из очереди и вставляете (\ $ d, c-1 \ $), если c = 0.

Цифры будут большими , а 64-разрядных целых чисел может быть недостаточно. Но если ваш результат не будет больше, чем say \ $ 2 ^ {112} \ $ (a, b до \ $ 2 ^ {28} \ $), вы можете вычислить результат один раз в двойной точности и один раз с использованием 64-битного знака. Двойная точность даст вам грубо правильный результат. 64-разрядный бит без знака даст вам точный последние 64 бит результата. Так как \ $ d, c \ $ будет больше, вычислим \ $ d ^ 4 - c ^ 4 \ $ при \ $ (d ^ 2 + c ^ 2) (d ^ 2-c ^ 2) \ $.

Это займет просто \ $ O (n) \ $ пространство, поэтому вы вряд ли сможете вычислить достаточно далеко, чтобы это стало проблемой.

ответил gnasher729 27 +03002016-10-27T11:45:24+03:00312016bEurope/MoscowThu, 27 Oct 2016 11:45:24 +0300 2016, 11:45:24
3

В духе фактической оптимизации самого алгоритма при игнорировании очень реальных проблем с переполнением, о которых справедливо упоминали другие, подумайте о том, что здесь происходит:

for (a = 1; a < 100000; a++) {
    for (b = 1; b < 300000; b++) {
        for (c = 1; c < 500000; c++) {
            for (d = 1; d < 500000; d++) {
                if (prop(a, b, c, d))
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

a, b, c и d все подняты до 4-й мощности на каждой итерации самого внутреннего цикла, даже если на d было изменено.

Кроме того, хотя добавление является менее дорогостоящей операцией, оно добавляет (не каламбур), когда вы делаете это миллионы раз без причины, поэтому нет необходимости суммировать левую часть уравнения на каждой итерации.

Это значительно сократило бы время вычисления:

for (a = 1; a < 100000; a++) {
    long int a4 = a * a * a * a;
    for (b = 1; b < 300000; b++) {
        long int a4b4 = a4 + b * b * b * b;
        for (c = 1; c < 500000; c++) {
            long int a4b4c4 = a4b4 + c * c * c * c;
            for (d = 1; d < 500000; d++) {
                if (a4b4c4 == d * d * d * d)
                    printf("FOUND IT!\na = %ld\nb = %ld\nc = %ld\nd = %ld\n", a, b, c, d);
            }
        }
    }
}

Я не нашел времени, чтобы оценить это, поэтому я не могу точно сказать, какие улучшения вы увидите, но я оставлю это как упражнение для читателя.

ответил Jaquez 31 +03002016-10-31T05:58:23+03:00312016bEurope/MoscowMon, 31 Oct 2016 05:58:23 +0300 2016, 05:58:23
2
  

Какие стратегии я мог бы предпринять, чтобы сделать вышеприведенный код быстрее?

Хранители времени

  1. Легко: узнайте, не более, комбинации, a = 1 to 95800, b = a to 217519, c = b to 414560 и d = c+1 or more. Использование комбинаций b = 1 to a-1 и c = 1 to b-1 является избыточным. d = 1 to c не может решить уравнение. @phoog

  2. После тестовых значений для a & b, значения c и d могут быть протестированы очень эффективно, увеличивая c или d, заметив нижнюю часть sum_ab + c**4 или d**4.

  3. Обратите внимание, что если a, b, c, d является решением, n*a, n*b, n*c, n*d также решение. Это означает, что код не должен проверять решение, где a,b,c,d все равно, как решение с нечетным значением было бы найдено с меньшими значениями.

  4. Каждая степень 4 целого числа имеет младшие значащие разряды 0000 или 0001. Добавление любых 3 из них должно быть равно 0 или 1 (наименьший существенный полубайт d). Это подразумевает, что все четные (которые мы можем игнорировать сверху) OR 1 из a,b,c должно быть нечетным, остальные 2 четные и d должны быть нечетное . Использование этого короткого замыкания линейно сокращает время выполнения.

  5. Задача состоит в использовании целочисленной математики, превышающей 64 бит. C не указывает стандартный тип шире, чем 64-битный. Различные платформы имеют более широкие типы, такие как __uint128 на gcc.

O (n 3 ). O (n 1 ).
Приведенный ниже код демонстрирует четкий тренд O (n 3 ) во времени, поскольку все больше и больше диапазона a,b,c,d предпринимаются.

Оцените завершение на моей машине 100% -ного прохода: 5,900,000 секунд (69 дней).

// Ref 
a = 95800; b = 217519; c = 414560; d = 422481


#include <assert.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define a (95800 + 1)
#define b (217519 + 1)
#define c (414560 + 1)
#define d (422481 + 1)

// Adjust to some 128 bit type on your platform
// Can leave it at 64-bit if `num/den < 65535/d` or about 0.155 for initial testing
typedef uintmax_t uint128t;

#define pow2(x) ((uint128t)(x)*(x))
#define pow4(x) (pow2((uint_fast64_t)(x)*(x)))
void doit_times(unsigned num, unsigned den) {
  time_t t0;
  time(&t0);
  printf("Portion:%6g%%\t", 100.0 * num / den);
  uint_fast32_t a_max = (uint_fast32_t) round(1.0 * a * num / den);
  uint_fast32_t b_max = (uint_fast32_t) round(1.0 * b * num / den);
  uint_fast32_t c_max = (uint_fast32_t) round(1.0 * c * num / den);
  uint_fast32_t d_max = (uint_fast32_t) round(1.0 * d * num / den);
  printf("a_max:%6" PRIuFAST32 " b_max:%6" PRIuFAST32 //
      " c_max:%6" PRIuFAST32 " d_max:%6" PRIuFAST32 "\t", //
      a_max, b_max, c_max, d_max);
  uint_fast32_t q0 = 1;
  for (uint_fast32_t qa = q0; qa <= a_max; qa++) {
    assert(qa <= 0xFFFF || sizeof(uint128t) > 4);
    uint128t qa4 = pow4(qa);
    for (uint_fast32_t qb = qa; qb <= b_max; qb++) {
      uint128t sum_ab4 = qa4 + pow4(qb);
      uint_fast32_t qc = qb;
      uint128t sum_abc4 = sum_ab4 + pow4(qc);
      uint_fast32_t qd = qc + 1;
      uint128t qd4 = pow4(qd);

      // Loop until `c` to too large.
      for (;;) {
        if (sum_abc4 < qd4) {
          if (qc >= c_max) break;
          qc++;
          sum_abc4 = sum_ab4 + pow4(qc);
        } else if (sum_abc4 > qd4) {
          //if (qd >= d_max) break;
          qd++;
          qd4 = pow4(qd);
        } else {
          puts("Success!");
          printf("a %" PRIuFAST32 "\n", qa);
          printf("b %" PRIuFAST32 "\n", qb);
          printf("c %" PRIuFAST32 "\n", qc);
          printf("d %" PRIuFAST32 "\n", qd);
          fflush(stdout);
          exit(0);
        }
      }

    }
  }
  time_t t1;
  time(&t1);
  double dt = difftime(t1, t0);
  printf("t:%g\n", dt);
  fflush(stdout);
}

int main(void) {
  for (unsigned num = 1; num <= 1000; num++) {
    doit_times(num, 1000);
  }
  puts("Done");
  return -1;
}

Выход

Portion:   0.1% a_max:    96 b_max:   218 c_max:   415 d_max:   422 t:0
Portion:   0.2% a_max:   192 b_max:   435 c_max:   829 d_max:   845 t:0
Portion:   0.3% a_max:   287 b_max:   653 c_max:  1244 d_max:  1267 t:1
Portion:   0.4% a_max:   383 b_max:   870 c_max:  1658 d_max:  1690 t:1
Portion:   0.5% a_max:   479 b_max:  1088 c_max:  2073 d_max:  2112 t:2
Portion:   0.6% a_max:   575 b_max:  1305 c_max:  2487 d_max:  2535 t:3
Portion:   0.7% a_max:   671 b_max:  1523 c_max:  2902 d_max:  2957 t:5
Portion:   0.8% a_max:   766 b_max:  1740 c_max:  3316 d_max:  3380 t:6
Portion:   0.9% a_max:   862 b_max:  1958 c_max:  3731 d_max:  3802 t:10
Portion:     1% a_max:   958 b_max:  2175 c_max:  4146 d_max:  4225 t:12
Portion:   1.1% a_max:  1054 b_max:  2393 c_max:  4560 d_max:  4647 t:16
Portion:   1.2% a_max:  1150 b_max:  2610 c_max:  4975 d_max:  5070 t:21
Portion:   1.3% a_max:  1245 b_max:  2828 c_max:  5389 d_max:  5492 t:25
Portion:   1.4% a_max:  1341 b_max:  3045 c_max:  5804 d_max:  5915 t:32
Portion:   1.5% a_max:  1437 b_max:  3263 c_max:  6218 d_max:  6337 t:39
Portion:   1.6% a_max:  1533 b_max:  3480 c_max:  6633 d_max:  6760 t:48
Portion:   1.7% a_max:  1629 b_max:  3698 c_max:  7048 d_max:  7182 t:56
Portion:   1.8% a_max:  1724 b_max:  3915 c_max:  7462 d_max:  7605 t:67
Portion:   1.9% a_max:  1820 b_max:  4133 c_max:  7877 d_max:  8027 t:77
Portion:     2% a_max:  1916 b_max:  4350 c_max:  8291 d_max:  8450 t:91
Portion:   2.1% a_max:  2012 b_max:  4568 c_max:  8706 d_max:  8872 t:104
Portion:   2.2% a_max:  2108 b_max:  4785 c_max:  9120 d_max:  9295 t:120

Обратите внимание, что тестирование d в пределах лимита не требуется, поскольку c превышает лимит сначала или вскоре после этого.

ответил chux 3 42016vEurope/Moscow11bEurope/MoscowThu, 03 Nov 2016 05:21:02 +0300 2016, 05:21:02
0

Что касается проверки того, является ли \ $ D \ $ четвертой степенью, не могли бы вы вычислить, скажем, все целые 4-й степени \ $ D \ $ до предела и сохранить их в хеше Таблица? Это изменило бы время проверки любого потенциального \ $ D ^ 4 \ $ до \ $ O (1) \ $ при увеличенной стоимости префактора в начале. Очевидно, вы не могли реализовать это для произвольно большого \ $ D \ $, но OP не проверяет это (жестко установленные пределы для \ A, B, C \ $), но хэш-таблица размером ~ 0.5M хорошо в пределах разума.

И если вы проверяете ограниченные \ $ A, B, C \ $, верхний предел для \ $ D \ $ можно найти из \ $ \ lceil {\ sqrt [4] {A ^ 4 + B ^ 4 + C ^ 4}} \ rceil \ $, который будет включать только одну операцию 4-го корня для поиска (исключение проверки).

ответил M. Withnall 28 +03002016-10-28T17:20:54+03:00312016bEurope/MoscowFri, 28 Oct 2016 17:20:54 +0300 2016, 17:20:54
-1

Как насчет поиска полномочий сначала, а затем вычисления эквивалентности? Это устранит перерасчет четвертой степени того же числа (из-за соображений пространства).

Что-то вроде этого можно запрограммировать:

%from numbers 1 to 500000
for k=1:50000:500000
%for indexing
i=0;
%compute the forth powers from k to k+49999 and store the powers in a 
for j=k:k+50000-1
a[i++]=j*j*j*j;
end
for x=1:50000
for y=1:50000
for z=1:50000
for t=1:50000
  if (a[x]+a[y]+a[z]==a[t])
   return;
end
end
end
end

Также можно использовать другие методы оптимизации.

ответил tempx 27 +03002016-10-27T14:38:57+03:00312016bEurope/MoscowThu, 27 Oct 2016 14:38:57 +0300 2016, 14:38:57

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132