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, 233-254, 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 Jean-Jacques 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
> 54321-12345
   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
a11 ... a1n
...
am1 ... amn
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 cdeg(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

  1. z=absmod(a,b)
  2. euclid(a,b,&q[ ],&r[ ],&s[ ],&t[ ],&n)
  3. z=gcd(x,y)
  4. z=gcdv(x,y,&u,&v)
  5. z=gcda(a[ ])
  6. z=gcdav(a[ ], &b[ ])
  7. egcd( )
  8. sgcd(N)
  9. lllgcd( )
  10. lllgcd0( )
  11. jacobigcd( )
  12. z=lcm(x,y)
  13. z=lcma(a[ ])
  14. z=length(n)
  15. z=pollard(x)
  16. z=nprime(x)
  17. z=nprimeap(a,b,m)
  18. z=jacobi(x,y)
  19. z=peralta(a,p)
  20. x=congr(a,b,m,&n)
  21. x=chinese(a,b,m,n,&l)
  22. x=chinesea(a[ ],m[ ],&l)
  23. z=mthroot(x,m)
  24. mthrootr(x,y,m,r)
  25. z=fund(d,&x,&y)
  26. z=pell(d,e,&x,&y)
  27. z=mpower(a,b,c)
  28. z=inv(a,m)
  29. z=elliptic(n,m,p)
  30. z=factor(n)
  31. z=tau(n)
  32. z=sigma(n)
  33. z=mobius(n)
  34. z=euler(n)
  35. z=lprimroot(n)
  36. z=orderm(a,n)
  37. z=lucas(n)
  38. serret(p,&x,&y)
  39. collatz(x,n)
  40. miller(m,b)
  41. lllhermite( )
  42. shermite(N)
  43. smith( )
  44. mlll( )
  45. fp( )
  46. slv( )
  47. cycle( )
  48. inhomfp( )
  49. e=rsae(p,q)
  50. encode(e,n)
  51. decode(e,p,q)
  52. addcubicr( )
  53. powercubicr( )
  54. ordercubicr( )
  55. convergents(a[ ],&p[ ],&q[ ])
  56. lagrange(f(X),&a[ ],m)
  57. z=perfectpower(n)
  58. axb( )
  59. addcubicm( )
  60. powercubicm( )
  61. ordercubicm( )
  62. leastqnr(p)
  63. sturm(f(X))
  64. rootexp(f(X), m)
  65. content(f(X))
  66. primitive(f(X))
  67. sqroot(a,n,&s[ ],&m)
  68. cornacchia(a,b,m)
  69. z=surd(d,t,u,v,&a[ ],&u[ ],&v[ ],&p[ ],&q[ ])
  70. patz(d,n)
  71. congq(a,b,c,n,&s[ ])
  72. binform(a,b,c,n,e)
  73. ceil(a,b)
  74. log(a,b,d,&a[ ],&l)
  75. logx(a,b,m,n)
  76. r=resultant(p(X),q(X))
  77. r=discriminant(p(X))
  78. q=deriv(p(X))
  79. c=primes(m,n)
  80. c=sturmsequence(f(X),b,e)
  81. p=cyclotomic(n)
  82. h=classnop(d)
  83. h=classnon(d)
  84. z=nearint(a,b)
  85. h=reduceneg(a,b,c)
  86. h=reducepos(a,b,c)
  87. h=classnop0(d)
  88. h=tableneg(m,n)
  89. h=tablepos(m,n)
  90. h=davison(l,m,n)
  91. h=raney(p,q,r,s)
  92. h=unimodular(p,q,r,s)
  93. twoadicsqrt(b,n,&a[])
  94. padicsqrt(b,n,p,&a[])
  95. ramanujan(n)
  96. repdefinite(a,b,c,m,print_flag)
  97. powerd(a,b,d,n,&aa,&bb)
  98. z=euclid1(a,b)
  99. z=rcfperiod(d)
  100. z=nscfperiod(d)
  101. z=nicfperiod(d)
  102. algcycle()
  103. spiral(n,&x,&y)
  104. n=spiralinverse(x,y)
  105. z=rcfperiod0(d,e,&x,&y)
  106. z=nscfperiod0(d,e,&x,&y)
  107. z=nicfperiod0(d,e,&x,&y)
  108. carmichael(n)
  109. carnielli()
  110. lupei()
  111. tangent(n)
  112. tangent(n,&x,&y)
  113. partition(n)
  114. twocycle( )
  115. mcycle( )
  116. pqcycle( )
  117. h=fg(p,q)
  118. nagell(d,n)
  119. frattini(d,n)
  120. branch(n,type,flag)
  121. pell4(d,e,&x,&y)
  122. stolt(a,b,c,n)
  123. z=lprimefactor(n)
  124. z=conwaycycles(a,b)
  125. conwayrangetest1(m1,m2,n1,n2)
  126. divisorlist(n)
  127. binaryform(a,b,c,n,&a[],&b[])
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[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.
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[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(√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[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]. 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][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. 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 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 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, RANGE/2 ≤ p < RANGE/2, where 0 < RANGE ≤ 232.
Trajectories containing an iterate whose magnitude exceeds a prescribed value INFINITY are deemed not to be ultimately cycling.
Here T(x)=int(mix/d) + Xi if x ≡ i (mod d), where the nonzero mi and Xi are entered from the keyboard. (See 3x+1 page.)
If when x ≡ i (mod d), T(x) is given in the form T(x) = (mix + ri)/d, where d divides mii + ri,
then Xi = (mii + ri)/d - int(mii/d) = ri/d + {mii/d}, where {t} denotes the fractional part of t.
In particular, if k is odd, the "3x+k" mapping is given by m0 = 1, m1 = 3, X0 = 0, X1 = (k + 1)/2.
The output is sent to cycle.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 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, each > 355142 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 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, 1307-1317,
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 = 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+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 786-787 of Cantor, Galyean and Zimmer, each interval is eventually converted to the interval (1, ∞). At each stage, a partial quotient an-1 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 fn(X) is calculated recursively by fn(X) = Xmfn-1(X)(X-1 + an-1), 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 fN(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^2-2
>content(-3X^2+6)
   -3

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 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+by2=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=±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 ax2+bx+c ≡ 0 (mod n), where gcd(a,n)=1, n > 0. The solutions are s[0],...,s[z-1].
z=binform(a,b,c,n,e)
This finds the primitive solution classes (if any) of the diophantine equation ax2+bxy+cy2=N, N non-zero, where D=b2-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 n2 ≡ D (mod 4|N|), -|N| < n ≤ |N|.
e=1 is verbose, e=0 is terse. The algorithm goes back to Lagrange (1770) and is described in a paper by the author.
z=absmod(a,b)
z=r=a(mod b) if r ≤ b/2, otherwise z=r-b.
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 logba, 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 logba. Output is sent to logx.out. (See paper where we are using Algorithm 1 with cutoff bb > c + b2√c.)
r=resultant(p,q)
This returns the resultant of non-constant polynomials p and q.
r=discriminant(p)
This returns the discriminant of a non-linear 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 < 1010. 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 sign-changes. 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 < 106 is squarefree and D is the field discriminant.
We locate all reduced irrationals of the form (b+√D)/(2|c|), where c is negative and 4c divides d-b2. 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=(b2-D)/(4c).
We are able to also determine if the Pell equation x2-Dy2=-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| < 106 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 class-number 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)=(a-absmod(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 ax2+bxy+cy2, 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=b2-4ac > 0, d is not a perfect square and we assume d < 106.
(See Henri Cohen's A course in computational number theory, Edition 1, pp. 257-261.)
h=classnop0(d)
Here 1 < d < 106 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)/(2|c|), where c is negative and 4*c divides d-b2. 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=(b2-d)/(4c). We are able to also determine if the Pell equation x2-d*y2=-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, 80-81.)
(Also see Henri Cohen's A course in computational number theory, page 260, First Edition.)
h=tableneg(m,n)
Here 1 ≤ m ≤ n < 106. 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 ≤ 106. 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, 105 ≥ n ≥ 0. We find h partial quotients a[i] of el/m. We cannot predict the value of h. The a[i] are also sent to davison.out. The program stops if 106 partial quotients a[i] are found.
h=raney(p,q,r,s)
Input: a non-singular matrix A=[p,q;r,s], p,q,r,s>=0, A!=I2, 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 non-negative powers of L and R, (at least one is positive) followed by a row-balanced B.
B=[a,b;c,d] is row-balanced 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 106.
Also BCMATH version.
(See G.N. Raney, On continued fractions and finite automata, Math, Annalen 206 (1973) 265-283.)
h=unimodular(p,q,r,s)
This program expresses a unimodular matrix A = [p,q;r,s] ≠ I2 or U = [0,1;1,0] with non-negative coefficients, as a product of the type P, UP, PU, or UPU, where P is a product of matrices of the form Ua=[a,1;1,0], a>0.
The representation is unique.
(See Kjell Kolden, Continued fractions and linear substitutions, Arch. Math. Naturvid. 50 (1949), 141-196.
Also see the corresponding BCMATH program.)
twoadicsqrt(b,n,&a[])
Here b > 0, b=8k+1.
Finds n terms a[0]...a[n-1]of the 2-adic square root x of b, x ≡ 1 (mod 4).
Output also sent to 2-adic.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[n-1]of the p-adic square root x of b, x ≡ u (mod p).
Output also sent to p-adic.out.
See lectures and solutions.
u=ramanujan(n)
Returns the value of Ramanujan's tau function, 0 < n < 216.
z=repdefinite(a,b,c,m,print_flag)
This solves the diophantine equation m=ax2+bxy+cy2, where d=b2-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 d-branched 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 non-negative 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 x2-dy2=±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 x2-dy2=±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 x2-dy2=±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 3x-1 mapping. See BCMath description. The output is sent to lupei_cycle.out.
tangent(n).
This computes the n-th tangent number t(n) for 1 ≤ n ≤ 2000. Also see BCMATH program.
bernoulli(n,&x,&y).
This computes the n-th Bernoulli number B(n) for 0 ≤ n ≤ 4000. Also see BCMATH program.
partition(n).
This computes the n-th partition number p(n) for 1 ≤ n ≤ 65535. Also see BCMATH program.
twocycle( ).
This tests for cycles of the 2-branched 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 x2 - dy2 = n, where d > 0 is not a perfect square and n is non-zero; it also gets a basis for the non-negative 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 non-negative solutions u+v√d of x2 - dy2 = n, where d > 0 is not a perfect square and n is non-zero. These solutions generate all non-negative solutions by multiplication by εn, n≥ 0, where ε=x1+y1√d is the fundamental solution of Pell's equation x2-dy2=n. The algorithm works only for small d and n and is based on an upper estimate v < y1√n due to Frattini 1891-92.
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 3-adic expansion of N, where the root node is (t,t,0) if type=0, (t,t2-t+1,t-1) if type=-1, (t,t2+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 x2 - dy2 = 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 ax2 + bxy + cy2 = n, where a > 0, n ≠ 0 and b2 - 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.
divisorlist(n).
This lists the divisors of n and sends the output to divisorlist.out.
binaryform(a,b,c,n,&a[],&b[]).
This finds representatives of the equivalence classes of solutions of ax2 + bxy + cy2 = n, where b2 - 4ac > 0 and is nonsquare, and sends the output to binaryform.out.
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

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 non-commercial purposes.

Email
http://www.numbertheory.org/keith.html

Last modified 2nd September 2022