CALC
CALC is a number theory calculator program which uses arbitrary precision integer arithmetic. It would be useful in a first course in number theory.
It is written in ANSI C and Yacc, along the lines of the calculator programs
hoc1,2,3, Kernighan and Pike, The Unix Programming Environment,
233254, 1984.
The source comes with two makefiles, one for Unix and one for Windows XP.
The Windows XP version and C source can be downloaded here.
I am grateful to JeanJacques Delgove for making changes to my code, so as to allow compilation under Visual C++ 2005. He also pointed out that CALC now can utilise AWK scripts such as the following, invoked by typing awk f calc.awk in a Unix window:
==== calc.awk ===========
BEGIN {
for(i=1;i < 8;i++){
cmd=sprintf("./calc \"(X+1)^%d\"",i)
while (((cmd)  getline a) > 0){
print a
}
close(cmd)
}
}
with output:
X + 1
X^2 + 2X + 1
X^3 + 3X^2 + 3X + 1
X^4 + 4X^3 + 6X^2 + 4X + 1
X^5 + 5X^4 + 10X^3 + 10X^2 + 5X + 1
X^6 + 6X^5 + 15X^4 + 20X^3 + 15X^2 + 6X + 1
X^7 + 7X^6 + 21X^5 + 35X^4 + 35X^3 + 21X^2 + 7X + 1
Also see AWK (Wikipedia).
To run CALC, type calc. One can then perform standard arithmetical calculations such as the following (the answer is written underneath the corresponding expression in red):
> 54321+12345
66666
> 5432112345
41976
> 54321*12345
670592745
> 54321^10
223705964292712579004593047803650384459784511201
> 54321/12345
4
> 54321%12345
4941
> x=12;y=13
> x*y
156
Arrays can be entered in two main ways:
(i) with an array identifer and a subscript:
> a[0] = 1;a[1] = 2;a[2] = 3;a[3] = 4;a[4] = 5
(ii) or in the following manner:
> a[ ] = {1, 2, 3, 4, 5, 6}
Note: Any previous values for the array will be erased and a new array will
be created with the values in the curly braces.
When one wishes to output the array a[ ], all that must be typed is:
> a[ ]
alone on a line.
If one uses the former method, the size of the array will be set to the
largest subscript entered. Any subscripts that have not been defined
and which are less than this subscript, are initialised to zero. eg.
> a[1]=2
> a[2]=3
> a[4]=10
> a[ ]
[0]:0
[1]:2
[2]:3
[3]:0
[4]:10
Some functions require an mxn integer matrix to be entered, either at the keyboard, or from a file. The format of such a file is as follows:
m n
a_{11} ... a_{1n}
...
a_{m1} ... a_{mn}
In other words, the first line contains the row and column size, separated by a space, while the entries on each row are separated by spaces and each row is terminated by a newline.
Calc is now capable of parsing polynomials. A polynomial can
be entered using the reserved symbol X (capital x) in the following manner:

> (X+2)^20

X^20 + 40X^19 + 760X^18 + 9120X^17 + 77520X^16 + 496128X^15 + 2480640X^14 + 9922560X^13 + 32248320X^12 + 85995520X^11 + 189190144X^10 + 343982080X^9 + 515973120X^8 + 635043840X^7 + 635043840X^6 + 508035072X^5 + 317521920X^4 + 149422080X^3 + 49807360X^2 + 10485760X + 1048576
Once a polynomial has been assigned to a variable, it is also possible
to evaluate that polynomial at a specified integer value.

> z=(X+2)^20

> z(2)

1099511627776
Polynomial division and remaindering are available as p/q and p%q.
Here p is initially multiplied by c^{deg(p)deg(q)} (if deg(p)deg(q) is not negative), where c is the leading coefficient of q.

> p=X^2+3X+7; q=2X+5

> p/q

> 2X+1

> p%q

> 23
Some more examples
To find z = gcd(4,6), type
> z = gcd(4,6)
Then z is stored and its value is printed:
> z
2
To find z = gcd(4,6) together with integers u,v satisfying z = 4u+6v,
type
>z = gcdv(4,6,&u,&v)
The values of u,v are stored and their values printed by typing
>u
1
>v
1
To find z=gcd(4,6,9), type
>x[0]=4;x[1]=6;x[2]=9
>z=gcda(x[ ])
To find z=gcd(4,6,9), together with integers b[0],b[1],b[2] satisfying
z=4*b[0]+6*b[1]+9*b[2], type
>x[0]=4;x[1]=6;x[2]=9
>z=gcdav(x[ ],&b[ ])
The values of b[0],b[1],b[2] are stored and their values printed by typing:
>b[ ]
[0]: 4
[1]: 4
[2]: 1
LIST OF FUNCTIONS
 z=absmod(a,b)
 euclid(a,b,&q[ ],&r[ ],&s[ ],&t[ ],&n)

z=gcd(x,y)

z=gcdv(x,y,&u,&v)

z=gcda(a[ ])

z=gcdav(a[ ], &b[ ])

egcd( )

sgcd(N)

lllgcd( )

lllgcd0( )

jacobigcd( )

z=lcm(x,y)

z=lcma(a[ ])

z=length(n)

z=pollard(x)

z=nprime(x)

z=nprimeap(a,b,m)

z=jacobi(x,y)

z=peralta(a,p)

x=congr(a,b,m,&n)

x=chinese(a,b,m,n,&l)

x=chinesea(a[ ],m[ ],&l)

z=mthroot(x,m)

mthrootr(x,y,m,r)

z=fund(d,&x,&y)

z=pell(d,e,&x,&y)

z=mpower(a,b,c)

z=inv(a,m)

z=elliptic(n,m,p)

z=factor(n)

z=tau(n)

z=sigma(n)

z=mobius(n)

z=euler(n)

z=lprimroot(n)

z=orderm(a,n)

z=lucas(n)

serret(p,&x,&y)

collatz(x,n)

miller(m,b)

lllhermite( )

shermite(N)

smith( )

mlll( )

fp( )

slv( )

cycle( )

inhomfp( )

e=rsae(p,q)

encode(e,n)

decode(e,p,q)

addcubicr( )

powercubicr( )

ordercubicr( )

convergents(a[ ],&p[ ],&q[ ])

lagrange(f(X),&a[ ],m)

z=perfectpower(n)

axb( )

addcubicm( )

powercubicm( )

ordercubicm( )

leastqnr(p)
 sturm(f(X))

rootexp(f(X), m)

content(f(X))
 primitive(f(X))
 sqroot(a,n,&s[ ],&m)
 cornacchia(a,b,m)
 z=surd(d,t,u,v,&a[ ],&u[ ],&v[ ],&p[ ],&q[ ])
 patz(d,n)
 congq(a,b,c,n,&s[ ])
 binform(a,b,c,n,e)
 ceil(a,b)
 log(a,b,d,&a[ ],&l)
 logx(a,b,m,n)
 r=resultant(p(X),q(X))
 r=discriminant(p(X))
 q=deriv(p(X))
 c=primes(m,n)
 c=sturmsequence(f(X),b,e)
 p=cyclotomic(n)
 h=classnop(d)
 h=classnon(d)
 z=nearint(a,b)
 h=reduceneg(a,b,c)
 h=reducepos(a,b,c)
 h=classnop0(d)
 h=tableneg(m,n)
 h=tablepos(m,n)
 h=davison(l,m,n)
 h=raney(p,q,r,s)
 h=unimodular(p,q,r,s)
 twoadicsqrt(b,n,&a[])
 padicsqrt(b,n,p,&a[])
 ramanujan(n)
 repdefinite(a,b,c,m,print_flag)
 powerd(a,b,d,n,&aa,&bb)
 z=euclid1(a,b)
 z=rcfperiod(d)
 z=nscfperiod(d)
 z=nicfperiod(d)
 algcycle()
 spiral(n,&x,&y)
 n=spiralinverse(x,y)
 z=rcfperiod0(d,e,&x,&y)
 z=nscfperiod0(d,e,&x,&y)
 z=nicfperiod0(d,e,&x,&y)
 carmichael(n)
 carnielli()
 lupei()
 tangent(n)
 tangent(n,&x,&y)
 partition(n)
 twocycle( )
 mcycle( )
 pqcycle( )
 h=fg(p,q)
 nagell(d,n)
 frattini(d,n)
 branch(n,type,flag)
 pell4(d,e,&x,&y)
 stolt(a,b,c,n)
 z=lprimefactor(n)
 z=conwaycycles(a,b)
 conwayrangetest1(m1,m2,n1,n2)
euclid(a,b,&q[],&r[],&s[],&t[],&m)
Here m=n+1 where we return four arrays arising from Euclid's algorithm are returned:
q[0]=NULL,...,q[n],q[n+1]=NULL (quotients)
r[0]=a,r[1]=b,...,r[n+1] (remainders)
s[0]=a,s[1]=0,...,s[n+1]
t[0]=a,t[1]=1,...,t[n+1]
r[k]=r[k+1]*q[k+1]+r[k+2], 0 < r[k+2] < r[k+1],
s[k]=q[k1]*s[k1]+s[k2],
t[k]=q[k1]*t[k1]+t[k2],
r[n]=gcd(a,b)=s[n]*a+t[n]*b.
These are printed in euclid.out.
z=gcd(x,y)
This returns the gcd of x and y.
z=gcdv(x,y,&u,&v)
As well as returning z=gcd(x,y), it gives numbers u and v arising from
Euclid's algorithm and satisfying the equation z = ux+vy.
z=gcda(a[])
z = gcd(a[0],...,a[n1]), where values for a[0],...,a[n1] having been
previously entered.
z=gcdav(a[],&b[])
z = gcd(a[0],...,a[n1]). Also gives integers b[0],...,b[n1] satisfying
z = b[0]a[0]+...+b[n1]a[n1].
egcd()
This is an implementation of Algorithm 1 of a recent paper by Havas, Majewski and Matthews.
Like lllgcd( ), it finds short multipliers for the gcd of m numbers,
using LLL ideas. It also finds all shortest vectors, unlike lllgcd( ),
which lists only one shortest multiplier. The file of m integers should
have as its first line m, then the integers should be listed on separate
lines.
An m x m matrix whose rows are X[1],...,X[m1],P is sent to an output
file egcdmat.out.
Here X[1],...,X[m1] form a LLL reduced basis for the lattice L defined
by the equation x_{1}d_{1}+ ··· +x_{m}d_{m}
= 0.
The multipliers are sent to an output file called egcdmult.out.
Sometimes the multipliers delivered by egcd() are shorter that those
of lllgcd( ). There is also an option to find all the shortest multipliers.
sgcd(N)
This performs the LLL algorithm on [I_{n}NA], where A is a
column vector of positive integers. This is Algorithm 2 of paper.
If N is sufficiently large, the last column will be reduced to ±NdE_{n},
where d = gcd(a[1],...,a[n]). Output is sent to sgcdbas.out.
lllgcd()
This performs a modification of LLL which is essentially a limiting
form of sgcd(N) for large N. It is superior to egcd() in that it avoids
inputting a large initial unimodular matrix and instead builds one from
the identity matrix at the outset. This is Algorithm 3 of a recent paper of Havas, Majewski and Matthews.
The file of m integers should have as its first line m, then the integers
should be listed on separate lines.
An m x m matrix whose rows are X[1],...,X[m] is sent to an output file
lllgcdmat.out.
Here X[1],...,X[m1] form a LLL reduced basis for the lattice L defined
by the equation x_{1}d_{1}+ ··· +x_{m}d_{m}
= 0.
The inhomogeneous version of the FinckePohst algorithm (see [Po2][191])
can then be used as an option to find a shortest multiplier vector by solving
the inequality
X[m]  x_{1}X[1] ··· x_{m1}X[m1]^{2}
≤ X[m]^{2}
in integers x_{1},...,x_{m1}.
Each time a shorter multiplier vector Q = X[m]x_{1}X[1] ··· x_{m1}X[m1]
is found, X[m] is replaced by Q, until the shortest Q is found. The multipliers
are sent to an output file called lllgcdmult.out. The unimodular matrix
P of a recent paper of Havas, Majewski and Matthews, is sent to lllgcdmat.out. In verbose mode, the intermediate steps are printed.
Note: if the shortest vector option is chosen, the last row of P has
been replaced by this vector.
lllgcd0()
This in general gives a better multiplier than lllgcd and is based
on an algorithm in a recent manuscript of
the author.
jacobigcd()
This performs Jacobi's extended gcd algorithm of 1869. In verbose mode
the intermediate steps are printed out.
z=lcm(x,y)
z=lcma(x)
z = lcm(x[0],...,x[n1]). (n is the size of the array)
z=length(n)
z is the number of decimal digits of n.
z=pollard(x)
This attempts to return a factor of a composite x using Pollard's p1
method.
z=nprime(x)
This finds the first integer after x which passes the strong base 2
pseudoprime test and the Lucas pseudoprime test. (See [Pom].)
This integer is likely to be prime.
z=nprimeap(a,b,m)
This finds the first p, p=b(mod a), m ≤ p, which passes the strong
base 2 pseudoprime test and the Lucas pseudoprime test. Here a must be
even, b odd, 1 ≤ b < a, gcd(a,b)=1, b ≤ m.
z=jacobi(x,y)
z is the value of the Jacobi symbol (x/y).
z=peralta(a,p)
Peralta's algorithm is used to return a square root z of a (mod p).
Here a is a quadratic residue mod p. (See [Per].)
x=congr(a,b,m,&n)
Returns the solution x of the congruence ax=b(mod m). Also n=m/gcd(a,m)
is returned.
x=chinese(a,b,m,n,&l)
Returns the solution x(mod l) of the system of congruences x=a(mod
m) and x=b(mod n). Also l=lcm(m,n) is returned.
x=chinesea(a[],m[],&l)
Returns the solution x(mod l) of the system of congruences x=a[i](mod
m[i]), 0 ≤ i < n. Also l=lcm(m[0],...,m[n1]) is returned. (n is
the size of the array)
z=mthroot(x,m)
The integer part of the mth root of x is returned. (See [Mat].)
mthrootr(x,y,m,r)
The mth root of x/y is computed to r decimal places.
z=fund(d,&x,&y)
x and y are returned, where x+yω is the fundamental unit of Q(√d).
z=Norm(x+yω) is also returned.
z=pell(d,e,&x,&y)
The continued fraction expansion of √d is periodic after the first term:
a[0],a[1],...,a[n1],2a[0],a[1],...,a[n1],2a[0],.... Also the section
a[1],...,a[n1] is palindromic. We print a[0] and half the palindrome,
iff e is nonzero, sending the output to a file called pell.out.
The least solution x and y of Pell's equation x^{2}dy^{2}=epsilon is returned, where epsilon=1 or 1. z=epsilon is also returned. (See [Ros][382], [Ven][62], [Knu][359], [Po1][359], [Sie][296], [Dav][109]. Also see manuscript.)
z=surd(d,t,u,v,&a[],&u[],&v[],&p[],&q[])
The continued fraction of the quadratic irrationality (u+t√d)/v, d > 1, nonsquare, v > 0, is determined as far as the period. The period length z is returned. a[] is the array of partial quotients, u[] and v[] give the complete convergents (u[i]+√d)/v[i], p[] and q[] give the convergents p[i]/q[i]. Output is sent to a file surd.out. (See [Ros][379381] and [Knu][359].)
z=mpower(a,b,c)
z=a^{b}(mod c) is returned.
z=inv(a,m)
z=a^{1}(mod m) is returned. (See [Knu][325].)
z=elliptic(n,m,p)
The elliptic curve method is used to try to find a factor z of a composite
number n.
Here m, p < 2^{32} and 1279 ≥ m > 10, p ≥ 1.
z=factor(n)
The factorization of n is performed using the multiple polynomial quadratic
sieve on numbers with < 55 digits. The elliptic curve algorithm is used
on numbers of ≥ 55 digits. If succesful, the primes are certified and z=ω(n), the number of distinct prime factors
of n. There is an option to choose the number p of elliptic curves used.
(p=7 suffices to factor 10^{100}+1. (See [Sil][325].)
z=tau(n)
The divisor function z=d(n) is returned.
z=sigma(n)
z=σ(n), the sum of the divisors of n, is returned.
z=mobius(n)
The Möbius function z=μ(n) is returned.
z=euler(n)
Euler's function z=φ(n) is returned.
z=lprimroot(n)
The least primitive root mod p, an odd prime, is returned.
z=orderm(a,n)
The order of a (mod n) is returned.
z=lucas(n)
n is subjected to a strong pseudoprime test to base 2, together with
a Lucas pseudoprime test. If z=0 is returned, n is composite, while if
z=1 is returned, then n is a Lucas probable prime, as well as a base 2
strong pseudoprime. (See [Pom].)
serret(p,&x,&y)
Here p is a prime of the form 4n+1. Serret's algorithm is used to return
integers x and y such that p=x^{2}+y^{2}. (See [Wag].)
collatz(x,n)
Collatz' 3x+1 conjecture is tested. (See [Lag]
and 3x+1 page.)
The iterates x,T(x),... are printed iff n is nonzero.
cycle()
This finds cycles for the generalised Collatz mapping T, arising from
starting numbers p, p ≤ RANGE/2. Trajectories containing an iterate
whose magnitude exceeds a prescribed value INFINITY are deemed not to be
ultimately cycling. T(x)=int(m[i]x/d)+X[i] if x=i (mod d). Here the d
nonzero moduli m[i] and d shifts X[i] are entered from the keyboard. (See
3x+1 page.)
If when x ≡ i (mod d), T(x) is given in the form T(x)=(m[i]x+r[i])/d, where d divides m[i]i+r[i], then X[i]=(m[i]i+r[i])/dint(m[i]i/d)=r[i]/d+{m[i]i/d}, where {t} denotes the fractional part of t.
In particular, if k is odd, the "3x+k" mapping is given by m[0]=1, m[1]=3, X[0]=0, X[1]=(k+1)/2.
The output is sent to collatz.out.
miller(m,b)
Here m > 1,b > 1 and m does not divide b. Miller's test to base b is
applied to m.
hermite( )
improvep()
The output unimodular matrix contained in the file hermitep.out is
improved using LLLBabai method, followed by Gauss lattice reduction.
The output is sent to a file improvep.out.
lllhermite()
This performs a recent LLL based Hermite normal form algorithm due
to Havas, Majewski and Matthews. The output unimodular matrix P of a recent paper is sent to lllhermitetrans.out. The HNF(A) is sent to lllhermitebas.out.
In verbose mode, the intermediate steps are printed out.
shermite(N)
This performs the LLL algorithm on [I_{m}N^{n}A_{1}...NA_{n}],
where A is an mxn matrix of integers. The columns are scaled down by the
corresponding power of N for viewing convenience. The output matrix, in
scaled form, is sent to shermitebas.out.
smith()
The Smith normal form SNF(A) of an integer matrix A is found using
a pivoting strategy due to G. Havas. A cutoff value is requested. When
coefficients grow above this value in size, the MLLL algorithm is used
to reduce coefficient explosion. MLLL is also used at the start of searching
for each invariant factor, as it often yields small vectors with potential
small invariant factor components. Unimodular matrices P and Q are found
such that PAQ = SNF(A). the invariant factors, P and Q are sent to files
smith.out, smithp.out and smithq.out, respectively, if desired.
mlll()
The MLLL algorithm of M. Pohst is applied to an integer matrix whose
first row is nonzero. (See [Po2][209210].)
The reduced matrix and transforming matrix are sent to files mlllbas.out
and mllltran.out, respectively.
axb()
This solves the linear system AX=B, where A,X,B have integer coefficients
and the first column of A is nonzero. In the case of solubility, file axb.out
is either the unique solution or a short solution or the shortest solution,
depending on choice; axbbas.out is a short basis for the nullspace of A.
(See preprint.)..
fp()
This is the simplest form of the FinckePohst algorithm (see [Po2][190]).
This takes an integer matrix with LI rows and a positive integer C as input
and finds all lattice vectors X with X^{2} ≤ C. The vectors
found are sent to a file called fp.out. It is a good idea to start with
a matrix whose rows are LLL reduced.
slv()
This finds the shortest vectors in the lattice spanned by the LI family
b[1],...,b[m]. It applies FinckePohst to examine all nonzero X in L satisfying
X^{2} < C = b[1]^{2}. If it finds an X shorter
than b[1], the new bound C=X^{2} is chosen. At the end, FinkePohst
has only the shortest vectors to enumerate.
Only the vectors with highest positive coefficient of b[i] are given.
The output is also sent to slv.out.
inhomfp()
The inhomogeneous version of the FinckePohst algorithm. (See [Po2][190].)
Input: An m x M integer matrix A whose first m1 rows are LI and which
are preferably in LLL reduced form. A positive integer C is also entered.
Let L be the lattice spanned by the first m1 rows of A and let P be the
last row of A.
Output: All lattice vectors X of L such that XP^{2} ≤ C.
The vectors found are sent to a file called inhomfp.out.
e=rsae(p,q)
Here p and q are distinct odd primes, each > 355142 and e is the least integer such that gcd(e, (p1)(q1))=1 and 32^{e}> pq.
e is for use as an RSA encryption modulus.
encode(e,n)
Here n=p*q, a product of primes each greater than 355142. (p*q > 126126126126.)
e is the RSA encryption modulus, found using rsae( ) above.
A string of noncontrol characters when entered from the keyboard or
from a file consisting of lines each containing less than 500 characters.
These characters have ascii values in the range 32126. The message string
is encoded using the RSA algorithm: every 4 characters are converted to
ascii, joined as strings and the resulting large number m is encoded as
n = m^{e}(mod p*q). The encoded numbers are sent to a file "encoded.out",
terminated by an entry 1.
decode(e,p,q)
This calculates the decryption modulus d, then decodes each number
n in the file "encode": m = n^{d}(mod pq). The ascii characters
are split off and the original string of message characters is reformed.
The encoded message is sent to a file "decoded.out", as well as to
the screen.
addcubicr()
The sum of two points on the elliptic curve y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
powercubicr()
The point nP, where P is on the elliptic curve y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
ordercubicr()
The order of point P on the elliptic curve: y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
addcubicm()
The sum of two points on the elliptic curve mod p: y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
powercubicm()
The point nP, where P is on the elliptic curve mod p: y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
ordercubicm()
The order of point P on the elliptic curve mod p: y^{2}+a_{1}xy+a_{3}y=x^{3}+a_{2}x^{2}+a_{4}x+a_{6}
is calculated. The discriminant is also calculated.
convergents(a[],&p[],&q[])
The convergents p[0]/q[0],...,p[n]/q[n] of the continued fraction [a[0];a[1],...,a[n]]
are returned as arrays p[] and q[]. Here a[0] is an integer, a[1],...,a[n]
are positive integers.
lagrange("poly",&a[],m)
"poly" 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.
(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) 112134; Addendum 267 (1974) ibid. 219220 and page 261, Number Theory with Applications, by R. Kumanduri and C. Romero.
P. Shiu, Computation of continued fractions without input values,
Math. Comp. 64 (1995), no. 211, 13071317,
and D.G. Cantor, P.H. Galyean, H.G. Zimmer [Can].
(If the root is rational, Lagrange's algorithm still finds the partial quotients
and we exit at the appropriate point.
We check that poly is squarefree, has leading coefficient positive, does not vanish at x=0 or x=1, is not linear and has exactly one positive real root t and that t > 1.
z=perfectpower(n)
Here n > 1. z= x if n = x^{k}, for some x, k > 1, NULL otherwise.
See E. Bach and J. Sorenson, "Sieve algorithms for perfect power testing",
Algorithmica 9 (1993) 313328.
z=leastqnr(p)
Returns n_{p}, the least quadratic nonresidue (mod p), if
n_{p} < 65536, otherwise returns 0.
sturm(f(X))
Prints rational open intervals that are guaranteed to each contain only
one real root of the polynomial f(X). It is assumed that f(X) has no rational roots. The output is sent to the file sturm.out and to the screen.
See ([Akr],page 341.)
rootexp(f(X),m)
Finds the first m+1 partial quotients of the continued fraction expansion of all real roots of a squarefree polynomial P(X) with no rational roots. The output is sent to a file rootexp.out and to
the screen. Here is an outline of the algorithm.
STURM(POLYI P) (p. 346 Akritas) returns a stack of open rational intervals that are guaranteed to enclose the real roots of P. It uses CAUCHY(POLYI P), which calculates an upperbound for the positive roots of P using using Cauchy's rule (Theorem 7.2.11, p. 350, Akritas).
Using the process described in pages 786787 of Cantor, Galyean and Zimmer, each interval is eventually converted to the interval (1, ∞). At each stage, a partial
quotient a_{n1} is calculated using rootInIntervalR(POLYI P, MPR *LEFT, MPR *RIGHT), using a bisection method (described on p. 789) to find the integer part of the unique root of P in the OPEN interval (LEFT, RIGHT); also a sequence of corresponding polynomials f_{n}(X) is calculated recursively by f_{n}(X) = X^{m}f_{n1}(X)(X^{1} + a_{n1}), 0 ≤ n ≤ N.
The method of Lagrange is then used to find the partial quotients arising from the unique root in the interval (1, ∞) of f_{N}(X). See Algorithm 11.3.5, p. 261 of Kumanduri and Romero.
Programmed by Sean Seefried in January 2000.
content(f(X))
Returns the content of a polynomial f(X) (the gcd of the coefficients of f(X) x sign of the leading coefficient of f(X)). eg.
> content(3X^2+6)
3
primitive(f(X))
Returns the primitive part pp(f(X)) of a polynomial f(X), where f(X)=content(f(X))*pp(f(X)). eg.
>primitive(3X^2+6)
X^22
>content(3X^2+6)
3
z=sqroot(a,n,&s[],&m)
This solves x^{2} ≡ a (mod n). The number z of solutions mod n is returned. The solutions have the form ±s[0],...,±s[r] (mod m).
If there is no solution, z=0 and s[0]=NULL are returned.
cornacchia(a,b,m)
See L'algorithme de Cornacchia, A. Nitaj, Expositiones Mathematicae, 358365. This returns the positive primitive solutions (x,y) of ax^{2}+by^{2}=m, x ≥ y, where a,b,m are positive integers, with m ≥ a+b, gcd(a,m)=1=gcd(a,b). (Actually Nitaj also requires gcd(b,m)=1, but this does not seem to be necessary.)
patz(d,n)
This returns the positive primitive fundamental solutions of the diophantine equation x^{2}d*y^{2}=±N, where d > 1 is not a perfect square. The algorithm is from a recent paper of the author and is in essence due to Lagrange 1770.
z=congq(a,b,c,n,&s[])
This solves the quadratic congruence ax^{2}+bx+c ≡ 0 (mod n), where gcd(a,n)=1, n > 0. The solutions are s[0],...,s[z1].
z=binform(a,b,c,n,e)
This solves the diophantine equation ax^{2}+bxy+cy^{2}=N, N nonzero, where D=b^{2}4ac > 0 and is not a perfect square. One solution from each equivalence class is printed, along with the corresponding solution n of the congruence n^{2} ≡ D (mod 4N), N < n ≤ N.
e=1 is verbose, e=0 is terse. Output is sent to binform.out. The algorithm goes back to Lagrange (1770) and described in a paper to be submitted by the author.
z=absmod(a,b)
z=r=a(mod b) if r ≤ b/2, otherwise z=rb.
z=ceil(a,b)
This finds the least integer z not less than a/b, where a and b are integers, b nonzero.
log(a,b,d,r,&m[],&l)
This delivers with a high degree of certainty, partial quotients a[s] of log_{b}a, where a > b > 1, d> 1, r > 0. l is the number of partial quotients returned. Output is sent to log.out. (See paper.)
logx(a,b,m,n)
We conjecture that this gives correct partial quotients of log_{b}a. Output is sent to logx.out. (See paper where we are using Algorithm 1 with cutoff bb > c + b^{2}√c.)
r=resultant(p,q)
This returns the resultant of nonconstant polynomials p and q.
r=discriminant(p)
This returns the discriminant of a nonlinear polynomial p.
q=deriv(p)
This returns the derivative of polynomial p.
c=primes(m,n)
This prints the c primes, (c ≥ 0), in the interval [m,n], where 1 < m,n < 10^{10}. Output is sent to primes.out.
c=sturmsequence(f,b,e)
This produces the Sturm sequence for the polynomial f, evaluates the sequence at x=b and calculates the number c of signchanges. No output if e=0, output if e=1.
f is checked to see if it is squarefree and nonlinear and that f(b) is nonzero.
p=cyclotomic(n)
This returns the nth cyclotomic polynomial, n < 65536. See H. Lüneburg Galoisfelder, Kreisteilungskörper und Schieberegisterfolgen, BI Mannheim 1979 or
MP473 notes.
h=classnop(d)
This returns the class number h=h(d) of a real quadratic field Q(√d). Here 1 < d < 10^{6} is squarefree and D is the field discriminant.
We locate all reduced irrationals of the form (b+√D)/(2c), where c is negative and 4c divides db^{2}. We use the PQa continued fraction algorithm of Lagrange to break the set into disjoint cycles, retaining one number from each cycle. Each reduced number then gives rise to a reduced form (a,b,c) of discriminant D, where a=(b^{2}D)/(4c).
We are able to also determine if the Pell equation x^{2}Dy^{2}=4 has a solution, thereby finding the norm of the fundamental unit.
(See Henri Cohen's A course in computational number theory, page 260, First Edition.)
h=classnon(d,e)
Here d < 0 and 1 < d < 10^{6} is squarefree and d=0 or 1(mod 4).
This is Henri Cohen's Algorithm 5.3.5, p. 228, for finding the class number h(d) of binary quadratic forms of discriminant d, when d < 0.
h(d) is returned in each case.
If e=1, we print only the primitive forms.
If e=0, we print primitive and imprimitive forms.
If d is the discriminant of an imaginary quadratic field K, then the primitive forms classnumber h(d) is also the class number of K.
Davenport's Higher Arithmetic has a table of forms, which lists the imprimitive ones with an asterisk.
z=nearint(a,b)
Here b > 0. Then z is the nearest integer to a/b and where z=t, if a/b=1/2+t, t an integer. In fact nearint(a,b)=(aabsmod(a,b))/b.
h=reduceneg(a,b,c)
This is Gauss's algorithm for reducing a positive definite binary quadratic form. See L.E. Dickson, Introduction to the theory of numbers, page 69.
The reduced form (A,B,C) satisfies A < B ≤ A, C ≥ A, with B ≥ 0 if C = A.
The number h of steps taken in the reduction is returned.
h=reducepos(a,b,c)
Given an indefinite binary quadratic form ax^{2}+bxy+cy^{2}, we use the PQa continued fraction algorithm to determine a reduced form and thence a cycle of reduced forms. A unimodular matrix is constructed which converts (a,b,c) to reduced form. The output is sent to reducepos.out.
Note: d=b^{2}4ac > 0, d is not a perfect square and we assume d < 10^{6}.
(See Henri Cohen's A course in computational number theory, Edition 1, pp. 257261.)
h=classnop0(d)
Here 1 < d < 10^{6} is not a perfect square and d is 0 or 1 (mod 4). Returns the number h of classes of binary quadratic forms of discriminant d
A complete set of reduced binary forms is given.
We locate all reduced irrationals of the form (b+√d)/(2c), where c is negative and 4*c divides db^{2}. We use the PQa continued fraction algorithm of Lagrange to break the set into disjoint cycles, retaining one number from each cycle. Each reduced number then gives rise to a reduced form (a,b,c) of discriminant d, where a=(b^{2}d)/(4c). We are able to also determine if the Pell equation x^{2}d*y^{2}=4 has a solution, by using the fact that the equation is soluble iff at least one of the above cycles is odd. If there is no solution, the reduced forms (a,b,c) have to be counted as well. (See G.B. Mathews, Theory of Numbers, 8081.)
(Also see Henri Cohen's A course in computational number theory, page 260, First Edition.)
h=tableneg(m,n)
Here 1 ≤ m ≤ n < 10^{6}. Calculates h(d) for all squarefree d with m ≤ d ≤ n. The number h of squarefree d in the range is returned. The output is also sent to tableneg.out.
h=tablepos(m,n)
Here 2 ≤ m ≤ n ≤ 10^{6}. Calculates h(d) for all squarefree d with m ≤ d ≤ n. The number h of squarefree d in the range is returned. The output is also sent to tablepos.out.
h=davison(l,m,n)
Here l,m > 0, 10^{5} ≥ n ≥ 0. We find h partial quotients a[i] of e^{l/m}. We cannot predict the value of h. The a[i] are also sent to davison.out. The program stops if 10^{6} partial quotients a[i] are found.
h=raney(p,q,r,s)
Input: a nonsingular matrix A=[p,q;r,s], p,q,r,s>=0, A!=I_{2}, A!=[0,1;1,0].
With L=[1,0;1,1] and R=[1,1;0,1], we express A uniquely as a product of nonnegative powers of L and R, (at least one is positive) followed by a rowbalanced B.
B=[a,b;c,d] is rowbalanced if (a < c & b > d) or (c < a & d > b) and a,b,c ≥ 0.
The number h of powers of L and R is returned.
The maximum number of partial quotients returned is 10^{6}.
Also BCMATH version.
(See G.N. Raney, On continued fractions and finite automata, Math, Annalen 206 (1973) 265283.)
h=unimodular(p,q,r,s)
This program expresses a unimodular matrix A = [p,q;r,s] ≠ I_{2} or U = [0,1;1,0] with nonnegative coefficients, as a product of the type P, UP, PU, or UPU, where P is a product of matrices of the form U_{a}=[a,1;1,0], a>0.
The representation is unique.
(See Kjell Kolden, Continued fractions and linear substitutions, Arch. Math. Naturvid. 50 (1949), 141196.
Also see the corresponding BCMATH program.)
twoadicsqrt(b,n,&a[])
Here b > 0, b=8k+1.
Finds n terms a[0]...a[n1]of the 2adic square root x of b, x ≡ 1 (mod 4).
Output also sent to 2adic.out.
padicsqrt(b,n,p,&a[])
Here b > 0, b a quadratic residue (mod p), p an odd prime.
First finds a square root u of b (mod p), 0 < u < p, and finds n terms a[0]...a[n1]of the padic square root x of b, x ≡ u (mod p).
Output also sent to padic.out.
See lectures and solutions.
u=ramanujan(n)
Returns the value of Ramanujan's tau function, 0 < n < 2^{16}.
z=repdefinite(a,b,c,m,print_flag)
This solves the diophantine equation m=ax^{2}+bxy+cy^{2}, where d=b^{2}4ac < 0, a > 0, c > 0. See BCMATH version for the algorithm due to Gauss. print_flag = 1 lists solutions and unimodular transformations, while print_flag = 0 lists only the solutions. The output is sent to repdefinite.out.
powerd(a,b,d,n,&aa,&bb).
Here (a + b√d)^{n} = aa + bb√d.
z=euclid1(a,b).
z is the length of Euclid's algorithm. Here a and b are positive integers.
z=rcfperiod(d).
z is the length of the period of the regular continued fraction expansion of √d.
z=nscfperiod(d).
z is the length of the period of the nearest square continued fraction expansion of √d.
z=nicfperiod(d).
z is the length of the period of the nearest integer continued fraction expansion of √d.
algcycle().
This looks for cycles for the dbranched generalised Collatz mapping T: ℤ[√d]→ℤ[√d]. The output is sent to alg_cycle.out.
spiral(n,&x,&y).
This produces the spiral point (x,y) corresponding to a nonnegative n. See page 99 of Concrete Mathematics, Graham, Knuth, Pataschnik.
n=spiralinverse(x,y).
This produces the n corresponding to spiral point (x,y). See page 99 of Concrete Mathematics, Graham, Knuth, Pataschnik.
z=nicfperiod0(d,e,&x,&y).
z is the length of the period of the nearest integer continued fraction expansion of √d. (x,y) is the fundamental solution of x^{2}dy^{2}=±1. e=1 prints the type of midpoint criteria.
z=nscfperiod0(d,e,&x,&y).
z is the length of the period of the nearest square continued fraction expansion of √d. (x,y) is the fundamental solution of x^{2}dy^{2}=±1. e=1 prints the type of midpoint criteria.
z=rcfperiod0(d,e,&x,&y).
z is the length of the period of the regular continued fraction expansion of √d. (x,y) is the fundamental solution of x^{2}dy^{2}=±1. e=1 prints the type of midpoint criteria.
carmichael(n).
This solves φ(x) = n. See BCMath description.
carnielli().
This finds cycles for Walter Carnielli's generalization of the 3x+1 mapping. See BCMath description. The output is sent to carnielli_cycle.out.
lupei().
This finds cycles for Lu Pei's generalization of the 3x1 mapping. See BCMath description. The output is sent to lupei_cycle.out.
tangent(n).
This computes the nth tangent number t(n) for 1 ≤ n ≤ 2000. Also see BCMATH program.
bernoulli(n,&x,&y).
This computes the nth Bernoulli number B(n) for 0 ≤ n ≤ 4000. Also see BCMATH program.
partition(n).
This computes the nth partition number p(n) for 1 ≤ n ≤ 65535. Also see BCMATH program.
twocycle( ).
This tests for cycles of the 2branched mapping T(x)=int(ai*x/mi)+xi, where i = mod(x,2), a0,a1 are nonzero multipliers, m0,m1 are nonzero positive divisors and x0,x1 are arbitrary shifts.
mcycle( ).
This tests for cycles of a 3x+1 type mapping due to Benoit Cloitre. Here f(x)=int((m+1)/m)x) and T(x)=f(x)/2 or (3f(x)+1)/2, according as f(x) is even or odd.
pqcycle( ).
This tests for cycles of a 3x+1 type mapping due to Benoit Cloitre. Here f(x)=int(p*x/q), p > q > 1, q not dividing p, and T(x)=f(x)/2 or (3f(x)+1)/2, according as f(x) is even or odd.
h=fg(p,q).
This calculates h=p(q(X)).
nagell(d,n).
This finds the fundamental solutions u+v√d of x^{2}  dy^{2} = n, where d > 0 is not a perfect square and n is nonzero; it also gets a basis for the nonnegative solutions. The algorithm works only for small d and is based on upper estimates for v due to Nagell, Tchebicheff and Frattini.
frattini(d,n).
This finds a set of nonnegative solutions u+v√d of x^{2}  dy^{2} = n, where d > 0 is not a perfect square and n is nonzero. These solutions generate all nonnegative solutions by multiplication by ε^{n}, n≥ 0, where ε=x_{1}+y_{1}√d is the fundamental solution of Pell's equation x^{2}dy^{2}=n. The algorithm works only for small d and n and is based on an upper estimate v < y_{1}√n due to Frattini 189192.
branch(n,type,flag).
This constructs the polynomial triple (k(t),x(t),y(t)) corresponding to a path given by the reverse digits of the negative 3adic expansion of
N, where the root node is (t,t,0) if type=0, (t,t^{2}t+1,t1) if type=1,
(t,t^{2}+t+1,t+1) if type=1. If flag is nonzero, we print (k(t),x(t),y(t)), otherwise only k(t). Dujella's unicity conjecture is equivalent to the statement that the values of k(t) as t ranges over t >= 2 for type =0, t >= 1 for type 1 and t>=2 for type 1, are distinct.
pell4(d,e,&x,&y).
This finds the least positive integer solution of x^{2}  dy^{2} = 4, where d > 1 is not a perfect square. This is a CALC version of the function fund in the BC program stolt.
stolt(a,b,c,n).
This finds the fundamental solutions of the diophantine equation ax^{2} + bxy + cy^{2} = n, where a > 0, n ≠ 0 and b^{2}  4ac is positive and not a perfect square. This is a CALC version of the function stolt in the BC program stolt.
lprimefactor(n).
This finds the least prime factor of n.
z=conwaycycles(a,b).
This performs Conway's least prime factor sequence and checks to see that it cycles. It returns the length of the cycle that is reached. See the BCMATH version.
conwayrangetest1(m1,m2,n1,n2).
Here m1 ≤ m2, n1 ≤ n2 are positive integers. We test Conway's conjecture for his least prime factor sequence in the range [m1,m2], [n1,n2] and believe that all sequences will eventually cycle. I believe that apart from cycles of length 1, there are just 6 cycles as of 31st May 2016. See the BCMATH version.
The author is grateful to Peter Adams for help in understanding the mysteries of Yacc and C programming.
Sean Vickery helped in getting the PC version up and running.
1999/2000 summer vacation scholar Sean Seefried vastly improved the error handling, added polynomials to the parser and coded sturm(), rootexp() and a much improved lagrange().
BIBLIOGRAPHY
 [Akr] A. Akritas.
 [Can] D.G. Cantor, P.H. Galyean, H.G. Zimmer, Continued Fraction Algorithm for Real Algebraic Numbers, Math. Comp. 119, 785791.
 [Dav] H. Davenport, The Higher Arithmetic, Cambridge University Press, 1992.

[Hav] G. Havas, B. Majewski, K.R. Matthews, Extended gcd and Hermite normal form algorithms via lattice basis reduction, Experimental Mathematics 7 No.2 (1998) 125136.

[Knu] D.E. Knuth, The Art of Computer Programming,
Volume 2. AddisonWesley, 1981.

[Lag] J. Lagarias, The 3x+1 problem and its generalizations, American Math. Monthly, 92 (1985) 323. Also see Jeff Lagarias' page and The Generalized 3x+1 Mapping by Keith Matthews.

[Mat] K.R. Matthews, Computing mth roots, College Mathematics Journal 19 (1988) 174176.
 [Mat_log] K.R. Matthews On a version of Shanks' algorithm for computing the continued fraction of log_{b}a (see published article)

[Per] R. Peralta, A simple and fast probabilistic
algorithm for computing square roots modulo a prime number, I.E.E.E.
Trans. Inform. Theory, IT32, (1986) 846847.

[Po1] M. Pohst and H. Zassenhaus, On unit computation
in real quadratic fields, EUROSAM '79, Lecture Notes in Computer Science Vol. 79,
(1972) 140152.

[Po2] M. Pohst and H. Zassenhaus, Algorithmic Algebraic
Number Theory, Cambridge University Press, 1989.

[Pom] C. Pomerance, J.L. Selfridge and S.S. Wagstaff
Jr., The Pseudoprimes to 25x10^{9}, Mathematics of Computation,
35 (1980) 10031026.

[Ros] K.H. Rosen Elementary Number Theory and Its Applications, (4th edition) AddisonWesley 2000
 [RosenShallit] D. Rosen, J. Shallit, A continued fraction algorithm for approximating all real polynomial roots, Math. Magazine, 51 (1978) 112116.

[Sie] W. Sierpinski, Theory of numbers, PWN,
1964.

[Sil] R.D. Silverman, The multiple polynomial quadratic
sieve, Mathematics of Computation, 48 (1987) 329339.

[Sim] C.C. Sims, Computing with finitely presented
groups, Cambridge University Press, 1994.

[Ven] B.A. Venkov, Elementary number theory,
WoltersNoordhoff, 1970.

[Wag] S. Wagon, The Euclidean algorithm strikes
again, American Math. Monthly, 97 (1990) 125134.
CALC is maintained by Keith Matthews.
The software is intended as a service to the number theory community,
but the author cannot be held responsible for any consequences, either
direct or indirect, which the use of this package may have. It can be freely
copied for noncommercial purposes.
Email
http://www.numbertheory.org/keith.html
Last modified 14th January 2019