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.
   • nearestsquare(d) returns the nearest square to a nonsquare d.
2. euclid0: euclid0(m,n) performs Euclid's algorithm and finds the partial quotients and convergents for the continued fraction of m/n, n> 0.
3. euclid: euclid(m,n) performs Euclid's algorithm.
4. euclid1: euclid(m,n) returns the length of Euclid's algorithm.
5. 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.
6. serret: serret(p) expresses a prime of the form 4n+1 as the sum of two squares using Serret's algorithm.
7. 3x+1: collatz(n) tests the 3x+1 conjecture.
8. 3x+371: s(n) tests the 3x+371 conjecture.
9. 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.
10. factors: factor(n) attempts to factor n using Brent-Pollard (slow).
11. lucas: lucas(n) performs the strong base 2 pseudoprime test and Lucas pseudoprime test on n. lucasnonverbose(n) is the non-verbose version of lucas.
    (Needs jacobi).
12. 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.
13. 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.
14. surd: surd(d,t,u,v) finds the continued fraction expansion of a quadratic irrational (u+t√d)/v, where d > 1 is not a square, t,u,v integers, v nonzero.
15. unit: unit(d) finds the fundamental unit of Q(√d).
16. fibonacci: f(m,n) prints the Fibonacci numbers F(m),...,F(n). l(m,n) prints the Lucas numbers L(m),...,L(n).
17. rootd: root(d) finds the continued fraction expansion of √d.
18. cfrac: cfrac(m,n) finds the continued fraction expansion of m/n.
19. proth: proth(h,m) investigates the primality of h*2^m+1, h < 2^m, using Proth's test.
    (Needs jacobi,lucas).
20. pollard: pollard(n) attempts to find a factor of n using the Pollard p-1 method.
21. challenge: s(n) tests another 3-branched generalized 3x+1 conjecture.
22. mordell: f(a,k) finds the integer solutions of y²=x³+a with x ≤ k. (Needs gcd.)
23. 4_branch: s(n) tests a 4-branched generalized 3x+1 function.
24. venturini1: s(n) tests a 6-branched generalized 3x+1 function of G. Venturini.
25. 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).
26. convergents: pn(a[],n) and qn(a[],n) compute the numerator and denominator of the continued fraction [a[0];a[1],...,a[n]].
27. 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.
28. lupei: s(n) tests Lu Pei's 3-branched generalized 3x+1 mapping.
29. 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.
30. tonelli: x=tonelli(a,p) returns a square root of a (mod p), deterministically.
31. 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.
32. 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).
33. leastqnr: leastqnr(p) returns n_p, the least quadratic nonresidue mod p. (Needs gcd).
34. rootd_modn: rootd_modn(d,n) finds all solutions of the congruence x² ≡ d (mod n) with 0 ≤ x ≤ n/2 for small n.
35. 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.
36. 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).
37. tomas1 and tomas2: These are generalised 3x+1 examples studied by Tomás Oliveira e Silva, where all trajectories are eventually periodic.
38. 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.
39. base: f(b,n), n > 0, b > 1, gives the base b expansion of n.
40. 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.
41. powerdd: powerdd(a,b,d,n) computes (a+b√d)ⁿ=z1+z2√d, with global variables z1 and z2. It also contains fund4(d,e) which finds the fundamental solution of x² – dy² = 4.
42. primes: primes(m,n) prints the primes in the interval [m,n], if m and n lie between 1 and 10¹⁰.
43. nprime: nprime(m) finds the least Lucas-base2 strong pseudoprime p satisfying p ≥ m.
44. 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.
45. 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.
46. 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.
47. 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.)
48. 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⁶.
49. 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).
50. 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.
51. 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.
52. 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 .
    
53. factorial: factorial_p(n,p) finds the power of a prime p dividing n!
    factorial(n) finds n!
    
54. 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.
    
55. 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₂.
56. 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.
57. 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.
58. 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.
59. patz: patz(d,n,e) finds fundamental solutions for the diophantine equation x²-dy²=n, d > 0, d not a perfect square, n ≠ 0 and n ≠ 1.
    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.
60. 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.
61. 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.
62. 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.
63. 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.
64. 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.
    
65. 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 √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
66. 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.
67. 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).
68. carmichael. carmichael(n,e) solves φ(x)=n. e=0 prints all solutions, if any, while e=1 tests Carmichael's conjecture. Needs files phi, lucas, jacobi.
69. tangent. Typing tangent(n) produces the n-th tangent number T_n. (8th June 2011)
70. bernoulli. Typing bernoulli(n) produces the n-th Bernoulli number. (9th June 2011)
71. kronecker. Typing kronecker(d,n) returns the Kronecker symbol (d,n), n ≥ 1.
72. 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)
73. 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.
74. 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.
75. 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.
76. zagier. zagier(d) gives the backward continued fraction of √d. zagier0(d,t,u,v) gives the least integer continued fraction of (u + t√d)/v.
77. 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.
78. pellab. Typing pellab(a,b,e) solves the diophantine equation ax²-by²=±1 and finds the least positive solution. We assume 1 < a < b, gcd(a,b)=1 and ab is not a perfect square.
79. squareroot2. Typing squares(n,flag) finds all representations of n as r²+s², r<s, gcd(r,s)=1. Here flag=1 is verbose and prints all parts of Euclid's algorithm as in the note. Typing squares1(n) only does that part of Euclid's algorithm as far as r_c-1 > √n ≥ r_c.
80. nagell_test. Typing nagell(d,n) finds the fundamental solutions, if any, of the diophantine equation x² - dy² = n, using the bounds given by T. Nagell in his book Introduction to Number Theory, 205-207.
81. markoff_triples. Typing test(n) prints all Markoff triples (x,y,z) with x < y < z and 5 ≤ z ≤ n. The recursion is due to Alan Offer.
82. patzpos. Typing patzpos(d,n,e), where e=0 or 1 finds the positive fundamental solutions of a generalised Pell equation x² - dy² = n, with e = 1 being prolix. The program needs squareroot2.
83. primepatterns. Typing primepatterns(a[],b[],s,u,v) will determine if a[0]m+b[0],...,a[s-1]m+b[s-1] form an admissible sequence and if so, will test the absolute values of these numbers for m in the range u ≤ m ≤ v for simultaneous primality. The program needs gcd, lucas, jacobi and a file containing the a[i] and b[i]. Here s ≤ 20.
84. posformrep. This solves the diophantine equation bt²+ctu+du²=n, where c²-4bd < 0, gcd(b,c,d)=1 and b > 0. Type bc posformrep squareroot reduceneg and then posrep(b,c,d,a,e), where e = 1 is verbose, e = 0 nonverbose. We use a transformation of variables due to Lagrange (1770) and then reduction to reduced form. See note.
85. posformrep2. This solves the diophantine equation bt²+ctu+du²=n, where c²-4bd < 0, gcd(b,c,d)=1 and b > 0. Type bc posformrep2 squareroot and then posrep(b,c,d,a,e), where e = 1 is verbose, e = 0 nonverbose. We use a transformation of variables due to Lagrange (1770) to get an equation Pu²+Quy+Ry²=1 and then perform Euclid's algorithm on -Q/2P. See paper.
86. stolt. This finds the fundamental solutions of the diophantine equation ax²+bxy+cy²=n, where d=b²-4ac > 0 is not a perfect square, gcd(a,b,c)=1 and a > 0. Type stolt(a,b,c,n). This is only useful for small d and an. Also typing fund(d,e), e = 0 or 1, finds the least positive solution of Pell's equation x² - dy²=4, with e=1 being verbose.
87. dujella_minus. The function negdujella1(k) finds all solutions of X²-(k²+1)Y,sup>2= -k² that lie in the range 1 < Y < k. The Dujella unicity conjecture is equivalent to there being at most one solution. If 2 ≤ m ≤ n, the function test(m,n) performs negdujella1(k) over the range k = m,...,n. We print the solutions (X,Y) and the corresponding solutions (x,y) under the mapping x = (k²+1)Y - kX, y = kY - X. To run the program one needs gcd and squareroot2.
88. phi now contains lprimrootneg(p), which finds the least negative primitive root mod p, lprimrootposmn(m,n), which lists the least primitive roots mod p for m ≤ p ≤ n, lprimrootnegmn(m,n), which lists the least negative primitive roots mod p for m ≤ p ≤ n. This need lucas and jacobi..
89. wildberger1. This solves Pell's equation x² - dy² = ±1 using an algorithm of Norman Wildberger. Type wild(d,e). e = 1 is verbose, e = 0 is nonverbose.
90. patz now contains patzgentest(d,n), which lists all fundamental solutions of x² - dy² = n and returns the their number. Type bc patz squareroot2 genfacs. Also patzgen creates a global array of these solutions. nagell_test contains ssw0(a,b,c,d,e,f), which solves the diophantine equation ax² + bxy + cy² + dx + ey +f = 0, where a ≠ 0, b² - 4ac > 0 is not a perfect square. These functions are accessed by typing bc patz squareroot2 genfacs nagell_test.
91. posformrep now contains posrepgen(a,b,c,n)(d,n), which finds all solutions (primitive and imprimitive) of ax²+bxy+cy²=n as a global array. These are listed by calling posrepgenlist(a,b,c,n). posrepgen(a,b,c,n) is used in ssw1(a,b,c,d,e,f) to solve the diophantine equation ax² + bxy + cy² + dx + ey +f = 0, where a ≠ 0, b² - 4ac < 0. These functions are accessed by typing bc posformrep squareroot2 genfacs reduceneg.
92. patz:. Type bc patz squareroot2 genfacs gcd. This now contains
    • generalizedpell(a,b,c,e). This solves the diophantine equation ax^2+bxy+cy^2=1, where a is nonzero, d=b^2-4ac>0 and is nonsquare. It produces the global solution class (generalizedpellfundx,generalizedpellfundy).
    • binary0(a,b,c,e). This finds the primitive solution classes of the diophantine equation ax² + bxy + cy² = n, where d = b² – 4ac > 0 is not a perfect square, gcd(a,b,c) = 1, gcd(a,n)=1 and n is nonzero. It produces the global solution classes (fund_x[i],fund_y[i]).
    • binary1(a,b,c,e). This finds the primitive solution classes of the diophantine equation ax² + bxy + cy² = n, where d = b² – 4ac > 0 is not a perfect square, gcd(a,b,c) = 1, a and n are nonzero. It produces the global solution classes (fund_x[i],fund_y[i]). binary(a,b,c,e) prints these solutions.
      binarygen(a,b,c,n,e) This finds all solutions (primitive and imprimitive) of ax^2+bxy+cy^2=n as a global array (globalbinarygenx[i],globalbinarygeny[i]). Here b² – 4ac > 0, s nonsquare, a and n are nonzero, gcd(a,b,c)=1. binarygenlist(a,b,c,n,e) prints these classes.
    • sswgeneral(a,b,c,d,e,f) . This function lives in patz and solves the diophantine equation ax² + bxy + cy² + dx + ey + f = 0, where (a,b,c) ≠ (0,0,0), using a method inspired by a lemma of John Robertson. One needs programs
      patz squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd ddzero-extra arithpartition.bc rob1
      This is one of the longest algorithms I have coded in BC and involved updating gcd, powerdd, genfacs, posformrep, reduceneg and constructing the files powertest3 and congruence. Files ddzero-extra and arithpartition.bc were added on December 10, 2019 and incorporate an improvement to the output when D = 0, due to Chi Chon Lei of Imperial College, London.
93. inversion. This contains inverse(m,n), which finds the inverse of m modulo n if it exists. The algorithm is from New algorithms for modular inversion and representation by binary quadratic forms arising from structure in the Euclidean algorithm, Christina Doran, Shen Lu and Barry R. Smith.
94. b2plus3bcplusc2. This contains rep(p), which expresses a prime of the form 5n+1 or 5n-1 as x² + xy - y² and x² + 3xy + y². The algorithm is from New algorithms for modular inversion and representation by binary quadratic forms arising from structure in the Euclidean algorithm, Christina Doran, Shen Lu and Barry R. Smith. Type bc lucas gcd b2plus3bcplusc2 jacobi squareroot2 then rep(p).
95. patz now contains aa1(a,b,c) and aa2(a,b,c,n). The former solves the diophantine equation ax²-by²=c, using the LMM method, while the latter solves the diophantine equation ax²+bxy+cy²=n, where b²-4ac > 0 and is not a perfect square.
96. stolt now contains stolt0(d,n,flag), which solves the diophantine equation x²-dy²=n, n ≠ 0 and a multiple of 4, d > 0 and nonsquare, using inequalities from a paper of Bengt Stolt.
97. phi now contains conwaycycles(a,b) and conwaysequence1(a,b,n). See BC program.
98. root3 contains cycle(x,y,n), which exhibits the trajectory (x,y), T(x,y), T²(x,y), Tⁿ(x,y), where T is the function defined in Example 10.2., page 101 of The Ultimate Challenge: The 3x+1 Problem, Ed. Jeffrey C. Lagarias, AMS January 2010. It seems likely that there are 103 cycles and that (x,y) = (-1,7) gives a divergent trajectory. The file gcd is needed.
99. partition.bc. Typing partition(m,1) prints the first m values of the partition function p(m).
    Typing partition(m,0) returns only p(m).
    Typing strict_partition(m,1) prints the first m values of the strict partition function q(m).
    Typing strict_partition(m,0) returns only q(m). Communicated by Thorsten Ehlers, May 3, 2019.
    See Computing partition numbers Q(N), an early version of a paper by Evangelos Georgiadis.
100. pell1 contains fund1(d,printflag), which finds the fundamental solution (globalr, globals) of Pell's equation x²-dy² = 1, using the midperiod approach and the continued fraction of √d.
101. pell4 contains fund4(d,printflag), which finds the fundamental solution (globalx, globaly) of Pell's equation x²-dy² = 4, using the midperiod approach and the continued fraction of √d.
102. sswgeneral(a,b,c,d,e,f) updated 27th January 2020. This function lives in patz and solves the diophantine equation ax² + bxy + cy² + dx + ey + f = 0, where (a,b,c) ≠ (0,0,0), using a method inspired by a lemma of John Robertson. One needs programs
     patz squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd ddzero-extra arithpartition.bc rob1 bigu stoltvialmm stoltarrange.bc
     .
103. stoltvialmm(d,n,printflag) lives in stoltvialmm and finds the Stolt fundamental solutions of x² – dy² = n, where n is divisible by 4.
     Needs stoltvialmm stoltarrange.bc squareroot patz genfacs
     binaryviasfs(a,b,c,n,printflag also lives in stoltvialmm and finds the stolt fundamental solutions of ax² +bxy + cy² = n. This has a wider applicability than the earlier version stolt(a,b,c,n).
104. bigu contains bigu(a,b,c,d,e,f,printflag), which solves the hyperbola case of the diophantine equation ax² + bxy + cy² + dx + ey + f = 0.
     Type bigu(a,b,c,d,e,f,printflag).
     Needs bc squareroot bigu stoltvialmm stoltarrange.bc patz genfacs.
105. squareroot now contains quadratic(a,b,c,n,printflag). This returns the number k of solutions of the congruence ax²+bx+c ≡ 0 (mod n), where a,n are positive . The solutions in the range 0 ≤ x < n are returned as global variables qglobalsol[0],…,qglobalsol[k-1]. printflag=1 prints the solutions. See note.
106. gcd now contains rationalization(a,b,c,d,dd), which rationalizes the denominator in (a + b√dd)/(c + d√ldd).
107. Added Ore's method of solving a system of linear congruences x ≡ a_i mod (m_i).
108. November 7, 2024 Added 3x+3korig, which performs the 3x+3^k algorithm. Type collatzk(y,k) where k ≥ 1 and y is an integer. The iterates of y will eventually lead to a multiple of 3^k, say 3^kz, at which point the trajectory becomes 3^k times the trajectory of z under the 3x+1 mapping.

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