## pefactorial.c

``` 1 /* The authors of this work have released all rights to it and placed it
2 in the public domain under the Creative Commons CC0 1.0 waiver
3 (http://creativecommons.org/publicdomain/zero/1.0/).
4
5 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
6 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
7 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
8 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
9 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12
14 */
15
16 #include <math.h>  /* for log() */
17 #include <stdlib.h>
18 #include <stdio.h>
19 #include "integer.h"
20
21 /* Return the power of the prime number p in the
22    factorization of n! */
23 int multiplicity(int n, int p) {
24     int q = n, m = 0;
25     if (p > n) return 0;
26     if (p > n/2) return 1;
27     while (q >= p) {
28         q /= p;
29         m += q;
30     }
31     return m;
32 }
33
34 unsigned char* prime_table(int n) {
35     int i, j;
36     unsigned char* sieve = (unsigned char*)calloc(n+1, sizeof(unsigned char));
37     sieve[0] = sieve[1] = 1;
38     for (i=2; i*i <= n; i++)
39         if (sieve[i] == 0)
40             for (j=i*i; j <= n; j+=i)
41                 sieve[j] = 1;
42     return sieve;
43 }
44 integer factorial(int n) {
45     unsigned char* primes = prime_table(n);
46     int num_primes;
47     int* p_list = malloc(sizeof(int)*n);
48     int* exp_list = malloc(sizeof(int)*n);
49     int p, i, bit;
50     integer result =
51       create_integer((int)ceil((n*log((double)n)/log(2) + 1)/COMPONENT_BITS));
52     integer temp     = create_integer(result.num_components*2);
53     set_zero_integer(result);
54     result.c[0] = 1;
55
56     num_primes = 0;
57     for(p = 2; p <= n; p++) {
58         if (primes[p] == 1) continue;
59         p_list[num_primes] = p;
60         exp_list[num_primes] = multiplicity(n,p);
61         num_primes++;
62     }
63
64     for(bit=sizeof(exp_list[0])*CHAR_BIT - 1; bit>=0; bit--) {
65         for(i=0; i<num_primes; i++) {
66             if ((exp_list[i] & (1 << bit)) != 0) {
67                 break;
68             }
69         }
70         if (i < num_primes) {
71             break;
72         }
73     }
74
75     for( ; bit>=0; bit--) {
76         multiply_integer(result, result, temp);
77         copy_integer(temp, result);
78         for(i=0; i<num_primes; i++) {
79             if ((exp_list[i] & (1 << bit)) != 0) {
80                 multiply_small_integer(result, p_list[i], result);
81             }
82         }
83     }
84
85     free_integer(temp);
86     free(primes);
87     return result;
88 }
89
90 int main(int argc, char* argv[]) {
91     int n = atoi(argv[1]);
92     puts(integer_to_string(factorial(n)));
93     return 0;
94 }
```

hijacker
hijacker
hijacker
hijacker

## integer.c

```  1 /* The authors of this work have released all rights to it and placed it
2 in the public domain under the Creative Commons CC0 1.0 waiver
3 (http://creativecommons.org/publicdomain/zero/1.0/).
4
5 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
6 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
7 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
8 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
9 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12
14 */
15
16 #include <stdlib.h>   /* for malloc() */
17
18 #include <string.h>   /* for memset(), memmove() */
19 #include <math.h>  /* for ceil */
20 #include "integer.h"
21
22 integer create_integer(int components) {
23     integer result;
24     result.num_components = components;
25     result.c = (component_t*)malloc(sizeof(component_t)*components);
26     return result;
27 }
28
29
30 void free_integer(integer i) {
31     free(i.c);
32 }
33
34
35 void set_zero_integer(integer i) {
36     memset(i.c, 0, sizeof(component_t)*i.num_components);
37 }
38
39
40 int is_zero_integer(integer x) {
41     int i;
42     for(i=0; i < x.num_components; i++) {
43         if (x.c[i] != 0) return 0;
44     }
45     return 1;
46 }
47
48 void copy_integer(integer source, integer target) {
49     memmove(target.c, source.c,
50             sizeof(component_t)*MIN(source.num_components, target.num_components));
51     if (target.num_components > source.num_components) {
52         memset(target.c + source.num_components, 0,
53                sizeof(component_t)*(target.num_components - source.num_components));
54     }
55 }
56
57 void add_integer(integer left, integer right, integer result) {
58     double_component_t carry = 0;
59     int i;
60     for(i=0; i<left.num_components || i<right.num_components || carry != 0; i++) {
61         double_component_t partial_sum = carry;
62         carry = 0;
63         if (i < left.num_components)  partial_sum += left.c[i];
64         if (i < right.num_components) partial_sum += right.c[i];
65         if (partial_sum > MAX_COMPONENT) {
66             partial_sum &= MAX_COMPONENT;
67             carry = 1;
68         }
69         result.c[i] = (component_t)partial_sum;
70     }
71     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
72 }
73
74 void subtract_integer(integer left, integer right, integer result) {
75     int borrow = 0;
76     int i;
77     for(i=0; i<left.num_components; i++) {
78         double_component_t lhs = left.c[i];
79         double_component_t rhs = (i < right.num_components) ? right.c[i] : 0;
80         if (borrow) {
81 	    if (lhs <= rhs) {
82 	        /* leave borrow set to 1 */
83 	        lhs += (MAX_COMPONENT + 1) - 1;
84 	    } else {
85 	        borrow = 0;
86 	        lhs--;
87 	    }
88 	}
89         if (lhs < rhs) {
90             borrow = 1;
91             lhs += MAX_COMPONENT + 1;
92         }
93         result.c[i] = lhs - rhs;
94     }
95     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
96 }
97
98 void multiply_small_integer(integer left, component_t right, integer result) {
99     double_component_t carry = 0;
100     int i;
101     for(i=0; i<left.num_components || carry != 0; i++) {
102         double_component_t partial_sum = carry;
103         carry = 0;
104         if (i < left.num_components)  partial_sum += left.c[i]*right;
105         carry = partial_sum >> COMPONENT_BITS;
106         result.c[i] = (component_t)(partial_sum & MAX_COMPONENT);
107     }
108     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
109 }
110
111 void multiply_integer(integer left, integer right, integer result) {
112     int i, lidx, ridx;
113     double_component_t carry = 0;
114     int max_size_no_carry;
115     int left_max_component  = left.num_components - 1;
116     int right_max_component = right.num_components - 1;
117     while(left.c[left_max_component]   == 0) left_max_component--;
118     while(right.c[right_max_component] == 0) right_max_component--;
119     max_size_no_carry = left_max_component + right_max_component;
120     for(i=0; i <= max_size_no_carry || carry != 0; i++) {
121         double_component_t partial_sum = carry;
122         carry = 0;
123         lidx = MIN(i, left_max_component);
124         ridx = i - lidx;
125         while(lidx >= 0 && ridx <= right_max_component) {
126             partial_sum += ((double_component_t)left.c[lidx])*right.c[ridx];
127             carry += partial_sum >> COMPONENT_BITS;
128             partial_sum &= MAX_COMPONENT;
129             lidx--; ridx++;
130         }
131         result.c[i] = partial_sum;
132     }
133     for ( ; i < result.num_components; i++) { result.c[i] = 0; }
134 }
135
136 int compare_integers(integer left, integer right) {
137     int i = MAX(left.num_components - 1, right.num_components - 1);
138     for ( ; i >= 0; i--) {
139         component_t left_comp =
140             (i < left.num_components) ? left.c[i] : 0;
141         component_t right_comp =
142             (i < right.num_components) ? right.c[i] : 0;
143         if (left_comp < right_comp)
144             return -1;
145         else if (left_comp > right_comp)
146             return  1;
147     }
148     return 0;
149 }
150
151 void shift_left_one_integer(integer arg) {
152     int i;
153     arg.c[arg.num_components - 1] <<= 1;
154     for (i = arg.num_components - 2; i >= 0; i--) {
155         arg.c[i + 1] |= arg.c[i] >> (COMPONENT_BITS - 1);
156         arg.c[i] <<= 1;
157     }
158 }
159
160 void shift_right_one_integer(integer arg) {
161     int i;
162     arg.c[0] >>= 1;
163     for (i = 1; i < arg.num_components; i++) {
164         arg.c[i - 1] |= (arg.c[i] & 1) << (COMPONENT_BITS - 1);
165         arg.c[i] >>= 1;
166     }
167 }
168
169 component_t mod_small_integer(integer left, component_t right) {
170     double_component_t mod_two_power = 1;
171     double_component_t result = 0;
172     int i, bit;
173     for(i=0; i<left.num_components; i++) {
174         for(bit=0; bit<COMPONENT_BITS; bit++) {
175             if ((left.c[i] & (1 << bit)) != 0) {
176                 result += mod_two_power;
177                 if (result >= right) {
178                     result -= right;
179                 }
180             }
181             mod_two_power <<= 1;
182             if (mod_two_power >= right) {
183                 mod_two_power -= right;
184             }
185         }
186     }
187     return (component_t)result;
188 }
189
190 void mod_integer(integer left, integer right, integer result) {
191     integer mod_two_power = create_integer(right.num_components + 1);
192     int i, bit;
193     set_zero_integer(result);
194     set_zero_integer(mod_two_power);
195     mod_two_power.c[0] = 1;
196     for(i=0; i<left.num_components; i++) {
197         for(bit=0; bit<COMPONENT_BITS; bit++) {
198             if ((left.c[i] & (1 << bit)) != 0) {
200                 if (compare_integers(result, right) >= 0) {
201                     subtract_integer(result, right, result);
202                 }
203             }
204             shift_left_one_integer(mod_two_power);
205             if (compare_integers(mod_two_power, right) >= 0) {
206                 subtract_integer(mod_two_power, right, mod_two_power);
207             }
208         }
209     }
210     free_integer(mod_two_power);
211 }
212
213 void divide_small_integer(integer left, component_t right, integer result) {
214     double_component_t dividend = 0;
215     int i;
216     for (i = left.num_components - 1; i >= 0; i--) {
217         dividend |= left.c[i];
218         result.c[i] = dividend/right;
219         dividend = (dividend % right) << COMPONENT_BITS;
220     }
221 }
222
223 integer string_to_integer(char* s) {
224     integer result = create_integer((int)ceil(LOG_2_10*strlen(s)/COMPONENT_BITS));
225     integer digit = create_integer(1);
226     int i;
227     for (i = 0; s[i] != '\0'; i++) {
228         multiply_small_integer(result, 10, result);
229         digit.c[0] = s[i] - '0';
231     }
232     free_integer(digit);
233     return result;
234 }
235
236
237 char* integer_to_string(integer x) {
238     int i, result_len;
239     char* result =
240         (char*)malloc((int)ceil(COMPONENT_BITS*x.num_components/LOG_2_10) + 2);
241     integer ten = create_integer(1);
242     ten.c[0] = 10;
243
244     if (is_zero_integer(x)) {
245         strcpy(result, "0");
246     } else {
247         for (i = 0; !is_zero_integer(x); i++) {
248             result[i] = (char)mod_small_integer(x, 10) + '0';
249             divide_small_integer(x, 10, x);
250         }
251         result[i] = '\0';
252     }
253
254     result_len = strlen(result);
255     for(i=0; i < result_len/2; i++) {
256         char temp = result[i];
257         result[i] = result[result_len - i - 1];
258         result[result_len - i - 1] = temp;
259     }
260
261     free_integer(ten);
262     return result;
263 }
264
```

hijacker
hijacker
hijacker
hijacker

## factorial.c

``` 1 /* The authors of this work have released all rights to it and placed it
2 in the public domain under the Creative Commons CC0 1.0 waiver
3 (http://creativecommons.org/publicdomain/zero/1.0/).
4
5 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
6 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
7 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
8 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
9 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12
14 */
15
16 #include <math.h>  /* for log() */
17 #include <stdlib.h>
18 #include <stdio.h>
19 #include "integer.h"
20
21 /* Return the power of the prime number p in the
22    factorization of n! */
23 int multiplicity(int n, int p) {
24     int q = n, m = 0;
25     if (p > n) return 0;
26     if (p > n/2) return 1;
27     while (q >= p) {
28         q /= p;
29         m += q;
30     }
31     return m;
32 }
33
34 unsigned char* prime_table(int n) {
35     int i, j;
36     unsigned char* sieve = (unsigned char*)calloc(n+1, sizeof(unsigned char));
37     sieve[0] = sieve[1] = 1;
38     for (i=2; i*i <= n; i++)
39         if (sieve[i] == 0)
40             for (j=i*i; j <= n; j+=i)
41                 sieve[j] = 1;
42     return sieve;
43 }
44 void exponent(unsigned int base, unsigned int power, integer result) {
45     int bit;
46     integer temp = create_integer(result.num_components*2);
47     set_zero_integer(result);
48     result.c[0] = 1;
49
50     bit=sizeof(power)*CHAR_BIT - 1;
51     while ((power & (1 << bit)) == 0) bit--;
52     for( ; bit>=0; bit--) {
53         multiply_integer(result, result, temp);
54         copy_integer(temp, result);
55         if ((power & (1 << bit)) != 0) {
56             multiply_small_integer(result, base, result);
57         }
58     }
59
60     free_integer(temp);
61 }
62 integer factorial(int n) {
63     unsigned char* primes = prime_table(n);
64     int p;
65     integer result =
66       create_integer((int)ceil((n*log((double)n)/log(2) + 1)/COMPONENT_BITS));
67     integer p_raised = create_integer(result.num_components);
68     integer temp     = create_integer(result.num_components*2);
69     set_zero_integer(result);
70     result.c[0] = 1;
71
72     for(p = 2; p <= n; p++) {
73         if (primes[p] == 1) continue;
74         exponent(p, multiplicity(n,p), p_raised);
75         multiply_integer(result, p_raised, temp);
76         copy_integer(temp, result);
77     }
78
79     free_integer(temp);
80     free_integer(p_raised);
81     free(primes);
82     return result;
83 }
84
85
86 int main(int argc, char* argv[]) {
87     int n = atoi(argv[1]);
88     puts(integer_to_string(factorial(n)));
89     return 0;
90 }
```

hijacker
hijacker
hijacker
hijacker

## integer.h

``` 1 /* The authors of this work have released all rights to it and placed it
2 in the public domain under the Creative Commons CC0 1.0 waiver
3 (http://creativecommons.org/publicdomain/zero/1.0/).
4
5 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
6 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
7 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
8 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
9 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
10 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
11 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
12
14 */
15
16 #ifndef _INTEGER_H_
17 #define _INTEGER_H_
18
19 #include <limits.h>  /* for CHAR_BIT */
20 typedef unsigned short component_t;
21 typedef unsigned long double_component_t;
22
23 #define MAX_COMPONENT  ((component_t)(-1))
24 #define COMPONENT_BITS  (sizeof(component_t)*CHAR_BIT)
25
26 #define LOG_2_10    3.3219280948873623478703194294894
27
28 #define MIN(x,y)  ((x)<(y) ? (x) : (y))
29 #define MAX(x,y)  ((x)>(y) ? (x) : (y))
30
31 typedef struct {
32     component_t* c;    /* least-significant word first */
33     int num_components;
34 } integer;
35
36
37 integer create_integer(int components);
38 void free_integer(integer i);
39 void set_zero_integer(integer i);
40 void copy_integer(integer source, integer target);
41 void add_integer(integer left, integer right, integer result);
42 void subtract_integer(integer left, integer right, integer result);
43 void multiply_small_integer(integer left, component_t right, integer result);
44 void multiply_integer(integer left, integer right, integer result);
45 int compare_integers(integer left, integer right);
46 void shift_left_one_integer(integer arg);
47 void shift_right_one_integer(integer arg);
48 component_t mod_small_integer(integer left, component_t right);
49 void mod_integer(integer left, integer right, integer result);
50 void divide_small_integer(integer left, component_t right, integer result);
51 integer string_to_integer(char* s);
52 char* integer_to_string(integer x);
53 int is_zero_integer(integer x);
54
55 #endif /* #ifndef _INTEGER_H_ */
56
```

hijacker
hijacker
hijacker
hijacker