---
tags: Template
title: Problem_name
author: Editorial-Slayers Team
license: Public domain
---

<style>
.markdown-body { max-width: 2048px; }
</style>

$\Huge \text{🌱 Code cho cái paper về Diophantine equation}$

-----
###### 🔗 Link: [paper](https://cs.uwaterloo.ca/journals/JIS/VOL23/Binner/binner4.pdf?fbclid=IwAR36ozR04cyipUzTOTvM9OhdzqfZTXBu5jDTj0CRLIx8QOf5n-bJAPWVOSo)

###### 📌 Tags: `Generating Function`, `Number Theory`
###### 👤 Writter: @matchamachiato
###### 👥 Contributor: [@](https://codeforces.com/profile/)
###### 📋 Content:

[TOC]

-----

Bài viết này chủ yếu sẽ nói đến phần cài đặt thuật toán tìm số lượng nghiệm của một phương trình $\text{Diophantine tuyến tính: ax + by + cz + d = n, với a, b, c, n cho trước}$.

Về phần chứng minh và lý thuyết sẽ nằm trong [LINK](https://cs.uwaterloo.ca/journals/JIS/VOL23/Binner/binner4.pdf?fbclid=IwAR36ozR04cyipUzTOTvM9OhdzqfZTXBu5jDTj0CRLIx8QOf5n-bJAPWVOSo).

Cụ thể, phần code sẽ là cách cài đặt phần **2. main theorem** trong paper:

### 2.1 Reduction to an equation with pairwise coprime coefficients

Phần này sẽ đưa phương trình ban đầu về phương trình có những hệ số là số nguyên tố cùng nhau.

Từ phương trình trên là $ax + by + cz = n$ ta sẽ đưa về $Ax + By + Cz = N$ với hệ số $\text{A, B, C, N}$ được tính như sau:

- Đặt $g_1, g_2, g_3$ lần lượt là $gcd(b, c), gcd(c, a), gcd(a, b)$. Chú ý rằng $gcd(g_1, g_2) = gcd(g_2, g_3) = gcd(g_1, g_3) = 1$
- Đặt $a_1, b_2, c_3$ là giá trị nghịch đảo modulo của $a, b, c$ với các modulo $g_1, g_2, g_3$. Hay nói cách khác:
$\begin{cases}
a_1 \equiv a^{-1} \pmod {g_1} \\
b_1 \equiv b^{-1} \pmod {g_2} \\
c_1 \equiv c^{-1} \pmod {g_3} \\
\end{cases}$

> **Pre declare:**
> [color=lightgreen]
:::success
:::spoiler code

    n /= gcd; a /= gcd; b /= gcd; c /= gcd;
	
    long long 
	g1 = __gcd(b, c), 
	g2 = __gcd(c, a), 
	g3 = __gcd(a, b); 

	long long
	a1 = binpow(a, phi(g1) - 1, g1),
	b2 = binpow(b, phi(g2) - 1, g2),
	c3 = binpow(c, phi(g3) - 1, g3);

	long long
	n1 = n*a1%g1,
	n2 = n*b2%g2,
	n3 = n*c3%g3;

	long long
	A = a/(g3*g2),
	B = b/(g3*g1),
	C = c/(g1*g2);

	ll N = (n - a*n1 - b*n2 - c*n3) / (g1*g2*g3);

:::

### 2.2 Statement of theorem and proof

Sau khi đã hạ hệ số thành những số nguyên tố cùng nhau. Ta chỉ cần đếm nghiệm của

Công thức và chứng minh. Nhưng ở đây mình chỉ đưa ra công thức.

Định nghĩa $N(a, b, c, n)$ là số nghiệm cho một phương trình $\text{ax + by + cz = n}$.

Một vài giá trị cần lưu ý $(*)$:
- Ta tìm $b_1'$ sao cho $b_1' \equiv -nb^{-1} \text{ (mod a)}$ với $1 \le b_1' \le a$. Ta lại tìm $c_1'$ sao cho $c_1' \equiv bc^{-1} \text{ (mod a)}$ với $1 \le c_1' \le a$
- Ta tìm $c_2'$ sao cho $c_2' \equiv -nc^{-1} \text{ (mod b)}$ với $1 \le c_2' \le b$. Ta lại tìm $a_2'$ sao cho $a_2' \equiv ca^{-1} \text{ (mod b)}$ với $1 \le a_2' \le b$
- Ta tìm $a_3'$ sao cho $a_3' \equiv -na^{-1} \text{ (mod c)}$ với $1 \le a_3' \le c$. Ta lại tìm $b_3'$ sao cho $b_3' \equiv ab^{-1} \text{ (mod a)}$ với $1 \le b_3' \le a$
- Ta định nghĩa: $N_1 = n(n + a + b + c) + cbb_1'(a + 1 - c_1'(b_1' - 1)) + acc_2'(b + 1 - a_2'(c_2 -1)) + baa_3'(c + 1 - b_3'(a_3'-1))$

Với $\text{a, b, c và n}$ là những số nguyên dương cho trước sao cho $gcd(a, b) = gcd(a, c) = gcd(b, c) = 1$. Cùng với những giá trị ở $(*)$ ta có công thức tính số nghiệm không âm của phương trình trên như sau:

$N(a, b, c, n) = \frac{N_1}{2abc} +
\overset{b_1'-1}{\underset{i = 1}{\large \Sigma}} \left \lfloor \frac{ic_1'}{a} \right \rfloor +
\overset{c_2'-1}{\underset{i = 1}{\Large \Sigma}} \left \lfloor \frac{ia_2'}{b} \right \rfloor +
\overset{a_3'-1}{\underset{i = 1}{\Large \Sigma}} \left \lfloor \frac{ib_3'}{c} \right \rfloor - 2$

# Code
> **Algo:** theo paper
> [color=lightgreen]

:::success
:::spoiler Code

#include <bits/stdc++.h>
#define ll long long
#define ull unsigned long long
#define pii pair<long long, long long>
#define st first
#define nd second
#define file "sweet3"
using namespace std;
const long long oo = 1e18;
const long long N = 2e5 + 5;

ll phi(ll n) {
    ll res = n;
    for (ll i = 2; i * i <= n; ++i) {
        if (n % i == 0) {
            while (n % i == 0) {
                n /= i;
            }
            res -= res / i;
        }
    }
    if (n != 1) {
        res -= res / n;
    }
    return res;
}

ll binpow(ll a, ll b, ll r){
	if (b == 0) return 1;
	ll c = binpow(a, b/2, r);
	c = c*c%r;
	if (b & 1) return c * a % r;
	return c;
}

ll cal_(ll n, ll b, ll a){
	ll ivb = binpow(b, phi(a) - 1, a);

	ll pw = a;

	while (pw > 1 && pw <= - n)
		pw *= a;

	n = (n + pw)%a;

	ll cur = n*ivb%a;

	// cout << pw << ' ';

	// cur = (cur + pw) % a;
	if (cur == 0) cur = a;
	
	// cout << "cal_(" << a << "): " << cur << '\n';
	
	return cur;
}

ll sum(ll u, ll tu, ll mau){
	// cout << u << ' ' << tu << ' ' << mau << '\n';
	// return 1;
	ll ans = 0;
	for (int i = 1; i < u; i ++){
		ans += floor(i*tu/mau);
	}
	return ans;
}

void number_of_solution(ll a, ll b, ll c, ll n){
	
	ll b_1 = cal_(-n, b, a), c_1 = cal_(b, c, a);
	ll c_2 = cal_(-n, c, b), a_2 = cal_(c, a, b);
	ll a_3 = cal_(-n, a, c), b_3 = cal_(a, b, c);

	ll N1 = n*(n + a + b + c) 
			+ c*b*b_1*(a + 1 - c_1*(b_1 - 1))
			+ a*c*c_2*(b + 1 - a_2*(c_2 - 1))
			+ b*a*a_3*(c + 1 - b_3*(a_3 - 1));

	ll ans = 	N1/(2*a*b*c) - 2 
				+ sum(b_1, c_1, a) 
				+ sum(c_2, a_2, b) 
				+ sum(a_3, b_3, c);

	cout << ans << '\n'; // << ' ' << N1/(2*a*b*c) - 2  << '\n';

}

ll n, a, b, c;

void solve(){
	// cout << '-' << '\n';
	cin >> n >> a >> b >> c;
	ll gcd = __gcd(c, __gcd(a, b));
	if (n % gcd) {
		cout << 0 << '\n';
		return;
	}

	n /= gcd; a /= gcd; b /= gcd; c /= gcd;

	ll 
	g1 = __gcd(b, c), 
	g2 = __gcd(c, a), 
	g3 = __gcd(a, b); 

	ll 
	a1 = binpow(a, phi(g1) - 1, g1),
	b2 = binpow(b, phi(g2) - 1, g2),
	c3 = binpow(c, phi(g3) - 1, g3);

	ll 
	n1 = n*a1%g1,
	n2 = n*b2%g2,
	n3 = n*c3%g3;

	ll 
	A = a/(g3*g2),
	B = b/(g3*g1),
	C = c/(g1*g2);

	ll N = (n - a*n1 - b*n2 - c*n3) / (g1*g2*g3);

	number_of_solution(A, B, C, N);

}

int main()
{
	ios_base::sync_with_stdio(0); cin.tie(0); cout.tie(0);
	#ifndef ONLINE_JUDGE
		freopen(file".inp","r",stdin); freopen(file".out","w",stdout);
	#endif
	
	ll t;
	cin >> t;
	while (t--)
		solve();

	return 0;
}	
/**  /\_/\
 *  (= ._.)
 *  / >0  \>1
**/


	// ll b_1 = cal_(-n, b, a), c_1 = cal_(b, c, a);
	// ll c_2 = cal_(-n, c, b), a_2 = cal_(c, a, b);
	// ll a_3 = cal_(-n, a, c), b_3 = cal_(a, b, c);

	// ll N1 = n*(n + a + b + c) 
	// 		+ c*b*b_1*(a + 1 - c_1*(b_1 - 1))
	// 		+ a*c*c_2*(b + 1 - a_2*(c_2 - 1))
	// 		+ b*a*a_3*(c + 1 - b_3*(a_3 - 1));

	// ll ans = 	N1/(2*a*b*c) - 2 
	// 			+ sum(b_1, c_1, a) 
	// 			+ sum(c_2, a_2, b) 
	// 			+ sum(a_3, b_3, c);

	// cout << ans << ' ' << N1/(2*a*b*c) - 2  << '\n';

:::

-----
## Bonus