Download code

Jump to: navigation, search

Back to Arbitrary-precision_integer_arithmetic_(C)

Download for Windows: zip

Download for UNIX: zip, tar.gz, tar.bz2

integer_sample.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 
13 Retrieved from: http://en.literateprograms.org/Arbitrary-precision_integer_arithmetic_(C)?oldid=19566
14 */
15 
16 #include <stdio.h>
17 #include <string.h>
18 #include <stdlib.h>  /* for atoi */
19 #include "integer.h"
20 
21 int main(int argc, char* argv[]) {
22     integer a1, a2, result;
23     a1.num_components = a2.num_components = 0; /* quelch compiler warnings */
24     if (argc > 2) a1 = string_to_integer(argv[2]);
25     if (argc > 3) a2 = string_to_integer(argv[3]);
26     if (strcmp(argv[1], "add") == 0) {
27         result = create_integer(a1.num_components + a2.num_components + 1);
28         add_integer(a1, a2, result);
29     } else if (strcmp(argv[1], "subtract") == 0) {
30         result = create_integer(a1.num_components + a2.num_components + 1);
31         subtract_integer(a1, a2, result);
32     } else if (strcmp(argv[1], "multiply") == 0) {
33         result = create_integer(a1.num_components + a2.num_components + 1);
34         multiply_integer(a1, a2, result);
35     } else if (strcmp(argv[1], "mod") == 0) {
36         result = create_integer(a2.num_components + 1);
37         mod_integer(a1, a2, result);
38     } else if (strcmp(argv[1], "factorial") == 0) {
39         int i, n = atoi(argv[2]);
40         result = create_integer(n*n + 1);
41         set_zero_integer(result);
42         result.c[0] = 1;
43         for(i=2; i<=n; i++) {
44             multiply_small_integer(result, i, result);
45         }
46     }
47     puts(integer_to_string(result));
48     return 0;
49 }


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 
 13 Retrieved from: http://en.literateprograms.org/Arbitrary-precision_integer_arithmetic_(C)?oldid=19566
 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) {
199                 add_integer(result, mod_two_power, result);
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     set_zero_integer(result);
226     integer digit = create_integer(1);
227     int i;
228     for (i = 0; s[i] != '\0'; i++) {
229         multiply_small_integer(result, 10, result);
230         digit.c[0] = s[i] - '0';
231         add_integer(result, digit, result);
232     }
233     free_integer(digit);
234     return result;
235 }
236 
237 
238 char* integer_to_string(integer x) {
239     int i, result_len;
240     char* result =
241         (char*)malloc((int)ceil(COMPONENT_BITS*x.num_components/LOG_2_10) + 2);
242     integer ten = create_integer(1);
243     ten.c[0] = 10;
244 
245     if (is_zero_integer(x)) {
246         strcpy(result, "0");
247     } else {
248         for (i = 0; !is_zero_integer(x); i++) {
249             result[i] = (char)mod_small_integer(x, 10) + '0';
250             divide_small_integer(x, 10, x);
251         }
252         result[i] = '\0';
253     }
254 
255     result_len = strlen(result);
256     for(i=0; i < result_len/2; i++) {
257         char temp = result[i];
258         result[i] = result[result_len - i - 1];
259         result[result_len - i - 1] = temp;
260     }
261 
262     free_integer(ten);
263     return result;
264 }
265 


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 
13 Retrieved from: http://en.literateprograms.org/Arbitrary-precision_integer_arithmetic_(C)?oldid=19566
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_ */


hijacker
hijacker
hijacker
hijacker