Совершенно случайно наткнулся вот на это вот равнество:
hermite's identity (в твитте с прикольным доказательством.
Зацените).
C помощью него наконец-то вывел способ подсчитать эти суммы за O(log n)! Я знал, что это можно сделать!
Итак, во-первых, можно переключатся между суммами floor(i*k/m) и (i*k)%m через вот это равенство:
floor(i*k/m) = (i*k - (i*k)%m) / m
Это нам позже поможет. В сумме через floor можно сократить GCD(m, k) в них. В сумме через модуль можно сделать k %= m, если оно больше m, да и сократить полные "обороты" в n.
Итак, можно допустить, что m > k, m >n, GCD(m, k) = 1, иначе преобразуем сумму к нужной форме и все упростим.
Дальше, применим равенство hermite's, взяв x = i/m, n = k.
Потом поменяем местами суммы. Под знаками суммы будет floor(i/m + t/k) (где t - новая переменная суммирования от 0 до k-1). Присмотритесь к этому выражению - там 2 числа меньших 1. Т.е. весь этот floor даст 1, только если t/k >= 1-i/m. Отсюда можно решить это неравнество, сдвинуть нижнюю границу суммирования и получить сумму единиц во внутренней сумме. Заменив ее всю на количество слагаемых там вылезает floor (t*m/k). Т.е. мы выразили сумму i*k/m через сумму t*m/k. Но мы же помним, что можно перейти к сумме модулей и там сократить множитель k в числителе. Таким образом мы вычисляем сумму для (m, k) через сумму для (k, m%k). Это точно такой же переход, как и в GCD, поэтому суммарно будет всего O(log n) итераций. Вообще, выкладки довольно нудные, ибо сумма для t*m/k будет не от 0 до какого-то n', а от n' до k-1. Но можно взять известное значение для суммы полного оборота (от 0 до k-1) и из нее вычесть сумму от 0 до n'-1. Эта сумма известна, потому что при GCD=1, она пройдется по всем остаткам в каком-то порядке.
Формула примерно такая получается:
FloorSum(n, k, m) = (m-1)*(k-1)/2 - (n1+1)*n1*(m/k)/2 + (n - m + 1)*(k-n1-1) - FloorSum(n1, m%k, k)
n' = floor(((m-n)*k-1)/m)
Вот код. Значения совпадают с тупым решением для всех чисел до 1000 и для миллиардов случайных чисел до миллиона.
// sum i = 0.. n floor(i * k / m)
// GCD(k, m) must be 1.
// n < m
// k < m
long long FloorSumInternal(long long n, long long k, long long m) {
if (k == 0 || n <= 0) return 0;
if (m == 1) return k*n*(n+1)/2;
const long long n1 = ((m-n)*k - 1)/m;
long long ans = (m-1)*(k-1)/2 - (n1+1)*n1*(m/k)/2 + (n - m + 1)*(k-n1-1);
ans -= FloorSumInternal(n1, m%k, k);
return ans;
}
// sum i = 0.. n floor(i * k / m)
long long FloorSum(long long n, long long k, long long m) {
if (k == 0 || n <= 0) return 0;
if (m == 1) return k*n*(n+1)/2;
const long long d = Gcd(m, k);
m /= d;
k /= d;
if (k >= m || n >= m) {
// (n*k*(n+1)/2 - ModSum(n, k, m, d))/m;
const long long nn = (n+1)%m-1;
const long long num_full = (n+1) / m;
const long long kk = k % m;
long long ans = 0;
ans = (k*n*(n+1) - kk*(nn+1)*nn)/m - num_full * (m-1);
ans /= 2;
ans += FloorSumInternal(nn, kk, m);
return ans;
}
return FloorSumInternal(n, k, m);
}
// sum i = 0.. n (i * k) % m;
long long ModSum(long long n, long long k, long long m)
{
long long d = Gcd(k, m);
if (k == 0 || m == 1) return 0;
// (i*k) % m = i*k-floor(i*k/m)*m
k %= m;
long long num_full = (n+1) / m;
long long ans = num_full * (m-d) * m/2;
n = (n+1)%m-1;
if (n > 0) {
ans += ((long long) k)*(n+1)*n/2 - m*FloorSum(n, k/d, m/d);
}
return ans;
}
Правда тут есть проблема, в процессе вычисления получаются весьма большие числа. Поэтому даже если финальный ответ в long long влезает, ответ может переполнится. Но все-равно, там числа больше куба от входных значений не используются, поэтому можно в оценке сложности это не учитывать.
Edit: прилизал код немного.
