Ticket #601: jmp.c

#define SPOOF_GMP

#include <stdlib.h>

#include <stdio.h>
#define dprintf printf

#ifdef SPOOF_GMP
#define jmpn_cmp __gmpn_cmp
#define jmpz_divexact __gmpz_divexact
#define jmpn_gcd_1 __gmpn_gcd_1
#define jmpz_tdiv_qr __gmpz_tdiv_qr
#define jmp_set_memory_functions __gmp_set_memory_functions
#define jmpz_tdiv_q __gmpz_tdiv_q
#define jmpz_tdiv_r __gmpz_tdiv_r
#define jmpz_fdiv_qr __gmpz_fdiv_qr
#define jmpz_add __gmpz_add
#define jmpz_and __gmpz_and
#define jmpz_com __gmpz_com
#define jmpz_gcd __gmpz_gcd
#define jmpz_ior __gmpz_ior
#define jmpz_init __gmpz_init
#define jmpz_mul __gmpz_mul
#define jmpz_sub __gmpz_sub
#define jmpz_xor __gmpz_xor
#endif


typedef unsigned int limb;

typedef struct mpz {
        int alloc; // zero indicates mp zero
        int size;  // negative indicates mp negativity
        limb *limbs;
} mpz;

typedef void *(*alloc_func ) (size_t); 
typedef void *(*realloc_func ) (void *, size_t, size_t);
typedef void (*free_func ) (void *, size_t);

static void default_free( void *a, size_t s ) { free( a ); }
static void *default_realloc( void *a, size_t olds, size_t news ) { return realloc( a, news ); }

static alloc_func gAlloc = malloc;
static realloc_func gRealloc = default_realloc;
static free_func gFree = default_free;

void jmp_set_memory_functions( alloc_func my_alloc, realloc_func my_realloc, free_func my_free )
{
        dprintf( "mp: set_memory_functions( %p, %p, %p)\n", my_alloc, my_realloc, my_free );
        if( my_alloc ) gAlloc = my_alloc;
        if( my_realloc ) gRealloc = my_realloc;
        if( my_free ) gFree = my_free;
}

void jmpz_init( mpz *m )
{
        dprintf( "mp: mpz_init( %p )\n", m );
        m->alloc = 1;
        m->size = 0;
        m->limbs = (limb*) gAlloc( sizeof(limb) );
}

void jmpz_add( mpz *retv, mpz *a, mpz *b )
{
        dprintf( "mp: mpz_add( %p, %p, %p )\n", retv, a, b );
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        int sx, sy;
        int overflow;
        int n = zmax( sa * a->size , sb * b->size ) + 1;
        int m;
        limb *x, *y, *z, *last_nzero;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb) * retv->alloc, sizeof(limb) * n );
        overflow = 0;
        x  = sa*a->size < sb*b->size ? a->limbs : b->limbs;
        sx = sa*a->size < sb*b->size ? sa : sb;
        y  = sa*a->size < sb*b->size ? b->limbs : a->limbs;
        sy = sa*a->size < sb*b->size ? sb : sa;
        z = retv->limbs;
        m = zmin( sa * a->size, sb * b->size );
        last_nzero = retv->limbs;
        for( ; m; m-- ) {
                *z = (unsigned int)(sx * *x) + (unsigned int)(sy * *y) + overflow;
                overflow = *z < *x ? 1 : 0;
                if( sx * sy < 0 ) overflow = 1 - overflow;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
                n--;
        }
        for( ; n; n-- ) {
                *z = (unsigned int)(sy * *y) + overflow;
                overflow = *z < *y  ? 1 : 0;
                if( sy < 0 ) overflow = 1 - overflow;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
        }
        
        retv->size = (last_nzero - retv->limbs)/sizeof(limb) + 1;

}

void jmpz_sub( mpz *retv, mpz *a, mpz *b )
{
        dprintf( "mp: mpz_sub( %p, %p, %p )\n", retv, a, b );
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        int sx, sy;
        int overflow;
        int n = zmax( sa * a->size , sb * b->size ) + 1;
        int m;
        limb *x, *y, *z, *last_nzero;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb) * retv->alloc, sizeof(limb) * n );
        overflow = 0;
        x  = sa*a->size < sb*b->size ? a->limbs : b->limbs;
        sx = sa*a->size < sb*b->size ? sa : -sb;
        y  = sa*a->size < sb*b->size ? b->limbs : a->limbs;
        sy = sa*a->size < sb*b->size ? -sb : sa;
        z = retv->limbs;
        m = zmin( sa * a->size, sb * b->size );
        last_nzero = retv->limbs;
        for( ; m; m-- ) {
                *z = (unsigned int)(sx * *x) + (unsigned int)(sy * *y) + overflow;
                overflow = *z < *x ? 1 : 0;
                if( sx * sy < 0 ) overflow = 1 - overflow;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
                n--;
        }
        for( ; n; n-- ) {
                *z = (unsigned int)(sy * *y) + overflow;
                overflow = *z < *y  ? 1 : 0;
                if( sy < 0 ) overflow = 1 - overflow;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
        }
        
        retv->size = (last_nzero - retv->limbs)/sizeof(limb) + 1;

}

// binary xor
// TODO: What to do with negatives? 
void jmpz_xor( mpz *retv, mpz *a, mpz *b )
{
        dprintf( "mp: mpz_xor( %p, %p, %p )\n", retv, a, b );
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        int n = zmax( sa*a->size, sb*b->size );
        int m = zmin( sa*a->size, sb*b->size );
        limb *x, *y, *z, *last_nzero;
        int sx, sy;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb)*retv->alloc, sizeof(limb)*n );
        x  = sa*a->size < sb*b->size ? a->limbs : b->limbs;
        sx = sa*a->size < sb*b->size ? sa : sb;
        y  = sa*a->size < sb*b->size ? b->limbs : a->limbs;
        sy = sa*a->size < sb*b->size ? sb : sa;
        z = retv->limbs;
        last_nzero = retv->limbs;
        for( ;m; m-- ) {
                *z = *x ^ *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
                n--;
        }
        for( ;n; n-- ) {
                *z = 0 ^ *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
        }
        
        retv->size = (last_nzero - retv->limbs)/sizeof(limb) + 1;

}

// binary or
// TODO: What to do with negatives? 
void jmpz_ior( mpz *retv, mpz *a, mpz *b )
{
        dprintf( "mp: mpz_xor( %p, %p, %p )\n", retv, a, b );
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        int n = zmax( sa*a->size, sb*b->size );
        int m = zmin( sa*a->size, sb*b->size );
        limb *x, *y, *z, *last_nzero;
        int sx, sy;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb)*retv->alloc, sizeof(limb)*n );
        x  = sa*a->size < sb*b->size ? a->limbs : b->limbs;
        sx = sa*a->size < sb*b->size ? sa : sb;
        y  = sa*a->size < sb*b->size ? b->limbs : a->limbs;
        sy = sa*a->size < sb*b->size ? sb : sa;
        z = retv->limbs;
        last_nzero = retv->limbs;
        for( ;m; m-- ) {
                *z = *x | *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
                n--;
        }
        for( ;n; n-- ) {
                *z = 0 | *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
        }
        
        retv->size = (last_nzero - retv->limbs)/sizeof(limb) + 1;

}

// binary or
// TODO: What to do with negatives? 
void jmpz_and( mpz *retv, mpz *a, mpz *b )
{
        dprintf( "mp: mpz_xor( %p, %p, %p )\n", retv, a, b );
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        int n = zmax( sa*a->size, sb*b->size );
        int m = zmin( sa*a->size, sb*b->size );
        limb *x, *y, *z, *last_nzero;
        int sx, sy;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb)*retv->alloc, sizeof(limb)*n );
        x  = sa*a->size < sb*b->size ? a->limbs : b->limbs;
        sx = sa*a->size < sb*b->size ? sa : sb;
        y  = sa*a->size < sb*b->size ? b->limbs : a->limbs;
        sy = sa*a->size < sb*b->size ? sb : sa;
        z = retv->limbs;
        last_nzero = retv->limbs;
        for( ;m; m-- ) {
                *z = *x & *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
                n--;
        }
        for( ;n; n-- ) {
                *z = 0 & *y;
                if( *z != 0 ) { last_nzero = z; }
                z++;
                x++;
                y++;
        }
        
        retv->size = (last_nzero - retv->limbs)/sizeof(limb) + 1;

}

// binary complement
// TODO: What to do with negatives? 
void jmpz_com( mpz *retv, mpz *a )
{
        dprintf( "mp: mpz_com( %p, %p )\n", retv, a );
        int sa = a->size < 0 ? -1 : 1;
        int n = sa * a->size;
        int i;
        if( retv->alloc < n ) retv->limbs = gRealloc( retv->limbs, sizeof(limb)*retv->alloc, sizeof(limb)*n );
        for( i=0; i< n; i++ ) {
                retv->limbs[i] = ~ a->limbs[i];
        }
        retv->size = a->size;
}

/*
int jmpz_cmp( mpz *a, mpz *b )
{
        int sa = a->size < 0 ? -1 : 1;
        int sb = b->size < 0 ? -1 : 1;
        if( sa < sb ) return -1;
        if( sa > sb ) return 1;

        limb r = 0;
        int n = sa * a->size;
        int i;
        for( i=n-1 ; i>=0; i-- ) {
                r = a->limbs[i]  - b->limbs[i];
                if( r ) break;
        }
        return sa * r;
}
*/

int jmpn_cmp (limb *a, limb *b, size_t n )
{
        dprintf( "mp: mpn_cmp(%p,%p, %zu)\n", a, b, n );
        limb r = 0;
        int i;
        for( i=n-1 ; i>=0; i-- ) {
                r = a[i]  - b[i];
                if( r ) break;
        }
        return (int) r;
}


void jmpz_set( mpz *t, mpz *b )
{
        t->alloc = b->alloc;
        t->size = b->size;
        t->limbs = gAlloc( t->alloc );
        memcpy( t->limbs, b->limbs, t->alloc );
}

// let a = b/c  where c divides b, but this algorithm is ridiculously slow (it repeatedly adds one to a and subtracts c from b)
int jmpz_divexact( mpz *a, mpz *b, mpz *c )
{
        mpz one;
        mpz t;
        int sgn_b, sgn_c;
        sgn_b = b->size > 0 ? 1 : -1;
        sgn_c = c->size > 0 ? 1 : -1;

        jmpz_init( &one )
        one.size = 1;
        one.limbs[0] = 1;

        jmpz_set( &t, b );
        t.size *= sgn_b;

        // Set a = 0 
        a->limbs = gRealloc( a->limbs, a->alloc, 1 );
        a->size = 0;
        a->alloc = 1;

        if( c->size == 0 ) { 
                // TODO: division by zero!
                return;
        }
        c->size *= sgn_c;

        while( t->size>0 ) {
                mpz tmp;
                jmpz_set( &tmp, a );
                jmpz_add( a, &tmp, &one );
                jmpz_set( &tmp, &t );
                jmpz_sub( &t, &t, c );
        }

        c->size *= sgn_c;

}

/*
TODO:
        ___gmpn_gcd_1
        ___gmpz_tdiv_qr
        ___gmpz_tdiv_q
        ___gmpz_tdiv_r
        ___gmpz_fdiv_qr
        ___gmpz_gcd
        ___gmpz_mul
*/