Implementação, com API parecida com a gmp (curiosamente a performance da adição e subtração é equiparável, a multiplicação e' bastante mais lenta), de operações aritméticas com inteiros de tamanho arbitrário. Os algoritmos implementados são os aprendidos na escola primária e o karatsuba para multiplicação, que só é chamado quando os números são de facto grandes.
Suporta inteiros negativos (mas não como exponente no int_pow()).
Este código depende de xalloc.h, que está também nesta wiki.
bigint.h:
/* * GPL * * Written by Diogo Sousa aka orium * Copyright (C) 2008, 2009, 2010, 2011 Diogo Sousa * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ #ifndef BIGINTS_H_ #define BIGINTS_H_ #include <stdbool.h> struct integer_str; typedef struct integer_str * integer; extern integer int_init(integer *); extern void int_free(integer); extern integer int_set(integer, integer); extern integer int_set_int(integer, int); extern integer int_set_string(integer, char *); extern integer int_set_string_radix(integer, char *, int); extern char *int_value(integer); extern char *int_value_radix(integer, int); extern int int_cmp(integer, integer); extern integer int_add(integer, integer, integer); extern integer int_sub(integer, integer, integer); extern integer int_mul(integer, integer, integer); extern integer int_div(integer, integer, integer); extern integer int_mod(integer, integer, integer); extern integer int_pow(integer, integer, integer); extern integer int_abs(integer); extern integer int_inc(integer); extern integer int_dec(integer); #endif
bigint.c:
/* * GPL * * Written by Diogo Sousa aka orium * Copyright (C) 2008, 2009, 2010, 2011 Diogo Sousa * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program. If not, see <http://www.gnu.org/licenses/>. */ /* Most of the algorithms used here were took from * The Art Of Computer Programming, Volume 2, Seminumerical Algorithms, * pages 265..280 (Multiple-Precision Arithmetic, The Classical Algorithms) */ #ifndef NDEBUG # define NDEBUG #endif #ifndef NDEBUG # include <stdio.h> #endif #include <stdlib.h> #include <stdbool.h> #include <string.h> #include <math.h> #include <assert.h> #include <stdint.h> #include "../xalloc/xalloc.h" #include "bigints.h" #define MUL_KARATSUBA_THRESHOLD (70*70) #define DEFAULT_DIGITS_ALLOCED 8 #define DEFAULT_DIGIT_EXPAND_EXTRA 8 #define STRING_INIT_SIZE 128 #define STRING_EXPAND_SIZE 128 #define NOMALIZE_SHRINK_DELTA 64 /* Radix should always be 10...0 */ #define RADIX 0x80000000 #define INT_END RADIX #define MOD_MASK 0x7FFFFFFF #define DIV_ROT 31 #define max(q,p) (((q) > (p)) ? (q) : (p)) #define min(q,p) (((q) < (p)) ? (q) : (p)) #define is_zero(i) (!(i)->size || ((i)->size == 1 && (i)->n[0] == 0)) #define is_one(i) ((i)->size == 1 && (i)->n[0] == 1) #define is_odd(i) ((i)->n[0]&1) #define dirty(i) ((i)->string == NULL) #define mark_dirty(i) ({ free((i)->string); (i)->string=NULL; }) #define tg_abs(x) (((x) >= 0) ? (x) : -(x)) /* type generic abs */ #define mod(x,m,mm) (((x)+(m))&(mm)) #define string_put(str,size,i,ch) ({ if ((i) >= *(size)) (str)=xrealloc(str, *(size)=(i)+STRING_EXPAND_SIZE); str[i]=(ch); }) typedef uint64_t double_digit; typedef int64_t double_digit_signed; typedef uint32_t digit; typedef int64_t carry_t; struct integer_str { digit *n; /* least significant digit first */ size_t size; size_t alloc_size; bool negative; char *string; /* if non-NULL it is not dirty */ short str_radix; }; static const char *char_radix="0123456789abcdefghijklmnopqrstuvwxyz"; static digit one_digits[]={1,INT_END}; static struct integer_str one={ .n=one_digits, .size=1, .alloc_size=0, .negative=false, .string=NULL }; static inline void expand(integer, size_t); static inline void move(integer, integer); static inline void normalize(integer); static inline void divide(integer, integer, integer, integer); static inline int digit_value(char); static inline void str_reverse(char *); static inline void double_int(integer); static inline void inc(integer); static inline int divide_int(integer, int); static inline integer int_mul_basecase(integer, integer, integer); static inline integer int_mul_karatsuba(integer, integer, integer); static inline integer int_init_dummy(integer *); #ifndef NDEBUG static void debug_print_integer(char *, integer); /* Prints the integer in a format python understands */ static void debug_print_integer(char *msg, integer v) { int pow=0; int i; printf("%s: ",msg); for (i=0; v->n[i] != INT_END; i++) printf("%llu*r**%d%s",(unsigned long long)v->n[i], pow++,(v->n[i+1] != INT_END) ? "+" : ""); putchar('n'); } #endif static inline void expand(integer i, size_t size) { assert(size); if (size >= i->alloc_size) { i->alloc_size=size+1+DEFAULT_DIGIT_EXPAND_EXTRA; i->n=xrealloc(i->n,sizeof(*i->n)*(i->alloc_size)); } } /* This make sure: * There are no leading zeros * Zero is not negative * Memory are not being wasted */ static inline void normalize(integer o) { int i; int new_size; digit *on; on=o->n; for (i=o->size-1; i > 0 && !on[i]; i--) ; on[i+1]=INT_END; new_size=i+1; if (o->size-new_size >= NOMALIZE_SHRINK_DELTA) expand(o,new_size); o->size=new_size; if (is_zero(o)) o->negative=false; } static inline void move(integer r, integer value) { mark_dirty(r); if (r == value) return; free(r->n); free(r->string); mark_dirty(value); *r=*value; value->n=NULL; } static inline void str_reverse(char *str) { int i; int j; char tmp; for (i=0, j=strlen(str)-1; i < j; i++, j--) { tmp=str[i]; str[i]=str[j]; str[j]=tmp; } } static inline int digit_value(char ch) { if (ch >= '0' && ch <= '9') return ch-'0'; return ch-'a'+10; } integer int_set(integer r, integer value) { if (r == value) return r; expand(r,value->size); memcpy(r->n,value->n,sizeof(*r->n)*(value->size+1)); r->size=value->size; r->negative=value->negative; mark_dirty(r); return r; } integer int_set_int(integer r, int value) { int i; double_digit ddvalue; if (value < 0) { r->negative=true; ddvalue=(double_digit)abs(value); } else ddvalue=(double_digit)value; expand(r,sizeof(ddvalue)); /* Wastes some memory but it is quickly normalized */ for (i=0; ddvalue; ddvalue/=RADIX, i++) r->n[i]=ddvalue%(double_digit)RADIX; if (!i) { r->n[0]=0; i=1; } r->n[i]=INT_END; r->size=i; normalize(r); mark_dirty(r); return r; } integer int_set_string(integer r, char *str) { return int_set_string_radix(r,str,10); } integer int_set_string_radix(integer r, char *str, int radix) { integer power; integer d; integer rad; int i; int_init(&power); int_init(&d); int_init(&rad); int_set_int(r,0); int_set_int(power,1); int_set_int(rad,radix); for (i=strlen(str)-1; i >= 0; i--) { if (!i && str[i] == '-') break; int_set_int(d,digit_value(str[i])); int_mul(d,d,power); int_add(r,r,d); int_mul(power,power,rad); } r->negative=str[0] == '-'; mark_dirty(r); normalize(r); int_free(power); int_free(d); int_free(rad); return r; } /* Returns the remainder */ static inline int divide_int(integer r, int d) { int i; double_digit v; carry_t carry=0; for (i=r->size-1; i >= 0; i--) { v=r->n[i]+carry*RADIX; r->n[i]=v/d; carry=v%d; } if (r->size-1 && !r->n[r->size-1]) { r->n[r->size-1]=INT_END; r->size--; } /* WARNING: This is not a public routine so the cleanup is not made: * mark_dirty(r); * normalize(r); */ return carry; } char * int_value(integer v) { return int_value_radix(v,10); } char * int_value_radix(integer v, int radix) { integer k; int i; int r; int str_size; if (v->str_radix != radix) mark_dirty(v); if (!dirty(v)) return v->string; int_init(&k); int_set(k,v); v->string=xmalloc(str_size=STRING_INIT_SIZE); for (i=0; !is_zero(k); i++) { r=divide_int(k,radix); assert(r < 16); string_put(v->string,&str_size,i,char_radix[r]); } if (v->negative) { string_put(v->string,&str_size,i,'-'); i++; } if (is_zero(v)) { string_put(v->string,&str_size,0,'0'); i=1; } string_put(v->string,&str_size,i,'0'); str_reverse(v->string); int_free(k); v->str_radix=radix; return v->string; } integer int_init(integer *i) { int_init_dummy(i); (*i)->n=xmalloc(sizeof(*(*i)->n)*DEFAULT_DIGITS_ALLOCED); (*i)->alloc_size=DEFAULT_DIGITS_ALLOCED; (*i)->size=1; (*i)->n[0]=0; (*i)->n[1]=INT_END; return *i; } static inline integer int_init_dummy(integer *i) { *i=xmalloc(sizeof(**i)); (*i)->n=NULL; (*i)->alloc_size=0; (*i)->size=0; (*i)->negative=false; (*i)->string=NULL; (*i)->str_radix=0; return *i; } void int_free(integer i) { free(i->n); free(i->string); free(i); } /* Returns > 0, 0 or < 0 if a > b, a == b or a < b */ int int_cmp(integer a, integer b) { digit *ad; digit *bd; digit *base; if (a == b) return 0; if (a->negative && !b->negative) return -1; if (!a->negative && b->negative) return 1; if (a->negative && b->negative) { int r; a->negative=false; b->negative=false; r=int_cmp(b,a); a->negative=true; b->negative=true; return r; } if (a->size > b->size) return 1; if (a->size < b->size) return -1; assert(a->size == b->size); base=a->n; for (ad=a->n+a->size-1, bd=b->n+b->size-1; ad >= base; ad--, bd--) if (*ad != *bd) return *ad < *bd ? -1 : 1; return 0; } integer int_add(integer res, integer a, integer b) { double_digit v; integer r; int i; int min_size; int max_size; digit *smaller; digit *larger; digit *rn; if (a->negative^b->negative) { if (!a->negative) { b->negative=false; int_sub(res,a,b); if (res != b) b->negative=true; } else { a->negative=false; int_sub(res,b,a); if (res != a) a->negative=true; } return res; } if (a != res && b != res) r=res; else int_init(&r); expand(r,max(a->size,b->size)+1); min_size=min(a->size,b->size); max_size=max(a->size,b->size); rn=r->n; if (min_size == a->size) { memcpy(rn+min_size,b->n+min_size, sizeof(*rn)*(b->size-min_size)); smaller=a->n; rn[b->size]=0; r->size=b->size; larger=b->n; } else { memcpy(rn+min_size,a->n+min_size, sizeof(*rn)*(a->size-min_size)); smaller=b->n; rn[a->size]=0; r->size=a->size; larger=a->n; } v=0; for (i=0; i < min_size; i++) { v+=(double_digit)larger[i]+(double_digit)smaller[i]; rn[i]=v&MOD_MASK; v=v >= RADIX; } for (; i < max_size && v; i++) { v+=rn[i]; rn[i]=v&MOD_MASK; v=v >= RADIX; } if (v) rn[i++]=v; r->size=max(i,r->size); rn[r->size]=INT_END; normalize(r); move(res,r); res->negative=a->negative && b->negative; if (a == res || b == res) int_free(r); return res; } integer int_sub(integer res, integer a, integer b) { integer r; int i; double_digit_signed v; digit *rn; digit *bn; digit *an; int bs; int as; /* If a and b are both negative this will work as well */ if (b->negative) { /* If b is negative: * * a-b=a+(-b) * */ b->negative=false; int_add(res,a,b); if (b != res) b->negative=true; return res; } assert(!b->negative); if (a->negative) { /* If a is negative: * * a-b=-((-a)+b) * */ a->negative=false; int_add(res,a,b); a->negative=true; res->negative=true; return res; } assert(!b->negative && !a->negative); if (int_cmp(b,a) > 0) { /* The result is negative * * a-b=-(b-a) * */ int_sub(res,b,a); res->negative=true; return res; } assert(!b->negative && !a->negative && int_cmp(a,b) >= 0); if (res != a && res != b) r=res; else int_init(&r); expand(r,a->size); assert(a->n[a->size] == INT_END); assert(b->n[b->size] == INT_END); /* we know a->size >= b->size */ rn=r->n; bn=b->n; an=a->n; bs=b->size; as=a->size; memcpy(rn+bs,an+bs,(as-bs)*sizeof(*r->n)); v=0; for (i=0; i < bs; i++) { v+=(double_digit_signed)an[i] -(double_digit_signed)bn[i]; rn[i]=mod(v,RADIX,MOD_MASK); v=-(v < 0); } for (; i < as; i++) { v+=rn[i]; rn[i]=mod(v,RADIX,MOD_MASK); v=-(v < 0); } if (v) rn[i]=mod(v,RADIX,MOD_MASK); r->n[a->size]=INT_END; r->size=i; r->negative=false; normalize(r); move(res,r); if (a == res || b == res) int_free(r); return res; } static inline integer int_mul_basecase(integer res, integer a, integer b) { int i; int j=0; integer r; double_digit v; digit *an; digit *bn; digit *rn; if (a->size > b->size) { integer tmp; tmp=a; a=b; b=tmp; } if (res != a && res != b) r=res; else int_init(&r); expand(r,a->size+b->size+1); an=a->n; bn=b->n; rn=r->n; memset(rn,0,sizeof(*rn)*(a->size+b->size+1)); for (i=0; bn[i] != INT_END; i++) { v=0; for (j=0; an[j] != INT_END; j++) { v+=(double_digit)an[j] *(double_digit)bn[i] +(double_digit)rn[i+j]; rn[i+j]=v&MOD_MASK; v>>=DIV_ROT; } rn[i+j]=v; } rn[i+j]=INT_END; r->size=i+j; r->negative=a->negative^b->negative; normalize(r); move(res,r); if (res == a || res == b) int_free(r); return res; } static inline integer int_mul_karatsuba(integer res, integer a, integer b) { int m=max(a->size,b->size)/2+max(a->size,b->size)%2; integer a0; integer a1; integer b0; integer b1; integer z1; integer z2; integer r; register double_digit v; int i; digit *rn; digit *z1n; digit *z2n; int z1s; int z2s; int mpi; int m2pi; int_init(&a0); int_init_dummy(&a1); int_init(&b0); int_init_dummy(&b1); int_init(&z1); int_init(&z2); if (res != a && res != b) r=res; else int_init(&r); /* set a0 and b0, the least significant digits of a and b */ a0->size=min(m,a->size); expand(a0,a0->size); memcpy(a0->n,a->n,a0->size*sizeof(*a0->n)); a0->n[a0->size]=INT_END; b0->size=min(m,b->size); expand(b0,b0->size); memcpy(b0->n,b->n,b0->size*sizeof(*b0->n)); b0->n[b0->size]=INT_END; /* set a1 and b1, the most significant digits of a and b */ a1->size=a->size-a0->size; a1->n=a->n+a0->size; b1->size=b->size-b0->size; b1->n=b->n+b0->size; normalize(a0); normalize(b0); int_mul(z2,a1,b1); int_mul(r,a0,b0); int_add(z1,a1,a0); int_add(a0,b1,b0); int_mul(b0,z1,a0); int_sub(a0,b0,z2); int_sub(z1,a0,r); expand(r,a->size+b->size+1); memset(r->n+r->size,0,sizeof(*r->n)*(a->size+b->size+1-r->size)); r->size=a->size+b->size+1; r->n[r->size]=INT_END; rn=r->n; z1n=z1->n; z2n=z2->n; z1s=z1->size; z2s=z2->size; v=0; for (i=0, mpi=m; i < z1s; i++, mpi++) { v+=(double_digit)rn[mpi]+(double_digit)z1n[i]; rn[mpi]=v&MOD_MASK; v=v >= RADIX; } for (; v; i++, mpi++) { v+=(double_digit)rn[mpi]; rn[mpi]=v&MOD_MASK; v=v >= RADIX; } for (i=0, m2pi=2*m; i < z2s; i++, m2pi++) { v+=(double_digit)rn[m2pi]+(double_digit)z2n[i]; rn[m2pi]=v&MOD_MASK; v=v >= RADIX; } for (; v; i++, m2pi++) { v+=(double_digit)rn[m2pi]; rn[m2pi]=v&MOD_MASK; v=v >= RADIX; } normalize(r); r->negative=a->negative^b->negative; int_free(a0); free(a1); int_free(b0); free(b1); int_free(z1); int_free(z2); move(res,r); if (res == a || res == b) int_free(r); return res; } integer int_mul(integer res, integer a, integer b) { if (a->size*b->size >= MUL_KARATSUBA_THRESHOLD) int_mul_karatsuba(res,a,b); else int_mul_basecase(res,a,b); return res; } integer int_abs(integer i) { i->negative=false; mark_dirty(i); return i; } integer int_inc(integer i) { return int_add(i,i,&one); } integer int_dec(integer i) { return int_sub(i,i,&one); } static inline void inc(integer r) { int i; double_digit v; carry_t carry=1; for (i=0; carry; i++) { if (i >= r->size) { expand(r,++r->size); r->n[i]=0; r->n[i+1]=INT_END; } v=(double_digit)r->n[i]+(double_digit)carry; r->n[i]=v&MOD_MASK; carry=v>>DIV_ROT; } /* WARNING: This is not a public routine so the cleanup is not made: * mark_dirty(r); * normalize(r); */ } static inline void double_int(integer r) { int i; double_digit v; carry_t carry=0; size_t size=r->size; for (i=0; i < size || carry; i++) { if (i >= r->size) { expand(r,++r->size); r->n[i]=0; } v=((double_digit)r->n[i]<<1)+(double_digit)carry; r->n[i]=v&MOD_MASK; carry=v>>DIV_ROT; } r->n[i]=INT_END; /* WARNING: This is not a public routine so the cleanup is not made: * mark_dirty(r); * normalize(r); */ } /* This division algorithms was token from Structured Programming, * by Dahl, Dijkstra and Hoare. */ static inline void divide(integer rq, integer rr, integer a, integer b) { integer t; integer r; integer q; bool an; bool bn; assert(!is_zero(b)); if (!int_cmp(a,b)) { if (rq != NULL) int_set_int(rq,1); if (rr != NULL) int_set_int(rr,0); return; } an=a->negative; bn=b->negative; a->negative=false; b->negative=false; int_init(&t); int_init(&r); int_init(&q); int_set(t,b); int_set(r,a); while (int_cmp(t,a) <= 0) { double_int(t); assert(!is_zero(t)); } while (int_cmp(t,b)) { divide_int(t,2); double_int(q); if (int_cmp(t,r) <= 0) { int_sub(r,r,t); inc(q); } } a->negative=an; b->negative=bn; if (rq != NULL) { move(rq,q); rq->negative=an^bn; } if (rr != NULL) move(rr,r); int_free(t); int_free(r); int_free(q); } integer int_div(integer r, integer a, integer b) { divide(r,NULL,a,b); return r; } integer int_mod(integer r, integer a, integer b) { divide(NULL,r,a,b); return r; } /* This computes b^p by the formula: b^p=b**(p/2)*b**(p/2)*b**(p%2) */ integer int_pow(integer r, integer b, integer p) { integer b_pow_halve_p; integer halve_p; integer res; assert(!p->negative); if (is_zero(p)) return int_set_int(r,1); if (is_one(p)) return int_set(r,b); int_init(&res); int_init(&b_pow_halve_p); int_init(&halve_p); int_set(halve_p,p); divide_int(halve_p,2); int_pow(b_pow_halve_p,b,halve_p); int_set(res,b_pow_halve_p); int_mul(res,b_pow_halve_p,b_pow_halve_p); if (is_odd(p)) int_mul(res,res,b); move(r,res); int_free(res); int_free(halve_p); int_free(b_pow_halve_p); r->negative=b->negative && is_odd(p); return r; }