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