Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/misc/volume40/fgmp/part01

Text File | 1993-11-01 | 51.1 KB | 2,119 lines

Newsgroups: comp.sources.misc
From: [e-mail] (Mark C. Henderson)
Subject: v40i074:  fgmp - Public domain multiple precision integer routines, Part01/01
Message-ID: <[e-mail]>
X-Md4-Signature: 87bfc6652ee9237438b9af0a463e74a9
Sender: [e-mail] (Kent Landfield)
Organization: Wimsey Information Services
Date: Mon, 1 Nov 1993 22:16:14 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Mark C. Henderson)
Posting-number: Volume 40, Issue 74
Archive-name: fgmp/part01
Environment: UNIX, MS-DOS

FGMP is a collection of routines for doing multiple precision integer
arithmetic with the same API as the GNU MP library (Version 1.3.2).
(In fact, only a subset of the functionality of GNU MP is implemented).

For instance, you can link the following trivial program with either
this code, or GNU MP and get the same results.
 ------------
 #include <stdio.h>
 #include "gmp.h"
 main()
 {
     MP_INT a; MP_INT b; MP_INT c;

     mpz_init_set_ui(&a,1); mpz_init_set_ui(&b,2); mpz_init(&c);
     mpz_add(&c,&a,&b);
     printf("\n%s\n", mpz_get_str(NULL,10,&c));
 }
 ------------

FGMP is considerably slower than GNU MP. The main advantage FGMP has over
GNU MP is that FGMP is in the public domain, while GNU MP is covered by
the GPL.

FGMP has been ported to the following platforms.:
Linux 0.99 (gcc 2.3, i486)
IBM RS6000/AIX 3.2 (IBM cc, gcc 2.4)
HP/UX 8.0 (HP cc, gcc 2.3)
Sun OS 4.1, Sun 3/4 (sun cc, gcc 2.2)
SGI Irix 4.0.5 (SGI cc, gcc 2.3)
DEC Alpha OSF/1 (DEC cc)
    (only lightly tested, 64 bit longs do make a difference,
     thanks to DEC for providing access on axposf.pa.dec.com).
     Define B64 for this platform.
MSDOS 286 C compiler (see credits above)

-------------cut here---------------
#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#        "End of shell archive."
# Contents:  notes gmp.h gmp.c
# Wrapped by mch@squid on Fri Oct 22 13:39:21 1993
PATH=/bin:/usr/bin:/usr/ucb:/usr/local/bin:/usr/lbin:$PATH ; export PATH
if test -f 'notes' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'notes'\"
else
echo shar: Extracting \"'notes'\" \(3666 characters\)
sed "s/^X//" >'notes' <<'END_OF_FILE'
XWELCOME TO FGMP.
X
XFGMP is a public domain implementation of a subset of the GNU gmp library
Xwith the same API.
X
XFor instance, you can link the following trivial program with either
Xthis code, or libgmp.a and get the same results.
X------------
X#include <stdio.h>
X#include "gmp.h"
Xmain()
X{
X    MP_INT a; MP_INT b; MP_INT c;
X
X    mpz_init_set_ui(&a,1); mpz_init_set_ui(&b,2); mpz_init(&c);
X    mpz_add(&c,&a,&b);
X    printf("\n%s\n", mpz_get_str(NULL,10,&c));
X}
X
X------------
X
XFGMP is really in the public domain. You can do whatever you want with
Xit.
X
XI wrote FGMP so that we would all have access to a (truly free)
Ximplementation of this subset of the API of GNU libgmp. I encourage
Xeveryone to distribute this as widely as possible.
X
XIf you need more documentation, I suggest you look at the file
Xgmp.texi which is included with the GNU gmp library.
X
XYou can send me bug reports, implementations of missing functions, flames
Xand rants by Email.
X
XAny submissions of new code to be integrated into fgmp must also be
Xplaced in the public domain (For the particularly dense, you can
Xrelease a new fgmp yourself under different licensing terms. This
Xis a condition for including a submission in a release of FGMP that
XI personally prepare).
X
XMark Henderson <[e-mail]>
X
X---
XThis is FGMP 1.0
X
XI hearby place this file and all of FGMP in the public domain.
X
XThanks to Paul Rouse <[e-mail]> for changes to get fgmp
Xto work on a 286 MSDOS compiler, the functions mpz_sqrt and
Xmpz_sqrtrem, plus other general bug fixes.
X
XThanks also to Erick Gallesio <[e-mail]> for a fix
Xto mpz_init_set_str
X
XDefine B64 if your "long" type is 64 bits. Otherwise we assume 32
Xbit longs. (The 64 bit version hasn't been tested enough)
X
X
X
XPlatforms:
X
XFGMP has been ported to the following platforms.:
X
XLinux 0.99 (gcc 2.3, i486)
XIBM RS6000/AIX 3.2 (IBM cc, gcc 2.4)
XHP/UX 8.0 (HP cc, gcc 2.3)
XSun OS 4.1, Sun 3/4 (sun cc, gcc 2.2)
XSGI Irix 4.0.5 (SGI cc, gcc 2.3)
XDEC Alpha OSF/1 (DEC cc)
X    (only lightly tested, 64 bit longs do make a difference,
X     thanks to DEC for providing access on axposf.pa.dec.com).
X     Define B64 for this platform.
XMSDOS 286 C compiler (see credits above)
X
X---
XSome differences between gmp and fgmp
X
X1. fgmp is considerably slower than gmp
X2. fgmp does not implement the following:
X    all mpq_*
X    internal mpn_* functions
X    mpz_perfect_square_p
X    mpz_inp_raw, mpz_out_raw
X    mp_set_memory_functions, mpz_out_str, mpz_inp_str
X3. fgmp implements the following in addition to the routines in GNU gmp.
X    int mpz_jacobi(MP_INT *a, MP_INT *b)
X    - finds the jacobi symbol (a/b)
X4. mpz_sizeinbase often overestimates the exact value
X
X5. To convert your gmp based program to fgmp (subject to the
Xabove)
X
X- recompile your source. Make sure to include the gmp.h file included
X  with fgmp rather than that included with gmp. (The point is to recompile
X  all files which include gmp.h)
X- link with gmp.o instead of libgmp.a
X
XHere's a complete sorted list of function implemented in fgmp:
X
X_mpz_realloc
Xmpz_abs
Xmpz_add
Xmpz_add_ui
Xmpz_and
Xmpz_clear
Xmpz_cmp
Xmpz_cmp_si
Xmpz_cmp_ui
Xmpz_div
Xmpz_div_2exp
Xmpz_div_ui
Xmpz_divmod
Xmpz_divmod_ui
Xmpz_fac_ui
Xmpz_gcd
Xmpz_gcdext
Xmpz_get_si
Xmpz_get_str
Xmpz_get_ui
Xmpz_init
Xmpz_init_set
Xmpz_init_set_si
Xmpz_init_set_str
Xmpz_init_set_ui
Xmpz_jacobi
Xmpz_mdiv
Xmpz_mdiv_ui
Xmpz_mdivmod
Xmpz_mdivmod_ui
Xmpz_mmod
Xmpz_mmod_ui
Xmpz_mod
Xmpz_mod_2exp
Xmpz_mod_ui
Xmpz_mul
Xmpz_mul_2exp
Xmpz_mul_ui
Xmpz_neg
Xmpz_or
Xmpz_pow_ui
Xmpz_powm
Xmpz_powm_ui
Xmpz_probab_prime_p
Xmpz_random
Xmpz_random2
Xmpz_set
Xmpz_set_si
Xmpz_set_str
Xmpz_set_ui
Xmpz_size
Xmpz_sizeinbase
Xmpz_sqrt
Xmpz_sqrtrem
Xmpz_sub
Xmpz_sub_ui
Xmpz_xor
END_OF_FILE
if test 3666 -ne `wc -c <'notes'`; then
    echo shar: \"'notes'\" unpacked with wrong size!
fi
# end of 'notes'
fi
if test -f 'gmp.h' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'gmp.h'\"
else
echo shar: Extracting \"'gmp.h'\" \(4201 characters\)
sed "s/^X//" >'gmp.h' <<'END_OF_FILE'
X/*
X * FREE GMP - a public domain implementation of a subset of the
X *           gmp library
X *
X * I hearby place the file in the public domain.
X *
X * Do whatever you want with this code. Change it. Sell it. Claim you
X *  wrote it.
X * Bugs, complaints, flames, rants: please send email to
X *    Mark Henderson <[e-mail]>
X * I'm already aware that fgmp is considerably slower than gmp
X *
X * CREDITS:
X *  Paul Rouse <[e-mail]> - generic bug fixes, mpz_sqrt and
X *    mpz_sqrtrem, and modifications to get fgmp to compile on a system
X *    with int and long of different sizes (specifically MSDOS,286 compiler)
X *  Also see the file "notes" included with the fgmp distribution, for
X *    more credits.
X *
X * VERSION 1.0
X */
X
X#include <stdio.h>
X#include <sys/types.h>
X
X/* for malloc and free */
X#include <stdlib.h>
X
X#ifndef NULL
X#define NULL ((void *)0)
X#endif
X
Xtypedef long mp_limb;
Xtypedef unsigned mp_size;
X
Xtypedef struct mp_int {
X    mp_limb *p;
X    short sn;
X    mp_size sz;
X} MP_INT;
X
X#ifdef __STDC__
X#define PROTO(x) x
X#else
X#define PROTO(x) ()
X#endif
X
Xvoid mpz_init PROTO((MP_INT *s));
Xvoid mpz_init_set PROTO((MP_INT *s, MP_INT *t));
Xvoid mpz_init_set_ui PROTO((MP_INT *s, unsigned long v));
Xvoid mpz_init_set_si PROTO((MP_INT *y, long v));
Xvoid mpz_clear PROTO((MP_INT *s));
Xvoid _mpz_realloc PROTO((MP_INT *x, mp_size size));
Xvoid mpz_set PROTO((MP_INT *y, MP_INT *x));
Xvoid mpz_set_ui PROTO((MP_INT *y, unsigned long v));
Xunsigned long mpz_get_ui PROTO((MP_INT *y));
Xlong mpz_get_si PROTO((MP_INT *y));
Xvoid mpz_set_si PROTO((MP_INT *y, long v));
Xint mpz_cmp PROTO((MP_INT *x, MP_INT *y));
Xvoid mpz_mul PROTO((MP_INT *ww,MP_INT *u, MP_INT *v));
Xvoid mpz_mul_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_div_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_mod_2exp PROTO((MP_INT *z, MP_INT *x, unsigned long e));
Xvoid mpz_add PROTO((MP_INT *zz,MP_INT *x,MP_INT *y));
Xvoid mpz_add_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mul_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_sub PROTO((MP_INT *z,MP_INT *x, MP_INT *y));
Xvoid mpz_sub_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_div PROTO((MP_INT *q, MP_INT *x, MP_INT *y));
Xvoid mpz_mdiv PROTO((MP_INT *q, MP_INT *x, MP_INT *y));
Xvoid mpz_mod PROTO((MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_divmod PROTO((MP_INT *q, MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mmod PROTO((MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mdivmod PROTO((MP_INT *q,MP_INT *r, MP_INT *x, MP_INT *y));
Xvoid mpz_mod_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mmod_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_div_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mdiv_ui PROTO((MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_divmod_ui PROTO((MP_INT *q,MP_INT *x,MP_INT *y, unsigned long n));
Xvoid mpz_mdivmod_ui PROTO((MP_INT *q,MP_INT *x,MP_INT *y, unsigned long n));
Xunsigned int mpz_sizeinbase PROTO((MP_INT *x, int base));
Xchar *mpz_get_str PROTO((char *s,  int base, MP_INT *x));
Xint mpz_set_str PROTO((MP_INT *x, char *s, int base));
Xint mpz_init_set_str PROTO((MP_INT *x, char *s, int base));
Xvoid mpz_random PROTO((MP_INT *x, mp_size size));
Xvoid mpz_random2 PROTO((MP_INT *x, mp_size size));
Xsize_t mpz_size PROTO((MP_INT *x));
Xvoid mpz_abs PROTO((MP_INT *, MP_INT *));
Xvoid mpz_neg PROTO((MP_INT *, MP_INT *));
Xvoid mpz_fac_ui PROTO((MP_INT *, unsigned long));
Xvoid mpz_gcd PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_gcdext PROTO((MP_INT *, MP_INT *, MP_INT *, MP_INT *, MP_INT *));
Xint mpz_jacobi PROTO((MP_INT *, MP_INT *));
Xint mpz_cmp_ui PROTO((MP_INT *, unsigned long));
Xint mpz_cmp_si PROTO((MP_INT *, long));
Xvoid mpz_and PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_or PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_xor PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_pow_ui PROTO((MP_INT *, MP_INT *, unsigned long));
Xvoid mpz_powm PROTO((MP_INT *, MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_powm_ui PROTO((MP_INT *, MP_INT *, unsigned long, MP_INT *));
Xint mpz_probab_prime_p  PROTO((MP_INT *, int));
Xvoid mpz_sqrtrem PROTO((MP_INT *, MP_INT *, MP_INT *));
Xvoid mpz_sqrt PROTO((MP_INT *, MP_INT *));
END_OF_FILE
if test 4201 -ne `wc -c <'gmp.h'`; then
    echo shar: \"'gmp.h'\" unpacked with wrong size!
fi
# end of 'gmp.h'
fi
if test -f 'gmp.c' -a "${1}" != "-c" ; then
  echo shar: Will not clobber existing file \"'gmp.c'\"
else
echo shar: Extracting \"'gmp.c'\" \(38872 characters\)
sed "s/^X//" >'gmp.c' <<'END_OF_FILE'
X/*
X * FREE GMP - a public domain implementation of a subset of the
X *           gmp library
X *
X * I hearby place the file in the public domain.
X *
X * Do whatever you want with this code. Change it. Sell it. Claim you
X *  wrote it.
X * Bugs, complaints, flames, rants: please send email to
X *    Mark Henderson <[e-mail]>
X * I'm already aware that fgmp is considerably slower than gmp
X *
X * CREDITS:
X *  Paul Rouse <[e-mail]> - generic bug fixes, mpz_sqrt and
X *    mpz_sqrtrem, and modifications to get fgmp to compile on a system
X *    with int and long of different sizes (specifically MSDOS,286 compiler)
X *  Also see the file "notes" included with the fgmp distribution, for
X *    more credits.
X *
X * VERSION 1.0
X */
X
X#include "gmp.h"
X#ifndef NULL
X#define NULL ((void *)0)
X#endif
X
X#define iabs(x) ((x>0) ? (x) : (-x))
X#define imax(x,y) ((x>y)?x:y)
X#define LONGBITS (sizeof(mp_limb) * 8)
X#define DIGITBITS (LONGBITS - 2)
X#define HALFDIGITBITS ((LONGBITS -2)/2)
X#ifndef init
X#define init(g) { g = (MP_INT *)malloc(sizeof(MP_INT));  mpz_init(g); }
X#endif
X#ifndef clear
X#define clear(g) { mpz_clear(g); free(g); }
X#endif
X#ifndef even
X#define even(a) (!((a)->p[0] & 1))
X#endif
X#ifndef odd
X#define odd(a) (((a)->p[0] & 1))
X#endif
X
X
X#ifndef B64
X/*
X * The values below are for 32 bit machines (i.e. machines with a
X *  32 bit long type)
X * You'll need to change them, if you're using something else
X * If DIGITBITS is odd, see the comment at the top of mpz_sqrtrem
X */
X#define LMAX 0x3fffffffL
X#define LC   0xc0000000L
X#define OVMASK 0x2
X#define CMASK (LMAX+1)
X#define FC ((double)CMASK)
X#define HLMAX 0x7fffL
X#define HCMASK (HLMAX + 1)
X#define HIGH(x) (((x) & 0x3fff8000L) >> 15)
X#define LOW(x)  ((x) & 0x7fffL)
X#else
X
X/* 64 bit long type */
X#define LMAX 0x3fffffffffffffffL
X#define LC 0xc000000000000000L
X#define OVMASK 0x2
X#define CMASK (LMAX+1)
X#define FC ((double)CMASK)
X#define HLMAX 0x7fffffffL
X#define HCMASK (HLMAX + 1)
X#define HIGH(x) (((x) & 0x3fffffff80000000L) >> 31)
X#define LOW(x) ((x) & 0x7fffffffL)
X
X#endif
X
X#define hd(x,i)  (((i)>=2*(x->sz))? 0:(((i)%2) ? HIGH((x)->p[(i)/2]) \
X    : LOW((x)->p[(i)/2])))
X#define dg(x,i) ((i < x->sz) ? (x->p)[i] : 0)
X
X#ifdef __GNUC__
X#define INLINE inline
X#else
X#define INLINE
X#endif
X
X#define lowdigit(x) (((x)->p)[0])
X
Xstruct is {
X    mp_limb v;
X    struct is *next;
X};
X
X
XINLINE static void push(i,sp)
Xmp_limb i;struct is **sp;
X{
X    struct is *tmp;
X    tmp = *sp;
X    *sp = (struct is *) malloc(sizeof(struct is));
X    (*sp)->v = i;
X    (*sp)->next=tmp;
X}
X
XINLINE static mp_limb pop(sp)
Xstruct is **sp;
X{
X    struct is *tmp;
X    mp_limb i;
X    if (!(*sp))
X        return (-1);
X    tmp = *sp;
X    *sp = (*sp)->next;
X    i = tmp->v;
X    tmp->v = 0;
X    free(tmp);
X    return i;
X}
X
Xvoid fatal(s)
Xchar *s;
X{
X    fprintf(stderr,"%s\n",s);
X    exit(123);
X}
X
Xvoid mpz_init(s)
XMP_INT *s;
X{
X    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X#ifdef DEBUG
X    printf("malloc: 8 bytes, %08lx\n", (long)(s->p));
X#endif
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    (s->p)[0] = 0;
X    (s->p)[1] = 0;
X    s->sn=0;
X    s->sz=2;
X}
X
Xvoid mpz_init_set(s,t)
XMP_INT *s,*t;
X{
X    int i;
X    s->p = (mp_limb *)malloc(sizeof(mp_limb) * t->sz);
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    for (i=0;i < t->sz ; i++)
X        (s->p)[i] = (t->p)[i];
X
X    s->sn = t->sn;
X    s->sz = t->sz;
X}
X
Xvoid mpz_init_set_ui(s,v)
XMP_INT *s;
Xunsigned long v;
X{
X    s->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X    if (!(s->p))
X        fatal("mpz_init: cannot allocate memory");
X    s->p[0] = (v & LMAX);
X    s->p[1] = (v & LC) >> DIGITBITS;
X    if (v)
X        s->sn = 1;
X    else
X        s->sn = 0;
X    s->sz = 2;
X}
X
Xvoid mpz_init_set_si(y,v)
XMP_INT *y;
Xlong v;
X{
X    y->p = (mp_limb *)malloc(sizeof(mp_limb)*2);
X    if (!(y->p))
X        fatal("mpz_init: cannot allocate memory");
X    if (v < 0) {
X        y->sn = -1;
X        y->p[0] = (-v) & LMAX;
X        y->p[1] = ((-v) & LC) >> DIGITBITS;
X    }
X    else if (v > 0) {
X        y->sn = 1;
X        y->p[0] = v & LMAX;
X        y->p[1] = (v & LC) >> DIGITBITS;
X    }
X    else {
X        y->sn=0;
X        y->p[0] = 0;
X        y->p[1] = 0;
X    }
X    y -> sz = 2;
X}
X
X
X
Xvoid mpz_clear(s)
XMP_INT *s;
X{
X#ifdef DEBUG
X    printf("free: %08lx\n", (long)(s->p));
X#endif
X    if (s->p)
X        free(s->p);
X    s->p=NULL;
X    s->sn=0;
X    s->sz=0;
X}
X
X/* number of leading zero bits in digit */
Xstatic int lzb(a)
Xmp_limb a;
X{
X    mp_limb i; int j;
X    j = 0;
X    for (i = ((mp_limb)1 << (DIGITBITS-1)); i && !(a&i) ; j++,i>>=1)
X        ;
X    return j;
X}
X
Xvoid _mpz_realloc(x,size)
XMP_INT *x; mp_size size;
X{
X    if (size > 1 && x->sz < size) {
X        int i;
X#ifdef DEBUG
X    printf("realloc %08lx to size = %ld ", (long)(x->p),(long)(size));
X#endif
X
X        if (x->p)
X            x->p=(mp_limb *)realloc(x->p,size * sizeof(mp_limb));
X        else
X            x->p=(mp_limb *)malloc(size * sizeof(mp_limb));
X#ifdef DEBUG
X    printf("becomes %08lx\n", (long)(x->p));
X#endif
X
X        if (!(x->p))
X            fatal("_mpz_realloc: cannot allocate memory");
X        for (i=x->sz;i<size;i++)
X            (x->p)[i] = 0;
X        x->sz = size;
X    }
X}
X
Xvoid dgset(x,i,n)
XMP_INT *x; unsigned int i; mp_limb n;
X{
X    if (n) {
X        if (i >= x->sz) {
X            _mpz_realloc(x,(mp_size)(i+1));
X            x->sz = i+1;
X        }
X        (x->p)[i] = n;
X    }
X}
X
Xstatic unsigned int digits(x)
XMP_INT *x;
X{
X    int i;
X    for (i = (x->sz) - 1; i>=0 && (x->p)[i] == 0 ; i--)
X        ;
X    return i+1;
X}
X
X/* y = x */
Xvoid mpz_set(y,x)
XMP_INT *y; MP_INT *x;
X{
X    unsigned int i,k = x->sz;
X    if (y->sz < k) {
X        k=digits(x);
X        _mpz_realloc(y,(mp_size)k);
X    }
X    if (y->sz > x->sz) {
X        mpz_clear(y); mpz_init(y); _mpz_realloc(y,(mp_size)(x->sz));
X    }
X
X    for (i=0;i < k; i++)
X        (y->p)[i] = (x->p)[i];
X
X    for (;i<y->sz;i++)
X        (y->p)[i] = 0;
X
X    y->sn = x->sn;
X}
X
Xvoid mpz_set_ui(y,v)
XMP_INT *y; unsigned long v; {
X    int i;
X    for (i=1;i<y->sz;i++)
X        y->p[i] = 0;
X    y->p[0] = (v & LMAX);
X    y->p[1] = (v & LC) >> (DIGITBITS);
X    if (v)
X        y->sn=1;
X    else
X        y->sn=0;
X}
X
Xunsigned long mpz_get_ui(y)
XMP_INT *y; {
X    return (y->p[0] | (y->p[1] << DIGITBITS));
X}
X
Xlong mpz_get_si(y)
XMP_INT *y; {
X    if (y->sn == 0)
X        return 0;
X    else
X        return (y->sn * (y->p[0] | (y->p[1] & 1) << DIGITBITS));
X}
X
Xvoid mpz_set_si(y,v)
XMP_INT *y; long v; {
X    int i;
X    for (i=1;i<y->sz;i++)
X        y->p[i] = 0;
X    if (v < 0) {
X        y->sn = -1;
X        y->p[0] = (-v) & LMAX;
X        y->p[1] = ((-v) & LC) >> DIGITBITS;
X    }
X    else if (v > 0) {
X        y->sn = 1;
X        y->p[0] = v & LMAX;
X        y->p[1] = (v & LC) >> DIGITBITS;
X    }
X    else {
X        y->sn=0;
X        y->p[0] = 0;
X        y->p[1] = 0;
X    }
X}
X
X/* z = x + y, without regard for sign */
Xstatic void uadd(z,x,y)
XMP_INT *z; MP_INT *x; MP_INT *y;
X{
X    mp_limb c;
X    int i;
X    MP_INT *t;
X
X    if (y->sz < x->sz) {
X        t=x; x=y; y=t;
X    }
X
X    /* now y->sz >= x->sz */
X
X    _mpz_realloc(z,(mp_size)((y->sz)+1));
X
X    c=0;
X    for (i=0;i<x->sz;i++) {
X        if (( z->p[i] = y->p[i] + x->p[i] + c ) & CMASK) {
X            c=1;
X            (z->p[i]) &=LMAX;
X        }
X        else
X            c=0;
X    }
X    for (;i<y->sz;i++) {
X        if ((z->p[i] = (y->p[i] + c)) & CMASK)
X            z->p[i]=0;
X        else
X            c=0;
X    }
X    (z->p)[y->sz]=c;
X}
X
X/* z = y - x, ignoring sign */
X/* precondition: abs(y) >= abs(x) */
Xstatic void usub(z,y,x)
XMP_INT *z; MP_INT *y; MP_INT *x;
X{
X    int i;
X    mp_limb b,m;
X    _mpz_realloc(z,(mp_size)(y->sz));
X    b=0;
X    for (i=0;i<y->sz;i++) {
X        m=((y->p)[i]-b)-dg(x,i);
X        if (m < 0) {
X            b = 1;
X            m = LMAX + 1 + m;
X        }
X        else
X            b = 0;
X        z->p[i] = m;
X    }
X}
X
X/* compare abs(x) and abs(y) */
Xstatic int ucmp(y,x)
XMP_INT *y;MP_INT *x;
X{
X    int i;
X    for (i=imax(x->sz,y->sz)-1;i>=0;i--) {
X        if (dg(y,i) < dg(x,i))
X            return (-1);
X        else if (dg(y,i) > dg(x,i))
X            return 1;
X    }
X    return 0;
X}
X
Xint mpz_cmp(x,y)
XMP_INT *x; MP_INT *y;
X{
X    int abscmp;
X    if (x->sn < 0 && y->sn > 0)
X        return (-1);
X    if (x->sn > 0 && y->sn < 0)
X        return 1;
X    abscmp=ucmp(x,y);
X    if (x->sn >=0 && y->sn >=0)
X        return abscmp;
X    if (x->sn <=0 && y->sn <=0)
X        return (-abscmp);
X}
X
X/* is abs(x) == 0 */
Xstatic int uzero(x)
XMP_INT *x;
X{
X    int i;
X    for (i=0; i < x->sz; i++)
X        if ((x->p)[i] != 0)
X            return 0;
X    return 1;
X}
X
Xstatic void zero(x)
XMP_INT *x;
X{
X    int i;
X    x->sn=0;
X    for (i=0;i<x->sz;i++)
X        (x->p)[i] = 0;
X}
X
X
X/* w = u * v */
Xvoid mpz_mul(ww,u,v)
XMP_INT *ww;MP_INT *u; MP_INT *v; {
X    int i,j;
X    mp_limb t0,t1,t2,t3;
X    mp_limb cc;
X    MP_INT *w;
X
X    w = (MP_INT *)malloc(sizeof(MP_INT));
X    mpz_init(w);
X    _mpz_realloc(w,(mp_size)(u->sz + v->sz));
X    for (j=0; j < 2*u->sz; j++) {
X        cc = (mp_limb)0;
X        t3 = hd(u,j);
X            for (i=0; i < 2*v->sz; i++) {
X                t0 = t3 * hd(v,i);
X                t1 = HIGH(t0); t0 = LOW(t0);
X                if ((i+j)%2)
X                    t2 = HIGH(w->p[(i+j)/2]);
X                else
X                    t2 = LOW(w->p[(i+j)/2]);
X                t2 += cc;
X                if (t2 & HCMASK) {
X                    cc = 1; t2&=HLMAX;
X                }
X                else
X                    cc = 0;
X                t2 += t0;
X                if (t2 & HCMASK) {
X                    cc++ ; t2&=HLMAX;
X                }
X                cc+=t1;
X                if ((i+j)%2)
X                    w->p[(i+j)/2] = LOW(w->p[(i+j)/2]) |
X                        (t2 << HALFDIGITBITS);
X                else
X                    w->p[(i+j)/2] = (HIGH(w->p[(i+j)/2]) << HALFDIGITBITS) |
X                        t2;
X
X            }
X        if (cc) {
X            if ((j+i)%2)
X                w->p[(i+j)/2] += cc << HALFDIGITBITS;
X            else
X                w->p[(i+j)/2] += cc;
X        }
X    }
X    w->sn = (u->sn) * (v->sn);
X    if (w != ww) {
X        mpz_set(ww,w);
X        mpz_clear(w);
X        free(w);
X    }
X}
X
X/* n must be < DIGITBITS */
Xstatic void urshift(c1,a,n)
XMP_INT *c1;MP_INT *a;int n;
X{
X    mp_limb cc = 0;
X    if (n >= DIGITBITS)
X        fatal("urshift: n >= DIGITBITS");
X    if (n == 0)
X        mpz_set(c1,a);
X    else {
X        MP_INT c; int i;
X        mp_limb rm = (((mp_limb)1<<n) - 1);
X        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz));
X        for (i=a->sz-1;i >= 0; i--) {
X            c.p[i] = ((a->p[i] >> n) | cc) & LMAX;
X            cc = (a->p[i] & rm) << (DIGITBITS - n);
X        }
X        mpz_set(c1,&c);
X        mpz_clear(&c);
X    }
X}
X
X/* n must be < DIGITBITS */
Xstatic void ulshift(c1,a,n)
XMP_INT *c1;MP_INT *a;int n;
X{
X    mp_limb cc = 0;
X    if (n >= DIGITBITS)
X        fatal("ulshift: n >= DIGITBITS");
X    if (n == 0)
X        mpz_set(c1,a);
X    else {
X        MP_INT c; int i;
X        mp_limb rm = (((mp_limb)1<<n) - 1) << (DIGITBITS -n);
X        mpz_init(&c); _mpz_realloc(&c,(mp_size)(a->sz + 1));
X        for (i=0;i < a->sz; i++) {
X            c.p[i] = ((a->p[i] << n) | cc) & LMAX;
X            cc = (a->p[i] & rm) >> (DIGITBITS -n);
X        }
X        c.p[i] = cc;
X        mpz_set(c1,&c);
X        mpz_clear(&c);
X    }
X}
X
Xvoid mpz_div_2exp(z,x,e)
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set(z,x);
X    else {
X        unsigned long i;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        MP_INT y; mpz_init(&y);
X        _mpz_realloc(&y,(mp_size)((x->sz) - digs));
X        for (i=0; i < (x->sz - digs); i++)
X            (y.p)[i] = (x->p)[i+digs];
X        if (bs) {
X            urshift(z,&y,bs);
X        }
X        else {
X            mpz_set(z,&y);
X        }
X        if (uzero(z))
X            z->sn = 0;
X        else
X            z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X
Xvoid mpz_mul_2exp(z,x,e)
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set(z,x);
X    else {
X        unsigned long i;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        MP_INT y; mpz_init(&y);
X        _mpz_realloc(&y,(mp_size)((x->sz)+digs));
X        for (i=digs;i<((x->sz) + digs);i++)
X            (y.p)[i] = (x->p)[i - digs];
X        if (bs) {
X            ulshift(z,&y,bs);
X        }
X        else {
X            mpz_set(z,&y);
X        }
X        z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X
Xvoid mpz_mod_2exp(z,x,e)
XMP_INT *z; MP_INT *x; unsigned long e;
X{
X    short sn = x->sn;
X    if (e==0)
X        mpz_set_ui(z,0L);
X    else {
X        unsigned long i;
X        MP_INT y;
X        long digs = (e / DIGITBITS);
X        unsigned long bs = (e % (DIGITBITS));
X        if (digs > x->sz || (digs == (x->sz) && bs)) {
X            mpz_set(z,x);
X            return;
X        }
X
X        mpz_init(&y);
X        if (bs)
X            _mpz_realloc(&y,(mp_size)(digs+1));
X        else
X            _mpz_realloc(&y,(mp_size)digs);
X        for (i=0; i<digs; i++)
X            (y.p)[i] = (x->p)[i];
X        if (bs) {
X            (y.p)[digs] = ((x->p)[digs]) & (((mp_limb)1 << bs) - 1);
X        }
X        mpz_set(z,&y);
X        if (uzero(z))
X            z->sn = 0;
X        else
X            z->sn = sn;
X        mpz_clear(&y);
X    }
X}
X
X
X/* internal routine to compute x/y and x%y ignoring signs */
Xstatic void udiv(qq,rr,xx,yy)
XMP_INT *qq; MP_INT *rr; MP_INT *xx; MP_INT *yy;
X{
X    MP_INT *q, *x, *y, *r;
X    int ns,xd,yd,j,f,i,ccc;
X    mp_limb zz,z,qhat,b,u,m;
X    ccc=0;
X
X    if (uzero(yy))
X        return;
X    q = (MP_INT *)malloc(sizeof(MP_INT));
X    r = (MP_INT *)malloc(sizeof(MP_INT));
X    x = (MP_INT *)malloc(sizeof(MP_INT));
X    y = (MP_INT *)malloc(sizeof(MP_INT));
X    if (!x || !y || !q || !r)
X        fatal("udiv: cannot allocate memory");
X    mpz_init(q); mpz_init(x);mpz_init(y);mpz_init(r);
X    _mpz_realloc(x,(mp_size)((xx->sz)+1));
X    yd = digits(yy);
X    ns = lzb(yy->p[yd-1]);
X    ulshift(x,xx,ns);
X    ulshift(y,yy,ns);
X    xd = digits(x);
X    _mpz_realloc(q,(mp_size)xd);
X    xd*=2; yd*=2;
X    z = hd(y,yd-1);;
X    for (j=(xd-yd);j>=0;j--) {
X#ifdef DEBUG
X    printf("y=");
X    for (i=yd-1;i>=0;i--)
X        printf(" %04lx", (long)hd(y,i));
X    printf("\n");
X    printf("x=");
X    for (i=xd-1;i>=0;i--)
X        printf(" %04lx", (long)hd(x,i));
X    printf("\n");
X    printf("z=%08lx\n",(long)z);
X#endif
X
X        if (z == LMAX)
X            qhat = hd(x,j+yd);
X        else {
X            qhat = ((hd(x,j+yd)<< HALFDIGITBITS) + hd(x,j+yd-1)) / (z+1);
X        }
X#ifdef DEBUG
X    printf("qhat=%08lx\n",(long)qhat);
X    printf("hd=%04lx/%04lx\n",(long)hd(x,j+yd),(long)hd(x,j+yd-1));
X#endif
X        b = 0; zz=0;
X        if (qhat) {
X            for (i=0; i < yd; i++) {
X                zz = qhat * hd(y,i);
X                u = hd(x,i+j);
X                u-=b;
X                if (u<0) {
X                    b=1; u+=HLMAX+1;
X                }
X                else
X                    b=0;
X                u-=LOW(zz);
X                if (u < 0) {
X                    b++;
X                    u+=HLMAX+1;
X                }
X                b+=HIGH(zz);
X                if ((i+j)%2)
X                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (u << HALFDIGITBITS);
X                else
X                    x->p[(i+j)/2] = (HIGH(x->p[(i+j)/2])
X                        << HALFDIGITBITS) | u;
X            }
X            if (b) {
X                if ((j+i)%2)
X                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
X                else
X                    x->p[(i+j)/2] -= b;
X            }
X        }
X#ifdef DEBUG
X        printf("x after sub=");
X        for (i=xd-1;i>=0;i--)
X            printf(" %04lx", (long)hd(x,i));
X        printf("\n");
X#endif
X        for(;;zz++) {
X            f=1;
X            if (!hd(x,j+yd)) {
X                for(i=yd-1; i>=0; i--) {
X                    if (hd(x,j+i) > hd(y,i)) {
X                        f=1;
X                        break;
X                    }
X                    if (hd(x,j+i) < hd(y,i)) {
X                        f=0;
X                        break;
X                    }
X                }
X            }
X            if (!f)
X                break;
X            qhat++;
X            ccc++;
X#ifdef DEBUG
X            printf("corrected qhat = %08lx\n", (long)qhat);
X#endif
X            b=0;
X            for (i=0;i<yd;i++) {
X                m = hd(x,i+j)-hd(y,i)-b;
X                if (m < 0) {
X                    b = 1;
X                    m = HLMAX + 1 + m;
X                }
X                else
X                    b = 0;
X                if ((i+j)%2)
X                    x->p[(i+j)/2] = LOW(x->p[(i+j)/2]) | (m << HALFDIGITBITS);
X                else
X                    x->p[(i+j)/2]
X                        = (HIGH(x->p[(i+j)/2]) << HALFDIGITBITS) | m;
X            }
X            if (b) {
X                if ((j+i)%2)
X                    x->p[(i+j)/2] -= b << HALFDIGITBITS;
X                else
X                    x->p[(i+j)/2] -= b;
X            }
X#ifdef DEBUG
X            printf("x after cor=");
X            for (i=2*(x->sz)-1;i>=0;i--)
X                printf(" %04lx", (long)hd(x,i));
X            printf("\n");
X#endif
X        }
X        if (j%2)
X            q->p[j/2] |= qhat << HALFDIGITBITS;
X        else
X            q->p[j/2] |= qhat;
X#ifdef DEBUG
X            printf("x after cor=");
X            for (i=xd - 1;i>=0;i--)
X                printf(" %04lx", (long)hd(q,i));
X            printf("\n");
X#endif
X    }
X    _mpz_realloc(r,(mp_size)(yy->sz));
X    zero(r);
X    urshift(r,x,ns);
X    mpz_set(rr,r);
X    mpz_set(qq,q);
X    mpz_clear(x); mpz_clear(y);
X    mpz_clear(q);  mpz_clear(r);
X    free(q); free(x); free(y); free(r);
X}
X
Xvoid mpz_add(zz,x,y)
XMP_INT *zz;MP_INT *x;MP_INT *y;
X{
X    int mg;
X    MP_INT *z;
X    if (x->sn == 0) {
X        mpz_set(zz,y);
X        return;
X    }
X    if (y->sn == 0) {
X        mpz_set(zz,x);
X        return;
X    }
X    z = (MP_INT *)malloc(sizeof(MP_INT));
X    mpz_init(z);
X
X    if (x->sn > 0 && y->sn > 0) {
X        uadd(z,x,y); z->sn = 1;
X    }
X    else if (x->sn < 0 && y->sn < 0) {
X        uadd(z,x,y); z->sn = -1;
X    }
X    else {
X        /* signs differ */
X        if ((mg = ucmp(x,y)) == 0) {
X            zero(z);
X        }
X        else if (mg > 0) {  /* abs(y) < abs(x) */
X            usub(z,x,y);
X            z->sn = (x->sn > 0 && y->sn < 0) ? 1 : (-1);
X        }
X        else { /* abs(y) > abs(x) */
X            usub(z,y,x);
X            z->sn = (x->sn < 0 && y->sn > 0) ? 1 : (-1);
X        }
X    }
X    mpz_set(zz,z);
X    mpz_clear(z);
X    free(z);
X}
X
Xvoid mpz_add_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_add(x,y,&z);
X    mpz_clear(&z);
X}
X
Xint mpz_cmp_ui(x,n)
XMP_INT *x; unsigned long n;
X{
X    MP_INT z; int ret;
X    mpz_init_set_ui(&z,n);
X    ret=mpz_cmp(x,&z);
X    mpz_clear(&z);
X    return ret;
X}
X
Xint mpz_cmp_si(x,n)
XMP_INT *x; long n;
X{
X    MP_INT z; int ret;
X    mpz_init(&z);
X    mpz_set_si(&z,n);
X    ret=mpz_cmp(x,&z);
X    mpz_clear(&z);
X    return ret;
X}
X
X
Xvoid mpz_mul_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mul(x,y,&z);
X    mpz_clear(&z);
X}
X
X
X/* z = x - y  -- just use mpz_add - I'm lazy */
Xvoid mpz_sub(z,x,y)
XMP_INT *z;MP_INT *x; MP_INT *y;
X{
X    MP_INT u;
X    mpz_init(&u);
X    mpz_set(&u,y);
X    u.sn = -(u.sn);
X    mpz_add(z,x,&u);
X    mpz_clear(&u);
X}
X
Xvoid mpz_sub_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_sub(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_div(q,x,y)
XMP_INT *q; MP_INT *x; MP_INT *y;
X{
X    MP_INT r;
X    short sn1 = x->sn, sn2 = y->sn;
X    mpz_init(&r);
X    udiv(q,&r,x,y);
X    q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    mpz_clear(&r);
X}
X
Xvoid mpz_mdiv(q,x,y)
XMP_INT *q; MP_INT *x; MP_INT *y;
X{
X    MP_INT r;
X    short sn1 = x->sn, sn2 = y->sn, qsign;
X    mpz_init(&r);
X    udiv(q,&r,x,y);
X    qsign = q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    /* now if r != 0 and q < 0 we need to round q towards -inf */
X    if (!uzero(&r) && qsign < 0)
X        mpz_sub_ui(q,q,1L);
X    mpz_clear(&r);
X}
X
Xvoid mpz_mod(r,x,y)
XMP_INT *r; MP_INT *x; MP_INT *y;
X{
X    MP_INT q;
X    short sn = x->sn;
X    mpz_init(&q);
X    if (x->sn == 0) {
X        zero(r);
X        return;
X    }
X    udiv(&q,r,x,y);
X    r->sn = sn;
X    if (uzero(r))
X        r->sn = 0;
X    mpz_clear(&q);
X}
X
Xvoid mpz_divmod(q,r,x,y)
XMP_INT *q; MP_INT *r; MP_INT *x; MP_INT *y;
X{
X    short sn1 = x->sn, sn2 = y->sn;
X    if (x->sn == 0) {
X        zero(q);
X        zero(r);
X        return;
X    }
X    udiv(q,r,x,y);
X    q->sn = sn1*sn2;
X    if (uzero(q))
X        q->sn = 0;
X    r->sn = sn1;
X    if (uzero(r))
X        r->sn = 0;
X}
X
Xvoid mpz_mmod(r,x,y)
XMP_INT *r; MP_INT *x; MP_INT *y;
X{
X    MP_INT q;
X    short sn1 = x->sn, sn2 = y->sn;
X    mpz_init(&q);
X    if (sn1 == 0) {
X        zero(r);
X        return;
X    }
X    udiv(&q,r,x,y);
X    if (uzero(r)) {
X        r->sn = 0;
X        return;
X    }
X    q.sn = sn1*sn2;
X    if (q.sn > 0)
X        r->sn = sn1;
X    else if (sn1 < 0 && sn2 > 0) {
X        r->sn = 1;
X        mpz_sub(r,y,r);
X    }
X    else {
X        r->sn = 1;
X        mpz_add(r,y,r);
X    }
X}
X
Xvoid mpz_mdivmod(q,r,x,y)
XMP_INT *q;MP_INT *r; MP_INT *x; MP_INT *y;
X{
X    short sn1 = x->sn, sn2 = y->sn, qsign;
X    if (sn1 == 0) {
X        zero(q);
X        zero(r);
X        return;
X    }
X    udiv(q,r,x,y);
X    qsign = q->sn = sn1*sn2;
X    if (uzero(r)) {
X        /* q != 0, since q=r=0 would mean x=0, which was tested above */
X        r->sn = 0;
X        return;
X    }
X    if (q->sn > 0)
X        r->sn = sn1;
X    else if (sn1 < 0 && sn2 > 0) {
X        r->sn = 1;
X        mpz_sub(r,y,r);
X    }
X    else {
X        r->sn = 1;
X        mpz_add(r,y,r);
X    }
X    if (uzero(q))
X        q->sn = 0;
X    /* now if r != 0 and q < 0 we need to round q towards -inf */
X    if (!uzero(r) && qsign < 0)
X        mpz_sub_ui(q,q,1L);
X}
X
Xvoid mpz_mod_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init(&z);
X    mpz_set_ui(&z,n);
X    mpz_mod(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mmod_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init(&z);
X    mpz_set_ui(&z,n);
X    mpz_mmod(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_div_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_div(x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mdiv_ui(x,y,n)
XMP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mdiv(x,y,&z);
X    mpz_clear(&z);
X}
Xvoid mpz_divmod_ui(q,x,y,n)
XMP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_divmod(q,x,y,&z);
X    mpz_clear(&z);
X}
X
Xvoid mpz_mdivmod_ui(q,x,y,n)
XMP_INT *q;MP_INT *x;MP_INT *y; unsigned long n;
X{
X    MP_INT z;
X    mpz_init_set_ui(&z,n);
X    mpz_mdivmod(q,x,y,&z);
X    mpz_clear(&z);
X}
X
X
X/* 2<=base <=36 - this overestimates the optimal value, which is OK */
Xunsigned int mpz_sizeinbase(x,base)
XMP_INT *x; int base;
X{
X    int i,j;
X    int bits = digits(x) * DIGITBITS;
X    if (base < 2 || base > 36)
X        fatal("mpz_sizeinbase: invalid base");
X    for (j=0,i=1; i<=base;i*=2,j++)
X        ;
X    return ((bits)/(j-1)+1);
X}
X
Xchar *mpz_get_str(s,base,x)
Xchar *s;  int base; MP_INT *x; {
X    MP_INT xx,q,r,bb;
X    char *p,*t,*ps;
X    int sz = mpz_sizeinbase(x,base);
X    mp_limb d;
X    if (base < 2 || base > 36)
X        return s;
X    t = (char *)malloc(sz+2);
X    if (!t)
X        fatal("cannot allocate memory in mpz_get_str");
X    if (!s) {
X        s=(char *)malloc(sz+2);
X        if (!s)
X            fatal("cannot allocate memory in mpz_get_str");
X    }
X    if (uzero(x)) {
X        *s='0';
X        *(s+1)='\0';
X        return s;
X    }
X    mpz_init(&xx); mpz_init(&q); mpz_init(&r); mpz_init(&bb);
X    mpz_set(&xx,x);
X    mpz_set_ui(&bb,(unsigned long)base);
X    ps = s;
X    if (x->sn < 0) {
X        *ps++= '-';
X        xx.sn = 1;
X    }
X    p = t;
X    while (!uzero(&xx)) {
X        udiv(&xx,&r,&xx,&bb);
X        d = r.p[0];
X        if (d < 10)
X            *p++  = (char) (r.p[0] + '0');
X        else
X            *p++  = (char) (r.p[0] + -10 + 'a');
X    }
X
X    p--;
X    for (;p>=t;p--,ps++)
X        *ps = *p;
X    *ps='\0';
X
X    mpz_clear(&q); mpz_clear(&r); free(t);
X    return s;
X}
X
X
Xint mpz_set_str(x,s,base)
XMP_INT *x; char *s; int base;
X{
X    int l,i;
X    int retval = 0;
X    MP_INT t,m,bb;
X    short sn;
X    unsigned int k;
X    mpz_init(&m);
X    mpz_init(&t);
X    mpz_init(&bb);
X    mpz_set_ui(&m,1L);
X    zero(x);
X    while (*s==' ' || *s=='\t' || *s=='\n')
X        s++;
X    if (*s == '-') {
X        sn = -1; s++;
X    }
X    else
X        sn = 1;
X    if (base == 0) {
X        if (*s == '0') {
X            if (*(s+1) == 'x' || *(s+1) == 'X') {
X                base = 16;
X                s+=2;   /* toss 0x */
X            }
X            else {
X                base = 8;
X                s++;    /* toss 0 */
X            }
X        }
X        else
X            base=10;
X    }
X    if (base < 2 || base > 36)
X        fatal("mpz_set_str: invalid base");
X    mpz_set_ui(&bb,(unsigned long)base);
X    l=strlen(s);
X    for (i = l-1; i>=0; i--) {
X        if (s[i]==' ' || s[i]=='\t' || s[i]=='\n')
X            continue;
X        if (s[i] >= '0' && s[i] <= '9')
X            k = (unsigned int)s[i] - (unsigned int)'0';
X        else if (s[i] >= 'A' && s[i] <= 'Z')
X            k = (unsigned int)s[i] - (unsigned int)'A'+10;
X        else if (s[i] >= 'a' && s[i] <= 'z')
X            k = (unsigned int)s[i] - (unsigned int)'a'+10;
X        else {
X            retval = (-1);
X            break;
X        }
X        if (k >= base) {
X            retval = (-1);
X            break;
X        }
X        mpz_mul_ui(&t,&m,(unsigned long)k);
X        mpz_add(x,x,&t);
X        mpz_mul(&m,&m,&bb);
X#ifdef DEBUG
X        printf("k=%d\n",k);
X        printf("t=%s\n",mpz_get_str(NULL,10,&t));
X        printf("x=%s\n",mpz_get_str(NULL,10,x));
X        printf("m=%s\n",mpz_get_str(NULL,10,&m));
X#endif
X    }
X    if (x->sn)
X        x->sn = sn;
X    mpz_clear(&m);
X    mpz_clear(&bb);
X    mpz_clear(&t);
X    return retval;
X}
X
Xint mpz_init_set_str(x,s,base)
XMP_INT *x; char *s; int base;
X{
X    mpz_init(x); return mpz_set_str(x,s,base);
X}
X
X#define mpz_randombyte (rand()& 0xff)
X
Xvoid mpz_random(x,size)
XMP_INT *x; unsigned int size;
X{
X    unsigned int bits = size * LONGBITS;
X    unsigned int digits = bits/DIGITBITS;
X    unsigned int oflow = bits % DIGITBITS;
X    unsigned int i,j;
X    long t;
X    if (oflow)
X        digits++;
X    _mpz_realloc(x,(mp_size)digits);
X
X    for (i=0; i<digits; i++) {
X        t = 0;
X        for (j=0; j < DIGITBITS; j+=8)
X            t = (t << 8) | mpz_randombyte;
X        (x->p)[i] = t & LMAX;
X    }
X    if (oflow)
X        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
X}
Xvoid mpz_random2(x,size)
XMP_INT *x; unsigned int size;
X{
X    unsigned int bits = size * LONGBITS;
X    unsigned int digits = bits/DIGITBITS;
X    unsigned int oflow = bits % DIGITBITS;
X    unsigned int i,j;
X    long t;
X    if (oflow)
X        digits++;
X    _mpz_realloc(x,(mp_size)digits);
X
X    for (i=0; i<digits; i++) {
X        t = 0;
X        for (j=0; j < DIGITBITS; j+=8)
X            t = (t << 8) | mpz_randombyte;
X        (x->p)[i] = (t & LMAX) % 2;
X    }
X    if (oflow)
X        (x->p)[digits-1] &= (((mp_limb)1 << oflow) - 1);
X}
X
Xsize_t mpz_size(x)
XMP_INT *x;
X{
X    int bits = (x->sz)*DIGITBITS;
X    size_t r;
X
X    r = bits/LONGBITS;
X    if (bits % LONGBITS)
X        r++;
X    return r;
X}
X
Xvoid mpz_abs(x,y)
XMP_INT *x; MP_INT *y;
X{
X    if (x!=y)
X        mpz_set(x,y);
X    if (uzero(x))
X        x->sn = 0;
X    else
X        x->sn = 1;
X}
X
Xvoid mpz_neg(x,y)
XMP_INT *x; MP_INT *y;
X{
X    if (x!=y)
X        mpz_set(x,y);
X    x->sn = -(y->sn);
X}
X
Xvoid mpz_fac_ui(x,n)
XMP_INT *x; unsigned long n;
X{
X    mpz_set_ui(x,1L);
X    if (n==0 || n == 1)
X        return;
X    for (;n>1;n--)
X        mpz_mul_ui(x,x,n);
X}
X
Xvoid mpz_gcd(gg,aa,bb)
XMP_INT *gg;
XMP_INT *aa;
XMP_INT *bb;
X{
X    MP_INT *g,*t,*a,*b;
X    int freeg;
X
X    t = (MP_INT *)malloc(sizeof(MP_INT));
X    g = (MP_INT *)malloc(sizeof(MP_INT));
X    a = (MP_INT *)malloc(sizeof(MP_INT));
X    b = (MP_INT *)malloc(sizeof(MP_INT));
X    if (!a || !b || !g || !t)
X        fatal("mpz_gcd: cannot allocate memory");
X    mpz_init(g); mpz_init(t); mpz_init(a); mpz_init(b);
X    mpz_abs(a,aa); mpz_abs(b,bb);
X
X    while (!uzero(b)) {
X        mpz_mod(t,a,b);
X        mpz_set(a,b);
X        mpz_set(b,t);
X    }
X
X    mpz_set(gg,a);
X    mpz_clear(g); mpz_clear(t); mpz_clear(a); mpz_clear(b);
X    free(g); free(t); free(a); free(b);
X}
X
Xvoid mpz_gcdext(gg,xx,yy,aa,bb)
XMP_INT *gg,*xx,*yy,*aa,*bb;
X{
X    MP_INT *x, *y, *g, *u, *q;
X
X    if (uzero(bb)) {
X        mpz_set(gg,aa);
X        mpz_set_ui(xx,1L);
X        if (yy)
X            mpz_set_ui(yy,0L);
X        return;
X    }
X
X    g = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(g);
X    q = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(q);
X    y = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(y);
X    x = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(x);
X    u = (MP_INT *)malloc(sizeof(MP_INT)); mpz_init(u);
X
X    if (!g || !q || !y || !x || !u)
X        fatal("mpz_gcdext: cannot allocate memory");
X
X    mpz_divmod(q,u,aa,bb);
X    mpz_gcdext(g,x,y,bb,u);
X    if (yy) {
X        mpz_mul(u,q,y);
X        mpz_sub(yy,x,u);
X    }
X    mpz_set(xx,y);
X    mpz_set(gg,g);
X    mpz_clear(g); mpz_clear(q); mpz_clear(y); mpz_clear(x); mpz_clear(u);
X    free(g); free(q); free(y); free(x); free(u);
X}
X
X
X/*
X *    a is a quadratic residue mod b if
X *    x^2 = a mod b      has an integer solution x.
X *
X *    J(a,b) = if a==1 then 1 else
X *             if a is even then J(a/2,b) * ((-1)^(b^2-1)/8))
X *             else J(b mod a, a) * ((-1)^((a-1)*(b-1)/4)))
X *
X *  b>0  b odd
X *
X *
X */
X
Xint mpz_jacobi(ac,bc)
XMP_INT *ac, *bc;
X{
X    int sgn = 1;
X    unsigned long c;
X    MP_INT *t,*a,*b;
X    if (bc->sn <=0)
X        fatal("mpz_jacobi call with b <= 0");
X    if (even(bc))
X        fatal("mpz_jacobi call with b even");
X
X    init(t); init(a); init(b);
X
X    /* if ac < 0, then we use the fact that
X     *  (-1/b)= -1 if b mod 4 == 3
X     *          +1 if b mod 4 == 1
X     */
X
X    if (mpz_cmp_ui(ac,0L) < 0 && (lowdigit(bc) % 4) == 3)
X        sgn=-sgn;
X
X    mpz_abs(a,ac); mpz_set(b,bc);
X
X    /* while (a > 1) { */
X    while(mpz_cmp_ui(a,1L) > 0) {
X
X        if (even(a)) {
X
X            /* c = b % 8 */
X            c = lowdigit(b) & 0x07;
X
X            /* b odd, then (b^2-1)/8 is even iff (b%8 == 3,5) */
X
X            /* if b % 8 == 3 or 5 */
X            if (c == 3 || c == 5)
X                sgn = -sgn;
X
X            /* a = a / 2 */
X            mpz_div_2exp(a,a,1L);
X
X        }
X        else {
X            /* note: (a-1)(b-1)/4 odd iff a mod 4==b mod 4==3 */
X
X            /* if a mod 4 == 3 and b mod 4 == 3 */
X            if (((lowdigit(a) & 3) == 3) && ((lowdigit(b) & 3) == 3))
X                sgn = -sgn;
X
X            /* set (a,b) = (b mod a,a) */
X            mpz_set(t,a); mpz_mmod(a,b,t); mpz_set(b,t);
X        }
X    }
X
X    /* if a == 0 then sgn = 0 */
X    if (uzero(a))
X        sgn=0;
X
X    clear(t); clear(a); clear(b);
X    return (sgn);
X}
X
Xvoid mpz_and(z,x,y) /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) & dg(y,i);
X    if (x->sn < 0 && y->sn < 0)
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
X
Xvoid mpz_or(z,x,y)  /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) | dg(y,i);
X    if (x->sn < 0 || y->sn < 0)
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
X
Xvoid mpz_xor(z,x,y)  /* not the most efficient way to do this */
XMP_INT *z,*x,*y;
X{
X    int i,sz;
X    sz = imax(x->sz, y->sz);
X    _mpz_realloc(z,(mp_size)sz);
X    for (i=0; i < sz; i++)
X        (z->p)[i] = dg(x,i) ^ dg(y,i);
X    if ((x->sn <= 0 && y->sn > 0) || (x->sn > 0 && y->sn <=0))
X        z->sn = (-1);
X    else
X        z->sn = 1;
X    if (uzero(z))
X        z->sn = 0;
X}
Xvoid mpz_pow_ui(zz,x,e)
XMP_INT *zz, *x;
Xunsigned long e;
X{
X    MP_INT *t;
X    unsigned long mask = (1L<< (LONGBITS-1));
X
X    if (e==0) {
X        mpz_set_ui(zz,1L);
X        return;
X    }
X
X    init(t);
X    mpz_set(t,x);
X    for (;!(mask &e); mask>>=1)
X        ;
X    mask>>=1;
X    for (;mask!=0; mask>>=1) {
X        mpz_mul(t,t,t);
X        if (e & mask)
X            mpz_mul(t,t,x);
X    }
X    mpz_set(zz,t);
X    clear(t);
X}
X
Xvoid mpz_powm(zz,x,ee,n)
XMP_INT *zz,*x,*ee,*n;
X{
X    MP_INT *t,*e;
X    struct is *stack = NULL;
X    int k,i;
X
X    if (uzero(ee)) {
X        mpz_set_ui(zz,1L);
X        return;
X    }
X
X    if (ee->sn < 0) {
X        return;
X    }
X    init(e); init(t); mpz_set(e,ee);
X
X    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
X        push(lowdigit(e) & 1,&stack);
X    k--;
X    i=pop(&stack);
X
X    mpz_mod(t,x,n);  /* t=x%n */
X
X    for (i=k-1;i>=0;i--) {
X        mpz_mul(t,t,t);
X        mpz_mod(t,t,n);
X        if (pop(&stack)) {
X            mpz_mul(t,t,x);
X            mpz_mod(t,t,n);
X        }
X    }
X    mpz_set(zz,t);
X    clear(t);
X}
X
Xvoid mpz_powm_ui(z,x,e,n)
XMP_INT *z,*x,*n; unsigned long e;
X{
X    MP_INT f;
X    mpz_init(&f);
X    mpz_set_ui(&f,e);
X    mpz_powm(z,x,&f,n);
X    mpz_clear(&f);
X}
X
X
X
X/* Miller-Rabin */
Xstatic int witness(x,n)
XMP_INT *x, *n;
X{
X    MP_INT *t,*e, *e1;
X    struct is *stack = NULL;
X    int trivflag = 0;
X    int k,i;
X
X    init(e1); init(e); init(t); mpz_sub_ui(e,n,1L); mpz_set(e1,e);
X
X    for (k=0;!uzero(e);k++,mpz_div_2exp(e,e,1L))
X        push(lowdigit(e) & 1,&stack);
X    k--;
X    i=pop(&stack);
X
X    mpz_mod(t,x,n);  /* t=x%n */
X
X    for (i=k-1;i>=0;i--) {
X        trivflag = !mpz_cmp_ui(t,1L) || !mpz_cmp(t,e1);
X        mpz_mul(t,t,t); mpz_mod(t,t,n);
X        if (!trivflag && !mpz_cmp_ui(t,1L)) {
X            clear(t); clear(e); clear(e1);
X            return 1;
X        }
X
X        if (pop(&stack)) {
X            mpz_mul(t,t,x);
X            mpz_mod(t,t,n);
X        }
X    }
X    if (mpz_cmp_ui(t,1L))  { /* t!=1 */
X        clear(t); clear(e); clear(e1);
X        return 1;
X    }
X    else {
X        clear(t); clear(e); clear(e1);
X        return 0;
X    }
X}
X
Xunsigned long smallp[] = {2,3,5,7,11,13,17};
Xint mpz_probab_prime_p(nn,s)
XMP_INT *nn; int s;
X{
X    MP_INT *a,*n;
X    int j,k,i;
X    long t;
X
X    if (uzero(nn))
X        return 0;
X    init(n); mpz_abs(n,nn);
X    if (!mpz_cmp_ui(n,1L)) {
X        clear(n);
X        return 0;
X    }
X    init(a);
X    for (i=0;i<sizeof(smallp)/sizeof(unsigned long); i++) {
X        mpz_mod_ui(a,n,smallp[i]);
X        if (uzero(a)) {
X            clear(a); clear(n); return 0;
X        }
X    }
X    _mpz_realloc(a,(mp_size)(n->sz));
X    for (k=0; k<s; k++) {
X
X        /* generate a "random" a */
X            for (i=0; i<n->sz; i++) {
X                t = 0;
X                for (j=0; j < DIGITBITS; j+=8)
X                    t = (t << 8) | mpz_randombyte;
X                (a->p)[i] = t & LMAX;
X            }
X            a->sn = 1;
X            mpz_mod(a,a,n);
X
X        if (witness(a,n)) {
X            clear(a); clear(n);
X            return 0;
X        }
X    }
X    clear(a);clear(n);
X    return 1;
X}
X
X
X/* Square roots are found by Newton's method, but since we are working
X * with integers we have to consider carefully the termination conditions.
X * If we are seeking x = floor(sqrt(a)), the iteration is
X *     x_{n+1} = floor ((floor (a/x_n) + x_n)/2) == floor ((a + x_n^2)/(2*x))
X * If eps_n represents the error (exactly, eps_n and sqrt(a) real) in the
X *  form:
X *     x_n = (1 + eps_n) sqrt(a)
X * then it is easy to show that for a >= 4
X *     if 0 <= eps_n, then either 0 <= eps_{n+1} <= (eps_n^2)/2
X *                         or x_{n+1} = floor (sqrt(a))
X *     if x_n = floor (sqrt(a)), then either x_{n+1} = floor (sqrt(a))
X *                               or x_{n+1} = floor (sqrt(a)) + 1
X * Therefore, if the initial approximation x_0 is chosen so that
X *     0 <= eps_0 <= 1/2
X * then the error remains postive and monotonically decreasing with
X *     eps_n <= 2 ^ (-(2^(n+1) - 1))
X * unless we reach the correct floor(sqrt(a)), after which we may oscillate
X * between it and the value one higher.
X * We choose the number of iterations, k, to guarantee
X *     eps_k sqrt(a) < 1,  so x_k <= floor (sqrt (a)) + 1
X * so finally x_k is either the required answer or one too big (which we
X * check by a simple multiplication and comparison).
X *
X * The calculation of the initial approximation *assumes* DIGITBITS is even.
X * It probably always will be, so for now let's leave the code simple,
X * clear and all reachable in testing!
X */
Xvoid mpz_sqrtrem (root, rem, a)
X    MP_INT *root;
X    MP_INT *rem;
X    MP_INT *a;
X{
X    MP_INT tmp;
X    MP_INT r;
X    mp_limb z;
X    unsigned long j, h;
X    int k;
X
X    if (a->sn < 0)
X        /* Negative operand, nothing correct we can do */
X        return;
X
X    /* Now, a good enough approximation to the root is obtained by writing
X     *     a = z * 2^(2j) + y,  4 <= z <= 15 and 0 <= y < 2^(2j)
X     * then taking
X     *     root = ciel (sqrt(z+1)) * 2^j
X     */
X
X    for (j = a->sz - 1; j != 0 && (a->p)[j] == 0; j--);
X    z = (a->p)[j];
X    if (z < 4) {
X        if (j == 0) {
X            /* Special case for small operand (main interation not valid) */
X            mpz_set_ui (root, (z>0) ? 1L : 0L);
X            if (rem)
X                mpz_set_ui (rem, (z>1) ? z-1L : 0L);
X            return;
X        } else {
X            z = (z << 2) + ((a->p)[j-1] >> (DIGITBITS-2));
X            j = j * DIGITBITS - 2;
X        }
X    } else {
X        j = j * DIGITBITS;
X        while (z & ~0xf) {
X            z >>= 2;
X            j += 2;
X        }
X    }
X    j >>= 1;                            /* z and j now as described above */
X    for (k=1, h=4; h < j+3; k++, h*=2);
X        /* 2^(k+1) >= j+3, since a < 2^(2j+4) */
X    mpz_init_set_ui (&r, (z>8) ? 4L : 3L);
X    mpz_mul_2exp (&r, &r, j);
X
X#ifdef DEBUG
X    printf ("mpz_sqrtrem: z = %lu, j = %lu, k = %d\n", (long)z, j, k);
X#endif
X
X    /* Main iteration */
X    mpz_init (&tmp);
X    for ( ; k > 0; k--) {
X        mpz_div (&tmp, a, &r);
X        mpz_add (&r, &tmp, &r);
X        mpz_div_2exp (&r, &r, 1L);
X    }
X
X    /* Adjust result (which may be one too big) and set remainder */
X    mpz_mul (&tmp, &r, &r);
X    if (mpz_cmp (&tmp, a) > 0) {
X        mpz_sub_ui (root, &r, 1L);
X        if (rem) {
X            mpz_sub (rem, a, &tmp);
X            mpz_add (rem, rem, &r);
X            mpz_add (rem, rem, &r);
X            mpz_sub_ui (rem, rem, 1L);
X        }
X    } else {
X        mpz_set (root, &r);
X        if (rem)
X            mpz_sub (rem, a, &tmp);
X    }
X
X    mpz_clear (&tmp);
X    mpz_clear (&r);
X}
X
X
Xvoid mpz_sqrt (root, a)
X    MP_INT *root;
X    MP_INT *a;
X{
X    mpz_sqrtrem (root, NULL, a);
X}
END_OF_FILE
if test 38872 -ne `wc -c <'gmp.c'`; then
    echo shar: \"'gmp.c'\" unpacked with wrong size!
fi
# end of 'gmp.c'
fi
echo shar: End of shell archive.
exit 0
--
Mark Henderson      [e-mail] (personal account)
RIPEM key available by key server/finger/E-mail
  MD5OfPublicKey: F1F5F0C3984CBEAF3889ADAFA2437433

exit 0 # Just in case...