Расчет 1^X + 2^X + ... + N^X мод 1000000007

Question

Есть ли алгоритм для расчета (1^x + 2^x + 3^x + ... + n^x) mod 1000000007?
Замечания: a^b является b-й степенью а.

Ограничения 1 <= n <= 10^16, 1 <= x <= 1000, Таким образом, значение N очень велико.

Я могу решить только за O(m log m) если m = 1000000007, Это очень медленно, потому что ограничение по времени составляет 2 секунды.

Есть ли у вас эффективный алгоритм?

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

Answer 1

Вы можете подвести итог серии

1**X + 2**X + ... + N**X

с помощью формулы Фолхабера, и вы получите многочлен с X + 1 сила для вычисления произвольного N,

Если вы не хотите вычислять числа Бернулли, вы можете найти многочлен, решив X + 2 линейные уравнения (для N = 1, N = 2, N = 3, ..., N = X + 2), который медленнее, но проще в реализации.

Давайте иметь пример для X = 2, В этом случае у нас есть X + 1 = 3 порядок полинома:

    A*N**3 + B*N**2 + C*N + D

Линейные уравнения

      A +    B +   C + D = 1              =  1 
    A*8 +  B*4 + C*2 + D = 1 + 4          =  5
   A*27 +  B*9 + C*3 + D = 1 + 4 + 9      = 14
   A*64 + B*16 + C*4 + D = 1 + 4 + 9 + 16 = 30 

Решив уравнения, получим

  A = 1/3
  B = 1/2
  C = 1/6
  D =   0 

Окончательная формула

  1**2 + 2**2 + ... + N**2 == N**3 / 3 + N**2 / 2 + N / 6

Теперь все, что вам нужно сделать, это поставить произвольно большое N в формулу. Пока что алгоритм имеет O(X**2) сложность (так как это не зависит от N).

Answer 2

Есть несколько способов ускорить возведение в моду. С этого момента я буду использовать ** обозначать "возвести в степень" и % для обозначения "модуль".

Сначала несколько замечаний. Это всегда так (a * b) % m является ((a % m) * (b % m)) % m, Это также всегда так a ** n такой же как (a ** floor(n / 2)) * (a ** (n - floor(n/2)), Это означает, что для показателя степени <= 1000 мы всегда можем завершить возведение в степень не более чем в 20 умножениях (и 21 моде).

Мы также можем пропустить довольно много расчетов, так как (a ** b) % m такой же как ((a % m) ** b) % m и если m значительно меньше, чем n, мы просто имеем несколько повторяющихся сумм с "хвостом" частичного повторения.

Answer 3

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

У меня сегодня нет времени для подробного ответа, но учти это.1^2 + 2^2 + 3^2 + ... + n^2 будет принимать O(N) шагов для непосредственного вычисления. Тем не менее, это эквивалентно (n(n+1)(2n+1)/6), который может быть вычислен за O(1) раз. Аналогичная эквивалентность существует для любой более высокой интегральной степениx.

Там может быть общее решение таких проблем; Я не знаю ни одного, и Вольфрам Альфа, кажется, тоже не знает ни одного. Однако в общем случае эквивалентное выражение имеет степень x+1 и может быть получено путем вычисления некоторых значений выборки и решения ряда линейных уравнений для коэффициентов.

Тем не менее, это будет трудно для большого x, такого как 1000, как в вашей проблеме, и, вероятно, не может быть сделано в течение 2 секунд.

Возможно, кто-то, кто знает больше математики, чем я, может превратить это в работоспособное решение?

Изменить: Ой, я вижу, Фабиан Пийке уже опубликовал полезную ссылку о формуле Фолхабера, прежде чем я написал.

Answer 4

Если вы хотите что-то простое и быстрое в реализации, попробуйте это:

Function Sum(x: Number, n: Integer) -> Number
  P := PolySum(:x, n)
  return P(x)
End

Function PolySum(x: Variable, n: Integer) -> Polynomial
  C := Sum-Coefficients(n)
  P := 0
  For i from 1 to n + 1
    P += C[i] * x^i
  End
  return P
End

Function Sum-Coefficients(n: Integer) -> Vector of Rationals
  A := Create-Matrix(n)
  R := Reduced-Row-Echelon-Form(A)
  return last column of R
End

Function Create-Matrix(n: Integer) -> Matrix of Integers
  A := New (n + 1) x (n + 2) Matrix of Integers
  Fill A with 0s
  Fill first row of A with 1s

  For i from 2 to n + 1
    For j from i to n + 1
      A[i, j] := A[i-1, j] * (j - i + 2)
    End
    A[i, n+2] := A[i, n]
  End
  A[n+1, n+2] := A[n, n+2]
  return A
End

объяснение

Наша цель - получить многочлен Q такой, что Q(x) = sum i^n for i from 1 to x, Знаю это Q(x) = Q(x - 1) + x^n => Q(x) - Q(x - 1) = x^nЗатем мы можем составить систему уравнений следующим образом:

d^0/dx^0( Q(x) - Q(x - 1) ) = d^0/dx^0( x^n )
d^1/dx^1( Q(x) - Q(x - 1) ) = d^1/dx^1( x^n )
d^2/dx^2( Q(x) - Q(x - 1) ) = d^2/dx^2( x^n )
                           ...                            .
d^n/dx^n( Q(x) - Q(x - 1) ) = d^n/dx^n( x^n )

При условии, что Q(x) = a_1*x + a_2*x^2 + ... + a_(n+1)*x^(n+1)тогда у нас будет n+1 линейные уравнения с неизвестными a1, ..., a_(n+1)и получается коэффициент cj умножение неизвестного aj в уравнении i следует шаблону (где (k)_p знак равно (k!)/(k - p)!):

• если j < i, cj знак равно 0
• иначе, cj знак равно (j)_(i - 1)

и независимая ценность iэто уравнение (n)_(i - 1), Объяснение, почему становится немного грязно, но вы можете проверить доказательство здесь.

Вышеуказанный алгоритм эквивалентен решению этой системы линейных уравнений. Множество реализаций и дальнейших объяснений можно найти в https://github.com/fcard/PolySum. Основным недостатком этого алгоритма является то, что он потребляет много памяти, даже моя самая эффективная версия памяти использует почти 1gb за n=3000, Но это быстрее, чем SymPy и Mathematica, поэтому я предполагаю, что все в порядке. Сравните с методом Шульца, который использует альтернативный набор уравнений.

Примеры

Легко применить этот метод вручную для небольших n, Матрица для n=1 является

| (1)_0   (2)_0   (1)_0 |     |   1       1       1   |
|   0     (2)_1   (1)_1 |  =  |   0       2       1   |

Применяя исключение Гаусса-Джордана, мы получаем

|   1       0      1/2  |
|   0       1      1/2  |

Результат = {a1 = 1/2, a2 = 1/2} => Q(x) = x/2 + (x^2)/2

Обратите внимание, что матрица всегда находится в форме ряда строк, нам просто нужно уменьшить ее.

За n=2:

| (1)_0   (2)_0   (3)_0   (2)_0 |     |   1       1       1       1   |
|   0     (2)_1   (3)_1   (2)_1 |  =  |   0       2       3       2   |
|   0       0     (3)_2   (2)_2 |     |   0       0       6       2   |

Применяя исключение Гаусса-Джордана, мы получаем

|   1       1       0       2/3 |      |   1       0       0       1/6 |
|   0       2       0         1 |  =>  |   0       1       0       1/2 |
|   0       0       1       1/3 |      |   0       0       1       1/3 |

Результат = {a1 = 1/6, a2 = 1/2, a3 = 1/3} => Q(x) = x/6 + (x^2)/2 + (x^3)/3

Ключом к скорости алгоритма является то, что он не вычисляет факториал для каждого элемента матрицы, а знает, что (k)_p знак равно (k)_(p-1) * (k - (p - 1)), следовательно A[i,j] знак равно (j)_(i-1) знак равно (j)_(i-2) * (j - (i - 2)) знак равно A[i-1, j] * (j - (i - 2)), поэтому он использует предыдущую строку для расчета текущей.