Once a polynomial has been assigned to a variable, it is also possible
to evaluate that polynomial at a specified integer value.
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,3)
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
- 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)
- hermite( )
- lllhermite( )
- shermite(N)
- smith( )
- mlll( )
- fp( )
- slv( )
- cycle( )
- inhomfp( )
- printa(b,n)
- 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)) (20/1/2000)
- rootexp(f(X), m) (20/1/2000)
- content(f(X))
- primitive(f(X))
- log(a,b,d,u,v,e) (
- log1(a,b,r,e) (
- sqroot(a,n,&s[ ],&m) (
- cornacchia(a,b,m) (
- z=surd(d,t,u,v,&a[],&u[],&v[],&p[],&q[]) (
4/5/2000)
- patz(d,n) (
27/3/2000)
27/3/2000)
14/4/2000)
18/4/2000)
15/5/2000)
euclid(a,b,&q[ ],&r[ ],&s[ ],&t[ ],&m)
Here m=n+1 and the arrays
q[0]=NULL,...,q[n],q[n+1]=NULL,
r[0]=a,r[1]=b,...,r[n+1],
s[0]=a,s[1]=[,...,s[n+1],
t[0]=a,t[1]=[,...,t[n+1],
arising from Euclid's algorithm are returned and printed in euclid.out.
r[k]=r[k+1]*q[k+1]+r[k+2], 0 < r[k+2] < r[k+1],
s[k]=-q[k-1]*s[k-1]+s[k-2],
t[k]=-q[k-1]*t[k-1]+t[k-2],
r[n]=gcd(a,b)=s[n]*a+t[n]*b.
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[n-1]), where values for a[0],...,a[n-1] having been previously entered.
z=gcdav(a[ ],&b[ ])
z = gcd(a[0],...,a[n-1]). Also gives integers b[0],...,b[n-1] satisfying z = b[0]a[0]+...+b[n-1]a[n-1].
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[m-1],P is sent to an output file egcdmat.out.
Here X[1],...,X[m-1] form a LLL reduced basis for the lattice L defined by the equation x1d1+...+xmdm = 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 [In|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 +-NdEn, 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[m-1] form a LLL reduced basis for the lattice L defined by the equation x1d1+...+xmdm = 0.
The inhomogeneous version of the Fincke-Pohst algorithm (see [Po2][191]) can then be used as an option to find a shortest multiplier vector by solving the inequality
||X[m] - x1X[1]-...-xm-1X[m-1]||2 <= ||X[m]||2
in integers x1,...,xm-1.
Each time a shorter multiplier vector Q = X[m]-x1X[1]-...-xm-1X[m-1]
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[n-1]). (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 p-1 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[n-1]) is returned. (n is the size of the array)
z=mthroot(x,m)
The integer part of the m-th root of x is returned. (See [Mat].)
mthrootr(x,y,m,r)
The m-th 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(
).
z=Norm(x+y ) is also returned.
z=pell(d,e,&x,&y)
The continued fraction expansion of
is periodic after the first
term:
a[0],a[1],...,a[n-1],2a[0],a[1],...,a[n-1],2a[0],.... Also the section a[1],...,a[n-1] 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 x2-dy2=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].)
z=surd(d,t,u,v,&a[],&u[],&v[],&p[],&q[])
The continued fraction of the quadratic irrationality (u+t
)/v,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]+ )/v[i], p[] and q[] give the convergents p[i]/q[i]. Output is sent to a file surd.out. (See [Ros][379-381] and [Knu][359].)
z=mpower(a,b,c)
z=ab(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 < 232 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. 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 10100+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=x2+y2. (See [Wag].)
collatz(x,n)
Collatz' 3x+1 conjecture is tested. (See [Lag] and my survey .) 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 survey (ps).)
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( )
The Hermite normal form HNF(A) of an integer matrix A is found using the Kannan-Bachem algorithm. (See [Sim][349-357].)
The matrix can be entered either from the keyboard or from a file such as
2 3 1 -4 -7 3 2 4in the case of a 2 x 3 matrix. The Hermite normal form is sent to a file called hermite.out. A unimodular integer matrix P such that PA = HNF(A) is sent to a file hermitep.out, if desired. The last line of this file contains the value of rank(A).
improvep( )
The output unimodular matrix contained in the file hermitep.out is improved using LLL-Babai 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 [Im|NnA1|...|NAn], 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 non-zero. (See [Po2][209-210].)
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 Fincke-Pohst 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 Fincke-Pohst 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, Finke-Pohst 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 Fincke-Pohst algorithm. (See [Po2][190].)
Input: An m x M integer matrix A whose first m-1 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 m-1 rows of A and let P be the last row of A.
Output: All lattice vectors X of L such that ||X-P||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 and e is the least integer such that gcd(e, (p-1)(q-1))=1 and 32e> 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 non-control 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 32-126. 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 = me(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 = nd(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 y2+a1xy+a3y=x3+a2x2+a4x+a6 is calculated. The discriminant is also calculated.
powercubicr( )
The point nP, where P is on the elliptic curve y2+a1xy+a3y=x3+a2x2+a4x+a6 is calculated. The discriminant is also calculated.
ordercubicr( )
The order of point P on the elliptic curve: y2+a1xy+a3y=x3+a2x2+a4x+a6 is calculated. The discriminant is also calculated.
addcubicm( )
The sum of two points on the elliptic curve mod p: y2+a1xy+a3y=x3+a2x2+a4x+a6 is calculated. The discriminant is also calculated.
powercubicm( )
The point nP, where P is on the elliptic curve mod p: y2+a1xy+a3y=x3+a2x2+a4x+a6 is calculated. The discriminant is also calculated.
ordercubicm( )
The order of point P on the elliptic curve mod p: y2+a1xy+a3y=x3+a2x2+a4x+a6 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) 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. and D.G. Cantor, P.H. Galyean, H.G. Zimmer [Can].
z=perfectpower(n)
Here n > 1. z= x if n = xk, for some x, k > 1, NULL otherwise.
See E. Bach and J. Sorenson, "Sieve algorithms for perfect power testing", Algorithmica 9 (1993) 313-328.
z=leastqnr(p)
Returns np, the least quadratic non-residue (mod p), if np < 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 partial quotients of the continued fraction expansion of all real roots of a squarefree polynomial f(X) with no rational roots, using Lagrange's method and a variation presented in D.G. Cantor, P.H. Galyean, H.G. Zimmer [Can], especially pp. 785-787. The output is sent to a file rootexp.out and to the screen.
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^2-2
>content(-3X^2+6)
-3
log(a,b,d,u,v,e)
This performs a discrete variant of Shank's logba algorithm. Here d>1, 1<=u<v are integers. The larger is v-u, the better chance of correctness of output. eg. v=2u is fine, but slow if u is say 500. v=u+10 is adequate. We do not guarantee the correctness of the output.
Roughly r partial quotients are outputted if d=10. If e=1, the convergents and decimal expansion of logba are printed, the latter truncated correct to as many decimal places as possible: the r-1th convergent is compared with the rth convergent and the decimal expansions are truncated from where they differ. e=0 prints only the partial quotients, Output is sent to log.out. (See Algorithm 3 of [Mat_log]
log1(a,b,d,r,e)
This performs a discrete variant of Shank's logba algorithm. Here 1< r is an integer. Also d > 1. We do not guarantee the correctness of the output. d=2 is fastest, but bigger d give more partial quotients.
Roughly r partial quotients are outputted. If e=1, the convergents, the integers A[i] and decimal expansion of logba are printed, the latter truncated correct to as many decimal places as possible: the r-1th convergent is compared with the rth convergent and the decimal expansions are truncated from where they differ. e=0 prints only the partial quotients, Output is sent to log.out. (See Algorithm 1 of [Mat_log]
z=sqroot(a,n,&s[],&m)
This solves x2
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=m are returned.
cornacchia(a,b,m)
See L'algorithme de Cornacchia, A. Nitaj, Expositiones Mathematicae, 358-365. This returns the positive primitive solutions (x,y) of ax2+b2=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 x2-d*y2=±1, 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.
The author is grateful to Peter Adams for help in understanding the mysteries of Yacc. Sean Vickery helped in getting the PC version up and running. 1999/2000 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, Elements of Computer Algebra with Applications, Wiley 1989.
- [Can] D.G. Cantor, P.H. Galyean, H.G. Zimmer, Continued Fraction Algorithm for Real Algebraic Numbers, Math. Comp. 119, 785-791,
- [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) 125-136.
- [Knu] D.E. Knuth, The Art of Computer Programming, Volume 2. Addison-Wesley, 1981.
- [Lag] J. Lagarias, The 3x+1 problem and its generalizations, American Math. Monthly, 92 (1985) 3-23. Also see Jeff Lagarias' Organic Mathematics Article on the 3x+1 Problem and its Generalizations and The Generalized 3x+1 Mapping by Keith Matthews.
- [Mat] K.R. Matthews, Computing mth roots, College Mathematics Journal 19 (1988) 174-176.
- [Mat_log] K.R. Matthews On a version of Shanks' algorithm for computing the continued fraction of logba (pdf 152K)
- [Per] R. Peralta, 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.
- [Po1] M. Pohst and H. Zassenhaus, On unit computation in real quadratic fields, Lecture Notes in Computer Science Vol. 72, (1972) 140-152.
- [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 25x109, Mathematics of Computation, 35 (1980) 1003-1026.
- [Ros] K. Rosen, Elementary number theory and its applications, Addison-Wesley, 1984.
- [Sie] W. Sierpinski, Theory of numbers, PWN, 1964.
- [Sil] R.D. Silverman, The multiple polynomial quadratic sieve, Mathematics of Computation, 48 (1987) 329-339.
- [Sim] C.C. Sims, Computing with finitely presented groups, Cambridge University Press, 1994.
- [Ven] B.A. Venkov, Elementary number theory, Wolters-Noordhoff, 1970.
- [Wag] S. Wagon, The Euclidean algorithm strikes again, American Math. Monthly, 97 (1990) 125-134.
The software is intended as a service to the scientific 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 non-commercial purposes.
[email protected]
http://www.maths.uq.edu.au/~krm/
Last modified 7th August 2000