• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# Untitled

a guest Jan 16th, 2013 82 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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.     }
RAW Paste Data
Top