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. OK, I Understand
 
Top