Виконання обчислень слід здійснювати в розширеному полі:
Правила операцій:
Додавання
Множення
Нехай x = a + b \sqrt{c}. Тоді шума сума є кінцевою геометричною прогресією:
Залишається реалізувати ділення a = x^{n + 1} - 1 на b = x - 1 в полі Z_p (\sqrt{c}):
Приклад
У першому прикладі
Обчислимо цей приклад за формулою:
У другому прикладі
Реалізація алгоритму
Функція modpow обчислює значення a^k~mod~p.
long long modpow(long long a, long long k, long long p)
{
long long res = 1;
a %= p;
while (k)
{
if (k & 1) res = res * a % p;
a = a * a % p;
k >>= 1;
}
return res;
}Оголосимо структуру для зберігання числа виду x + y \sqrt{c}.
struct Num
{
long long x, y; // x + y*sqrt(c)
};Реалізуємо додавання в полі Z_p (\sqrt{c}). Функція add повертає суму чисел та .
Num add(Num& a, Num& b)
{
Num res;
res.x = (a.x + b.x) % p;
res.y = (a.y + b.y) % p;
return res;
}Реалізуємо віднімання в полі Z_p (\sqrt{c}). Функція sub повертає різницю чисел та .
Num sub(Num& a, Num& b)
{
Num res;
res.x = (a.x - b.x + p) % p;
res.y = (a.y - b.y + p) % p;
return res;
}Реалізуємо множення в полі Z_p (\sqrt{c}). Функція mult повертає добуток чисел та .
Num mult(Num& a, Num& b)
{
Num res;
res.x = (a.x * b.x + c * a.y % p * b.y) % p;
res.y = (a.x * b.y + a.y * b.x) % p;
return res;
}Функція pow реалізує піднесення до степеня base^n, де — число в полі Z_p (\sqrt{c}).
Num pow(Num base, long long n)
{
Num res = { 1, 0 }; // одиниця поля
while (n > 0)
{
if (n & 1)
res = mult(res, base);
base = mult(base, base);
n >>= 1;
}
return res;
}Функція divide реалізує ділення в полі Z_p (\sqrt{c}). Вона повертає частку чисел та . Ділення реалізуємо через множення на спряжене:
Num divide(Num& a, Num& b)
{Нехай conj = b_x - b_y \sqrt{c}.
Num conj = { b.x, (p - b.y) % p };Обчислюємо n = (a_x + a_y \sqrt{c})(b_x - b_y \sqrt{c}).
Num n = mult(a, conj);Обчислюємо denom = b_x^2 - b_y^2 c.
long long denom =
((b.x * b.x) % p - ((b.y * b.y) % p * c) % p + p) % p;Обчислюємо обернений елемент і повертаємо n / denom = n \cdot denom^{-1}.
long long inv_denom = modpow(denom, p - 2, p); // обернений елемент
return Num((n.x * inv_denom) % p, (n.y * inv_denom) % p);
}Основна частина програми. Читаємо вхідні дані.
scanf("%lld %lld %lld %lld %lld", &a, &b, &c, &n, &p);Ініціалізуємо число x = a + b \sqrt{c} та одиницю поля one = 1 + 0 \sqrt{c}.
Num x = { a % p, b % p };
Num one = { 1, 0 };Обчислюємо num = x^{n+1} - 1, denom = x - 1.
Num num = pow(x, n + 1);
num = sub(num, one);
Num denom = sub(x, one);Обчислюємо і виводимо відповідь res = num / denom = \frac{x^{n+1} - 1}{x - 1}.
Num res = divide(num, denom);
printf("%lld %lld\n", res.x, res.y);