[C] Operações aritméticas com inteiros grandes

[orium]

Implementação, com API parecida com a gmp (mas com a performance monstruosamente inferior, claro), de operações aritméticas com inteiros de tamanho arbitrário. Os algoritmos implementados são os aprendidos na escola primária. Suporta inteiros negativos (mas não como exponente no int_pow()).

Este código depende de xalloc.h, que ja' ca' foi postado.

bigints.h:

/*
*  GPL
*
*  Written by Diogo Sousa aka orium
*  Copyright © 2008 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 <stdint.h>
#include <stdbool.h>

/* For best performace double_digit should be the size of one word
* and digit should be halve word
*/
typedef uint32_t double_digit;
typedef int32_t double_digit_signed;
typedef uint16_t digit;
typedef int32_t carry_t;

/* RADIX must be even
* STRING_RADIX must not exceed RADIX
* STRING_RADIX must not exceed 16
* INT_END must be greater (or equal) than RADIX
* RADIX must be greater than 16
*/
#define RADIX ((int32_t)((digit)~0)-1)
#define INT_END ((digit)~0)

#define STRING_RADIX 10

#if STRING_RADIX > 16
#error STRING_RADIX must not exceed 16
#endif

struct integer_str
{
digit *n; /* least significant digit first */
size_t size;

bool negative;

char *string; /* if non-NULL it is not dirty */
};

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 char *int_value(integer);

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

bigints.c:

/*
*  GPL
*
*  Written by Diogo Sousa aka orium
*  Copyright © 2008 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/>.
*/

#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include "xalloc.h"
#include "bigints.h"

#define STRING_INIT_SIZE 128
#define STRING_EXPAND_SIZE 512

#define max(q,p) (((q) > (p)) ? (q) : (p))
#define is_zero(i) ((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 mod(x,m) (((x) >= 0) ? (x)%(m) : (m)-((abs(x))%(m)))
#define string_put(str,size,i,ch) ({					\
		if ((i) >= *(size))				\
			(str)=xrealloc(str,			\
					   *(size)=(i)+STRING_EXPAND_SIZE);	\
		str[i]=(ch);					\
	})

static const char *char_radix="0123456789abcdef";

static digit one_digits[]={1,INT_END};
static digit two_digits[]={2,INT_END};
static digit ten_digits[]={10,INT_END};

static struct integer_str one={
.n=one_digits,
.size=1,
.negative=false,
.string=NULL
};

static struct integer_str two={
.n=two_digits,
.size=1,
.negative=false,
.string=NULL
};

static struct integer_str ten={
.n=ten_digits,
.size=1,
.negative=false,
.string=NULL
};

static void expand(integer, size_t);
static void move(integer, integer);
static void normalize(integer);
static void halve(integer);
static void divide(integer, integer, integer, integer);
static int digit_value(char);

static void
expand(integer i, size_t size)
{
assert(size);

i->n=xrealloc(i->n,sizeof(*i->n)*(size+1));
}

/* This make sure:
*   There are no leading zeros
*   Zero is not negative
*   Memory are not being wasted
*/
static void
normalize(integer o)
{
int i;

for (i=o->size-1; i && !o->n[i]; i--)
	;

o->n[i+1]=INT_END;
o->size=i+1;

expand(o,o->size);

if (is_zero(o))
	o->negative=false;
}

static void
move(integer r, integer value)
{
free(r->n);
free(r->string);

mark_dirty(value);
*r=*value;
value->n=NULL;
}

static 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 int
digit_value(char ch)
{
if (ch >= '0' && ch <= '9')
	return ch-'0';

assert(0);

return -1;
}

integer
int_set(integer r, integer value)
{
int i;

if (r == value)
	return r;

expand(r,value->size);

for (i=0; value->n[i] != INT_END; i++)
	r->n[i]=value->n[i];

r->n[i]=INT_END;

r->size=value->size;
r->negative=value->negative;

mark_dirty(r);

return r;
}

integer
int_set_int(integer r, int value)
{
int i;

if (value < 0)
{
	r->negative=true;
	value=abs(value);
}

expand(r,((int)ceil(log(value)/log(RADIX)))+4);

for (i=0; value; value/=RADIX, i++)
	r->n[i]=value%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)
{
integer power;
integer d;
int i;

int_init(&power);
int_init(&d);

int_set_int(r,0);
int_set_int(power,1);

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,&ten);
}

if (str[0] == '-')
	r->negative=true;

mark_dirty(r);
normalize(r);

int_free(power);
int_free(d);

return r;
}

char *
int_value(integer v)
{
integer k;
integer m;
integer str_radix;
int i;
int str_size;

if (!dirty(v))
	return v->string;

int_init(&k);
int_init(&m);
int_init(&str_radix);

int_set(k,v);
int_set_int(str_radix,STRING_RADIX);

v->string=xmalloc(str_size=STRING_INIT_SIZE);

for (i=0; !is_zero(k); i++)
{
	divide(k,m,k,str_radix);

	assert(m->size == 1);
	assert(m->n[0] < 16);

	string_put(v->string,&str_size,i,char_radix[m->n[0]]);
}

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');

v->string=xrealloc(v->string,i+1);

str_reverse(v->string);

int_free(k);
int_free(m);
int_free(str_radix);

return v->string;
}

integer
int_init(integer *i)
{
*i=xmalloc(sizeof(**i));

(*i)->n=xmalloc(sizeof(*(*i)->n)*2);
(*i)->size=1;
(*i)->negative=false;
(*i)->string=NULL;

(*i)->n[0]=0;
(*i)->n[1]=INT_END;

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)
{
int i;

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);

for (i=a->size-1; i >= 0; i--)
	if (a->n[i] != b->n[i])
		return (a->n[i]-b->n[i] > 0) ? 1 : -1;

return 0;
}

integer
int_add(integer res, integer a, integer b)
{
integer r;
int i;
double_digit v;
carry_t carry;

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;
}

int_init(&r);

expand(r,max(a->size,b->size)+1);

carry=0;

for (i=0; i < max(a->size,b->size) || carry; i++)
{
	v=carry;

	if (i < a->size)
		v+=a->n[i];

	if (i < b->size)
		v+=b->n[i];

	r->n[i]=v%RADIX;
	carry=v >= RADIX;
}

r->n[i]=INT_END;
r->size=i;

if (a->negative && b->negative)
	r->negative=true;

normalize(r);
move(res,r);
int_free(r);

return res;
}

integer
int_sub(integer res, integer a, integer b)
{
integer r;	
int i;
double_digit_signed v;
carry_t carry;

	/* If a and b are both negative this will work too */
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);

int_init(&r);

expand(r,a->size);

carry=0;

assert(a->n[a->size] == INT_END);
assert(b->n[b->size] == INT_END);

for (i=0; i < a->size || carry; i++)
{
	v=(i < a->size) ? a->n[i] : 0;

	if (i < b->size)
		v-=b->n[i];

	v+=carry;

	r->n[i]=mod(v,RADIX);

	carry=(v < 0) ? -1 : 0;
}

r->n[i]=INT_END;
r->size=i;

normalize(r);
move(res,r);
int_free(r);

return res;
}

integer
int_mul(integer res, integer a, integer b)
{
integer r;
int i;
int j;
carry_t carry;
double_digit v;
int_init(&r);

expand(r,a->size+b->size+1);

for (i=0; i < a->size+b->size+1; i++)
	r->n[i]=0;

for (i=0; b->n[i] != INT_END; i++)
{
	carry=0;

	for (j=0; a->n[j] != INT_END; j++)
	{
		v=a->n[j]*b->n[i]+r->n[i+j]+carry;
		r->n[i+j]=v%RADIX;
		carry=v/RADIX;
	}

	r->n[i+j]=carry;
}

r->n[i+j]=INT_END;
r->size=i+j;

r->negative=a->negative^b->negative;

normalize(r);
move(res,r);
int_free(r);

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 void
halve(integer r)
{
int i;
double_digit v;
carry_t carry;

carry=0;

for (i=r->size-1; i >= 0; i--)
{
	v=r->n[i]+carry*RADIX;
	r->n[i]=v>>1;
	carry=v&1;
}

mark_dirty(r);
normalize(r);
}

/* This division algorithms was token from Structured Programming,
* by Dahl, Dijkstra and Hoare.
*/
static 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)
	int_mul(t,t,&two);

while (int_cmp(t,b))
{
	halve(t);
	int_mul(q,q,&two);

	if (int_cmp(t,r) <= 0)
	{
		int_sub(r,r,t);
		int_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);

halve(halve_p);

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;
}

Um pequeno exemplo:

#include <stdio.h>
#include "bigints.h"

int main(void)
{
integer i1;
integer i2;

int_init(&i1);
int_init(&i2);

int_set_int(i1,100);
int_set_int(i2,4);

int_mul(i1,i1,i2);

printf("%s\n",int_value(i1));

int_free(i1);
int_free(i2);

return 0;
}