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