Jump to content

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


Recommended Posts

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;
}
Link to comment
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
 Share

×
×
  • Create New...

Important Information

By using this site you accept our Terms of Use and Privacy Policy. We have placed cookies on your device to help make this website better. You can adjust your cookie settings, otherwise we'll assume you're okay to continue.