This week only. Pastebin PRO Accounts Christmas Special! Don't miss out!Want more features on Pastebin? Sign Up, it's FREE!
Guest

Untitled

By: a guest on Jan 16th, 2013  |  syntax: None  |  size: 4.50 KB  |  views: 59  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. m = ∏ (p_i ^ e_i)
  2.        
  3. F(p, n) =   ∏ k   = p^(n/p) * (n/p)!
  4.          1<=k<=n
  5.            p | k
  6.        
  7. G(p, n) =    ∏ k
  8.           1<=k<=n
  9.          gcd(k,n)=1
  10.        
  11. n! = F(p, n) * G(p, n)
  12.        
  13. m
  14. n! = p^K * ∏ G(p, n/(p^j))
  15.           j=0
  16.        
  17. n! / (r! * (n-r)!) = p^M * (∏ G(p, (n/p^j)) * [ ∏ G(p, r/(p^j)) * ∏ G(p, (n-r)/(p^j)) ]^(-1)
  18.        
  19. H(p, e, n) ≡ (-1)^(n/(p^e)) * H(p, e, n % (p^e)) (mod p^e)
  20.        
  21. // invert k modulo p, k and p are supposed coprime
  22. unsigned long long invertMod(unsigned long long k, unsigned long long p) {
  23.     unsigned long long q, pn = 1, po = 0, r = p, s = k;
  24.     unsigned odd = 1;
  25.     do {
  26.         q = r/s;
  27.         q = pn*q + po;
  28.         po = pn;
  29.         pn = q;
  30.         q = r%s;
  31.         r = s;
  32.         s = q;
  33.         odd ^= 1;
  34.     }while(pn < p);
  35.     return odd ? p-po : po;
  36. }
  37.  
  38. // Calculate the binomial coefficient (n choose k) modulo (prime^exponent)
  39. // requires prime to be a prime, exponent > 0, and 0 <= k <= n,
  40. // furthermore supposes prime^exponent < 2^32, otherwise intermediate
  41. // computations could have mathematical results out of range.
  42. // If k or (n-k) is small, a direct computation would be more efficient.
  43. unsigned long long binmod(unsigned long long prime, unsigned exponent,
  44.                           unsigned long long n, unsigned long long k) {
  45.     // The modulus, prime^exponent
  46.     unsigned long long ppow = 1;
  47.     // We suppose exponent is small, so that exponentiation by repeated
  48.     // squaring wouldn't gain much.
  49.     for(unsigned i = 0; i < exponent; ++i) {
  50.         ppow *= prime;
  51.     }
  52.     // array of remainders of products
  53.     unsigned long long *remainders = malloc(ppow * sizeof *remainders);
  54.     if (!remainders) {
  55.         fprintf(stderr, "Allocation failuren");
  56.         exit(EXIT_FAILURE);
  57.     }
  58.     for(unsigned long long i = 1; i < ppow; ++i) {
  59.         remainders[i] = i;
  60.     }
  61.     for(unsigned long long i = 0; i < ppow; i += prime) {
  62.         remainders[i] = 1;
  63.     }
  64.     for(unsigned long long i = 2; i < ppow; ++i) {
  65.         remainders[i] *= remainders[i-1];
  66.         remainders[i] %= ppow;
  67.     }
  68.  
  69.     // Now to business.
  70.     unsigned long long pmult = 0, ntemp = n, ktemp = k, mtemp = n-k,
  71.                        numer = 1, denom = 1, q, r, f;
  72.     if (prime == 2 && exponent > 2) {
  73.         f = 0;
  74.     } else {
  75.         f = 1;
  76.     }
  77.     while(ntemp) {
  78.         r = ntemp % ppow;
  79.         q = ntemp / ppow;
  80.         numer *= remainders[r];
  81.         numer %= ppow;
  82.         if (q & f) {
  83.             numer = ppow - numer;
  84.         }
  85.         ntemp /= prime;
  86.         pmult += ntemp;
  87.     }
  88.     while(ktemp) {
  89.         r = ktemp % ppow;
  90.         q = ktemp / ppow;
  91.         denom *= remainders[r];
  92.         denom %= ppow;
  93.         if (q & f) {
  94.             denom = ppow - denom;
  95.         }
  96.         ktemp /= prime;
  97.         pmult -= ktemp;
  98.     }
  99.     while(mtemp) {
  100.         r = mtemp % ppow;
  101.         q = mtemp / ppow;
  102.         denom *= remainders[r];
  103.         denom %= ppow;
  104.         if (q & f) {
  105.             denom = ppow - denom;
  106.         }
  107.         mtemp /= prime;
  108.         pmult -= mtemp;
  109.     }
  110.     // free memory before returning, we don't use it anymore
  111.     free(remainders);
  112.     if (pmult >= exponent) {
  113.         return 0;
  114.     }
  115.     while(pmult > 0) {
  116.         numer = (numer * prime) % ppow;
  117.         --pmult;
  118.     }
  119.     return (numer * invertMod(denom, ppow)) % ppow;
  120. }
  121.        
  122. int factorial_prime_power(int f, int p) {
  123.     // When p is prime, returns k such that p^k divides f!
  124.     int k = 0;
  125.     while (f > p) {
  126.         f = f / p;
  127.         k = k + f;
  128.     }
  129.     return k;
  130. }
  131. int binomial_prime_power(int n, int r, int p) {
  132.     // when p is prime, returns k such that p^k divides nCr
  133.     return factorial_prime_power(n,p) - factorial_prime_power(r,p) - factorial_prime_power(n-r,p);
  134. }
  135. int powmod(int p, int k, int m) {
  136.     // quickly calculates p^k mod m
  137.     int res = 1;
  138.     int q = p;
  139.     int j = k;
  140.     while (j > 0) {
  141.         // invariant:  p^k is congruent to res * q^j
  142.         if (j is odd) {
  143.             res = (res * q) % m;
  144.             j = (j-1)/2;
  145.         } else {
  146.             j = j / 2;
  147.         }
  148.         q = (q * q) % m;
  149.     }
  150.     return res;
  151. }
  152. int big_binomial(int n, int r, int m) {
  153.     if (n < r or r < 0) {
  154.         return 0;
  155.     }
  156.     int res = 1;
  157.     for(p in all primes from 2 to n) {
  158.         k = binomial_prime_power(n,r,p);
  159.         res = (res * powmod(p,k,m)) % m;
  160.     }
  161.     return res;
  162. }
  163.        
  164. int factorial(int x) {
  165.         int i;
  166.         for(i=1; i<x; i++)
  167.         {
  168.             x *= i;
  169.         }
  170.         return x;
  171.     }
  172.  
  173.     int main()
  174.     {
  175.       //use like this
  176.       factorial(n)/( factorial(r) - factorial(n-r) );
  177.     //do your stuffs
  178.     }
clone this paste RAW Paste Data