NAME
mpsetminbits, mpnew, mpfree, mpbits, mpnorm, mpcopy, mpassign, mprand, mpnrand, strtomp, mpfmt, mptoa, betomp, mptobe, mptober, letomp, mptole, mptolel, mptoui, uitomp, mptoi, itomp, uvtomp, mptouv, vtomp, mptov, mptod, dtomp, mpdigdiv, mpadd, mpsub, mpleft, mpright, mpmul, mpexp, mpmod, mpmodadd, mpmodsub, mpmodmul, mpdiv, mpcmp, mpsel, mpfactorial, mpextendedgcd, mpinvert, mpsignif, mplowbits0, mpvecdigmuladd, mpvecdigmulsub, mpvecadd, mpvecsub, mpveccmp, mpvecmul, mpmagcmp, mpmagadd, mpmagsub, crtpre, crtin, crtout, crtprefree, crtresfree – extended precision arithmetic

SYNOPSIS
#include <u.h>
#include <libc.h>
#include <mp.h>

mpint*    mpnew(int n)

void      mpfree(mpint *b)

void      mpsetminbits(int n)

void      mpbits(mpint *b, int n)

mpint*    mpnorm(mpint *b)

mpint*    mpcopy(mpint *b)

void      mpassign(mpint *old, mpint *new)

mpint*    mprand(int bits, void (*gen)(uchar*, int), mpint *b)

mpint*    mpnrand(mpint *n, void (*gen)(uchar*, int), mpint *b)

mpint*    strtomp(char *buf, char **rptr, int base, mpint *b)

char*     mptoa(mpint *b, int base, char *buf, int blen)

int       mpfmt(Fmt*)

mpint*    betomp(uchar *buf, uint blen, mpint *b)

int       mptobe(mpint *b, uchar *buf, uint blen, uchar **bufp)

void      mptober(mpint *b, uchar *buf, int blen)

mpint*    letomp(uchar *buf, uint blen, mpint *b)

int       mptole(mpint *b, uchar *buf, uint blen, uchar **bufp)

void      mptolel(mpint *b, uchar *buf, int blen)

uint      mptoui(mpint*)

mpint*    uitomp(uint, mpint*)

int       mptoi(mpint*)

mpint*    itomp(int, mpint*)

mpint*    vtomp(vlong, mpint*)

vlong     mptov(mpint*)

mpint*    uvtomp(uvlong, mpint*)

uvlong    mptouv(mpint*)

mpint*    dtomp(double, mpint*)

double    mptod(mpint*)

void      mpadd(mpint *b1, mpint *b2, mpint *sum)

void      mpmagadd(mpint *b1, mpint *b2, mpint *sum)

void      mpsub(mpint *b1, mpint *b2, mpint *diff)

void      mpmagsub(mpint *b1, mpint *b2, mpint *diff)

void      mpleft(mpint *b, int shift, mpint *res)

void      mpright(mpint *b, int shift, mpint *res)

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      mptrunc(mpint *b, int n, mpint *res)

void      mpxtend(mpint *b, int n, mpint *res)

void      mpasr(mpint *b, int n, mpint *res)

void      mpmul(mpint *b1, mpint *b2, mpint *prod)

void      mpexp(mpint *b, mpint *e, mpint *m, mpint *res)

void      mpmod(mpint *b, mpint *m, mpint *remainder)

void      mpdiv(mpint *dividend, mpint *divisor,    mpint *quotient,
mpint *remainder)

void      mpmodadd(mpint *b1, mpint *b2, mpint *m, mpint *sum)

void      mpmodsub(mpint *b1, mpint *b2, mpint *m, mpint *diff)

void      mpmodmul(mpint *b1, mpint *b2, mpint *m, mpint *prod)

int       mpcmp(mpint *b1, mpint *b2)

int       mpmagcmp(mpint *b1, mpint *b2)

void      mpsel(int s, mpint *b1, mpint *b2, mpint *res)

mpint*    mpfactorial(ulong n)

void      mpextendedgcd(mpint *a, mpint *b, mpint *d, mpint *x,
mpint *y)

void      mpinvert(mpint *b, mpint *m, mpint *res)

int       mpsignif(mpint *b)

int       mplowbits0(mpint *b)

void      mpdigdiv(mpdigit *dividend, mpdigit divisor,
mpdigit *quotient)

void      mpvecadd(mpdigit *a, int alen, mpdigit *b, int blen,
mpdigit *sum)

void      mpvecsub(mpdigit *a, int alen, mpdigit *b, int blen,
mpdigit *diff)

void      mpvecdigmuladd(mpdigit *b, int n, mpdigit m, mpdigit *p)

int       mpvecdigmulsub(mpdigit *b, int n, mpdigit m, mpdigit *p)

void      mpvecmul(mpdigit *a, int alen, mpdigit *b, int blen,
mpdigit *p)

int       mpveccmp(mpdigit *a, int alen, mpdigit *b, int blen)

CRTpre* crtpre(int nfactors, mpint **factors)

CRTres* crtin(CRTpre *crt, mpint *x)

void      crtout(CRTpre *crt, CRTres *r, mpint *x)

void      crtprefree(CRTpre *cre)

void      crtresfree(CRTres *res)

mpint     *mpzero, *mpone, *mptwo

DESCRIPTION
These routines perform extended precision integer arithmetic. The basic type is mpint, which points to an array of mpdigits, stored in little–endian order:
typedef struct mpint mpint;
struct mpint
{
int    sign;     /* +1 or –1 */
int    size;     /* allocated digits */
int    top;      /* significant digits */
mpdigit     *p;
char flags;
};

The sign of 0 is +1.

The size of mpdigit is architecture–dependent and defined in /$cputype/include/u.h. Mpints are dynamically allocated and must be explicitly freed. Operations grow the array of digits as needed.

In general, the result parameters are last in the argument list.

Routines that return an mpint will allocate the mpint if the result parameter is nil. This includes strtomp, itomp, uitomp, btomp, and dtomp. These functions, in addition to mpnew and mpcopy, will call sysfatal (see perror(2)) if the allocation fails.

Input and result parameters may point to the same mpint. The routines check and copy where necessary.

Mpnew creates an mpint with an initial allocation of n bits. If n is zero, the allocation will be whatever was specified in the last call to mpsetminbits or to the initial value, 1056. Mpfree frees an mpint. Mpbits grows the allocation of b to fit at least n bits. If b–>top doesn't cover n bits, mpbits increases it to do so. Unless you are writing new basic operations, you can restrict yourself to mpnew(0) and mpfree(b).

Mpnorm normalizes the representation by trimming any high order zero digits. All routines except mpbits return normalized results.

Mpcopy creates a new mpint with the same value as b while mpassign sets the value of new to be that of old.

Mprand creates an n bit random number using the generator gen. Gen takes a pointer to a string of uchar's and the number to fill in.

Mpnrand uses gen to generate a uniform random number x, 0 ≤ x < n.

Strtomp and mptoa convert between ASCII and mpint representations using the base indicated. Only the bases 2, 4, 8, 10, 16, 32, and 64 are supported. Strtomp skips any leading spaces or tabs. Strtomp's scan stops when encountering a digit not valid in the base. If base is zero then C–style prefixes are interpreted to find the base: 0x for hexadecimal, 0b for binary and 0 for octal. Otherwise decimal is assumed. rptr is not zero, *rptr is set to point to the character immediately after the string converted. If the parse terminates before any digits are found, strtomp return nil. Mptoa returns a pointer to the ASCII filled buffer. If the parameter buf is nil, the buffer is allocated. Setting base to zero uses hexadecimal default. Mpfmt can be used with fmtinstall(2) and print(2) to print ASCII representations of mpints. The conventional verb is B, for which mp.h provides a pragma. The precision in the format string changes the base, defaulting to hexadecimal when omited.

Mptobe and mptole convert an mpint to a byte array. The former creates a big endian representation, the latter a little endian one. If the destination buf is not nil, it specifies the buffer of length blen for the result. If the representation is less than blen bytes, the rest of the buffer is zero filled. If buf is nil, then a buffer is allocated and a pointer to it is deposited in the location pointed to by bufp. Sign is ignored in these conversions, i.e., the byte array version is always positive.

Mptober and mptolel fill blen lower bytes of an mpint into a fixed length byte array. Mptober fills the bytes right adjusted in big endian order so that the least significant byte is at buf[blen–1] while mptolel fills in little endian order; left adjusted; so that the least significat byte is filled into buf[0].

Betomp, and letomp convert from a big or little endian byte array at buf of length blen to an mpint. If b is not nil, it refers to a preallocated mpint for the result. If b is nil, a new integer is allocated and returned as the result.

The integer (and floating point) conversions are:
mptoui    mpint–>unsigned int
uitomp    unsigned int–>mpint
mptoi     mpint–>int
itomp     int–>mpint
mptouv    mpint–>unsigned vlong
uvtomp    unsigned vlong–>mpint
mptov     mpint–>vlong
vtomp     vlong–>mpint
mptod     mpint–>double
dtomp     double–>mpint

When converting to the base integer types, if the integer is too large, the largest integer of the appropriate sign and size is returned.

When converting to and from floating point, results are rounded using IEEE 754 "round to nearest". If the integer is too large in magnitude, mptod returns infinity of the appropriate sign.

The mathematical functions are:
mpadd         sum = b1 + b2.
mpmagadd      sum = abs(b1) + abs(b2).
mpsub         diff = b1 – b2.
mpmagsub      diff = abs(b1) – abs(b2).
mpleft         res = b<<shift.
mpright        res = b>>shift.
mpmul        prod = b1*b2.
mpexp         if m is nil, res = b**e. Otherwise, res = b**e mod m.
mpmod        remainder = b % m.
mpdiv         quotient = dividend/divisor. remainder = dividend % divisor.
mpcmp         returns –1, 0, or +1 as b1 is less than, equal to, or greater than b2.
mpmagcmp      the same as mpcmp but ignores the sign and just compares magnitudes.
mpsel         assigns b1 to res when s is not zero, otherwise b2 is assigned to res.
mpfactorial      returns n!.

Logical operations (treating negative numbers using two's complement):
mpand      res = b1 & b2.
mpbic      res = b1 & ~b2.
mpor       res = b1 | b2.
mpxor      res = b1 ^ b2.
mpnot      res = ~b1.
mpasr      res = b>>shift (mpasr, unlike mpright, uses two's complement).
mptrunc     truncates b to n bits and stores the result in res. The result is never negative.
mpxtend     truncates b to n bits, sign extends the MSB and stores the result in res.

Modular arithmetic:
mpmodadd     sum = b1+b2 mod m.
mpmodsub     diff = b1–b2 mod m.
mpmodmul     prod = b1*b2 mod m.

Mpextendedgcd computes the greatest common denominator, d, of a and b. It also computes x and y such that a*x + b*y = d. Both a and b are required to be positive. If called with negative arguments, it will return a gcd of 0.

Mpinvert computes the multiplicative inverse of b mod m.

Mpsignif returns the number of significant bits in b. Mplowbits0 returns the number of consecutive zero bits at the low end of the significant bits. For example, for 0x14, mpsignif returns 5 and mplowbits0 returns 2. For 0, mpsignif and mplowbits0 both return 0.

The remaining routines all work on arrays of mpdigit rather than mpint's. They are the basis of all the other routines. They are separated out to allow them to be rewritten in assembler for each architecture. There is also a portable C version for each one.
mpdigdiv          quotient = dividend[0:1] / divisor.
mpvecadd          sum[0:alen] = a[0:alen–1] + b[0:blen–1]. We assume alen >= blen and that sum has room for alen+1 digits.
mpvecsub          diff[0:alen–1] = a[0:alen–1] – b[0:blen–1]. We assume that alen >= blen and that diff has room for alen digits.
mpvecdigmuladd     p[0:n] += m * b[0:n–1]. This multiplies a an array of digits times a scalar and adds it to another array. We assume p has room for n+1 digits.
mpvecdigmulsub     p[0:n] –= m * b[0:n–1]. This multiplies a an array of digits times a scalar and subtracts it from another array. We assume p has room for n+1 digits. It returns +1 is the result is positive and –1 if negative.
mpvecmul         p[0:alen+blen] = a[0:alen–1] * b[0:blen–1]. We assume that p has room for alen+blen+1 digits.
mpveccmp         This returns –1, 0, or +1 as a – b is negative, 0, or positive.

mptwo, mpone and mpzero are the constants 2, 1 and 0. These cannot be freed.

Time invariant computation

In the field of cryptography, it is sometimes neccesary to implement algorithms such that the runtime of the algorithm is not depdenent on the input data. This library provides partial support for time invariant computation with the MPtimesafe flag that can be set on input or destination operands to request timing safe operation. The result of a timing safe operation will also have the MPtimesafe flag set and is not normalized.

Chinese remainder theorem

When computing in a non–prime modulus, n, it is possible to perform the computations on the residues modulo the prime factors of n instead. Since these numbers are smaller, multiplication and exponentiation can be much faster.

Crtin computes the residues of x and returns them in a newly allocated structure:
typedef struct CRTres      CRTres;
{
int    n;     /* number of residues */
mpint       *r[n];      /* residues */
};

Crtout takes a residue representation of a number and converts it back into the number. It also frees the residue structure.

Crepre saves a copy of the factors and precomputes the constants necessary for converting the residue form back into a number modulo the product of the factors. It returns a newly allocated structure containing values.

Crtprefree and crtresfree free CRTpre and CRTres structures respectively.

SOURCE
/sys/src/libmp