BC PROGRAMS WRITTEN BY KEITH MATTHEWS

1. gcd:
   • sign(n).
   • abs(n).
   • mod(a,b), (b > 0): returns a(mod b).
   • int(a,b): returns the integer part of a/b, b nonzero.
   • gcd(m,n).
   • gcd1(m,n): gcd(m,n)=gcd1(m,n)·m+gcd2(m,n)·n.
   • gcd2(m,n).
   • lcm(m,n).
   • inv(a,n): returns the inverse of a mod m.
   • cong(a,b,m): solves ax ≡ b (mod m).
   • chinese(a,b,m,n): solves x ≡ a(mod m) and x ≡ b(mod n).
   • gcda(m[],n): finds gcd(m[0],...,m[n-1]) and expresses it as a linear combination of m[0],...,m[n-1].
   • lcma(m[],n): finds lcm(m[0],...,m[n-1]).
   • chinesea(a[],m[ ],n): solves x ≡ a[i](mod m[i]), i=1,...,n-1.
   • chineseb(a[],b[],m[],n): solves a[i]x ≡ b[i](mod m[i]), i=1,...,n-1.
   • mpower(a,b,c): returns a^b (mod c).
   • exp(a,b): returns a^b.
   • mthroot(a,b,m): returns the integer part of the mth-root of a/b. Here a,b and m are positive integers, m > 1. (See K.R. Matthews, Computing mth roots, College Mathematics Journal 19 (1988) 174-176.)
   • mthrootr(a,b,m,r): this gives the mth-root of a/b truncated to r places.
   • binomial(n,m): returns the binomial coefficient.
   • gcd3(a,b,c): returns gcd(a,b,c).
   • ceiling(a,b): returns least integer x>= a/b.
   • bezout(a,b), when a and b are non-negative, returns d = gcd(a,b) = globalu*a + globalv*b, where globalu and globalv are global variables. bezout1(a,b) deals with arbitrary integers a and b.
2. euclid: euclid(m,n) performs Euclid's algorithm.
3. euclid1: euclid(m,n) returns the length of Euclid's algorithm.
4. jacobi:
   • jacobi(m,n) calculates the Jacobi symbol
   • peralta(a,p) finds a square root of a quadratic residue a mod p, using an algorithm of Rene Peralta.
5. serret: serret(p) expresses a prime of the form 4n+1 as the sum of two squares using Serret's algorithm.
6. 3x+1: collatz(n) tests the 3x+1 conjecture.
7. 3x+371: s(n) tests the 3x+371 conjecture.
8. phi: (slow-uses Brent-Pollard only)
   • omega(n) returns the number of distinct prime factors of n.
   • phi(n) returns the value of Euler's function.
   • tau(n) returns the number of divisors of n.
   • sigma(n) returns the sum of the divisors of n.
   • mu(n) returns the value of the Mobius Function.
   • lprimroot(p) returns the least primitive root mod p.
   • orderm(a,m) returns the order of a mod m.
9. factors: factor(n) attempts to factor n using Brent-Pollard (slow).
10. lucas: lucas(n) performs the strong base 2 pseudoprime test and Lucas pseudoprime test on n.
    (Needs jacobi).
11. decimal: period(m,n,b) outputs the period digits of the base b expansion of m/n, where m,n,b (b > 1) are positive integers, 1≤m < n.
12. pell: pell(d,e) finds the least solution of Pell's equations x²-d*y²=±1,±2,±3 and least primitive solution of x²-d*y²=±4. If e=1, the complete and partial quotients are printed; if e=0, this detail is suppressed.
13. surd: surd(d,t,u,v) finds the continued fraction expansion of a quadratic irrational (u+t*sqrt(d))/v, where d > 1 is not a square, t,u,v integers, v nonzero.
14. unit: unit(d) finds the fundamental unit of Q(sqrt(d)).
15. fibonacci: f(m,n) prints the Fibonacci numbers F(m),...,F(n). l(m,n) prints the Lucas numbers L(m),...,L(n).
16. rootd: root(d) finds the continued fraction expansion of sqrt(d).
17. cfrac: cfrac(m,n) finds the continued fraction expansion of m/n.
18. proth: proth(h,m) investigates the primality of h*2^m+1, h < 2^m, using Proth's test.
    (Needs jacobi,lucas).
19. pollard: pollard(n) attempts to find a factor of n using the Pollard p-1 method.
20. challenge: s(n) tests another 3-branched generalized 3x+1 conjecture.
21. mordell: f(a,k) finds the integer solutions of y²=x³+a with x ≤ k. (Needs gcd.)
22. 4_branch: s(n) tests a 4-branched generalized 3x+1 function.
23. venturini1: s(n) tests a 6-branched generalized 3x+1 function of G. Venturini.
24. lra: (needs gcd)
    • lnearint(m,n) returns the (left) nearest integer to m/n.
    • rnearint(m,n) returns the (right) nearest integer to m/n.
    • lmodd(m,n) returns the (left) least remainder of sign(n)m mod |n|.
    • rmodd(m,n) returns the right) least remainder of sign(n)m mod |n|.
    • lra(m,n) performs the least remainder algorithm on m, n.
    • nicf(m,n) prints the nearest integer continued fraction of m/n: m/n = [a₀ - 1/a₁ - ··· - 1/a_n].
      We write this as (a₀,a₁,...,a_n).
25. convergents: pn(a[],n) and qn(a[],n) compute the numerator and denominator of the continued fraction [a[0];a[1],...,a[n]].
26. lagrange: lagrange(a[],n,m) uses the method of Lagrange (1797) to find the first m+1 partial quotients of t, where f(x)=a[n]xⁿ+···+a[0], a[n] > 0, is a polynomial with integer coefficients, having no rational roots and having exactly one real positive root t, this being > 1.
27. lupei: s(n) tests Lu Pei's 3-branched generalized 3x+1 mapping.
28. recursion (Recursive bc programs)
    • fib(n): returns the nth Fibonacci number, n ≥ 0.
    • luc(n): returns the nth Lucas number, n ≥ 0.
    • fac(n): returns factorial(n) if n ≥ 1.
29. tonelli: x=tonelli(a,p) returns a square root of a (mod p), deterministically.
30. discrete_log: r=shanks(n,g,p). Here g is a primitive root (mod p) and g^r ≡ n (mod p), 0 ≤ r < p-1.
    Note: p < 2³²-2¹⁶=4294901760, in order to satisy BC array upper bound length of 2¹⁶-1.
    If r does not exist, we return -1.
    We use Shanks' giant steps-baby steps approach as described in Algorithmische Zahlentheorie by Otto Forster, pp 65-66.
31. forster_log: r=shanks(n,g,p). Similar to discrete_log, except that p is no longer necessarily a prime. Now needs phi to provide orderm(a,m).
32. leastqnr: leastqnr(p) returns n_p, the least quadratic nonresidue mod p. (Needs gcd).
33. rootd_modn: rootd_modn(d,n) finds all solutions of the congruence x² ≡ d (mod n) with 0 ≤ x ≤ n/2 for small n.
34. thue: Here d>1, is not a square, n > 1 an odd integer not dividing d-1. u² ≡ d (mod n), 1 < u < n. Then thue(d,u,p) finds x,y such that x²-dy²=kn, with small k.
35. sqroot: sqroot(d,n), n > 1, finds all solutions of x² ≡ d (mod n). It returns -1 if there is no solution, otherwise returns the number of solutions (mod n).
36. tomas1 and tomas2: These are generalised 3x+1 examples studied by Tomás Oliveira e Silva, where all trajectories are eventually periodic.
37. log: log(a,b,d,r,e) performs a discrete variant of Shank's log_ba algorithm. Here 1 < r is an integer. Also d > 1. We do not guarantee the correctness of the output. Bigger (d,r) give more partial quotients. e=1 prints the A[i] and AA[i], while e=0 prints only the m[i] and mm[i]. See paper.
    I suggest the user runs test(a,b,d,m,n), over a range (m,n), where m ≤ n, to get an idea of the correct partial quotients. This runs log1(a,b,d,r) for r=m,...,n. To get an idea of the correct answer when m=n=r, we recommend taking m=r-t, n=r+t, with 1 ≤ t ≤ (say) 2.
38. base: f(b,n), n > 0, b > 1, gives the base b expansion of n.
39. perfect_power: perfect_power(n) produces 0 if n is not a perfect power; otherwise returns x and p, where n=x^p and p is the least prime with this property.
40. powerdd: powerdd(a,b,d,n) computes (a+b√d)ⁿ=z1+z2√d, with global variables z1 and z2.
41. primes: primes(m,n) prints the primes in the interval [m,n], if ma and n lie between 1 and 10¹⁰.
42. nprime: nprime(m) finds the least Lucas-base2 strong pseudoprime p satisfying p ≥ m.
43. nprimeap: nprimeap(m,a,b) finds the least Lucas-base2 strong pseudoprime p of the form p=ak+b and satisfying p ≥ m. Here 0 < b < a and gcd(a,b)=1.
44. sturm: sturm(a[],n,b,e) prints a sturm polynomial sequence for the (squarefree) polynomial a[n]xⁿ+···+a[0], evaluates the sequence at x=b and calculates its sign variation. e=0 suppresses printing.
45. john: johna(b[],n) takes an array of positive integers b[0],...,b[n-1] and replaces them by an array of positive integers a[0],...,a[n-1], where
    
    1. a[i] divides b[i] for 0 ≤ i < n;
    2. gcd(a[i],a[j])=1 for 0 ≤ i < j< n;
    3. lcm(b[0],...,b[n-1])=a[0]·a[1]···a[n-1].
    john(a,b) does the case n=2.
    The program is due to John Campbell. Needs the program gcd.
46. reducepos: reduce(a,b,c) takes as input an indefinite binary quadratic form ax²+bxy+cy² and uses the PQa algorithm to find a cycle of reduced forms. (See explanation.)
47. classnoneg: class_number(d,flag,table_flag) lists the reduced binary quadratic forms of negative discriminant d and returns their number h(d) if flag=0. If flag=1, only the primitive forms are counted. Only the class number is printed if table_flag=0.
    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.
    Algorithm 5.3.5 of Henri Cohen's A course in computational algebraic number theory is used.
    table(m,n) prints h(-d) for all -d ≡ 0 or 1 (mod 4) in the range 3 ≤ m ≤ d ≤ n, where n<10⁶.
48. reduceneg:
    • reduce(a,b,c) takes as input a positive definite binary quadratic form ax²+bxy+cy² and uses an algorithm of Gauss to find the equivalent reduced form and unimodular transforming matrix.
    • compose(a1,b1,c1,a2,b2,c2) composes two positive definite binary quadratic forms, with output reduced. Needs bezout from gcd.
    • compose1(a1,b1,c1,a2,b2,c2) nonverbose version of compose(). Needs bezout1 from gcd.
    • power_compose(a,b,c,n) outputs (a,b,c)ⁿ, where (a,b,c) represents a positive definite form. The output is reduced.
    • power_compose1(a,b,c,n) nonverbose version of tt>power_compose(a,b,c,n).
    • ordercomposite(a,b,c,flag) computes the order under composition, of the positive definite form (a,b,c). flag =1 is verbose, flag=0 nonvebose. We need phi to supply omega(h(d)) and gcd to supply bezout(a,b).
49. classnopos: class_number(d) (1 < d < 10⁶ and squarefree) finds the class number of the real quadratic field Q(√d) and the sign of the fundamental unit. A list of reduced binary quadratic forms corresponding to the ideal classes is also given.
    table(m,n) prints h(d) for all squarefree d in the range m ≤ d < n,where n < 10⁶.
    class_number0(d) (d > 0, not a perfect square and 0 or 1 (mod 4)) returns the class-number of binary quadratic forms of discriminant d. Also the solubility of x²-d*y²=-4 is determined.
    This is basically the same program as class_number(d), except that in the case of non-solubility of x²-d*y²=-4, we count the form (-a,b,c) as well as (a,b,c), a > 0 and this means we return twice the value that class_number(d) would otherwise have returned. Regarding this point, see G.B. Mathews, Theory of Numbers, pp. 80-81.
50. unimodular: unimodular(p,q,r,s) expresses a unimodular matrix A ≠ I₂ or U=[0,1,1,0] with non-negative coefficients, as a product of one of the following forms:
    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), 141-196.
    The number n of matrices in the product U₀ U_n-1 is returned.
51. binomial: binomial_p(n,k,p) finds the power of a prime p dividing the binomial coefficient .
    binomial(n,k) prints the prime power factorization of the binomial coefficient .
    
52. factorial: factorial_p(n,p) finds the power of a prime p dividing n!
    factorial(n) finds n!
    
53. p-adic: 2adic(a,n) returns the first n binary digits of a 2-adic sqroot x of a positive integer a=8k+1. Here x=1 or 5 (mod 8).
    padic(a,p,n) returns the first n p-adic digits of a p-adic sqroot x of a positive integer a which is a quadratic residue (mod p). Here x=b (mod p), where b²=a (mod p) and 0 < b < p.
    
54. raney: raney(p,q,r,s) expresses a nonsingular matrix A=[p,q;r,s] (≠ I_2 or U=[0,1;1,0]) as a product of positive powers of R=[1,1;0,1] and L=[1,0;1,1], followed by a row-balanced matrix D=[a,b;c,d]. ie. a < c & b > d or a > c & b < d. The number of terms L and R is returned.
    With U_a=[a,1;1,0], note that U_a0...U_a2n=R^a₀L^a₁...R^a_2nU₀ and that U_a0...U_a2n+1=R^a₀L^a₁...L^a_2n+1I₂.
55. davison: davison(l,m,n) performs the algorithm of J.L. Davison's paper An algorithm for the continued fraction of e^l/m, Proceedings of the Eighth Manitoba Conference on Numerical Mathematics and Computing (Univ. Manitoba, Winnipeg, 1978), 169--179, Congress. Numer., XXII, Utilitas Math.
    With n ≥ 0, we first find the n* of Davison's Proposition 4.1 and apply Raney's factorisation to A₀...A_k, for n* ≤ k ≤ n*+n.
    The number (count) of partial quotients of e^l/m found is returned. count becomes positive for all large n.
56. squareroot: This is an improved version of sqroot and contains cornacchia(a,b,m). This finds all positive primitive solutions of ax²+by²=m, where a > 0, b > 0, a+b > m > 0, gcd(a,b)=1=gcd(a,m). If a=b=1, we get solutions with y ≤ x.
    r=sqroot(d,n,e) returns the solutions of x²=d (mod n). r is the number of solutions (mod n). If e=1, we print the solutions (mod reduced_modulus) as reduced_solution[0],...,reduced_solution[count-1]. If e=0, solutions are not printed. Used eg. in cornacchia(). If omega(n) > 1, we use the Chinese remainder theorem after solving mod qglobal[i]^kglobal[i], i=0,...,omega(n)-1. The array solution[0],...,solution[numbr-1] consists of the solutions (mod n) in the range 0 ≤ x ≤ n/2 and is used in cornacchia().
    See A. Nitaj, L'algorithme de Cornacchia, Expositiones Mathematicae 13 (1995), 358-365.
57. phi now contains sigmak(k,n) and tau(n), where r=sigmak(k,n), k > 0,n > 0, returns the sum of the k-th powers of the divisors of n and u=tau(n) returns Ramanujan's tau function.
    We use the simplest formula for tau(n) given on page 140 of T.M. Apostol, Modular functions and Dirichlet series in number theory, 20-22, 140.
58. patz: patz(d,n,e) finds fundamental solutions for the diophantine equation x²-dy²=n, d > 0, d not a perfect square. We only find the fundamental solutions for classes P satisfying P² ≡ d (mod |n|) with 0 ≤ P ≤ |n|/2.
    e=1 is verbose and prints the partial quotients, complete quotients and convergents up to the period (resp. double period) according as the period-length of omega[j] and omega*[j] is even or odd. e=0 prints only the fundamental solutions.
    Needs squareroot.
59. squareroot now contains quadratic(a,b,c,n,flag). This returns the number k of solutions of the congruence ax²+bx+c ≡ 0 (mod n), where gcd(a,n)=1. The solutions in the range 0 ≤ x < n are returned as global variables quadratic_solution[0],...,quadratic_solution[k-1]. flag=1 prints the solutions.
60. patz now contains binary(a,b,c,n). This solves ax²+bxy+cy²=n, where n is non-zero and b²-4ac is positive and not a perfect square. The method is from the paper The Diophantine equation ax²+bxy+cy²=N, D=b²-4ac > 0, Journal de Théorie des Nombres de Bordeaux, 14 (2002) 257-270 by K.R. Matthews.
61. decimal2rational. Typing d2r(pre[],per[],k,r,b) will output the rational number with base b preperiod given by pre[k-1],...,pre[0] (if present) and period per[r-1],...,per[0] consisting of integers in the range [0,b-1].
    If there is no preperiod, we let pre[0]=0 and k=0.
    Example: Take pre[0]=0, per[0]=2, per[1]=1, k=0, r=2, b=10.
    Then typing d2r(pre[],per[],k,r,b) outputs .1212···=4/33.
62. sqrtd_period. Typing z=period(d) outputs the period-length z of the continued fraction of √d. We use the Pohst-Zassenhaus half-period trick. See paper.
63. nipell. (a) Typing nipell(61,0) finds the smallest solution of Pell equation x²-61y²=1 or -1, using the nearest integer continued fraction of √61. The period (or semi-period in the case that the negative Pell equation is soluble) is also printed.
    Typing nipell(61,1) prints partial quotients, complete convergents and convergents.
    (b) Typing nicf_pqa(d,t,u,v,e), e=0 or 1, finds the nearest integer continued fraction of (u + t√d)/v of Hurwitz.
    Typing nicf_pqa0(d,t,u,v,e), e=0 or 1, finds the nearest integer continued fraction of (u + t√d)/v in Perron's book. When e=1, the period partial numerators and denominators are in capitals.
    
64. nscf_pell. This program solves the Pell equations x² - dy² = ±1, using the NSCF algorithm - nearest square continued fraction algorithm. We follow the algorithm of A.A. Krisnaswami Ayyangar.
    (a) Typing nscf_pell(d,1) prints the partial quotients, convergents, complete convergents of the period of the nearest square continued fraction of sqrt(d), as well as the Pell equation solutions, whereas nscf_pell(d,0) prints only the continued fraction and solutions of the Pell equation.
    (b) Typing nscf_pqa(d,t,u,v,e), e=0 or 1, finds the nearest square continued fraction of (u + t√d)/v. P
65. spiral. This calculates the spiral of Concrete Mathematics, Graham, Knuth, Patashnik, Exercise 40, page 99. Needs gcd. Type spiral(n) and then x(n) and y(n) to get the coordinates of the n-th point of the spiral. Typing n=inverse_spiral(x,y) gives the inverse function.
66. equivalent. This produces a quadratic surd η=(A+√d)/C=(pξ+q)/(rξ+s), where ξ=(a+√d)/c and Δ=ps-qr=±1. Here A=Δ(prt+a(qr+ps)+qsc), C=Δ(r²t+2rsa+s^c) and t=(a²-d)/c. If ξ is in standard form, so is η. Type equiv(a,c,d,p,q,r,s).
67. carmichael. carmichael(n,e) solves φ(x)=n. e=0 prints all solutions, if any, while e=1 tests Carmichael's conjecture.
68. tangent. Typing tangent(n) produces the n-th tangent number T_n. (8th June 2011)
69. bernoulli. Typing bernoulli(n) produces the n-th Bernoulli number. (9th June 2011)
70. kronecker. Typing kronecker(d,n) returns the Kronecker symbol (d,n), n ≥ 1.
71. partition. Typing partition(m,1) prints the first m values of the partition function p(m) Typing partition(m,0) returns only p(m). (23rd June 2011)
72. cloitre. Typing cloitre(x) prints the trajectory of a mapping devised by Benoit Cloitre, who conjectured that all trajectories starting with a positive integer will eventually reach 0. Experimentally it also seems that all trajectories starting with a negative integer will eventually reach one of the cycles
    • -1,-2,-1;
    • -4,-17,-6,-4;
    • -19,-92,-31,-152,-51,-34,-167,-56,-19.
73. cloitrem. Typing cloitrem(m) if m is even or m=4n+1, produces the Markov matrix Q(2m) for the mapping F_m devised by Benoit Cloitre. The mapping is given by typing t(x,m). Finally, typing function(m) produces a listing of the d=2m branches of the mapping in two ways.
74. reducepos_zagier. Typing reduce(a,b,c) uses an algorithm from Don Zagier's book Zetafunktionen und quadratische Körper, Eine Einführung in die höhere Zahlentheorie, to find a reduced form equivalent to (a,b,c) and also prints the cycle this generates.
75. zagier. zagier(d) gives the backward continued fraction of √d. zagier0(d,t,u,v) gives the least integer continued fraction of (u + t\sqrt(d))/v.
76. genfacs. Given the canonical factorization p[0]^a[0]···p[t-1]^a[t-1]of n, typing exps(a[],t) lists the exponent vectors of the divisors of n, whereas typing divisors(p[],a[],t) lists the divisors of n.

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