• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

_no0B Apr 25th, 2019 117 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <ext/rope>
2. #include <bits/stdtr1c++.h>
3. #include <ext/pb_ds/assoc_container.hpp>
4. #include <ext/pb_ds/tree_policy.hpp>
5. using namespace __gnu_cxx;
6. using namespace __gnu_pbds;
7. using namespace namespace;
8. #define clr(ar) memset(ar, 0, sizeof(ar))
9. #define read() freopen("lol.txt", "r", stdin)
10. #define ran(a, b) ((((rand() << 15) ^ rand()) % ((b) - (a) + 1)) + (a))
11. #define dbg(x) cout << #x << " = " << x << endl
12. /// matrix
13. #define MOD 1000000007
14. struct Matrix{
15.     int row, col;
16.     int ar[101][101]; /// Change matrix size here, also change to long long for safety
17.     Matrix(){} ///Beware if matrix can contain negative numbers in matrix exponentiation problems
18.     Matrix(int n, int m, int diagonal = 0){
19.         clr(ar);
20.         row = n, col = m;
21.         for (int i = min(n, m) - 1; i >= 0; i--) ar[i][i] = diagonal;
22.     }
23.     /// To multiply two matrices A and B, the number of columns in A must equal the number of rows in B
24.     Matrix operator* (const Matrix& other) const{ /// hash = 709758
25.         int i, j, k;
26.         Matrix res(row, other.col);
27.         long long x, y = (long long)MOD * MOD;
28.         for(i = 0; i < row; i++){
29.             for(j = 0; j < other.col; j++){
30.                 for(k = 0, x = 0; k < col; k++){
31.                     x += ((long long)ar[i][k] * other.ar[k][j]); /// replace matrix other with its transpose matrix to reduce cache miss
32.                     if (x >= y) x -= y;
33.                 }
34.                 res.ar[i][j] = x % MOD;
35.             }
36.         }
37.         return res;
38.     }
39.     Matrix operator^ (long long n) const{
40.         Matrix x = *this, res = Matrix(row, col, 1);
41.         while (n){
42.             if (n & 1) res = res * x;
43.             n = n >> 1, x = x * x;
44.         }
45.         return res;
46.     }
47.     /// Transpose matrix, T[i][j] = ar[j][i]
48.     Matrix transpose(){
49.         Matrix res = Matrix(col, row);
50.         for (int i = 0; i < row; i++){
51.             for (int j = 0; j < col; j++){
52.                 res.ar[j][i] = ar[i][j];
53.             }
54.         }
55.         return res;
56.     }
57.     /// rotates the matrix 90 degrees clockwise
58.     Matrix rotate(){
59.         Matrix res = this->transpose();
60.         for (int i = 0; i < res.row; i++) reverse(res.ar[i], res.ar[i] + res.col);
61.         return res;
62.     }
63.     inline void print(){
64.         for (int i = 0; i < row; i++){
65.             for (int j = 0; j < col; j++){
66.                 printf("%d%c", ar[i][j], ((j + 1) == col) ? 10 : 32);
67.             }
68.         }
69.     }
70. };
71. int main(){
72.     Matrix a = Matrix(4, 5, 1);
73.     int k = 0;
74.     for (int i = 0; i < a.row; i++){
75.         for (int j = 0; j < a.col; j++){
76.             a.ar[i][j] = ++k;
77.         }
78.     }
79.     a.print();
80.     puts("");
81.     Matrix b = a.rotate();
82.     b.print();
83.     return 0;
84.     Matrix x = Matrix(5, 5, 5);
85.     Matrix y = x ^ 5;
86.     x.print();
87.     y.print();
88. }
89.
90.
91. /***
92.
93. Sometimes we may need to maintain more than one recurrence, where they are interrelated.
94. For example, let a recurrence relation be:
95. g(n) = 2g(n-1) + 2g(n-2) + f(n), where, f(n) = 2f(n-1) + 2f(n-2).
96. Here, recurrence g(n) is dependent upon f(n) and the can be calculated in the same matrix
97. but of increased dimensions. Lets design the matrices A, B then we'll try to find matrix M.
98.
99.
100.                         Matrix A                Matrix B
101.
102.                         |  g(n)  |              | g(n+1) |
103.                         | g(n-1) |              |  g(n)  |
104.                         | f(n+1) |              | f(n+2) |
105.                         |  f(n)  |              | f(n+1) |
106.
107.
108. Here, g(n+1) = 2g(n) + 2g(n-1) + f(n+1) and f(n+2) = 2f(n+1) + 2f(n).
109. Now, using the above process, we can generate the objective matrix M as follows:
110.
111. | 2  2  1  0 |
112. | 1  0  0  0 |
113. | 0  0  2  2 |
114. | 0  0  1  0 |
115.
116.
117. /// Matrix Rotations:
118.
119. Rotate by +90:
120. Transpose
121. Reverse each row
122.
123. Rotate by -90:
124. Transpose
125. Reverse each column
126.
127.
128. Rotate by +180:
129. Method 1: Rotate by +90 twice
130. Method 2: Reverse each row and then reverse each column
131.
132. Rotate by -180:
133. Method 1: Rotate by -90 twice
134. Method 2: Reverse each column and then reverse each row
135. Method 3: Reverse by +180 as they are same
136.
137. ***/
138.
139.
140.
141. /// Lucas theorem to calculate binomial co-efficients modulo a prime
142. #define MAXP 100010
143. namespace lc{
144.     int MOD = 1000000007;
145.     int fact[MAXP], inv[MAXP];
146.     /// Call once with the modulo prime
147.     void init(int prime){
148.         MOD = prime;
149.         fact[0] = 1, inv[MOD - 1] = MOD - 1;
150.         for (int i = 1; i < MOD; i++) fact[i] = ((long long)fact[i - 1] * i) % MOD;
151.         for (int i = MOD - 2; i >= 0; i--) inv[i] = ((long long)inv[i + 1] * (i + 1)) % MOD;
152.     }
153.     inline int count(int n, int k){
154.         if (k > n) return 0;
155.         int x = ((long long)inv[n - k] * inv[k]) % MOD;
156.         return ((long long)x * fact[n]) % MOD;
157.     }
158.     /// Lucas theorem, calculates binomial(n, k) modulo MOD, MOD must be a prime
159.     inline int binomial(long long n, long long k){
160.         if (k > n) return 0;
161.         int res = 1;
162.         k = min(k, n - k);
163.         while (k && res){
164.             res = ((long long)res * count(n % MOD, k % MOD)) % MOD;
165.             n /= MOD, k /= MOD;
166.         }
167.         return res;
168.     }
169.     /*** Alternate and extended functionalities ***/
170.     /// Must call init with prime before (Or set lc::MOD = prime)
171.     /// Computes (n! / (p ^ (n / p))) % p in O(p log(n)) time, p MUST be a prime
172.     /// That is, calculating n! without p's powers
173.     /// For instance factmod(9, 3) = (1 * 2 * 4 * 5 * 2 * 7 * 8 * 1) % 3 = 1
174.     inline int factmod(long long n, int p){
175.         int i, res = 1;
176.         while (n > 1) {
177.             if ((n / p) & 1) res = ((long long)res * (p - 1)) % p;
178.             for (i = n % p; i > 1; i--) res = ((long long)res * i) % p;
179.             n /= p;
180.         }
181.         return (res % p);
182.     }
183.     inline int expo(int a, int b){
184.         int res = 1;
185.         while (b){
186.             if (b & 1) res = (long long)res * a % MOD;
187.             a = (long long)a * a % MOD;
188.             b >>= 1;
189.         }
190.         return res;
191.     }
192.     /// Trailing zeros of n! in base p, p is a prime
193.     inline long long fact_ctz(long long n, long long p){
194.         long long x = p, res = 0;
195.         while (n >= x){
196.             res += (n / x);
197.             x *= p;
198.         }
199.         return res;
200.     }
201.     /// Calculates binomial(n, k) modulo MOD, MOD must be a prime
202.     inline int binomial2(long long n, long long k){
203.         if (k > n) return 0;
204.         if (fact_ctz(n, MOD) != (fact_ctz(k, MOD) + fact_ctz(n - k, MOD))) return 0;
205.         int a = factmod(n - k, MOD), b = factmod(k, MOD), c = factmod(n, MOD);
206.         int x = ((long long)expo(a, MOD - 2) * expo(b, MOD - 2)) % MOD;
207.         return ((long long)x * c) % MOD;
208.     }
209. }
210. int main(){
211.     lc::init(997);
212.     printf("%d\n", lc::binomial(10, 5));
213.     printf("%d\n", lc::binomial(1996, 998));
214.     lc::MOD = 10007;
215.     printf("%d\n", lc::binomial2(10, 5));
216.     printf("%d\n", lc::binomial2(1996, 998));
217.     return 0;
218. }
219.
220.
221. /// maximum divisors
222. unsigned long long n, res, idx;
223. int p, primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71};
224. unsigned long long mul(unsigned long long a, unsigned long long b){
225.     unsigned long long res = 0;
226.     while (b){
227.         if (b & 1LL) res = (res + a);
228.         if (res > n) return 0;
229.         a = (a << 1LL);
230.         b >>= 1LL;
231.     }
232.     return res;
233. }
234. void backtrack(int i, int lim, unsigned long long val, unsigned long long r){
235.     if ((r > res) || (r == res && val < idx)) res = r, idx = val;
236.     if (i == p) return;
237.     int d;
238.     unsigned long long x = val;
239.     for (d = 1; d <= lim; d++){
240.         x = mul(x, primes[i]);
241.         if (x == 0) return;
242.         backtrack(i + 1, d, x, r * (d + 1));
243.     }
244. }
245. int main(){
246.     /* Tested for n <= 10^18 */
247.     p = sizeof(primes) / sizeof(int);
248.     while (scanf("%llu", &n) != EOF){
249.         res = 0;
250.         backtrack(0, 100, 1, 1);
251.         printf("%llu = %llu\n", idx, res);
252.     }
253.     return 0;
254. }
255.
256.
257.
258. /// Miller Rabin
259. #define MAX 1000010
260. namespace prm{ /// hash = 130793
261.     bitset <MAX> flag;
262.     long double op = 0.0;
263.     int p = 0, prime[78777];
264.     const int base[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
265.     void Sieve(){
266.         int i, j, x;
267.         for (i = 3; i < MAX; i += 2) flag[i] = true;
268.         for (i = 3, flag[2] = true; (i * i) < MAX; i += 2){
269.             if (flag[i]){
270.                 for (j = (i * i), x = i << 1; j < MAX; j += x){
271.                     flag[j] = false;
272.                 }
273.             }
274.         }
275.         for (i = 2; i < MAX; i++){
276.             if (flag[i]) prime[p++] = i;
277.         }
278.     }
279.     void init(){
280.         if (!flag[2]) Sieve();
281.     }
282.     inline long long mul(long long a, long long b, long long MOD){
283.        if ((MOD < 3037000500LL)) return ((a * b) % MOD);
284.        long double res = a;
285.        res *= b;
286.        long long c = (long long)(res * op);
287.        a *= b;
288.        a -= c * MOD;
289.        if (a >= MOD) a -= MOD;
290.        if (a < 0) a += MOD;
291.        return a;
292.     }
293.     inline long long expo(long long x, long long n, long long m){
294.         long long res = 1;
295.         while (n){
296.             if (n & 1) res = mul(res, x, m);
297.             x = mul(x, x, m);
298.             n >>= 1;
299.         }
300.         return (res % m);
301.     }
302.     inline bool miller_rabin(long long p){
303.         if (p < MAX) return flag[p];
304.         if ((p + 1) & 1) return false;
305.         for (int i = 1; i < 10; i++){ /// basic iterations
306.             if (!(p % prime[i])) return false;
307.         }
308.         long long a, m, x, s = p - 1, y = p - 1;
309.         op = (long double)1 / p, s = s >> __builtin_ctzll(s);
310.         for (int i = 0; i < 7; i++) {
311.             x = s, a = (base[i] % y) + 1;
312.             m = expo(a, x, p);
313.             while ((x != y) && (m != 1) && (m != y)) m = mul(m, m, p), x <<= 1;
314.             if ((m != y) && !(x & 1)) return false;
315.         }
316.         return true;
317.     }
318.     inline long long countdiv(long long n){
319.         int i, j, c;
320.         long long x, res = 1;
321.         for (i = 0; i < p; i++){
322.             x = prime[i];
323.             if ((x * x * x) > n) break;
324.             c = 1;
325.             while (!(n % x)) c++, n /= x;
326.             res *= c;
327.         }
328.         if (miller_rabin(n)) res <<= 1;
329.         else if (n > 1) {
330.             x = sqrt((long double)0.95 + n); /// may be change to sqrtl() ?
331.             if ((x * x) == n && miller_rabin(x)) res *= 3;
332.             else res <<= 2;
333.         }
334.         return res;
335.     }
336. }
337. int main(){
338.     prm::init();
339. }
340.
341.
342.
343. /// mobius function
344. #define LEN 666666
345. #define MAX 10000010
346. int len = 0, prime[LEN];
347. char mu[MAX] = {0}, flag[MAX] = {0};
348. /// mu[1] = 1, mu[n] = 0 if n has a squared prime factor,
349. /// mu[n] = 1 if n is square-free with even number of prime factors
350. /// mu[n] = -1 if n is square-free with odd number of prime factors
351. void Mobius(){
352.     mu[1] = 1;
353.     int i, j, k;
354.     for (i = 2; i < MAX; i++){
355.         if (!flag[i]) mu[i] = -1, prime[len++] = i;
356.         for (j = 0; j < len && (k = i * prime[j]) < MAX; j++){
357.             flag[k] = true;
358.             if (!(i % prime[j])){
359.                 mu[k] = 0;
360.                 break;
361.             }
362.             else mu[k] -= mu[i];
363.         }
364.     }
365. }
366. void Mobius(){ /// Simplified NlogN version
367.     int i, j;
368.     mu[1] = 1;
369.     for (i = 1; i < MAX; i++){
370.         for (j = i + i; j < MAX; j += i){
371.             mu[j] -= mu[i];
372.         }
373.     }
374. }
375. void Mobius(){ /// Simplified version optimized
376.     int i, j;
377.     mu[1] = 1;
378.     for (i = 1; i < MAX; i++){
379.         if (mu[i]){
380.             for (j = i + i; j < MAX; j += i){
381.                 mu[j] -= mu[i];
382.             }
383.         }
384.     }
385. }
386. int main(){
387.     Mobius();
388. }
389.
390.
391.
392. /// phi
393. #define MAX 4000010
394. int phi[MAX];
395. void Totient(){
396.     clr(phi);
397.     int i, j, x;
398.     for (i = 1; i < MAX; i++){
399.         phi[i] += i;
400.         for (j = (i << 1); j < MAX; j += i){
401.             phi[j] -= phi[i];
402.         }
403.     }
404. }
405. int main(){
406.     Totient();
407. }
408.
409.
410.
411. /// Pick's Theorem
412. struct Point{
413.     long long x, y;
414.     Point(){}
415.     Point(long long x, long long y) : x(x), y(y){}
416. };
417. /// twice the area of polygon
418. long long double_area(Point poly[], int n){
419.     long long res = 0;
420.     for (int i = 0, j = n - 1; i < n; j = i++){
421.         res += ((poly[j].x + poly[i].x) * (poly[j].y - poly[i].y));
422.     }
423.     return abs(res);
424. }
425. /// number of lattice points strictly on polygon border
426. long long on_border(Point poly[], int n){
427.     long long res = 0;
428.     for (int i = 0, j = n - 1; i < n; j = i++){
429.         res += __gcd(abs(poly[i].x - poly[j].x), abs(poly[i].y - poly[j].y));
430.     }
431.     return res;
432. }
433. /// number of lattice points strictly inside polygon
434. long long interior(Point poly[], int n){
435.     long long res = 2 + double_area(poly, n) - on_border(poly, n);
436.     return res / 2;
437. }
438. int main(){
439.     Point ar[10];
440.     ar[0] = Point(0, 0);
441.     ar[1] = Point(3, 0);
442.     ar[2] = Point(3, 3);
443.     ar[3] = Point(0, 3);
444.     dbg(on_border(ar, 4)); /// 12
445.     dbg(interior(ar, 4)); /// 4
446.     return 0;
447. }
448.
449.
450.
451. /// Pollard Rho
452. const long long LIM = LLONG_MAX;
453. long long mul(long long a, long long b, long long m){
454.     long long x, res;
455.     if (a < b) swap(a, b);
456.     if (!b) return 0;
457.     if (a < (LIM / b)) return ((a * b) % m);
458.     res = 0, x = (a % m);
459.     while (b){
460.         if (b & 1){
461.             res = res + x;
462.             if (res >= m) res -= m;
463.         }
464.         b >>= 1;
465.         x <<= 1;
466.         if (x >= m) x -= m;
467.     }
468.     return res;
469. }
470. long long expo(long long x, long long n, long long m){
471.     long long res = 1;
472.
473.     while (n){
474.         if (n & 1) res = mul(res, x, m);
475.         x = mul(x, x, m);
476.         n >>= 1;
477.     }
478.     return (res % m);
479. }
480. const int small_primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 51, 53, 59, 61, 67, 71};
481. bool miller_rabin(long long p, int lim){
482.     long long a, s, m, x, y;
483.     s = p - 1, y = p - 1;
484.     while (!(s & 1)) s >>= 1;
485.     while (lim--){
486.         x = s;
487.         a = (rand() % y) + 1;
488.         m = expo(a, x, p);
489.         while ((x != y) && (m != 1) && (m != y)){
490.             m = mul(m, m, p);
491.             x <<= 1;
492.         }
493.         if ((m != y) && !(x & 1)) return false;
494.     }
495.     return true;
496. }
497. void brent_pollard_rho(uint64_t n, vector <uint64_t> &v){
498.     if (miller_rabin(n, 10)){
499.         v.push_back(n);
500.         return;
501.     }
502.     uint64_t a, g, x, y;
503.     y = 1;
504.     a = rand() % n;
505.     x = rand() % n;
506.     for (int i = 0; ((i * i) >> 1) < n; i++){
507.         x = mul(x, x, n);
508.         x += a;
509.         if (x < a) x += (ULLONG_MAX - n) + 1;
510.         x %= n;
511.         g = __gcd(n, y - x);
512.         if ((g != 1) && (g != n)){
513.             n /= g;
514.             brent_pollard_rho(g, v);
515.             if (n != g) brent_pollard_rho(n, v);
516.             else if (miller_rabin(n, 10)) v.push_back(n);
517.             return;
518.         }
519.         if (!(i & (i - 1))) y = x;
520.     }
521.     brent_pollard_rho(n, v);
522. }
523. void factorize(uint64_t n, vector <uint64_t> &v){
524.     srand(time(0));
525.     int i, j, x;
526.     for (i = 0; i < 21; i++){
527.         x = small_primes[i];
528.         while ((n % x) == 0){
529.             n /= x;
530.             v.push_back(x);
531.         }
532.     }
533.     if (n > 1) brent_pollard_rho(n, v);
534.     sort(v.begin(), v.end());
535. }
536. int main(){
537.     int i, t;
538.     uint64_t n, x;
539.     cin >> t;
540.     while (t--){
541.         cin >> n;
542.         vector <uint64_t> v;
543.         factorize(n, v);
544.         sort(v.begin(), v.end());
545.         v.push_back(-1);
546.         int len = v.size();
547.         int c = 0, counter = 0;
548.         printf("%llu = ", n);
549.         for (i = 0; (i + 1) < len; i++){
550.             if (v[i] == v[i + 1]) counter++;
551.             else{
552.                 if (c) printf(" * ");
553.                 if (counter) printf("%llu^%d", v[i], ++counter);
554.                 else printf("%llu", v[i]);
555.                 c++, counter = 0;
556.             }
557.         }
558.         puts("");
559.     }
560.     return 0;
561. }
562.
563.
564. /// Prime Counting Functions
565. #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))
566. #define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))
567. #define isprime(x) (( (x) && ((x)&1) && (!chkbit(ar, (x)))) || ((x) == 2))
568. namespace pcf{
569.     ///   Prime-Counting Function
570.     ///   initialize once by calling init()
571.     ///   Legendre(m) and Lehmer(m) returns the number of primes less than or equal to m
572.     #define MAXN 1000010
573.     #define MAX_PRIMES 1000010
574.     #define PHI_M 10010
575.     #define PHI_N 48
576.     unsigned int ar[(MAXN >> 6) + 5] = {0};
577.     int len = 0, primes[MAX_PRIMES], counter[MAXN], phi_dp[PHI_N][PHI_M];
578.     void Sieve(int n, unsigned int* ar, int* primes, int& len){
579.         n++;
580.         setbit(ar, 0), setbit(ar, 1);
581.         int i, j, k, lim = sqrt(n) + 1;
582.         for (i = 3; i < lim; i++, i++){
583.             if (!chkbit(ar, i)){
584.                 k = i << 1;
585.                 for (j = (i * i); j < n; j += k) setbit(ar, j);
586.             }
587.         }
588.         for (i = 1; i < n; i++){
589.             if (isprime(i)) primes[len++] = i;
590.             counter[i] = len;
591.         }
592.     }
593.     void Sieve(int n){
594.         Sieve(n, ar, primes, len);
595.     }
596.     void init(){
597.         Sieve(MAXN - 1);
598.         int m, n, res;
599.         for (n = 0; n < PHI_N; n++){
600.             for (m = 0; m < PHI_M; m++){
601.                 if (!n) res = m;
602.                 else res = phi_dp[n - 1][m] - phi_dp[n - 1][m / primes[n - 1]];
603.                 phi_dp[n][m] = res;
604.             }
605.         }
606.     }
607.     int phi(int m, int n){
608.         if (m < PHI_M && n < PHI_N) return phi_dp[n][m];
609.         if (n == 1) return ((++m) >> 1);
610.         if (primes[n - 1] >= m) return 1;
611.         return phi(m, n - 1) - phi(m / primes[n - 1], n - 1);
612.     }
613.     int Legendre(int m){
614.         if (m < MAXN) return counter[m];
615.         int lim = sqrt(m) + 1;
616.         int n = upper_bound(primes, primes + len, lim) - primes;
617.         return phi(m, n) + (n - 1);
618.     }
619.     int Lehmer(int m){
620.         if (m < MAXN) return counter[m];
621.         int i, j, a, b, c, w, lim, res;
622.         b = sqrt(m), c = Lehmer(cbrt(m)), a = Lehmer(sqrt(b)), b = Lehmer(b);
623.         res = phi(m, a) + (((b + a - 2) * (b - a + 1)) >> 1);
624.         for (i = a; i < b; i++){
625.             w = m / primes[i];
626.             lim = Lehmer(sqrt(w)), res -= Lehmer(w);
627.             if (i <= c){
628.                 for (j = i; j < lim; j++){
629.                     res += j;
630.                     res -= Lehmer(w / primes[j]);
631.                 }
632.             }
633.         }
634.         return res;
635.     }
636. }
637. int main(){
638.     pcf::init();
639.     return 0;
640. }
641.
642.
643.
644. /// Segmented Sieve
645. #define MAX 1000010
646. #define BASE_SQR 216
647. #define BASE_LEN 10010
648. #define BASE_MAX 46656
649. #define chkbit(ar, i) (((ar[(i) >> 6]) & (1 << (((i) >> 1) & 31))))
650. #define setbit(ar, i) (((ar[(i) >> 6]) |= (1 << (((i) >> 1) & 31))))
651. int p, primes[BASE_LEN];
652. unsigned int base[(BASE_MAX >> 6) + 5], isprime[(MAX >> 6) + 5];
653. void Sieve(){
654.     clr(base);
655.     int i, j, k;
656.     for (i = 3; i < BASE_SQR; i++, i++){
657.         if (!chkbit(base, i)){
658.             k = i << 1;
659.             for (j = (i * i); j < BASE_MAX; j += k){
660.                 setbit(base, j);
661.             }
662.         }
663.     }
664.     p = 0;
665.     for (i = 3; i < BASE_MAX; i++, i++){
666.         if (!chkbit(base, i)){
667.             primes[p++] = i;
668.         }
669.     }
670. }
671. int SegmentedSieve(long long a, long long b){
672.     long long j, k, x;
673.     int i, d, counter = 0;
674.     if (a <= 2 && 2 <= b) counter = 1; /// 2 is counted separately if in range
675.     if (!(a & 1)) a++;
676.     if (!(b & 1)) b--;
677.     if (a > b) return counter;
678.     clr(isprime);
679.     for (i = 0; i < p; i++){
680.         x = primes[i];
681.         if ((x * x) > b) break;
682.         k = x << 1;
683.         j = x * ((a + x - 1) / x);
684.         if (!(j & 1)) j += x;
685.         else if (j == x) j += k;
686.         while (j <= b){
687.             setbit(isprime, j - a);
688.             j += k;
689.         }
690.     }
691.     /// Other primes in the range except 2 are added here
692.     d = (b - a + 1);
693.     for (i = 0; i < d; i++, i++){
694.         if (!chkbit(isprime, i) && (a + i) != 1) counter++;
695.     }
696.     return counter;
697. }
698. int main(){
699.     Sieve();
700.     int T = 0, t, i, j, a, b;
701.     scanf("%d", &t);
702.     while (t--){
703.         scanf("%d %d", &a, &b);
704.         printf("Case %d: %d\n", ++T, SegmentedSieve(a, b));
705.     }
706.     return 0;
707. }
708.
709.
710.
711. /// simplex
712. #define MAXC 1010
713. #define MAXV 1010
714. #define EPS 1e-13
715. #define MINIMIZE -1
716. #define MAXIMIZE +1
717. #define LESSEQ -1
718. #define EQUAL 0
719. #define GREATEQ 1
720. #define INFEASIBLE -1
721. #define UNBOUNDED 666
722. namespace lp{
723.     long double val[MAXV], ar[MAXC][MAXV];
724.     int m, n, solution_flag, minmax_flag, basis[MAXC], index[MAXV];
725.     inline void init(int nvars, long double f[], int flag){
726.         solution_flag = 0;
727.         ar[0][nvars] = 0.0;
728.         m = 0, n = nvars, minmax_flag = flag;
729.         for (int i = 0; i < n; i++){
730.             ar[0][i] = f[i] * minmax_flag;
731.         }
732.     }
733.     inline void add_constraint(long double C[], long double lim, int cmp){
734.         m++, cmp *= -1;
735.         if (cmp == 0){
736.             for (int i = 0; i < n; i++) ar[m][i] = C[i];
737.             ar[m++][n] = lim;
738.             for (int i = 0; i < n; i++) ar[m][i] = -C[i];
739.             ar[m][n] = -lim;
740.         }
741.         else{
742.             for (int i = 0; i < n; i++) ar[m][i] = C[i] * cmp;
743.             ar[m][n] = lim * cmp;
744.         }
745.     }
746.     inline void init(){
747.         for (int i = 0; i <= m; i++) basis[i] = -i;
748.         for (int j = 0; j <= n; j++){
749.             ar[0][j] = -ar[0][j], index[j] = j, val[j] = 0;
750.         }
751.     }
752.     inline void pivot(int m, int n, int a, int b){
753.         for (int i = 0; i <= m; i++){
754.             if (i != a){
755.                 for (int j = 0; j <= n; j++){
756.                     if (j != b){
757.                         ar[i][j] -= (ar[i][b] * ar[a][j]) / ar[a][b];
758.                     }
759.                 }
760.             }
761.         }
762.         for (int j = 0; j <= n; j++){
763.             if (j != b) ar[a][j] /= ar[a][b];
764.         }
765.         for (int i = 0; i <= m; i++){
766.             if (i != a) ar[i][b] = -ar[i][b] / ar[a][b];
767.         }
768.         ar[a][b] = 1.0 / ar[a][b];
769.         swap(basis[a], index[b]);
770.     }
771.     inline long double solve(){
772.         init();
773.         int i, j, k, l;
774.         for (; ;){
775.             for (i = 1, k = 1; i <= m; i++){
776.                 if ((ar[i][n] < ar[k][n]) || (ar[i][n] == ar[k][n] && basis[i] < basis[k] && (rand() & 1))) k = i;
777.             }
778.             if (ar[k][n] >= -EPS) break;
779.
780.             for (j = 0, l = 0; j < n; j++){
781.                 if ((ar[k][j] < (ar[k][l] - EPS)) || (ar[k][j] < (ar[k][l] - EPS) && index[i] < index[j] && (rand() & 1))){
782.                     l = j;
783.                 }
784.             }
785.             if (ar[k][l] >= -EPS){
786.                 solution_flag = INFEASIBLE;
787.                 return -1.0;
788.             }
789.             pivot(m, n, k, l);
790.         }
791.         for (; ;){
792.             for (j = 0, l = 0; j < n; j++){
793.                 if ((ar[0][j] < ar[0][l]) || (ar[0][j] == ar[0][l] && index[j] < index[l] && (rand() & 1))) l = j;
794.             }
795.             if (ar[0][l] > -EPS) break;
796.             for (i = 1, k = 0; i <= m; i++){
797.                 if (ar[i][l] > EPS && (!k || ar[i][n] / ar[i][l] < ar[k][n] / ar[k][l] - EPS || (ar[i][n] / ar[i][l] < ar[k][n] / ar[k][l] + EPS && basis[i] < basis[k]))){
798.                     k = i;
799.                 }
800.             }
801.             if (ar[k][l] <= EPS){
802.                 solution_flag = UNBOUNDED;
803.                 return -666.0;
804.             }
805.             pivot(m, n, k, l);
806.         }
807.         for (i = 1; i <= m; i++){
808.             if (basis[i] >= 0) val[basis[i]] = ar[i][n];
809.         }
810.         solution_flag = 1;
811.         return (ar[0][n] * minmax_flag);
812.     }
813. }
814.
815.
816. /// Thomas algorithm to solve tri-digonal system of equations in O(n)
817. /// Equation of the form: (x_prev * l) + (x_cur * p) + (x_next * r) = rhs
818. struct equation{
819.     long double l, p, r, rhs;
820.     equation(){}
821.     equation(long double l, long double p, long double r, long double rhs = 0.0):
822.         l(l), p(p), r(r), rhs(rhs){}
823. };
824. vector <long double> thomas_algorithm(int n, vector <struct equation> ar){
825.     ar[0].r = ar[0].r / ar[0].p;
826.     ar[0].rhs = ar[0].rhs / ar[0].p;
827.     for (int i = 1; i < n; i++){
828.         long double v = 1.0 / (ar[i].p - ar[i].l * ar[i - 1].r);
829.         ar[i].r = ar[i].r * v;
830.         ar[i].rhs = (ar[i].rhs - ar[i].l * ar[i - 1].rhs) * v;
831.     }
832.     for (int i = n - 2; i >= 0; i--) ar[i].rhs = ar[i].rhs - ar[i].r * ar[i + 1].rhs;
833.     vector <long double> res;
834.     for (int i = 0; i < n; i++) res.push_back(ar[i].rhs);
835.     return res;
836. }
837.
838.
839.
840. /// Karatsuba
841. #define MAX 131072 /// Must be a power of 2
842. #define MOD 1000000007
843. unsigned long long temp[128];
844. int ptr = 0, buffer[MAX * 6];
845. /// n is a power of 2
846. void karatsuba(int n, int *a, int *b, int *res){ /// hash = 829512
847.     int i, j, h;
848.     if (n < 17){ /// Reduce recursive calls by setting a threshold
849.         for (i = 0; i < (n + n); i++) temp[i] = 0;
850.         for (i = 0; i < n; i++){
851.             if (a[i]){
852.                 for (j = 0; j < n; j++){
853.                     temp[i + j] += ((long long)a[i] * b[j]);
854.                 }
855.             }
856.         }
857.         for (i = 0; i < (n + n); i++) res[i] = temp[i] % MOD;
858.         return;
859.     }
860.     h = n >> 1;
861.     karatsuba(h, a, b, res);
862.     karatsuba(h, a + h, b + h, res + n);
863.     int *x = buffer + ptr, *y = buffer + ptr + h, *z = buffer + ptr + h + h;
864.
865.     ptr += (h + h + n);
866.     for (i = 0; i < h; i++){
867.         x[i] = a[i] + a[i + h], y[i] = b[i] + b[i + h];
868.         if (x[i] >= MOD) x[i] -= MOD;
869.         if (y[i] >= MOD) y[i] -= MOD;
870.     }
871.     karatsuba(h, x, y, z);
872.     for (i = 0; i < n; i++) z[i] -= (res[i] + res[i + n]);
873.     for (i = 0; i < n; i++){
874.         res[i + h] = (res[i + h] + z[i]) % MOD;
875.         if (res[i + h] < 0) res[i + h] += MOD;
876.     }
877.     ptr -= (h + h + n);
878. }
879. /// multiplies two polynomial a(degree n) and b(degree m) and returns the result modulo MOD in a
880. /// returns the degree of the multiplied polynomial
881. /// note that a and b are changed in the process
882. int mul(int n, int *a, int m, int *b){ /// hash = 903808
883.     int i, r, c = (n < m ? n : m), d = (n > m ? n : m), *res = buffer + ptr;
884.     r = 1 << (32 - __builtin_clz(d) - (__builtin_popcount(d) == 1));
885.     for (i = d; i < r; i++) a[i] = b[i] = 0;
886.     for (i = c; i < d && n < m; i++) a[i] = 0;
887.     for (i = c; i < d && m < n; i++) b[i] = 0;
888.
889.     ptr += (r << 1), karatsuba(r, a, b, res), ptr -= (r << 1);
890.     for (i = 0; i < (r << 1); i++) a[i] = res[i];
891.     return (n + m - 1);
892. }
893. int a[MAX * 2], b[MAX * 2];
894. int main(){
895.     int i, j, k, n = MAX - 10;
896.     for (i = 0; i < n; i++) a[i] = ran(1, 1000000000);
897.     for (i = 0; i < n; i++) b[i] = ran(1, 991929183);
898.     clock_t start = clock();
899.     mul(n, a, n, b);
900.     dbg(a[n / 2]);
901.     for (i = 0; i < (n << 1); i++){
902.         if (a[i] < 0) puts("YO");
903.     }
904.     printf("%0.5f\n", (clock() - start) / (1.0 * CLOCKS_PER_SEC));
905.     return 0;
906. }
907.
908.
909. /// Integral Determinant
910. /// Finds the determinant of a square matrix
911. /// Returns 0 if the matrix is singular or degenerate (hence no determinant exists)
912. /// Absolute value of final answer should be < MOD / 2
913. #define MAX 1010
914. const long long MOD = 4517409488245517117LL;
915. const long double OP = (long double)1 / 4517409488245517117LL;
916. long long mul(long long a, long long b){
917.    long double res = a;
918.    res *= b;
919.    long long c = (long long)(res * OP);
920.    a *= b;
921.    a -= c * MOD;
922.    if (a >= MOD) a -= MOD;
923.    if (a < 0) a += MOD;
924.    return a;
925. }
926. long long expo(long long x, long long n){
927.     long long res = 1;
928.     while (n){
929.         if (n & 1) res = mul(res, x);
930.         x = mul(x, x);
931.         n >>= 1;
932.     }
933.     return res;
934. }
935. int gauss(int n, long long ar[MAX][MAX]){
936.     long long x, y;
937.     int i, j, k, l, p, counter = 0;
938.     for (i = 0; i < n; i++){
939.         for (p = i, j = i + 1; j < n && !ar[p][i]; j++){
940.             p = j;
941.         }
942.         if (!ar[p][i]) return -1;
943.         for (j = i; j < n; j++){
944.             x = ar[p][j], ar[p][j] = ar[i][j], ar[i][j] = x;
945.         }
946.         if (p != i) counter++;
947.         for (j = i + 1; j < n; j++){
948.             x = expo(ar[i][i], MOD - 2);
949.             x = mul(x, MOD - ar[j][i]);
950.             for (k = i; k < n; k++){
951.                 ar[j][k] = ar[j][k] + mul(x, ar[i][k]);
952.                 if (ar[j][k] >= MOD) ar[j][k] -= MOD;
953.             }
954.         }
955.     }
956.     return counter;
957. }
958. long long determinant(int n, long long ar[MAX][MAX]){
959.     int i, j, free;
960.     long long res = 1;
961.     for (i = 0; i < n; i++){
962.         for (j = 0; j < n; j++){
963.             if (ar[i][j] < 0) ar[i][j] += MOD;
964.         }
965.     }
966.     free = gauss(n, ar);
967.     if (free == -1) return 0; /// Determinant is 0 so matrix is not invertible, singular or degenerate matrix
968.     for (i = 0; i < n; i++) res = mul(res, ar[i][i]);
969.     if (free & 1) res = MOD - res;
970.     if ((MOD - res) < res) res -= MOD; /// Determinant can be negative so if determinant is more close to MOD than 0, make it negative
971.     return res;
972. }
973. int n;
974. long long ar[MAX][MAX];
975. int main(){
976.     int t, i, j, k, l;
977.     while (scanf("%d", &n) != EOF){
978.         if (n == 0) break;
979.         for (i = 0; i < n; i++){
980.             for (j = 0; j < n; j++){
981.                 scanf("%lld", &ar[i][j]);
982.             }
983.         }
984.         printf("%lld\n", determinant(n, ar));
985.     }
986.     return 0;
987. }
988.
989.
990.
991. /// Gaussian Elimination in Modular Field
992. #define MAXROW 1010
993. #define MAXCOL 1010
994. int expo(int a, int b, int MOD){
995.     int res = 1;
996.     while (b){
997.         if (b & 1) res = (long long)res * a % MOD;
998.         a = (long long)a * a % MOD;
999.         b >>= 1;
1000.     }
1001.     return res;
1002. }
1003. /// Gaussian elimination in field MOD (MOD should be a prime)
1004. int gauss(int n, int m, int MOD, int ar[MAXROW][MAXCOL], vector<int>& res){
1005.     res.assign(m, 0);
1006.     vector <int> pos(m, -1);
1007.     int i, j, k, l, p, d, free_var = 0;
1008.     const long long MODSQ = (long long)MOD * MOD;
1009.     for (j = 0, i = 0; j < m && i < n; j++){
1010.         for (k = i, p = i; k < n; k++){
1011.             if (ar[k][j] > ar[p][j]) p = k;
1012.         }
1013.         if (ar[p][j]){
1014.             pos[j] = i;
1015.             for (l = j; l <= m; l++) swap(ar[p][l], ar[i][l]);
1016.             d = expo(ar[i][j], MOD - 2, MOD);
1017.             for (k = 0; k < n && d; k++){
1018.                 if (k != i && ar[k][j]){
1019.                     int x = ((long long)ar[k][j] * d) % MOD;
1020.                     for (l = j; l <= m && x; l++){
1021.                         if (ar[i][l]) ar[k][l] = (MODSQ + ar[k][l] - ((long long)ar[i][l] * x)) % MOD;
1022.                     }
1023.                 }
1024.             }
1025.             i++;
1026.         }
1027.     }
1028.     for (i = 0; i < m; i++){
1029.         if (pos[i] == -1) free_var++;
1030.         else res[i] = ((long long)ar[pos[i]][m] * expo(ar[pos[i]][i], MOD - 2, MOD)) % MOD;
1031.     }
1032.     for (i = 0; i < n; i++) {
1033.         long long val = 0;
1034.         for (j = 0; j < m; j++) val = (val + ((long long)res[j] * ar[i][j])) % MOD;
1035.         if (val != ar[i][m]) return -1;
1036.     }
1037.     return free_var;
1038. }
1039.
1040.
1041.
1042. /// Gaussian Elimination Extended
1043. #define EPS 1e-9
1044. #define MAXROW 512
1045. #define MAXCOL 512
1046.
1047. /***
1048.
1049. Gauss-Jordan Elimination
1050.
1051. n = number of linear equations
1052. m = number of variables
1053. ar[i][m] = right-hand side value of constants
1054.
1055. For instance, the system of linear equations becomes:
1056.
1057. 2x + y -z = 8      ----->  (i)
1058. -3x -y + 2z = -11  ----->  (ii)
1059. -2x + y + 2z = -3  ----->  (iii)
1060.
1061. n = 3 (x, y, z), m = 3 (i, ii, iii)
1062. ar[0] = {2, 1, -1, 8}    ----->  (i)
1063. ar[1] = {-3, -1, 2, -11} ----->  (ii)
1064. ar[2] = {-2, 1, 2, -3}   ----->  (iii)
1065.
1066.
1067. Returns -1 when there is no solution
1068. Otherwise returns the number of independent variables (0 for an unique solution)
1069. Contains a solution in the vector res on successful completion
1070. Note that the array is modified in the process
1071.
1072. Notes:
1073. For solving problems on graphs with probability/expectation, make sure the graph
1074. is connected and a single component. If not, then re-number the vertex and solve
1075. for each connected component separately.
1076.
1077. ***/
1078. int gauss(int n, int m, double ar[MAXROW][MAXCOL], vector<double>& res){ /// hash = 835176
1079.     res.assign(m, 0);
1080.     vector <int> pos(m, -1);
1081.     int i, j, k, l, p, free_var = 0;
1082.     for (j = 0, i = 0; j < m && i < n; j++){
1083.         for (k = i, p = i; k < n; k++){
1084.             if (abs(ar[k][j]) > abs(ar[p][j])) p = k;
1085.         }
1086.         if (abs(ar[p][j]) > EPS){
1087.             pos[j] = i;
1088.             for (l = j; l <= m; l++) swap(ar[p][l], ar[i][l]);
1089.             for (k = 0; k < n; k++){
1090.                 if (k != i){
1091.                     double x = ar[k][j] / ar[i][j];
1092.                     for (l = j; l <= m; l++) ar[k][l] -= (ar[i][l] * x);
1093.                 }
1094.             }
1095.             i++;
1096.         }
1097.     }
1098.     for (i = 0; i < m; i++){
1099.         if (pos[i] == -1) free_var++;
1100.         else res[i] = ar[pos[i]][m] / ar[pos[i]][i];
1101.     }
1102.     for (i = 0; i < n; i++) {
1103.         double val = 0.0;
1104.         for (j = 0; j < m; j++) val += (res[j] * ar[i][j]);
1105.         if (abs(val - ar[i][m]) > EPS) return -1;
1106.     }
1107.     return free_var;
1108. }
1109.
1110.
1111.
1112. /// FFT extended
1113. #define MAXN 1048576 /// 2 * MAX at least
1114. /// Change long double to double if not required
1115. namespace fft{
1116.     int len, last = -1, step = 0, rev[MAXN];
1117.     long long C[MAXN], D[MAXN], P[MAXN], Q[MAXN];
1118.     struct complx{
1119.         long double real, img;
1120.         inline complx(){
1121.             real = img = 0.0;
1122.         }
1123.         inline complx conjugate(){
1124.             return complx(real, -img);
1125.         }
1126.         inline complx(long double x){
1127.             real = x, img = 0.0;
1128.         }
1129.         inline complx(long double x, long double y){
1130.             real = x, img = y;
1131.         }
1132.         inline complx operator + (complx other){
1133.             return complx(real + other.real, img + other.img);
1134.         }
1135.         inline complx operator - (complx other){
1136.             return complx(real - other.real, img - other.img);
1137.         }
1138.         inline complx operator * (complx other){
1139.             return complx((real * other.real) - (img * other.img), (real * other.img) + (img * other.real));
1140.         }
1141.     } u[MAXN], v[MAXN], f[MAXN], g[MAXN], dp[MAXN], inv[MAXN];
1142.     /// Pre-process roots, inverse roots and fft leaf index
1143.     void build(int& a, long long* A, int& b, long long* B){
1144.         while (a > 1 && A[a - 1] == 0) a--;
1145.         while (b > 1 && B[b - 1] == 0) b--;
1146.         len = 1 << (32 - __builtin_clz(a + b) - (__builtin_popcount(a + b) == 1));
1147.         for (int i = a; i < len; i++) A[i] = 0;
1148.         for (int i = b; i < len; i++) B[i] = 0;
1149.         if (!step++){
1150.             dp[1] = inv[1] = complx(1);
1151.             for (int i = 1; (1 << i) < MAXN; i++){
1152.                 double theta = (2.0 * acos(0.0)) / (1 << i);
1153.                 complx mul = complx(cos(theta), sin(theta));
1154.                 complx inv_mul = complx(cos(-theta), sin(-theta));
1155.
1156.                 int lim = 1 << i;
1157.                 for (int j = lim >> 1; j < lim; j++){
1158.                     dp[2 * j] = dp[j], inv[2 * j] = inv[j];
1159.                     inv[2 * j + 1] = inv[j] * inv_mul;
1160.                     dp[2 * j + 1] = dp[j] * mul;
1161.                 }
1162.             }
1163.         }
1164.         if (last != len){
1165.             last = len;
1166.             int bit = (32 - __builtin_clz(len) - (__builtin_popcount(len) == 1));
1167.             for (int i = 0; i < len; i++) rev[i] = (rev[i >> 1] >> 1) + ((i & 1) << (bit - 1));
1168.         }
1169.     }
1170.     /// Fast Fourier Transformation, iterative divide and conquer
1171.     void transform(complx *in, complx *out, complx* ar){
1172.         for (int i = 0; i < len; i++) out[i] = in[rev[i]];
1173.         for (int k = 1; k < len; k <<= 1){
1174.             for (int i = 0; i < len; i += (k << 1)){
1175.                 for (int j = 0; j < k; j++){
1176.                     complx z = out[i + j + k] * ar[j + k];
1177.                     out[i + j + k] = out[i + j] - z;
1178.                     out[i + j] = out[i + j] + z;
1179.                 }
1180.             }
1181.         }
1182.     }
1183.     /// Fast Fourier Transformation, iterative divide and conquer unrolled and optimized
1184.     void transform_unrolled(complx *in, complx *out, complx* ar){
1185.         for (int i = 0; i < len; i++) out[i] = in[rev[i]];
1186.         for (int k = 1; k < len; k <<= 1){
1187.             for (int i = 0; i < len; i += (k << 1)){
1188.                 complx z, *a = out + i, *b = out + i + k, *c = ar + k;
1189.                 if (k == 1){
1190.                     z = (*b) * (*c);
1191.                     *b = *a - z, *a = *a + z;
1192.                 }
1193.                 for (int j = 0; j < k && k > 1; j += 2, a++, b++, c++){
1194.                     z = (*b) * (*c);
1195.                     *b = *a - z, *a = *a + z;
1196.                     a++, b++, c++;
1197.                     z = (*b) * (*c);
1198.                     *b = *a - z, *a = *a + z;
1199.                 }
1200.             }
1201.         }
1202.     }
1203.     bool equals(int a, long long* A, int b, long long* B){
1204.         if (a != b) return false;
1205.         for (a = 0; a < b && A[a] == B[a]; a++){}
1206.         return (a == b);
1207.     }
1208.     /// Square of a polynomial
1209.     int square(int a, long long* A){
1210.         build(a, A, a, A);
1211.         for (int i = 0; i < len; i++) u[i] = complx(A[i], 0);
1212.         transform_unrolled(u, f, dp);
1213.         for (int i = 0; i < len; i++) u[i] = f[i] * f[i];
1214.         transform_unrolled(u, f, inv);
1215.         for (int i = 0; i < len; i++) A[i] = (f[i].real / (long double)len) + 0.5;
1216.         return a + a - 1;
1217.     }
1218.     /// Multiplies two polynomials A and B and return the coefficients of their product in A
1219.     /// Function returns degree of the polynomial A * B
1220.     int multiply(int a, long long* A, int b, long long* B){
1221.         if (equals(a, A, b, B)) return square(a, A); /// Optimization
1222.         build(a, A, b, B);
1223.         for (int i = 0; i < len; i++) u[i] = complx(A[i], B[i]);
1224.         transform_unrolled(u, f, dp);
1225.         for (int i = 0; i < len; i++){
1226.             int j = (len - 1) & (len - i);
1227.             u[i] = (f[j] * f[j] - f[i].conjugate() * f[i].conjugate()) * complx(0, -0.25 / len);
1228.         }
1229.         transform_unrolled(u, f, dp);
1230.         for (int i = 0; i < len; i++) A[i] = f[i].real + 0.5;
1231.         return a + b - 1;
1232.     }
1233.     /// Modular multiplication
1234.     int mod_multiply(int a, long long* A, int b, long long* B, int mod){
1235.         build(a, A, b, B);
1236.         int flag = equals(a, A, b, B);
1237.         for (int i = 0; i < len; i++) A[i] %= mod, B[i] %= mod;
1238.         for (int i = 0; i < len; i++) u[i] = complx(A[i] & 32767, A[i] >> 15);
1239.         for (int i = 0; i < len; i++) v[i] = complx(B[i] & 32767, B[i] >> 15);
1240.         transform_unrolled(u, f, dp);
1241.         for (int i = 0; i < len; i++) g[i] = f[i];
1242.         if (!flag) transform_unrolled(v, g, dp);
1243.         for (int i = 0; i < len; i++){
1244.             int j = (len - 1) & (len - i);
1245.             complx c1 = f[j].conjugate(), c2 = g[j].conjugate();
1246.             complx a1 = (f[i] + c1) * complx(0.5, 0);
1247.             complx a2 = (f[i] - c1) * complx(0, -0.5);
1248.             complx b1 = (g[i] + c2) * complx(0.5 / len, 0);
1249.             complx b2 = (g[i] - c2) * complx(0, -0.5 / len);
1250.             v[j] = a1 * b2 + a2 * b1;
1251.             u[j] = a1 * b1 + a2 * b2 * complx(0, 1);
1252.         }
1253.         transform_unrolled(u, f, dp);
1254.         transform_unrolled(v, g, dp);
1255.         long long x, y, z;
1256.         for (int i = 0; i < len; i++){
1257.             x = f[i].real + 0.5, y = g[i].real + 0.5, z = f[i].img + 0.5;
1258.             A[i] = (x + ((y % mod) << 15) + ((z % mod) << 30)) % mod;
1259.         }
1260.         return a + b - 1;
1261.     }
1262.     /// Multiplies two polynomials where intermediate and final values fits in long long
1263.     int long_multiply(int a, long long* A, int b, long long* B){
1264.         int mod1 = 1.5e9;
1265.         int mod2 = mod1 + 1;
1266.         for (int i = 0; i < a; i++) C[i] = A[i];
1267.         for (int i = 0; i < b; i++) D[i] = B[i];
1268.         mod_multiply(a, A, b, B, mod1);
1269.         mod_multiply(a, C, b, D, mod2);
1270.         for (int i = 0; i < len; i++){
1271.             A[i] = A[i] + (C[i] - A[i] + (long long)mod2) * (long long)mod1 % mod2 * mod1;
1272.         }
1273.         return a + b - 1;
1274.     }
1275.     int build_convolution(int n, long long* A, long long* B){
1276.         int i, m, d = 0;
1277.         for (i = 0; i < n; i++) Q[i] = Q[i + n] = B[i];
1278.         for (i = 0; i < n; i++) P[i] = A[i], P[i + n] = 0;
1279.         n *= 2, m = 1 << (32 - __builtin_clz(n) - (__builtin_popcount(n) == 1));
1280.         for (i = n; i < m; i++) P[i] = Q[i] = 0;
1281.         return n;
1282.     }
1283.     /***
1284.         Computes the circular convolution of A and B, denoted A * B, in C
1285.         A and B must be of equal size, if not normalize before calling function
1286.         Example to demonstrate convolution for n = 5:
1287.
1288.         c0 = a0b0 + a1b4 + a2b3 + a3b2 + a4b1
1289.         c1 = a0b1 + a1b0 + a2b4 + a3b3 + a4b2
1290.         ...
1291.         c4 = a0b4 + a1b3 + a2b2 + a3b1 + a4b0
1292.
1293.
1294.         Note: If linear convolution is required, pad with zeros appropriately, as in multiplication
1295.
1296.     ***/
1297.     /// Returns the convolution of A and B in A
1298.     void convolution(int n, long long* A, long long* B){
1299.         int len = build_convolution(n, A, B);
1300.         multiply(len, P, len, Q);
1301.         for (int i = 0; i < n; i++) A[i] = P[i + n];
1302.     }
1303.     /// Modular convolution
1304.     void mod_convolution(int n, long long* A, long long* B, int mod){
1305.         int len = build_convolution(n, A, B);
1306.         mod_multiply(len, P, len, Q, mod);
1307.         for (int i = 0; i < n; i++) A[i] = P[i + n];
1308.     }
1309.     /// Convolution in long long
1310.     void long_convolution(int n, long long* A, long long* B){
1311.         int len = build_convolution(n, A, B);
1312.         long_multiply(len, P, len, Q);
1313.         for (int i = 0; i < n; i++) A[i] = P[i + n];
1314.     }
1315.     /// Hamming distance vector of B with every substring of length |pattern| in str
1316.     /// str and pattern consists of only '1' and '0'
1317.     /// str = "01111000010011111111110010001101000100011110101111"
1318.     /// pattern = "1001101001101110101101000"
1319.     /// Sum of values in hamming distance vector = 321
1320.     vector <int> hamming_distance(const char* str, const char* pattern){
1321.         int n = strlen(str), m = strlen(pattern);
1322.         for (int i = 0; i < n; i++) P[i] = Q[i] = 0;
1323.         for (int i = 0; i < n; i++) P[i] = str[i] == '1' ? 1 : -1;
1324.         for (int i = 0, j = m - 1; j >= 0; i++, j--) Q[i] = pattern[j] == '1' ? 1 : -1;
1325.         vector <int> res;
1326.         fft::multiply(n, P, m, Q);
1327.         for (int i = 0; (i + m) <= n; i++){
1328.             res.push_back(m - ((P[i + m - 1] + m) >> 1));
1329.         }
1330.         return res;
1331.     }
1332. }
1333.
1334.
1335.
1336. /// Fast Integer Cube and Square Root
1337. unsigned int fast_sqrt(unsigned int n){
1338.     unsigned int c, g;
1339.     c = g = 0x8000;
1340.     for (; ;){
1341.         if ((g * g) > n) g ^= c;
1342.         c >>= 1;
1343.         if (!c) return g;
1344.         g |= c;
1345.     }
1346. }
1347. int fast_cbrt(int n){
1348.     int x, r = 30, res = 0;
1349.     for (; r >= 0; r -= 3){
1350.         res <<= 1;
1351.         x = (3 * res * (res + 1)) + 1;
1352.         if ((n >> r) >= x){
1353.             res++;
1354.             n -= (x << r);
1355.         }
1356.     }
1357.     return res;
1358. }
1359. unsigned long long fast_sqrt(unsigned long long n){
1360.     unsigned long long c, g;
1361.     c = g = 0x80000000;
1362.     for (; ;){
1363.         if ((g * g) > n) g ^= c;
1364.         c >>= 1;
1365.         if (!c) return g;
1366.         g |= c;
1367.     }
1368. }
1369. unsigned long long fast_cbrt(unsigned long long n){
1370.     int r = 63;
1371.     unsigned long long x, res = 0;
1372.     for (; r >= 0; r -= 3){
1373.         res <<= 1;
1374.         x = (res * (res + 1) * 3) + 1;
1375.         if ((n >> r) >= x){
1376.             res++;
1377.             n -= (x << r);
1378.         }
1379.     }
1380.     return res;
1381. }
1382.
1383.
1384.
1385. /// Extended Euclid
1386. /// Bezout's identity, ax + by = gcd(a,b)
1387. int extended_gcd(int a, int b, int& x, int& y){
1388.     if (!b){
1389.         y = 0, x = 1;
1390.         return a;
1391.     }
1392.     int g = extended_gcd(b, a % b, y, x);
1393.     y -= ((a / b) * x);
1394.     return g;
1395. }
1396. /// Linear diophantine equation, ax + by = c
1397. void diophantine(int a, int b, int c, int& x, int& y){
1398.     int g = extended_gcd(a, b, x, y);
1399.     assert((c % g) == 0); /// c must be a multiply of g
1400.     c /= g;
1401.     x *= c, y *= c;
1402. }
1403.
1404.
1405.
1406. /// Eulerian Numbers
1407. #define MAX 5010
1408. #define MOD 1000000007
1409. /***
1410. Eulerian number A(n,k) is the number of permutations of 1 to n in which exactly k elements are greater than their previous element
1411. Eulerian triangle for n = 1 to 7 and k = 0 to n - 1 below
1412.
1413. 1
1414. 1 1
1415. 1 4 1
1416. 1 11 11 1
1417. 1 26 66 26 1
1418. 1 57 302 302 57 1
1419. 1 120 1191 2416 1191 120 1
1420.
1421. ***/
1422. int dp[MAX][MAX];
1423. void generate(){
1424.     int n, k;
1425.     for (n = 0; n < MAX; n++){
1426.         for (k = 1, dp[n][0] = 1; k <= n; k++){
1427.             dp[n][k] = ((long long)dp[n - 1][k] * (k + 1) + (long long)dp[n - 1][k - 1] * (n - k)) % MOD;
1428.         }
1429.     }
1430. }
1431. int main(){
1432.     int n, k;
1433.     generate();
1434.     for (n = 1; n <= 7; n++){
1435.         for (k = 0; k < n; k++){
1436.             printf("%d ", dp[n][k]);
1437.         }
1438.         puts("");
1439.     }
1440.     return 0;
1441. }
1442.
1443.
1444.
1445. /// Discrete Log (General)
1446. int extended_gcd(int a, int b, int& x, int& y){
1447.     /// Bezout's identity, ax + by = gcd(a,b)
1448.     if (!b){
1449.         y = 0, x = 1;
1450.         return a;
1451.     }
1452.     int g = extended_gcd(b, a % b, y, x);
1453.     y -= ((a / b) * x);
1454.     return g;
1455. }
1456. int discrete_log(int g, int h, int p){ /// hash = 167626
1457.     /***
1458.         * returns smallest x such that (g^x) % p = h, -1 if none exists
1459.         * function returns x, the discrete log of h with respect to g modulo p
1460.     ***/
1461.     if (h >= p) return -1;
1462.     if ((g % p) == 0){
1463.         if (h == 1) return 0; /// return -1 if strictly positive integer solution is required
1464.         else return -1;
1465.     }
1466.     int i, c, x, y, z, r, m, counter = 0;
1467.     long long v = 1, d = 1, mul = 1, temp = 1 % p;
1468.     for (int i = 0; i < 100; i++){
1469.         if (temp == h) return i;
1470.         temp = (temp * g) % p;
1471.     }
1472.     while ((v = __gcd(g, p)) > 1){
1473.         if (h % v) return -1;
1474.         h /= v, p /= v;
1475.         d = (d * (g / v)) % p;
1476.         counter++;
1477.     }
1478.     m = ceil(sqrt(p)); /// may be change to sqrtl() ?
1479.     tr1::unordered_map <int, int> mp;
1480.     for (i = 0; i < m; i++){
1481.         if (!mp[mul]) mp[mul] = i + 1;
1482.         mul = (mul * g) % p;
1483.     }
1484.     for (i = 0; i < m; i++){
1485.         z = extended_gcd(d, p, x, y);
1486.         c = p / z;
1487.         r = ((((long long)x * h) / z) % p + p) % p;
1488.         if (mp[r]) return ((i * m) + mp[r] + counter - 1);
1489.         d = (d * mul) % p;
1490.     }
1491.     return -1;
1492. }
1493. int main(){
1494.     int g, h, p, res;
1495.     for (; ;){
1496.         scanf("%d %d %d", &g, &p, &h);
1497.         if (!g && !p && !h) break;
1498.         res = discrete_log(g, h % p, p);
1499.         if (res == -1) puts("No Solution");
1500.         else printf("%d\n", res);
1501.     }
1502.     return 0;
1503. }
1504.
1505.
1506. /// Discrete Log
1507. int expo(long long x, int n, int m){
1508.     long long res = 1;
1509.
1510.     while (n){
1511.         if (n & 1) res = (res * x) % m;
1512.         x = (x * x) % m;
1513.         n >>= 1;
1514.     }
1515.     return (res % m);
1516. }
1517. int extended_gcd(int a, int b, int& x, int& y){
1518.     /// Bezout's identity, ax + by = gcd(a,b)
1519.     if (!b){
1520.         y = 0, x = 1;
1521.         return a;
1522.     }
1523.     int g = extended_gcd(b, a % b, y, x);
1524.     y -= ((a / b) * x);
1525.     return g;
1526. }
1527. int inverse_modulo(int a, int m){
1528.     /// inverse exists if and only if a and m are co-prime
1529.     int x, y, inv;
1530.     extended_gcd(a, m, x, y);
1531.     inv = (x + m) % m;
1532.     return inv;
1533. }
1534. int discrete_log(int g, int h, int p){
1535.     /***
1536.         * returns smallest x such that (g^x) % p = h, -1 if none exists
1537.         * p must be a PRIME
1538.         * function returns x, the discrete log of h with respect to g modulo p
1539.     ***/
1540.     if (h >= p) return -1;
1541.     if ((g % p) == 0){
1542.         if (h == 1) return 0; /// return -1 if strictly positive integer solution is required
1543.         else return -1;
1544.     }
1545.     unordered_map <int, int> mp;
1546.     int i, q, r, m = ceil(sqrt(p)); /// may be change to sqrtl() ?
1547.     long long d = 1, inv = expo(inverse_modulo(g, p), m, p);
1548.     for (i = 0; i <= m; i++){
1549.         if (!mp[d]) mp[d] = i + 1;
1550.         d *= g;
1551.         if (d >= p) d %= p;
1552.     }
1553.     d = h;
1554.     for (q = 0; q <= m; q++){
1555.         r = mp[d];
1556.         if (r) return ((m * q) + (--r));
1557.         d *= inv;
1558.         if (d >= p) d %= p;
1559.     }
1560.     return -1;
1561. }
1562. int main(){
1563.     int T = 0, t, g, h, p;
1564.     scanf("%d", &t);
1565.     while (t--){
1566.         scanf("%d %d %d", &g, &h, &p);
1567.         int x = discrete_log(g, h, p);
1568.         printf("Case %d: %d\n", ++T, x);
1569.     }
1570.     return 0;
1571. }
1572.
1573.
1574.
1575. /// Digits of Factorial
1576. const double A = 4.0 * acos(0.0);
1577. const double B = log10(exp(1.0));
1578. long long digfact(long long n){
1579.     if (n == 0 || n == 1) return 1;
1580.     double x = ((0.5 * log(A * n)) + (n * log(n)) - n) * B;
1581.     return ceil(x);
1582. }
1583. int main(){
1584.     int t;
1585.     long long n;
1586.     scanf("%d", &t);
1587.     while (t--){
1588.         scanf("%lld", &n);
1589.         printf("%lld\n", digfact(n));
1590.     }
1591.     return 0;
1592. }
1593.
1594.
1595.
1596. /// Chinese Remainder Theorem
1597. #define dbg(x) cout << #x << " = " << x << endl
1598. using namespace std;
1599. long long extended_gcd(long long a, long long b, long long& x, long long& y){
1600.     /// Bezout's identity, ax + by = gcd(a,b)
1601.     if (!b){
1602.         y = 0, x = 1;
1603.         return a;
1604.     }
1605.     long long g = extended_gcd(b, a % b, y, x);
1606.     y -= ((a / b) * x);
1607.     return g;
1608. }
1609. long long mod_inverse(long long a, long long m){
1610.     /// inverse exists if and only if a and m are co-prime
1611.     long long x, y, inv;
1612.     extended_gcd(a, m, x, y);
1613.     inv = (x + m) % m;
1614.     return inv;
1615. }
1616. long long CRT(int n, long long* ar, long long* mods){
1617.     /// finds the unique solution x modulo M (product of mods) for which x % mods[i] = ar[i]
1618.     /// mods must be pairwise co-prime
1619.     int i, j;
1620.     long long x, y, res = 0, M = 1;
1621.     for (i = 0; i < n; i++) M *= mods[i];
1622.     for (i = 0; i < n; i++){
1623.         x = M / mods[i];
1624.         y = mod_inverse(x, mods[i]);
1625.         res = (res + (((x * ar[i]) % M) * y)) % M;
1626.     }
1627.     return res;
1628. }
1629. int main(){
1630.     long long mods[] = {3, 5, 7};
1631.     long long ar[] = {2, 3, 2};
1632.     long long res = CRT(3, ar, mods);
1633.     dbg(res);
1634. }
1635.
1636.
1637.
1638. /// Binomial Coefficients Modulo Prime Powers
1639. #define MAXP 10000010
1640. #define ran(a, b) ((((rand() << 15) ^ rand()) % ((b) - (a) + 1)) + (a))
1641. using namespace std;
1642. namespace crt{
1643.     long long extended_gcd(long long a, long long b, long long& x, long long& y){
1644.         if (!b){
1645.             y = 0, x = 1;
1646.             return a;
1647.         }
1648.
1649.         long long g = extended_gcd(b, a % b, y, x);
1650.         y -= ((a / b) * x);
1651.         return g;
1652.     }
1653.     long long mod_inverse(long long a, long long m){
1654.         long long x, y, inv;
1655.         extended_gcd(a, m, x, y);
1656.         inv = (x + m) % m;
1657.         return inv;
1658.     }
1659.     long long chinese_remainder(vector <long long> ar, vector <long long> mods){
1660.         int i, j;
1661.         long long x, y, res = 0, M = 1;
1662.         for (i = 0; i < ar.size(); i++) M *= mods[i];
1663.         for (i = 0; i < ar.size(); i++){
1664.             x = M / mods[i];
1665.             y = mod_inverse(x, mods[i]);
1666.             res = (res + (((x * ar[i]) % M) * y)) % M;
1667.         }
1668.         return res;
1669.     }
1670. }
1671. namespace bin{
1672.     int dp[MAXP];
1673.     long long mod = 0;
1674.     long long trailing(long long x, long long p){
1675.         long long res = 0;
1676.         while (x){
1677.             x /= p;
1678.             res += x;
1679.         }
1680.         return res;
1681.     }
1682.     long long expo(long long x, long long n, long long m){
1683.         if (!n) return 1;
1684.         else if (n & 1) return ((expo(x, n - 1, m) * x) % m);
1685.         else{
1686.             long long r = expo(x, n >> 1, m);
1687.             return ((r * r) % m);
1688.         }
1689.     }
1690.     long long factorial(long long x, long long p){
1691.         long long res = expo(dp[mod - 1], x / mod, mod);
1692.         if (x >= p) res = res * factorial(x / p, p) % mod;
1693.         return res * dp[x % mod] % mod;
1694.     }
1695.     /// binomial co-efficients modulo p^q (p must be a prime)
1696.     long long binomial(long long n, long long k, long long p, long long q){
1697.         if (k > n) return 0;
1698.         if (n == k || k == 0) return 1;
1699.         int i, j;
1700.         for (i = 0, mod = 1; i < q; i++) mod *= p;
1701.         long long t = trailing(n, p) - trailing(k, p) - trailing(n - k, p);
1702.         if (t >= q) return 0;
1703.         assert(mod < MAXP);
1704.         for (dp[0] = 1, i = 1; i < mod; i++){
1705.             dp[i] = (long long)dp[i - 1] * (i % p ? i : 1) % mod;
1706.         }
1707.         long long res = factorial(n, p) * expo(factorial(k, p) * factorial(n - k, p) % mod, (mod / p) * (p - 1) - 1, mod) % mod;
1708.         while (t--) res = res * p % mod;
1709.         return res;
1710.     }
1711.     /// binomial co-efficients modulo m (p_i^q_i of m must be less than MAXP)
1712.     long long binomial(long long n, long long k, long long m){
1713.         if (k > n || m == 1) return 0;
1714.         if (n == k || k == 0) return 1;
1715.         vector < pair <int, int > > factors;
1716.         for (long long i = 2; i * i <= m; i++){
1717.             int c = 0;
1718.             while (m % i == 0){
1719.                 c++;
1720.                 m /= i;
1721.             }
1722.             if (c) factors.push_back(make_pair(i, c));
1723.         }
1724.         if (m > 1) factors.push_back(make_pair(m, 1));
1725.         vector <long long> ar, mods;
1726.         for (int i = 0; i < factors.size(); i++){
1727.             long long x = 1;
1728.             for (int j = 0; j < factors[i].second; j++) x *= factors[i].first;
1729.             mods.push_back(x), ar.push_back(binomial(n, k, factors[i].first, factors[i].second));
1730.         }
1731.         return crt::chinese_remainder(ar, mods);
1732.     }
1733. }
1734.
1735.
1736.
1737. ///All Divisors in Range
1738. #define MAX 1000010
1739. short len[MAX];
1740. int L[MAX], *divisors[MAX];
1741. void Generate(){
1742.     int i, j, k, c, d, x;
1743.     len[0] = len[1] = L[0] = L[1] = 1;
1744.     for (i = 4; i < MAX; i++, i++) len[i] = 2;
1745.     for (i = 3; (i * i) < MAX; i += 2){
1746.         if (!len[i]){
1747.             for (j = (i * i), d = i << 1; j < MAX; j += d){
1748.                 if (!len[j]) len[j] = i;
1749.             }
1750.         }
1751.     }
1752.     for (i = 2; i < MAX; i++){
1753.         if (!len[i]) L[i] = i;
1754.         else{
1755.             L[i] = len[i];
1756.             x = L[i /len[i]];
1757.             if (x > L[i]) L[i] = x;
1758.         }
1759.     }
1760.     divisors[1] = (int *) malloc(sizeof(int));
1761.     divisors[1][0] = 1, len[1] = 1;
1762.     for (i = 2; i < MAX; i++){
1763.         c = 0, k = i;
1764.         while (k > 1 && L[k] == L[i]){
1765.             c++;
1766.             k /= L[k];
1767.         }
1768.         len[i] = (c + 1) * len[k];
1769.         divisors[i] = (int *) malloc(sizeof(int) * len[i]);
1770.         for (x = 1, j = 0, len[i] = 0; j <= c; j++, x *= L[i]){
1771.             for (d = 0; d < len[k]; d++){
1772.                 divisors[i][len[i]++] = divisors[k][d] * x;
1773.             }
1774.         }
1775.     }
1776. }
1777. int main(){
1778.     Generate();
1779.     return 0;
1780. }
1781.
1782.
1783.
1784. /// Bernoulli Numbers
1785. #define MAX 2510
1786. #define MOD 1000000007
1787. int S[MAX][MAX], inv[MAX], factorial[MAX], bernoulli[MAX];
1788. int expo(long long x, int n){
1789.     long long res = 1;
1790.     while (n){
1791.         if (n & 1) res = (res * x) % MOD;
1792.         x = (x * x) % MOD;
1793.         n >>= 1;
1794.     }
1795.     return (res % MOD);
1796. }
1797. void Generate(){
1798.     int i, j;
1799.     long long x, y, z, lim = (long long)MOD * MOD;
1800.     for (i = 1, factorial[0] = 1; i < MAX; i++) factorial[i] = ((long long) factorial[i - 1] * i) % MOD;
1801.     for (i = 0; i < MAX; i++) inv[i] = expo(i, MOD - 2);
1802.     for (i = 1, S[0][0] = 1; i < MAX; i++){
1803.         for (j = 1, S[i][0] = 0; j <= i; j++) S[i][j] = ( ((long long)S[i - 1][j] * j) + S[i - 1][j - 1]) % MOD;
1804.     }
1805.     bernoulli[0] = 1;
1806.     for (i = 1; (i + 1) < MAX; i++){
1807.         if ((i & 1) && i > 1) bernoulli[i] = 0;
1808.         else{
1809.             for (j = 0, x = 0, y = 0; j <= i; j++){
1810.                 z = ((long long) factorial[j] * inv[j + 1]) % MOD;
1811.                 z = (z * S[i][j]) % MOD;
1812.                 if (j & 1) y += z;
1813.                 else x += z;
1814.             }
1815.             bernoulli[i] = (lim + x - y) % MOD;
1816.         }
1817.     }
1818. }
1819. int main(){
1820.     Generate();
1821.     printf("%d\n", bernoulli[10]);
1822.     return 0;
1823. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top