Skip to content

Commit fcab3b9

Browse files
committed
Add submission for SANDWICH
- Algorithms used - Extended Euclid's Algorithm - Multiplicative modular inverse - Generalized Lucas' Theorem - Chinese Remainder Theorem
1 parent 0a62f6e commit fcab3b9

File tree

1 file changed

+228
-0
lines changed

1 file changed

+228
-0
lines changed

LongChallenge/May16/may8.cpp

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
#include <bits/stdc++.h>
2+
#define ll long long
3+
#define ull unsigned long long
4+
#define vi vector<ll>
5+
#define pp pair<ll,ll>
6+
#define mp make_pair
7+
#define PI acos(-1.0)
8+
#define all(v) v.begin(),v.end()
9+
#define pb push_back
10+
#define FOR(i,a,b) for(i=a;i<b;i++)
11+
#define FREV(i,a,b) for(i=a;i>=b;i--)
12+
#define SULL(n) scanf("%llu", &n)
13+
#define INF 1e18
14+
15+
#ifndef ONLINE_JUDGE
16+
#define gc getchar
17+
#define pc putchar
18+
#else
19+
#define gc getchar_unlocked
20+
#define pc putchar_unlocked
21+
#endif
22+
23+
using namespace std;
24+
25+
int read_int() {
26+
char c = gc();
27+
while((c < '0' || c > '9') && c != '-') c = gc();
28+
int ret = 0, neg = 0;
29+
if (c == '-') neg = 1, c = gc();
30+
while(c >= '0' && c <= '9') {
31+
ret = 10 * ret + c - 48;
32+
c = gc();
33+
}
34+
return neg ? -ret : ret;
35+
}
36+
37+
ll read_ll() {
38+
char c = gc();
39+
while((c < '0' || c > '9') && c != '-') c = gc();
40+
ll ret = 0;
41+
int neg = 0;
42+
if (c == '-') neg = 1, c = gc();
43+
while(c >= '0' && c <= '9') {
44+
ret = 10 * ret + c - 48;
45+
c = gc();
46+
}
47+
return neg ? -ret : ret;
48+
}
49+
50+
vector<pp > prime_factors;
51+
vi prime_powers;
52+
53+
// Legendre's theorem to return highest power of a prime p in n!
54+
// http://www.cut-the-knot.org/blue/LegendresTheorem.shtml
55+
ll legendre(ll n, ll p) {
56+
ll ans = 0;
57+
while (n > 0) {
58+
n = n / p;
59+
ans = ans + n;
60+
}
61+
return ans;
62+
}
63+
64+
// Extended Euclid's Algorithm
65+
ll xGCD(ll a, ll b, ll &x, ll &y) {
66+
if(b == 0) {
67+
x = 1;
68+
y = 0;
69+
return a;
70+
}
71+
72+
ll x1, y1, gcd = xGCD(b, a % b, x1, y1);
73+
x = y1;
74+
y = x1 - (a / b) * y1;
75+
return gcd;
76+
}
77+
78+
// Multiplicative modular inverse using Extended Euclid's Algorithm
79+
ll find_inverse(ll n, ll mod) {
80+
ll x,y;
81+
xGCD(n,mod,x,y);
82+
83+
if (x < 0) {
84+
x = x + mod;
85+
}
86+
87+
return x;
88+
}
89+
90+
// factorize into prime factors
91+
void factorize(ll n) {
92+
ll power,i,val;
93+
prime_factors.clear();
94+
prime_powers.clear();
95+
for(i=2; i*i <= n; i++) {
96+
power = 0, val = 1;
97+
while (n%i == 0) {
98+
n = n / i;
99+
power++;
100+
val = val * i;
101+
}
102+
if(power != 0) {
103+
prime_factors.pb(mp(i,power));
104+
prime_powers.pb(val);
105+
}
106+
}
107+
if (n != 1) {
108+
prime_factors.pb(mp(n,1));
109+
prime_powers.pb(n);
110+
}
111+
}
112+
113+
// Generalized form of Lucas' Theorem
114+
// https://math.stackexchange.com/questions/60206/lucas-theorem-but-without-prime-numbers
115+
ll nCr_mod_prime_power(ll n, ll m, ll prime, ll exponent, ll prime_power) {
116+
ll highest_divisible_power = legendre(n,prime) - legendre(m,prime) - legendre(n-m,prime);
117+
118+
if (highest_divisible_power >= exponent) {
119+
return 0;
120+
}
121+
122+
exponent = exponent - highest_divisible_power;
123+
ll i, new_prime_power = 1;
124+
FOR(i,0,exponent) {
125+
new_prime_power = new_prime_power * prime;
126+
}
127+
128+
vi factorials(new_prime_power);
129+
factorials[0] = 1;
130+
FOR(i,1,new_prime_power) {
131+
if (i % prime == 0) {
132+
factorials[i] = factorials[i-1];
133+
}
134+
else {
135+
factorials[i] = (factorials[i-1] * i) % new_prime_power;
136+
}
137+
}
138+
139+
ll numerator = 1, denominator = 1, r = n-m;
140+
bool is_negative = false;
141+
int digits = 0;
142+
143+
while (n > 0) {
144+
if (digits >= exponent && factorials[new_prime_power - 1] != 1) {
145+
is_negative = is_negative ^ ((n&1) ^ (m&1) ^ (r&1));
146+
}
147+
148+
numerator = (numerator * factorials[n % new_prime_power]) % new_prime_power;
149+
denominator = (denominator * factorials[m % new_prime_power]) % new_prime_power;
150+
denominator = (denominator * factorials[r % new_prime_power]) % new_prime_power;
151+
152+
n = n / prime;
153+
m = m / prime;
154+
r = r / prime;
155+
156+
digits++;
157+
}
158+
159+
ll answer = (numerator * find_inverse(denominator, new_prime_power)) % new_prime_power;
160+
161+
if (is_negative && (prime != 2 || exponent < 3)) {
162+
answer = new_prime_power - answer;
163+
}
164+
165+
FOR(i,0,highest_divisible_power) {
166+
answer = (answer * prime) % prime_power;
167+
}
168+
169+
return answer;
170+
}
171+
172+
// Chinese Remainder Theorem
173+
// https://brilliant.org/wiki/chinese-remainder-theorem/
174+
ll chinese_remainder_theorem(vi &remainders, ll mod) {
175+
if (remainders.size() == 1) {
176+
return remainders[0];
177+
}
178+
179+
ll i, len = remainders.size();
180+
ll answer = 0;
181+
FOR(i,0,len) {
182+
answer = (answer + ((remainders[i] * (mod / prime_powers[i]))%mod * find_inverse(mod / prime_powers[i], prime_powers[i]))%mod)%mod;
183+
}
184+
return answer;
185+
}
186+
187+
188+
// nCr modulo m (where m may be prime or composite)
189+
ll nCr_modulo(ll n, ll m, ll mod) {
190+
if (mod == 1) {
191+
return 0;
192+
}
193+
ll i;
194+
factorize(mod);
195+
ll len = prime_powers.size();
196+
vi remainders;
197+
198+
FOR(i,0,len) {
199+
remainders.pb(nCr_mod_prime_power(n,m,prime_factors[i].first, prime_factors[i].second, prime_powers[i]));
200+
}
201+
return chinese_remainder_theorem(remainders,mod);
202+
}
203+
204+
int main() {
205+
ll i,j,t,n,k,m,r,mod;
206+
t = read_ll();
207+
while(t--) {
208+
n = read_ll();
209+
k = read_ll();
210+
mod = read_ll();
211+
212+
m = n / k;
213+
if (n % k != 0) {
214+
m++;
215+
}
216+
217+
if (n%k == 0) {
218+
printf("%lld 1\n", m);
219+
continue;
220+
}
221+
222+
n = m - n - 1 + m*k;
223+
r = min(m-1,m*k - n);
224+
225+
printf("%lld %lld\n", m, nCr_modulo(n,r,mod));
226+
}
227+
return 0;
228+
}

0 commit comments

Comments
 (0)