ref: 74761f99b3ab80b2df581c1c141d28f9681f6c3d
dir: /posix/mp.h/
#pragma once
typedef uint32_t mpdigit;
typedef union FPdbleword FPdbleword;
union FPdbleword
{
double x;
struct {
#if BYTE_ORDER == LITTLE_ENDIAN
uint32_t lo;
uint32_t hi;
#else
uint32_t hi;
uint32_t lo;
#endif
};
};
#define sysfatal(s) do{ fprintf(stderr, "%s\n", s); abort(); }while(0)
#define mpdighi (mpdigit)(1U<<(Dbits-1))
#define DIGITS(x) ((Dbits - 1 + (x))/Dbits)
// for converting between int's and mpint's
#define MAXUINT ((uint32_t)-1)
#define MAXINT (MAXUINT>>1)
#define MININT (MAXINT+1)
// for converting between vlongs's and mpint's
#define MAXUVLONG (~0ULL)
#define MAXVLONG (MAXUVLONG>>1)
#define MINVLONG (MAXVLONG+1ULL)
extern int dec16(uint8_t*, int, char*, int);
extern int enc16(char*, int, uint8_t*, int);
extern uint32_t dec16chr(int);
extern int enc16chr(int);
/*
* the code assumes mpdigit to be at least an int
* mpdigit must be an atomic type. mpdigit is defined
* in the architecture specific u.h
*/
typedef struct mpint mpint;
struct mpint
{
int sign; /* +1 or -1 */
uint32_t size; /* allocated digits */
uint32_t top; /* significant digits */
mpdigit *p;
char flags;
};
enum
{
MPstatic= 0x01, /* static constant */
MPnorm= 0x02, /* normalization status */
MPtimesafe= 0x04, /* request time invariant computation */
MPfield= 0x08, /* this mpint is a field modulus */
Dbytes= sizeof(mpdigit), /* bytes per digit */
Dbits= Dbytes*8 /* bits per digit */
};
/* allocation */
void mpsetminbits(int n); /* newly created mpint's get at least n bits */
mpint* mpnew(int n); /* create a new mpint with at least n bits */
void mpfree(mpint *b);
void mpbits(mpint *b, uint32_t n); /* ensure that b has at least n bits */
mpint* mpnorm(mpint *b); /* dump leading zeros */
mpint* mpcopy(mpint *b);
void mpassign(mpint *old, mpint *new);
/* random bits */
mpint* mprand(int bits, void (*gen)(uint8_t*, int), mpint *b);
/* return uniform random [0..n-1] */
mpint* mpnrand(mpint *n, void (*gen)(uint8_t*, int), mpint *b);
/* conversion */
mpint* strtomp(const char*, char**, int, mpint*); /* ascii */
char* mptoa(mpint*, int, char*, int);
mpint* letomp(uint8_t*, uint32_t, mpint*); /* byte array, little-endian */
int mptole(mpint*, uint8_t*, uint32_t, uint8_t**);
void mptolel(mpint *b, uint8_t *p, int n);
mpint* betomp(uint8_t*, uint32_t, mpint*); /* byte array, big-endian */
int mptobe(mpint*, uint8_t*, uint32_t, uint8_t**);
void mptober(mpint *b, uint8_t *p, int n);
uint32_t mptoui(mpint*); /* unsigned int */
mpint* uitomp(uint32_t, mpint*);
int mptoi(mpint*); /* int */
mpint* itomp(int, mpint*);
uint64_t mptouv(mpint*); /* unsigned int64_t */
mpint* uvtomp(uint64_t, mpint*);
int64_t mptov(mpint*); /* int64_t */
mpint* vtomp(int64_t, mpint*);
double mptod(mpint*); /* double */
mpint* dtomp(double, mpint*);
/* divide 2 digits by one */
void mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);
/* in the following, the result mpint may be */
/* the same as one of the inputs. */
void mpadd(mpint *b1, mpint *b2, mpint *sum); /* sum = b1+b2 */
void mpsub(mpint *b1, mpint *b2, mpint *diff); /* diff = b1-b2 */
void mpleft(mpint *b, int shift, mpint *res); /* res = b<<shift */
void mpright(mpint *b, int shift, mpint *res); /* res = b>>shift */
void mpmul(mpint *b1, mpint *b2, mpint *prod); /* prod = b1*b2 */
void mpexp(mpint *b, mpint *e, mpint *m, mpint *res); /* res = b**e mod m */
void mpmod(mpint *b, mpint *m, mpint *remainder); /* remainder = b mod m */
/* logical operations */
void mpand(mpint *b1, mpint *b2, mpint *res);
void mpbic(mpint *b1, mpint *b2, mpint *res);
void mpor(mpint *b1, mpint *b2, mpint *res);
void mpnot(mpint *b, mpint *res);
void mpxor(mpint *b1, mpint *b2, mpint *res);
void mpasr(mpint *b, int shift, mpint *res);
/* modular arithmetic, time invariant when 0≤b1≤m-1 and 0≤b2≤m-1 */
void mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum); /* sum = b1+b2 % m */
void mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff); /* diff = b1-b2 % m */
void mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod); /* prod = b1*b2 % m */
/* quotient = dividend/divisor, remainder = dividend % divisor */
void mpdiv(mpint *dividend, mpint *divisor, mpint *quotient, mpint *remainder);
/* return neg, 0, pos as b1-b2 is neg, 0, pos */
int mpcmp(mpint *b1, mpint *b2);
/* res = s != 0 ? b1 : b2 */
void mpsel(int s, mpint *b1, mpint *b2, mpint *res);
/* extended gcd return d, x, and y, s.t. d = gcd(a,b) and ax+by = d */
void mpextendedgcd(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y);
/* res = b**-1 mod m */
void mpinvert(mpint *b, mpint *m, mpint *res);
/* bit counting */
uint32_t mpsignif(mpint*); /* number of sigificant bits in mantissa */
uint32_t mplowbits0(mpint*); /* k, where n = 2**k * q for odd q */
/* well known constants */
extern mpint *mpzero, *mpone, *mptwo;
/* sum[0:alen] = a[0:alen-1] + b[0:blen-1] */
/* prereq: alen >= blen, sum has room for alen+1 digits */
void mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *sum);
/* diff[0:alen-1] = a[0:alen-1] - b[0:blen-1] */
/* prereq: alen >= blen, diff has room for alen digits */
void mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *diff);
/* p[0:n] += m * b[0:n-1] */
/* prereq: p has room for n+1 digits */
void mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p);
/* p[0:n] -= m * b[0:n-1] */
/* prereq: p has room for n+1 digits */
int mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p);
/* p[0:alen+blen-1] = a[0:alen-1] * b[0:blen-1] */
/* prereq: alen >= blen, p has room for m*n digits */
void mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
void mpvectsmul(mpdigit *a, int alen, mpdigit *b, int blen, mpdigit *p);
/* sign of a - b or zero if the same */
int mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen);
int mpvectscmp(mpdigit *a, int alen, mpdigit *b, int blen);
/* divide the 2 digit dividend by the one digit divisor and stick in quotient */
/* we assume that the result is one digit - overflow is all 1's */
void mpdigdiv(mpdigit *dividend, mpdigit divisor, mpdigit *quotient);
/* playing with magnitudes */
int mpmagcmp(mpint *b1, mpint *b2);
void mpmagadd(mpint *b1, mpint *b2, mpint *sum); /* sum = b1+b2 */
void mpmagsub(mpint *b1, mpint *b2, mpint *sum); /* sum = b1+b2 */
/* chinese remainder theorem */
typedef struct CRTpre CRTpre; /* precomputed values for converting */
/* twixt residues and mpint */
typedef struct CRTres CRTres; /* residue form of an mpint */
struct CRTres
{
int n; /* number of residues */
mpint *r[1]; /* residues */
};
CRTpre* crtpre(int, mpint**); /* precompute conversion values */
CRTres* crtin(CRTpre*, mpint*); /* convert mpint to residues */
void crtout(CRTpre*, CRTres*, mpint*); /* convert residues to mpint */
void crtprefree(CRTpre*);
void crtresfree(CRTres*);
/* fast field arithmetic */
typedef struct Mfield Mfield;
struct Mfield
{
mpint m;
int (*reduce)(Mfield*, mpint*, mpint*);
};
mpint *mpfield(mpint*);
Mfield *gmfield(mpint*);
Mfield *cnfield(mpint*);