Detailed description of the major CALC functions

R0=2¹⁶
USI stands for unsigned int
USL stands for unsigned long
An MPI is a multiprecision integer
An MPMATI is a matrix of multiprecision integers
An MPR is a multiprecision rational
An MPMATR is a matrix of multiprecision rationals
An MPIA is an array of MPI's
A POLYI is a polynomial with MPI coefficients
X is a reserved symbol, for use in polynomials
Stack is used in STURM and also in the wrappers to functions for later use in init.c. It is only in the latter context that I use Stack.

MPI *ABSI(MPI *Aptr): i5m.c

Returns |Aptr|.

void ADD_TO_MPIA(MPIA MA, MPI *V, USL n): i5I.c

Adds the supplied MPI at the subscript n. If slot already exists, the MPI at that slot is freed and the new one is added. If n is greater than the number of slots, then the array is correctly resized and the new MPI is added. Slots between the previous last slots and the new subscript n are initialised to zero.

MPI *ADD0I(MPI *Aptr, MPI *Bptr): i1.c

Returns Aptr+Bptr, Aptr >= 0, Bptr >= 0.

MPI *ADD0_I(MPI *Aptr, unsigned long b): i1.c

Returns Aptr+b, Aptr >= 0, b < R0.

MPI *ADDI(MPI *Aptr, MPI *Bptr): i1.c

Returns Aptr+Bptr.

MPI *ADDM(MPI *Aptr, MPI *Bptr, MPI *Mptr): i5m.c

Returns Aptr+Bptr (mod Mptr), where 0 <=Aptr,Bptr < Mptr.

MPMATI *ADDMATI(MPMATI *Mptr, MPMATI *Nptr): i6I.c

Returns Aptr+Bptr.

MPR *ADDR(MPR *Aptr, MPR *Bptr): i5R.c

Returns Aptr+Bptr

void ADD_CUBICR(MPR *X1, MPR *Y1, MPR *X2, MPR *Y2, MPR **Xptr, MPR **Yptr, MPR *A1, MPR *A2, MPR *A3, MPR *A4, MPR *A6): cubicr.c

(Xptr,Yptr) is the sum of the two points (X1,Y1) and (X2,Y2) on the elliptic curve y²+A1·xy+A3·y=X³+A2·X²+A4·x+A6.
See D. Husemoller, Elliptic curves, page 25.

void ADD_ELLIPTIC_Q(MPR *X1, MPR *Y1, MPR *X2, MPR *Y2, MPR **Xptr, MPR **Yptr, MPR *A, MPR *B): elliptic.c

(Xptr, Yptr) is the sum of the two points (X1,Y1) and (X2,Y2) on the rational elliptic curve y²=x³+AX+B, where 4A³+27b² is nonzero.

MPMATI *ADD_MULT_ROWI(USI p, USI q, MPI *Aptr, MPMATI *Mptr): i6I.c

Returns the matrix obtained by adding Aptr times row p to row q of Mptr.

MPMATI *ADD_MULT_ROWI0(USI p, USI q, MPI *Aptr, MPMATI *Mptr): i6I.c

Overwrites Mptr by adding Aptr times row p to row q of Mptr.

unsigned long ADDm(USL a, USL b, USL m): i5m.c

Returns a+b(mod m), where 0 <= a,b < m < 2³².

void AXB(): nfunc.c

This solves the linear system AX=B, where the coefficients of A,X,B are integers. A short solution is found in the case of solubility, if N(A) is nontrivial. The method is LLL-based.

MPMATI *BASIS_REDUCTION(MPMATI *Bptr, MPMATI **Eptr, USI rowstage, USI m1, USI n1): LLL.c

Input: Bptr, a matrix of MPI's, whose first row is not zero.
Output: an MPMATI whose rows form a LLL reduced basis for the lattice spanned by the rows of Bptr in the sense of the paper "Factoring polynomials with rational coefficients" by A. K. Lenstra, H. W. Lenstra and L. Lovász, Math. Ann. 261, 515-534 (1982).
We use the modified version in "Solving exponential Diophantine equations using lattice basis reduction algorithms" by B. M. M. De Weger, J. No. Theory 26, 325-367 (1987). A change of basis matrix Eptr is also returned.
De Weger's algorithm has been changed to cater for arbitrary matrices whose rows are now in general linearly dependent.
We use the fact that the Gram Schmidt process detects the first row which is a linear combination of the preceding rows. We employ a modification of the LLL algorithm outlined by M. Pohst in J. Symbolic Computation (1987)4, 123-127.
We call this the MLLL algorithm.
The last sigma rows of the matrix Eptr are relation vectors.
(m1, n1) is usually taken to be (3, 4) for a quick answer, but (1,1), while slower, usually provides shorter basis vectors and multiplier.

MPI *BIG_MTHROOT(MPI *Aptr, unsigned int m): i8.c

The integer part of the mth root of the positive MPI Aptr, 1 < m < R0, is obtained by Newton's method, using the integer part function. (See the article by [Matthews].)

unsigned int BINARYB(MPI *N): binary.c

Returns the number of binary digits of N.

MPI *BINOMIAL(USI n, USI m): i5m.c

returns n choose m, where n >= m are unsigned integers.

MPI *BRENT_POLLARD(MPI *Nptr): primes1.c

The Brent-Pollard method returns a proper factor of a composite MPI Nptr. (see R. Brent, BIT 20, 176 - 184).

MPMATI *BUILDMATI(unsigned int m, unsigned int n): i6I.c

Allocates space for an m x n matrix of MPI's.

MPMATR *BUILDMATR(unsigned int m, unsigned int n): i6R.c

Allocates space for an m x n matrix of MPR's.

MPI *BUILDMPI(unsigned int n): i5I.c

Mallocs space for an MPI of length n.
If there is an MPI of this size in the bank, then use it rather than malloc.

MPIA BUILDMPIA(): i5I.c

Allocates space for an array initially of size 11 (enough to hold a[0] to a[10] and sets these slots to contain the zero MPI. Extra MPI's are added using ADD_TO_MPIA.

MPR *BUILDMPR( ): i5R.c

Mallocs space for an MPR.

MPI *CHANGE(unsigned long n): i5I.c

Converts n, 0 <= n < (R0)² to an MPI.

MPI *CHANGEI(long n): i5I.c

Converts n, 0 <= |n| < R0 to an MPI.

MPI *CHINESE(MPI *A, MPI *B, MPI *M, MPI *N, MPI **Mptr): nfunc.c

Returns the solution mod Mptr=lcm[M,N] and Mptr) of the simultaneous congruences X = A (mod M) and X = B (mod N), if soluble; otherwise returns NULL.

MPI *CHINESE_ARRAY(MPI *A[ ], MPI *M[ ], MPI **Mptr, USI n): nfunc.c

Returns the solution mod Mptr=lcm[M[0],...,M[n-1] and Mptr) of the congruences X = A[i] (mod M[i]),0 <= i < n, if soluble; otherwise returns NULL.

MPMATR *CHOLESKYR(MPMATR *A): i6R.c

Input: The positive definite matrix A.
Output: The Cholesky decomposition of A.
(See U. Finke and M. Pohst, "Improved methods for calculating vectors of short length in a lattice, including a complexity analysis." Math. Comp. 44, 463-471, 1985.

MPI *COLLATZ(MPI *Dptr, *Eptr): nfunc.c

The Collatz 3x+1 function. The iterates x,T(x),.. are printed iff Eptr is nonzero.

MPMATI *COLSUBI(USI p, USI q, MPI *Aptr, MPMATI *Mptr): i6I.c

Returns the result of subtracting Aptr times the p-th column of Mptr from the q-th.
0 <= p, q <= Mprt->C - 1.

MPI *COLSUMI(MPMATI *Mptr, USI j): i9.c

Returns the sum of the elements of column j of Mptr.

int COMPAREI(MPI *Aptr, MPI *Bptr): i1.c

Compares MPI's: Returns 1 if Aptr > Bptr, 0 if Aptr = Bptr, -1 if Aptr < Bptr.

int COMPARER(MPR *Aptr, MPR *Bptr): i6R.c

Compares MPR's: Returns 1 if Aptr > Bptr, 0 if Aptr = Bptr, -1 if Aptr < Bptr.

MPI *CONGR(MPI *A, MPI *B, MPI *M, MPI **N): nfunc.c

Returns the least solution (mod N) of the congruence AX=B(mod M) (and N), where N = M / gcd(A, M); otherwise returns the null pointer.

MPI *CONTENTPI(POLYI Pptr): pI.c

Cptr is the content of the polynomial Pptr.

void CONVERGENTS(MPI *A[], MPI **P[], MPI **Q[], MPI *N): i5I.c

Returns the convergents P[0]/Q[0],...,P[N]/Q[N] to [A[0];A[1],...,A[n]] as arrays P[ ] and Q[ ].

unsigned long CONVERTI(MPI *N): i5I.c

Returns N as an unsigned long, providing 0 < N < 2³².

MPI *COPYI(MPI *Aptr): i5I.c

Returns a copy of the MPI Aptr.

MPR *COPYR(MPR *Aptr): i5R.c

Returns a copy of the MPR Aptr.

MPMATI *COPYMATI(MPMATI *Aptr): i6I.c

Returns a copy of the MPMATI Aptr.

MPMATR *COPYMATR(MPMATR *Aptr): i6R.c

Returns a copy of the MPMATR Aptr.

VOID CORNACCHIA(MPI *A, MPI *B, MPI *M): primes1.c

This prints the positive primitive solutions (x,y) of Ax²+By²=M, where A,B,M are positive integers, with gcd(A,M)=1=gcd(A,B).

void CYCLE(USL d, MPI *m[ ], MPI *X[ ], USL INFINITY, USL RANGE): collatz.c

This function searches all trajectories of the d-branched generalized Collatz function, which start from p, |p| <= RANGE/2 (RANGE an even integer). INFINITY is an upper bound for the size of a trajectory, above which the trajecory is deemed to be divergent. Floyd's cycle finding algorithm is used. Also see survey (pdf) by the author.

MPI *DIVM(MPI *Aptr, MPI *Bptr, MPI *Mptr): i5m.c

Returns (Aptr / Bptr) mod (Mptr).
Here 0 <= Aptr, Bptr < Mptr and gcd(Bptr, Mptr) = 1.

unsigned long DIVm(USL a, USL b, USL m): i5m.c

Returns a / b mod m if m > 0
Here 0 <= a, b < m < R0, gcd(b, m) = 1 if m > 0>

MPI *DOTRI(MPMATI *Mptr, USI i, USI j): i7I.c

Returns the dot product of rows i and j in Mptr.

MPI *EFACTOR(MPI *N, USI m, USI p): elliptic.c

The elliptic curve method is used to try to find a factor of a composite number N.
Here m, p < 2³² and 1279 >= m > 10, p >= 1.

unsigned int EQUALI (MPI *Aptr, MPI *Bptr): i5I.c

Returns 1 if Aptr = Bptr, otherwise = 0.

unsigned int EQUALR (MPR *Aptr, MPR *Bptr): i5R.c

Returns 1 if Aptr = Bptr, otherwise = 0.

unsigned int EQMINUSONECI(MPCI *Mptr): i3.c

Returns 1 if Mptr = -1, 0 otherwise.

unsigned int EQMINUSONECR(MPCR *Mptr): i3.c

Returns 1 if Mptr = -1, 0 otherwise.

unsigned int EQMINUSONEI(MPI *Mptr): i3.c

Returns 1 if Mptr = -1, 0 otherwise.

unsigned int EQMINUSONER(MPR *Mptr): i3.c

Returns 1 if Mptr = -1, 0 otherwise.

unsigned int EQONEI(MPI *Mptr): i3.c

Returns 1 if Mptr = 1, 0 otherwise.

unsigned int EQONER(MPR *Mptr): i3.c

Returns 1 if Mptr = 1, 0 otherwise.

unsigned int EQONECI(MPCI *Mptr): i3.c

Returns 1 if Mptr = 1, 0 otherwise.

unsigned int EQONECR(MPCR *Mptr): i3.c

Returns 1 if Mptr = 1, 0 otherwise.

unsigned int EQZEROI(MPI *Mptr): i3.c

Returns 1 if Mptr = 0, but 0 otherwise.

unsigned int EQZEROR(MPR *Mptr): i3.c

Returns 1 if Mptr = 0, but 0 otherwise.

unsigned int EQZEROCI(MPCI *Mptr): i3.c

Returns 1 if Mptr = 0, but 0 otherwise.

unsigned int EQZEROCR(MPCR *Mptr): i3.c

Returns 1 if Mptr = 0, but 0 otherwise.

MPI *EUCLIDI(MPI *Pptr, MPI *Qptr, MPI **Hptr, MPI **Kptr): nfunc.c

Returns gcd(Pptr, Qptr) = Hptr · Pptr + Kptr · Qptr.

void EUCLID(MPI *Aptr, MPI *Bptr, MPI **Q[ ], MPI **R[ ], MPI **S[ ], MPI **T[ ], MPI **Dptr)

Returns Q[0]=NULL,Q[1],...Q[n],Q[n+1]=NULL,
* R[0],...R[n + 1],
* S[0],...S[n + 1], T[0],...T[n + 1]
* for Euclid's algorithm on R[0]=Aptr, R[1]=Bptr.
* R[0]=R[1]*Q[1]+R[2]
* R[1]=R[2]*Q[2]+R[3]
* .....
* R[n-2]=R[n-1]*Q[n-1]+R[n]
* R[n-1]=R[n]*Q[n], R[n+1]=0.
* S[0]=1,S[1]=0, S[j]=s[j-2]-Q[j-1]*S[j-1],
* T[0]=0,T[1]=1, T[j]=T[j-2]-Q[j-1]*T[j-1], j=2,...,n+1
* Here *Dptr = n.

MPI *EULER(MPI *N): primes1.c

Returns Euler's function phi(N).

FACTOR(Mptr *Nptr): primes1.c

Attempts to factor Nptr using the Multiple Polynomial Quadratic sieve.
If that fails, it uses the elliptic curve method. (No Peter Montgomery fine-tuning.)
I suggest one uses it for 55 or less digit numbers.

unsigned long *FFI(USL N, USL *b, USL w, USL p): i5m.c

Fast Fourier Interpolation.
Here N is a power of 2, b is an array of N elements from Z_p,
w is a primitive N-th root of unity mod p and N | p - 1.
Outputs a(x)=a[0]x⁰+···+a[0]x^N-1, where a(w^k)=b[k], 0 <= k < N.
(See "Elements of Algebra and Algebraic Computing", J.D. Lipson, p.303.
I have been unable to free the arrays B and C, due to the way the recursion works.

unsigned long *FFP(USL N, USL *a, USL *b, USL m, USL n, USL w, USL p): i5m.c

Input: arrays a and b of USL's mod p, representing polynomials
of degrees m, n, respectively; N = 2^e > m + n, N | p - 1. w is a primitive N-th root of unity mod p.
Output: array c mod p, representing a(x)b(x).

MPI *CRA(USL n, USL *a, USL *m): i5m.c

Garner's Chinese Remainder Theorem.
Page 180, Algorithms for computer algebra, K.O. Geddes, S.R. Czapor, G. Labahn.
Here gcd(m[i],m[j])=1 if i != j.
Returns the least remainder mod(m[0].....m[n]) of u = a[i]mod(m[i]), 0 < =i < = n.

MPI *FFM(MPI *a, MPI *b, USL K): i5m.c

Returns the product of a=(a[0],...,a[m])_B and b=(b[0],...,b[n])_B, B=2¹⁶,
m = a->D, n = b->D, using the Discrete Fast Fourier Transform.
Let M = min(m,n). Then a(x)b(x)=c(x), where 0 <= c[k] < (M+1)B².
Using the CRA mod fp[i] for 0 <= i <= K-1, enables us to reconstruct c(B),
provided that fp[0]···fp[K-1] >= (M+1)B².
We also need m+n < N = 2^e, where N | fp[i] - 1.
If m <B and n < B, then M < B, then K=3 primes suffice as fp[i]>=B
and we take e>=17.
If m < 2²⁶ and n < 2²⁶, then M < 2²⁶, then K=4 primes suffice and we take e > = 27, N = 2^27.
External variables:
fp[0] = 2013265921, lprimroot[0] = 31, primitive Nth root = 440564289;
fp[1] = 2281701377, lprimroot[1] = 3, primitive Nth root = 129140163;
fp[2] = 3221225473, lprimroot[2] = 5, primitive Nth root = 229807484;
fp[3] = 3489660929, lprimroot[3] = 3, primitive Nth root = 1392672017.
See Elements of Algebra and Algebraic Computing, J.D. Lipson, p.310.
I do not use FFM, as my implementation seems slow. If the user wants to invoke it, go to i1.c and uncomment the relevant parts of MULTI( ).

MPI *FIBONACCI(USI n): functions.c

Returns the nth Fibonacci number.

void FINCKE_POHST(MPMATR *A, MPR *C): i6R.c

Input: A matrix of integers A with LI rows spanning a lattice L.
Output: The integer vectors X with ||X||² <= C, highest nonzero coord <= 0.
(See Improved methods for calculating vectors of short length
in a lattice including a complexity analysis,
U. Fincke and M. Pohst, Mathematics of Computation, 44, 1985, 463-471.

MPI *FINPUTI(FILE *f, unsigned int *uptr): i5I.c

Converts decimal input from stream into an MPI Mptr.
Ignores the combination '\' followed by '\n'.
If a rubbish character is met before decimal input, Mptr is set to 0
and 0 is returned. All characters up to and including the first newline
met are wiped.
If a rubbish character is met immediately after decimal input,
uptr = 0 is returned and all characters up to and including
the first newline met are wiped. Otherwise 1 is returned.
In any case Mptr is set equal to any inputted decimal.

MPMATI *FINPUTMATI(FILE *infile): i6I.c

Inputs a matrix of MPI's from infile.

MPR *FINPUTR(FILE *f, unsigned int *uptr): i5R.c

Converts the ratio of two decimal inputs from stream into an MPR.
uptr = 0 if input fails, 1 if successful.

void FPRINTI(FILE *outfile, MPI *Mptr): i5I.c

The MPI Mptr is printed in decimal notation to outfile.
No new-line is incorporated.

void FPRINTMATI(FILE *outfile, USI i1, USI i2, USI j1, USI j2, MPMATI *Mptr): i6I.c

Printing an MPMATI to outfile.

void FPRINTMATR(FILE *outfile, USI i1, USI i2, USI j1, USI j2, MPMATR *Mptr): i6I.c

Printing an MPMATR to outfile.

void FPRINTR(FILE *outfile, MPR *Aptr): i5R.c

prints the MPR Aptr as (Aptr->N)/(Aptr->D).

MPR *FRAC_PARTI(MPI *Aptr, MPI *Bptr): i5R.c

Returns the fractional part of Aptr/Bptr.

MPR *FRAC_PARTR(MPR *Aptr): i5R.c

Returns the fractional part of Aptr.

void FREEMATI(MPMATI *Mptr): i6I.c

Frees the storage allocated to the 2-dimensional array Mptr->V.

void FREEMATR(MPMATR *Mptr): i6R.c

Frees the storage allocated to the 2-dimensional array Mptr->V.

void FREEMPI(MPI *Mptr): i5I.c

Deallocates the space alloted for the MPI Mptr.

void FREEMPIA(MPI *Mptr): i5I.c

Frees an MPIA previously returned by BUILDMPIA. It will free all the the MPI's in the MPIA.

void FREEMPR(MPR *Mptr): i5R.c

Deallocates the space alloted for the MPR Mptr.

MPI *FUND_UNIT(MPI *D, MPI **Xptr, MPI **Yptr): nfunc.c

This is a program for finding the fundamental unit of Q(sqrt(D)).
The algorithm is based on:
K. Rosen, Elementary number theory and its applications, p382,
B.A. Venkov, Elementary Number theory, p.62
D. Knuth, Art of computer programming, Vol.2, p359,
with Pohst's trick of using half the period.
w=(1+sqrt(D))/2 if D=1 (mod 4), w=sqrt(D) otherwise.
The norm of the fundamental unit Xptr + Yptr·w is returned.

MPI *GCD(MPI *Aptr, MPI *Bptr): i5I.c

Returns GCD(|Aptr|, |Bptr|).

MPI *GCD_ARRAY(MPI *M[ ], unsigned int n): i5I.c

Returns GCD(M[0], ..., M[n - 1]).

MPI *GCD_ARRAYV(MPI *M[ ], MPI **Y[ ], USI n): nfunc.c

Returns d=gcd(M[0],...,M[n-1]) and an array Y[ ] of MPI's
such that d = M[0]Y[0]+···+M[n-1]Y[n-1]. Here n > 1.

unsigned long GCDm(USL m, USL n): i5m.c

Returns gcd(m,n)>

void GetReturn( ): menu.c

Waits for a return to be entered from the keyboard.

unsigned int GetYN( ): menu.c

Gets a character from the keyboard, making sure it's a
y or an n (either case). If at first the user doesn't succeed,
he/she tries, tries again. 0 is returned if n, 1 if y.

MPMATI *HERMITE1(MPMATI *Aptr, USI *T, USI nz): i6I.c

(Kannan-Bachem) Returns the Hermite normal form of Aptr.

MPMATI *HERMITE1P(MPMATI *Aptr, MPMATI *Pptr, MPMATI **Qptr, USI *T, USI nz): i6I.c

(Kannan-Bachem) Returns the Hermite normal form of Aptr and a transforming unimodular matrix Pptr.

MPMATI *IDENTITYI(USI n): i6I.c

Returns the identity matrix of size n.

MPI *INPUTI(unsigned int *uptr): i5I.c

Inputs an MPI from the keyboard.

MPR *INPUTR(unsigned int *uptr): i5I.c

Inputs an MPR from the keyboard.
uptr=1 if no corruption takes place, else 0.

MPMATI *INPUTMATI( ): i6I.c

Inputs a matrix of MPI's from stdin.

MPI *INT0(MPI *Aptr, MPI *Bptr): i2.c

Returns the integer part of Aptr/Bptr,
where Aptr,Bptr > 0.

MPI *INT0_(MPI *Aptr, unsigned long b): i2.c

Returns the integer part of Aptr/b,
where Aptr > 0 and 0 < b < R0.

MPI *INT(MPI *Aptr, MPI *Bptr): i2.c

Returns the integer part of Aptr/Bptr,
where Bptr > 0.

MPI *INT_(MPI *Aptr, USL b): i2.c

Returns the integer part of Aptr/b, 0 < b < R0.

MPI *INTI(MPI *Aptr, MPI *Bptr): i2.c

Returns the integer part of Aptr/Bptr.

MPR *INTR(MPR *Aptr): i5R.c

Returns the integer part of Aptr.

MPI *INVERSEM(MPI *Nptr, MPI *Mptr): i5m.c

Returns the inverse of Nptr (mod Mptr).
Here gcd(Nptr,Mptr) = 1, 1 <= Nptr < Mptr.
(See Knuth, p. 325.)

MPR *INVERSER(MPR *Aptr): i5R.c

Returns 1/ Aptr.

unsigned long INVERSEm(USL n, USL m): i5m.c

Returns the inverse of n (mod m).
Here 1 <= n < m < 2³², gcd(n, m) = 1.

MPI *JACOB(MPI *M, MPI *N): func.c

Returns the Jacobi symbol (M/N), N odd, N > 0.

int JACOBI(USL n, USL m): qres.c

Returns the Jacobi symbol (n/m), m odd, 0 < n < m.

MPMATI *JACOBIGCD(MPMATI *DD, MPI **Aptr, USI m): LLLGCD.c

Input: an m x 1 vector DD of positive MPI's.
Output: Aptr = gcd of the DD[i]. Also we return a set of
multipliers using a version of a method of Jacobi
A unimodular transforming matrix B is returned.

void LAGRANGE(POLY P, **AA[], MPI *M): i5I.c

f(x)=a[n]xⁿ+···+a[0], a[n] > 0, is a polynomial with integer coefficients, having no rational roots and having exactly one real positive root x, this being > 1. The method of Lagrange (1769) is used to find the the first m+1 partial quotients aa[0],···aa[m] of x. WARNING: the array a[] will be changed after lagrange is called.
Then a further call to lagrange will produce subsequent partial quotients. (See Knuth, Art of computer programming, volume 2, problem 13, 4.5.3.
Also S. Lang and H. Trotter,Continued fractions for some algebraic numbers J. für Math. 255 (1972) 112-134; Addendum 267 (1974) ibid. 219-220.
E. Bombieri and A. van der Poorten, Continued fractions of algebraic numbers, Computational algebra and number theory (Sydney, 1992), 137-152, Math. Appl., 325, Kluwer Acad. Publ.
P. Shiu, Computation of continued fractions without input values, Math. Comp. 64 (1995), no. 211, 1307-1317.

MPI *LCM(MPI *Aptr, MPI *Bptr): i5I.c

Returns lcm(Aptr,Bptr).

MPI *LCM_ARRAY(MPIA M): i5I.c

Returns lcm(M[0],...,M[n - 1]).

MPI *LEASTQNR(MPI *P): qres.c

Returns the least quadratic non-residue mod P.

unsigned long LENGTHI(MPI *Mptr): i5I.c

Returns the number of decimal digits in the MPI Mptr
increased by 1 if Mptr is negative.

unsigned int LENGTHm(USL n): i5m.c

Returns the number of decimal digits in the unsigned int n.

MPI *LENGTHSQRI(MPMATI *Mptr, USI i): LLL.c

Returns the square of the length of row i of matrix Mptr.

MPMATI *LLLGCD(MPMATI *DD, MPI **Aptr, USI m, USI m1, USI n1): LLLGCD.c

Input: an m x 1 vector of MPI's.
Output: gcd of the vector of DD[i]. We return a small set of
multipliers using the LLL method of Havas, Majewski and Matthews.
matrix B of the algorithm is returned.
(m1, n1) is usually taken to be (3, 4) for a quick answer,
but (1,1), while slower, usually provides a shorter basis vectors.

MPI *LPRIMROOT(MPI *P): primes1.c

Returns the least primitive root mod P, an odd prime;
returns NULL if factorization of P - 1 is unsuccessful.

MPI *LUCAS(MPI *N)

Here N is odd, N > 1.
If LUCAS(N) returns 1, then N is a strong base 2 pseudoprime
and a Lucas probable prime; if it returns 0, then N is composite.
(See The Pseudoprimes to 25·10⁹,
Mathematics of computation, 35 (1980) 1003-1026.
At the end of this paper it is conjectured that if N is a strong
base 2 pseudoprime and a Lucas probable prime, then N is in fact a prime.

int MAX(int i, int j): i1.c

int MIN(int i, int j): i1.c

MPI *MINUSI(MPI *Aptr): i5I.c

Returns -Aptr.

MPR *MINUSI(MPR *Aptr): i5R.c

Returns -Aptr.

MPI *MINUSM(MPI *Aptr, MPI *Mptr): i5m.c

Returns -Aptr (mod Mptr).
Here 0 <= Aptr < Mptr.

unsigned long MINUSm(USL a, USL m): i5m.c

Returns -a (mod m) if m > 0.
Here 0 <= a < m < 2³².

MPI *MINUS_ONEI( ): i5I.c

Returns -1 as an MPI.

MPR *MINUS_ONER( ): i5R.c

Returns -1 as an MPR.

MPI *MOBIUS(MPI *N): primes1.c

Returns the Mobius function mu(N),
returns NULL if factorization unsuccessful.

MPI *MOD(MPI *Aptr, MPI *Bptr): i2.c

Returns Aptr (mod Bptr), Bptr > 0.

MPI *MOD0(MPI *Aptr, MPI *Bptr): i2.c

Returns Aptr (mod Bptr), Aptr >= 0, Bptr > 0.

unsigned long MOD0_(MPI *Aptr, unsigned long b): i2.c

Returns Aptr (mod b), Aptr >= 0, 0 < b < R0.

unsigned long MOD_(MPI *Aptr, unsigned long b): i2.c

Returns Aptr (mod b), 0 < b < R0.

unsigned long MODINT0_(MPI *Aptr, unsigned long b, MPI **Qptr): i2.c

Returns Aptr (mod b) and Qptr = int(Aptr/b).
Here Aptr >= 0, b is a positive integer, b < R0.

MPI *MPOWER(MPI *Aptr, MPI *Bptr, MPI *Cptr): i5m.c

Returns (Aptr)^Bptr (mod Cptr).
Here Cptr > 0, Bptr >= 0.
We use an analogue of the Russian Peasant Multiplication algorithm and conserve the quantity w=zx^y(mod c).
Initially z=1,x=1,y=b.
If y is odd, y -> y-1, z -> z*x(mod c);
if y is even, y -> y/2, x -> x²(mod c).
Eventually y=0. Then w=zx⁰(mod c)=z and the final value of z gives a^b(mod c).

MPI *MPOWER_(long a, MPI *Bptr, MPI *Cptr) i5m.c

Returns a^Bptr (mod Cptr).
Here 0 < |a| < R0, Cptr > 0, Bptr >= 0.

MPI *MPOWER_M(MPI *Aptr, USL b, MPI *Cptr): i5m.c

Returns (Aptr)^b (mod Cptr).
Here Cptr > 0, 0 <= b < Ro.

void MTHROOT(MPI *Aptr, MPI *Bptr, unsigned int m, unsigned int r): i8.c

Aptr and Bptr are positive MPI'S. 0 < mr < R0².
The mthroot of Aptr/Bptr is computed to r decimal places, r >= 0.

MPR *MTHROOTR(MPR *Nptr, unsigned int m, unsigned int r): i8.c

The m-throot of the positive MPR Nptr is computed to r decimal places,
m, r are integers, r >= 0, 0 < mr < R0².

MPI *MULTABC(MPI *A, MPI *B, MPI *C): i1.c

Returns A + BC.

MPR *MULTABCR(MPR *A, MPR *B, MPR *C): i1.c

Returns A + BC.

MPI *MULTI(MPI *Aptr, MPI *Bptr): i1.c

Returns Aptr·Bptr;.

MPI *MULTR(MPR *Aptr, MPR *Bptr): i5R.c

Returns Aptr·Bptr;.

MPI *MULTI3(MPI *A, MPI *B, MPI *C): i1.c

Returns ABC.

MPR *MULTR3(MPR *A, MPR *B, MPR *C): cubicr.c

Returns ABC. Here 0 <= Aptr, Bptr < Mptr.

MPMATI *MULTMATI(MPMATI *Mptr, MPMATI *Nptr): i6I.c

Returns Mptr·Bptr.

MPMATR *MULTMATR(MPMATR *Mptr, MPMATR *Nptr): i6R.c

Returns Mptr·Bptr.

MPI *MULT_I(MPI *Aptr, long b): i1.c

Returns Aptr·b, where |b| < R0.

MPI *MULT_II(MPI *Aptr, USL b): i1.c

Returns Aptr·b, where b is an USL.

MPI *MULTM(MPI *Aptr, MPI *Bptr, MPI *Mptr): i5m.c

Returns Aptr·Bptr (mod Mptr).

unsigned long MULTm(USL a, USL b, USL m): i5m.c

Returns ab (mod m) if m > 0; here 0 <= a,b < m < 2³².

MPI *NEAREST_INTI(MPI *Aptr, MPI *Bptr): i5I.c

Returns the nearest integer to Aptr/Bptr,
choosing downwards if half an odd integer.

MPI *NEAREST_INTR(MPR *Aptr): i5R.c

Returns the nearest integer to Aptr,
choosing downwards if half an odd integer.

MPI *NEXTPRIMEAP(MPI *A, MPI *B, MPI *M): utility.c

Finds the first Lucas probable prime P, A | P - B, M <= P.
Here A is even, B is odd, A > 1 , A > B >= 1, gcd(A,B) = 1, M > 1.

MPI *NEXT_PRIME(MPI *M, USI *hptr): utility.c

Returns a probable prime Q with Q = M + hptr.

MPI *ONEI( ): i5I.c

Returns 1.

MPR *ONER( ): i5R.c

Returns 1.

unsigned int ORDERECP(MPI *X, MPI *Z, MPI *P, MPI *Q, MPI *N): elliptic.c

Calculates the order of the point (X,Y,Z) on the elliptic curve
Y²Z=X³+PXZ²+QZ³ (mod N), N a prime.

MPI *ORDERM(MPI *A, MPI *M): primes1.c

Returns the order of A (mod M).

MPI *ORDERP(MPI *A, MPI *P): primes1.c

Returns the order of A (mod P), P a prime.

void PATZ(MPI *D, MPI *N): lagrange.c

Returns the fundamental solutions (x,y) (with x > 0) of x²-Dy²=± N, where D > 1 is not a perfect square and N is non-zero.
The algorithm goes back to Lagrange and has been strangely forgotten by textbook writers. (See the preprint (pdf 173K) of Keith Matthews.)

MPI *PEL(MPI *D, MPI*E, MPI **Xptr, MPI **Yptr): nfunc.c

This finds the least solution of Pell's equation x² - Dy² = ±1.
The algorithm is based on K. Rosen,
Elementary number theory and its applications, p382,
B.A. Venkov, Elementary Number theory, p.62
and D. Knuth, Art of computer programming, Vol.2, p359,
with Pohst's trick of using half the period.
The partial quotients are printed iff E is nonzero.
The norm of the least solution is returned.

void PELL(MPI *Dptr, MPI *Eptr): nfunc.c

This finds the period of the regular continued fraction
of square-root d, as well as the least solution x,y
of the Pellian equation x²-dy²=±1.
The algorithm is from Sierpinski's Theory of Numbers, p.296.
and Davenport's The Higher Arithmetic, p.109.
Here sqrt(d)=a[0]+1/a[1]+···+1/a[n-1]+1/2*a[0]+1/···.
The partial quotients are printed iff Eptr is nonzero.
The length n of the period a[1],···,a[n-1],2*a[0] is printed.

MPI *PERFECT_POWER(MPI *N): primes1.c

If N > 1, this returns X if N=X^k for some X, k > 1, otherwise NULL.

MPI *POLLARD(MPI *Nptr): primes1.c

Pollard's p-1 method of factoring Nptr.

MPI *POWERI(MPI *Aptr, unsigned n): i1.c

Returns (Aptr)ⁿ, where 0 <= n < R0².

MPR *POWERR(MPR *Aptr, unsigned n): i5R.c

Returns (Aptr)ⁿ, where 0 <= n < R0².

void POWER_CUBICR(MPR *X1, MPR *Y1, MPR **Xptr, MPR **Yptr, MPR *A1, MPR *A2, MPR *A3, MPR *A4, MPR *A6, unsigned int n): cubicr.c

(Xptr,Yptr)= n(X1,Y1), where 0 <= n < R0² and (X1, Y1)
is on the elliptic curve y²+A1·xy+A3·y=X³+A2·X²+A4·x+A6.
(See D. Husemoller, Elliptic curves, page 25.)

MPI *POWER_I(long a, unsigned n): i1.c

Returns aⁿ, where 0 <=n < R0².
a is an integer, |a| < R0.

unsigned long POWER_m(USL a, USL y, USL m): i5m.c

Returns a^y (mod m).
Here 0 <= a < m < R0 and 0 <= y.

unsigned long POWERm(USL a, MPI *Bptr, USL m): i5m.c

Returns a^Bptr (mod m).
Here 0 <= a < m < R0 and 0 <= Bptr.

POLYI PRIMITIVEPI(POLYI Pptr): pI.c

Returns the primitive part of the polynomial Pptr.

void PRINTI(MPI *Mptr): i5I.c

Prints the MPI Mptr.

void PRINTR(MPR *Mptr): i5R.c

Prints the MPR Mptr.

void PRINTMATI(USI i1, USI i2, USI j1, USI j2, MPMATI *Mptr): i6I.c

Prints the matrix Mptr, rows i1-i2, cols j1-j2.

void PRINTMATR(USI i1, USI i2, USI j1, USI j2, MPMATR *Mptr): i5R.c

Prints the matrix Mptr, rows i1-i2, cols j1-j2.

unsigned long RANDOMm(USL x): i5m.c

input: seed x, output: a "random number" a x + c (mod m).
a = 1001, m = R0 = 65536, c = 65.
From H. Lüneburg, On the Rational Normal Form of Endomorphisms,
B.I. WissenSchaftsverlag, Mannheim/Wien/Zurich, 1987.
See Knuth Vol 2, Theorem A, p. 16.

MPR *RATIOI(MPI *Aptr, MPI *Bptr): i5R.c

Returns Aptr/Bptr.

MPR *RATIOR(MPR *Aptr, MPR *Bptr): i5R.c

Returns Aptr/Bptr.

void readme(): readme.c

Contains the readme manual for CALC.

MPR *RECIPROCAL(unsigned long n): i5R.c

Returns 1/n, where 0 < n < R0.

void ROOTEXPANSION(Stack s): roots.c

Here POLYI P = stackPop(s), MPI *M = stackPop(s). This finds the first M partial quotients of the continued fraction expansion of all real roots of the supplied polynomial P using Lagrange's method and methods presented in a paper by Cantor, Galyean and Zimmer.

MPMATI *ROWSUBI(USI p, USI q, MPI *Aptr, MPMATI *Mptr): i6I.c

Updates Mptr: subtracts Aptr x the p-th row of Mptr from the q-th.
0 <= p, q < Mprt->R.

MPI *ROWSUMI(MPMATI *Mptr, USI i): i9.c

Returns the sum of the elements in row i of Mptr.
0 <= i < Mptr->R.

int RSV(MPI *Aptr, MPI *Bptr): i1.c

Returns 1 if |Aptr| > |Bptr|,
0 if |Aptr| = |Bptr|, -1 if |Aptr| < |Bptr|.

unsigned long RUSSm(USL a, USL b, USL c, USL p): i5m.c

input: unsigned long ints a, b, c, p, with p > 0.
output: a + bc (mod p).
Russian Peasant multiplication algorithm. Uses the identities
RUSSm(a, b, 2c, p) = RUSSm(a, 2b, c, p),
RUSSm(a, b, c + 1, p) = RUSSm(a + b, b, c, p).
If a, b, c and p are less than 2³², so is RUSSm(a, b, c, p).
From H. Lüneburg, On the Rational Normal Form of Endomorphisms,
B.I. WissenSchaftsverlag, Mannheim/Wien/Zurich, 1987, pp 18-19.
Lüneburg's restriction to 2p<2³² removed by me on 18/4/94.

void SCHNORRGCD(MPI *N): nfunc.c

Applies LLL to [I_m|ND] with N large
to get small multipliers for the extended gcd problem.
See [Matthews].)

void SERRET(MPI *P, MPI **Xptr, MPI **Yptr): nfunc.c

This program finds positive integers X, Y such that
X²+Y²=P, where P=4n=1 is a prime.
The algorithm goes back to Serret and is from the book
Computational methods in number theory, part 1,
edited by H.W. Lenstra and R. Tijdeman.

MPI *SIGMA(MPI *N): primes1.c

Returns sigma(N), the sum of the divisors of N,
Returns NULL if factorization unsuccessful.

MPR *SLVECTOR(MPMATR *A, MPR *C, MPR **VALUE): i6R.c

Input: A matrix of integers A with LI LLL reduced rows spanning a lattice L. C = length-squared of a small lattice vector. Output: if 0 is returned, this means that VALUE is the shortest length. Otherwise a shorter length is returned and VALUE will be NULL.
Used in nfunc.c, in SLVECTORX(), starting with a LLL-reduced matrix and C=||A₁||², eventually 0 is returned. Then FINCKE-POHST is applied to get all the shortest vectors in L with highest nonzero coord > 0 are printed.
(See "Improved methods for calculating vectors of short length in a lattice, including a complexity analysis, U. Fincke and M. Pohst, Mathematics of Computation, 44, 1985, 463-471.)

MPI **SMITHI(MPMATI *Mptr, MPMATI **Pptr, MPMATI **Qptr, USI *ptr, MPI *Eptr): i9.c

Returns the invariant factors of Mptr.
Pptr and Qptr are invertible matrices such that Pptr·Mptr·Qptr
is the Smith normal form Nptr. ptr is the number of invariant factors.
See Rings, Modules and Linear Algebra, B. Hartley and T.O. Hawkes,
Chapman and Hall, 1970. We use a pivoting strategy suggested by George Havas. Eptr is the cutoff above which we bring in MLLL.

MPI **SMITHI1(MPMATI *Mptr, USI *ptr, MPI *Eptr): i9.c

Same as SMITHI, but with no transforming matrices.

void SPRINTI(char *buffer, MPI *Mptr): i5I.c

The decimal digits of the MPI Mptr are placed in the string buffer.
No new-line is incorporated.

void SPRINTR(char *buffer, MPR *Aptr): i5R.c

Prints the MPR Aptr as (Aptr->N)/(Aptr->D).

MPI *SQROOT1(MPI *A, MPI *P, USL n): primes1.c

This returns a solution of x² A (mod Pⁿ), where P is an odd prime not dividing A. Returns -1 otherwise.
Let d=gcd(a,pⁿ).
Case 1: d=1. Use the standard recursive approach - two solutions if soluble.
                If P=2, care is needed, as in the case of solubility, (i) there is 1 solution if n=1, (ii) 2 solutions if n=2 and (iii) 4 solutions if n 3.
Case 2: d=Pⁿ. Then x=0 is the solution.
Case 3: d=P^r, 1 < r < n. Write a=P^rA, where P does not divide A.

1. If r is odd, there is no solution.
2. If r is even, say r=2R. Put x=P^RX.
   Then X² A (mod P^n-2R), which reduces to Case 1.

MPI *SQROOT2(MPI *A, USL n): primes1.c

This returns a solution of x² A (mod 2ⁿ), A odd. Returns -1 otherwise. Also returns a global variable sqroot2_number=1,2,4 if n=1,2,or n > 2. This variable is used in SQROOT

MPI *SQROOT3(MPI *A, MPI *P, USL n, MPI**EXPONENT): primes1.c

This returns a solution of x² A (mod Pⁿ), where P is an odd prime possibly dividing A. Returns -1 otherwise. Also returns a "reduced moduli" explained as follows:
If P does not divide A, the story is well-known.
If P divides A, there are two cases:
(i) Pⁿ divides A,
(ii) Pⁿ does not divide A.
In case (i) there are two cases:
(a) n=2m+1, (b) n=2m.
In case (a), the solution is x 0 (mod P^(m+1)).
In case (b), the solution is x 0 (mod P^m).
(These are called reduced moduli and are returned as EXPONENT. This variable is used in SQROOT.)
In case (i) gcd(A,Pⁿ)=P^r, r < n. If r is odd, no solution.
If r=2m, write x=(P^m)*X, A=(P^2m)*A1, P not dividing A1.
Then x² A (mod Pⁿ) becomes X² A1 (mod P^n-2m).
If P is odd, this has 2 solutions X mod P^n-2m and we get two solutions x (mod P^n-m). If P=2, we get 1 solution mod (2^n-m) if n-2m=1, 2 solutions mod (2^n-m) if n-2m=2, 4 solutions mod (2^n-m) if n-2m > 2.

MPI *SQROOT(MPI *A, MPI *N, MPIA *SOL, MPI **MODULUS, USI *lptr): primes1.c

This finds the solutions of x²A (mod N) as residue classes ±SOL[0],...,±SOL[lptr-1] (mod MODULUS), where 0 <= x < MODULUS/2. The number of solutions mod N is returned. If there is no solution, -1 is returned, SOL[0]=NULL and lptr=0.

MPI *SQRTM(MPI *x, MPI *p): qres.c

Calculates sqrt(x) (mod p) using A simple and fast probabilistic
algorithm for computing square roots modulo a prime number,
I.E.E.E. Trans. Inform. Theory, IT-32, 1986, 846-847, R. Peralta.
Here x is a quadratic residue mod p. x can be negative.

unsigned long SQRTm(USL x, USL p): qres.c

Returns sqrt(x) (mod p), as above. Here 1 <=x < p.
x is a quadratic residue mod p.

Stack STURM(POLYI P): roots.c

Returns a stack of intervals about each of the real roots of P. Each interval contains only one root. If there are no roots, it returns an empty stack. P is supposed to have no rational roots, but is not necessarily squarefree.

MPI *SUB0I(MPI *Aptr, MPI *Bptr): i1.c

Returns Aptr-Bptr, where Aptr >= Bptr >= 0.

MPI *SUB0_I(MPI *Aptr, unsigned long b): i1.c

Returns Aptr-b, where Aptr >= b>= 0.

MPI *SUBI(MPI *Aptr, MPI *Bptr): i1.c

Returns Aptr-Bptr.

MPR *SUBR(MPR *Aptr, MPR *Bptr): i5R.c

Returns Aptr-Bptr.

MPI *SUBM(MPI *Aptr, MPI *Bptr, MPI *Mptr): i5m.c

returns Aptr - Bptr (mod Mptr).
Here 0 <= Aptr, Bptr < Mptr.

unsigned long SUBm(USL a, USL b, USL m): i5m.c

Returns a - b (mod m) if m > 0; here 0 <= a, b < m < 2³².

unsigned int SURD(MPI *D, MPI *T, MPI *U, MPI *V): nfunc.c

This function uses the continued fraction algorithm expansion
in K. Rosen, Elementary Number theory and its applications,
p.379-381 and Knuth's The art of computer programming,
Vol.2, p. 359. It locates the first complete quotient that is reduced
and then uses the function REDUCED(D,U,V,i) above to locate and return
the period of the continued fraction expansion of (U+T*sqrt(D))/V.

MPMATI *SWAP_COLSI(USI p, USI q, MPMATI *Mptr): i6I.c

Interchanges cols p and q (C notation) of Mptr.

MPMATI *SWAP_COLSI1(USI p, USI q, MPMATI *Mptr): i6I.c

Interchanges cols p and q (C notation) of Mptr.
Updates Mptr.

MPMATI *SWAP_ROWSSI(USI p, USI q, MPMATI *Mptr): i6I.c

Interchanges rows p and q (C notation) of Mptr.

MPMATI *SWAP_ROWSSI1(USI p, USI q, MPMATI *Mptr): i6I.c

Interchanges rows p and q (C notation) of Mptr.

void SelOpt( ): menu.c

This function simply prints the "SELECT OPTION: " prompt.

MPMATR *TRANSPOSER(MPMATR *Mptr): i6R.c

Returns the transpose of Mptr.

MPI *TWOI( ): i5I.c

Returns 2.

void TryAgain( ): menu.c

Prints "Try again"

MPI *ZEROI( ): i5I.c

Returns 0.

MPR *ZEROR( ): i5R.c

Returns 0.

MPMATI *ZEROMNI(USI m, USI n): i6I.c

Returns the zero max n matrix.

MPMATR *ZEROMNR(USI m, USI n): i6R.c

Returns the zero max n matrix.

void init( ): init.c

Install built-in function names in table.

void initv( ): init.c

Install built-in void function names in table.

Last modified 6th July 2000