API tools faq
paste
Advertisement
_no0B

Untitled

_no0B
Apr 25th, 2019
223
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 54.57 KB | None | 0 0
raw download clone embed print report
  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. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!