wataru
@wataru
Разработчик на С++, экс-олимпиадник.

Можно ли быстрее чем за O(N) подсчитать сумму S(N,K,M) = sum i=0..N K*i%M?

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

Мои наблюдения:
Можно найти GCD K и M, потом можно и то и другое поделить и дальше можно допустить, что k и m взаимнопросты.

Дальше, можно считать, что N < M. Потому что полные "круги" из-за взаимнопростоты K и M соберут все остатки по модулю M в перемешку и их можно подсчитать формулой M(M-1)/2. Останется только часть в конце длины N%M.

Можно перейти от модуля к делениям с округлениями вниз, воспользуясь формулой floor(a/b)*b = a - a%b, но как это использовать я не придумал.

Казалось бы, тут что-то вроде решения задачи Иосифа напрашивается - как-то хитро считать арифметическую прогрессию, пока k*i не перевалит за m. Есть сколько-то начальных позиций, куда k*i приземляется, делая полный оборот через m, и вот от этих точек арифметические прогрессии с шагом +k можно сложить за O(1). Проблема только в том, что часть этих прогрессий будет на 1 шаг длинее остальных. Вроде тут напрашивается какая-то рекурсия, но я не смог это формализовать.
  • Вопрос задан
  • 836 просмотров
Решения вопроса 1
wataru
@wataru Автор вопроса, куратор тега Математика
Разработчик на С++, экс-олимпиадник.
Совершенно случайно наткнулся вот на это вот равнество: 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: прилизал код немного.
Ответ написан
Комментировать
Пригласить эксперта
Ответы на вопрос 3
@Sumor
Если я правильно понимаю, то это просто арифметическая прогрессия на кольце вычетов по модулю M.
Так что нужно рассчитать сумму членов арифметической прогрессии по формуле:
((K + K*N)/2 * N) % M
Ответ написан
Alexandroppolus
@Alexandroppolus
кодир
Мысли вслух.

К наблюдениям ещё можно добавить довольно кэпский поинт K < M, думаю понятно почему.

Если К "маленькое", например меньше корня из N, то у нас тут мало прогрессий и они довольно длинные. Границы прогрессий и "точки приземления" можно искать с шагом не менее sqrt(N) - делаем такой шаг, потом в его окрестности за О(1) отыскиваем границу, - итого сложность не более O(sqrt(N)).

Вот что делать для "больших" К, пока не совсем понятно...
Ответ написан
Комментировать
ALLIGATOR
@ALLIGATOR
def calc(n,k,m):
  l = [i*k%m for i in range(m)]
  return sum(l)*((n+1)//m) + sum(l[:(n+1)%m])
  

print(calc(7, 3, 5))


https://www.online-ide.com/5RpHe9U6kS
Ответ написан
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы