Project Euler #211 - проблема эффективности

Вот решение, учитывающее дух [Проекта] Эйлера. [Внимание: спойлер. Я старался делать подсказки медленными, чтобы вы могли прочитать только часть ответа и подумать самостоятельно, если хотите. :)]

Когда вы сталкиваетесь с проблемой, связанной с числами, одна из хороших стратегий (как вы, вероятно, уже знаете из решения 210 задач проекта Эйлера) состоит в том, чтобы рассмотреть небольшие примеры, найти закономерность и доказать ее. [Последняя часть может быть необязательной в зависимости от вашего отношения к математике ;-)]

Однако в этой задаче просмотр небольших примеров -- для n=1,2,3,4,..., вероятно, не даст вам никаких подсказок. Но есть еще один смысл «небольших примеров» при решении задач теории чисел, который вы, вероятно, уже знаете: простые числа — это строительные блоки натуральных чисел, так что начните с простых.

Делителями простого числа p являются только 1 и p, поэтому сумма квадратов его делителей равна 1+p².
Для степени простого числа p^{k< /sup>, его единственными делителями являются 1, p, p2, pk, поэтому сумма квадратов его делителей равна 1+p+p 2++pk=(pk+1-1)/(p-1).
Это был самый простой случай: вы Решил задачу для всех чисел только с одним простым множителем.}

Пока ничего особенного. Теперь предположим, что у вас есть число n, имеющее два простых делителя, скажем, n=pq. Тогда его делители равны 1, p, q и pq, поэтому сумма квадратов его делителей равна 1+p²+q²+p² q²=(1+p²)(1+q²).
А как насчет n=p^aq^b? Чему равна сумма квадратов его множителей?

[.................................Опасно читать ниже этой строки............... ....]

Это _{0ca, 0db}(p^cq^d)² = ((p^a+1 -1)/(p-1))((q^b+1-1)/(q-1)).

Это должно дать вам подсказку как о том, каков ответ, так и о том, как его доказать: сумма делителей n — это просто произведение (ответа) для каждой из степеней простых чисел в его факторизации, поэтому все, что вам нужно нужно разложить 64000000 на множители (что очень легко сделать даже в уме :-)) и умножить ответ на каждую (=оба, потому что единственные простые числа 2 и 5) его степени простого числа.

Это решает проблему проекта Эйлера; теперь мораль отнять от него.

Более общий факт здесь касается мультипликативных функций — функций на натуральных числах, таких что f(mn) = f(m)f(n), когда gcd(m,n)=1, т. е. m и n не имеют общих простых множителей. Если у вас есть такая функция, то значение функции в конкретном числе полностью определяется ее значениями в простых степенях (вы можете это доказать?)

Немного более сложный факт, который вы можете попытаться доказать[это не это сложно], заключается в следующем: если у вас есть мультипликативная функция f [здесь, f(n)=n^{2] и вы определяете функцию F как F(n) = _{d делит n}f(d), (как это было сделано здесь), тогда F(n) также является мультипликативной функцией.}

[На самом деле что-то очень красивое верно, но не смотрите на это просто так. тем не менее, и вам, вероятно, никогда не понадобится. :-)]

Я думаю, что ваш алгоритм не самый эффективный из возможных. Подсказка: возможно, вы начинаете не с той стороны.

редактировать: я хотел бы добавить, что выбор 64000000 в качестве верхнего предела, вероятно, является способом автора проблемы сказать вам, чтобы вы подумали о чем-то лучшем.

edit: Несколько советов по эффективности:

• вместо

(setf l (append l (...)))

вы можете использовать

(push (...) l)

который деструктивно изменяет ваш список, добавляя новую ячейку с вашим значением как car и прежним l как cdr, а затем указывает l на эту ячейку. Это намного быстрее, чем добавление, которое должно проходить по списку один раз. Если вам нужен список в другом порядке, вы можете перевернуть его после его завершения (но здесь это не требуется).

• почему вы сортируете l?

• вы можете сделать (> current (/ num current)) более эффективным, сравнив его с квадратным корнем из числа (которое нужно вычислять только один раз для каждого числа).

• возможно ли более эффективно находить множители числа?

И подсказка по стилю: вы можете указать область видимости l в объявлении do:

(do ((l ())
     (current 1 (+ current 1)))
    ((> current (/ num current))
     l)
  ...)

Я бы атаковал это, выполнив простую факторизацию числа (например: 300 = 2 ^ 2 * 3 ^ 1 * 5 ^ 2), что относительно быстро, особенно если вы генерируете это с помощью сита. Исходя из этого, относительно просто сгенерировать коэффициенты путем повторения i=0..2; j=0..1; k=0..2 и выполнение 2^i * 3^j * 5^k.

5 3 2
-----
0 0 0 = 1
0 0 1 = 2
0 0 2 = 4
0 1 0 = 3
0 1 1 = 6
0 1 2 = 12
1 0 0 = 5
1 0 1 = 10
1 0 2 = 20
1 1 0 = 15
1 1 1 = 30
1 1 2 = 60
2 0 0 = 25
2 0 1 = 50
2 0 2 = 100
2 1 0 = 75
2 1 1 = 150
2 1 2 = 300

Это может быть недостаточно быстро

Хитрый трюк, который вам не хватает, заключается в том, что вам вообще не нужно разлагать числа на множители. Сколько чисел из 1..N кратны 1? N Сколько чисел из 1..N кратны 2? Н/2

Хитрость заключается в суммировании факторов каждого числа в списке. Чтобы получить 1, добавьте 1^2 к каждому числу в списке. Чтобы получить 2, прибавьте 2^2 ко всем остальным числам. Чтобы получить число 3, прибавьте 3^2 к каждому третьему числу.

Не проверяйте делимость вообще. В конце нужно проверить, является ли сумма полным квадратом, вот и все. В C++ у меня это сработало за 58 секунд.

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

открывающая скобка

sqrt(k), чтобы найти делители k (путем пробного деления), возвести каждый из них в квадрат (постоянное время на множитель) и суммировать их (постоянное время на множитель). Это ₂(k), который я буду называть x.

плюс

не уверен, какова сложность хорошего алгоритма целочисленного квадратного корня, но, конечно, не хуже, чем sqrt (x) (тупое пробное умножение). x вполне может быть большим-O больше, чем k, поэтому я воздержусь от суждения здесь, но x, очевидно, ограничен сверху k ^ 3, потому что k имеет не более k делителей, каждый из которых не больше k, и, следовательно, его квадрат не больше, чем k ^ 2. Прошло так много времени с тех пор, как я получил степень по математике, и я понятия не имею, как быстро сходится Ньютон-Рафсон, но я подозреваю, что это быстрее, чем sqrt (x), и если все остальное не удается, двоичная отбивка - это log (x).

закрывающая скобка

умножается на n (поскольку k находится в диапазоне 1 .. n).

Итак, если ваш алгоритм хуже, чем O(n * sqrt(n^3)) = O(n ^ (5/2)), в случае немого sqrt, или O(n * (sqrt(n) + log( n ^ 3)) = O (n ^ 3/2) в случае smart-sqrt, я думаю, что что-то пошло не так, что должно быть идентифицировано в алгоритме. На этом этапе я застрял, потому что не могу отладить ваш LISP .

О, я предположил, что арифметика использует постоянное время для используемых чисел. Это чертовски хорошо должно быть для чисел размером всего 64 миллиона, и куб этого едва помещается в 64-битное целое число без знака. Но даже если ваша реализация LISP делает арифметику хуже, чем O(1), она не должна быть хуже, чем O(log n), так что это не сильно повлияет на сложность. Конечно, это не сделает его суперполиномиальным.

Тут кто-то приходит и говорит мне, насколько я неправ.

Ой, я только что посмотрел на ваши фактические цифры времени. Они не хуже экспоненциальных. Игнорируя первое и последнее значения (поскольку маленькие времена не поддаются точному измерению, и вы, соответственно, не закончили), умножение n на 10 умножает время не более чем на 30. 30 составляет около 10 ^ 1,5, что примерно подходит для грубой силы, как описано выше.

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

Я переработал программу с некоторыми примечаниями, сделанными из комментариев здесь. Функция «факторы» теперь немного более эффективна, и мне также пришлось изменить функцию σ_(2)(n), чтобы принять новый вывод.

«факторы» пошли от вывода, например:

$ (factors 10) => (1 2 5 10)

иметь лайк

$ (factors 10) => ((2 5) (1 10))

Пересмотренная функция выглядит так:

(defun o_2 (n)
"sum of squares of divisors"
  (reduce #'+ (mapcar (lambda (x) (* x x)) (reduce #'append (factors n)))))

После скромных перезаписей, которые я сделал, я сэкономил всего около 7 секунд в расчете на 100 000.

Похоже, мне придется поднять свою задницу и написать более прямой подход.