Operações aritméticas com Inteiros Grandes
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
, 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;
}