BC number theory programs

BC (version 1.06) (see the manual) is a language that supports arbitrary precision integer arithmetic calculations with interactive execution of statements. It is an idiosyncratic programming language which closely resembles parts of the C language. While it is somewhat slow for some purposes (all code is executed as it is read), the author has found bc very useful in teaching and research. The code is very easy to read and students can see directly how an algorithm has been implemented. It is also useful as a template for harder-to-write, but faster C code.
bc-1.06 is a superior version of the basic bc program which comes with the UNIX operating system.

To run a bc program such as gcd below, type bc gcd. This loads the bc program and the bc program gcd.

As a calculator, bc has a number of standard functions eg.

5+3, 5*3, 5/3, (=integer part), 5%3, (=least remainder), 5^3, sqrt(n).

global array variables can be entered: m[0]=5;m[1]=3;m[2]=7

Typing z=lcma(m[ ],3) and then z, prints and stores z, the lcm of 5, 3 and 7.

Consult the bc manual for more information.
Bugs and idiosyncrasies:

  1. When using the if-else construction, insert a "\", as follows: 
              if(expression){statement1}\
          else{ statement2}
    Alternatively, as pointed out by     Anton Stiglic, use the construction 
              if(expression){
               statement1
          }else{
               statement2
          }
  2. There is a bug in bc-1.06 under some Solaris platforms, which causes my program squareroot below to give a bus error, though not under linux. Hopefully the bug will not be present when bc-1.07 is released.
  3. The maximum number of auto variables allowed in a bc function seems to be 24!
  4. The continue statement only works for for loops.

BC PROGRAMS WRITTEN BY KEITH MATTHEWS

These are downloadable as a gzipped tar file.
They may contain bugs, both on the algorithmic and programming side.

  1. gcd:
  2. euclid: euclid(m,n) performs Euclid's algorithm.
  3. euclid1: euclid(m,n) returns the length of Euclid's algorithm.
  4. jacobi:
  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)
  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 x2-d*y2=±1,±2,±3 and least primitive solution of x2-d*y2=±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*2m+1, h < 2m, 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. 3branch: s(n) tests a 3-branched generalized 3x+1 conjecture.
  21. challenge: s(n) tests another 3-branched generalized 3x+1 conjecture.
  22. mordell: f(a,k) finds the integer solutions of y2=x3+a with x ≤ k. (Needs gcd.)
  23. venturini1: s(n) tests a 6-branched generalized 3x+1 function of G. Venturini.
  24. lra: (needs gcd)
  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]xn+···+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)
  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 gr ≡ n (mod p), 0 ≤ r < p-1.
    Note: p < 232-216=4294901760, in order to satisy BC array upper bound length of 216-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 np, the least quadratic nonresidue mod p. (Needs gcd).
  33. rootd_modn: rootd_modn(d,n) finds all solutions of the congruence x2 ≡ 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. u2 ≡ d (mod n), 1 < u < n. Then thue(d,u,p) finds x,y such that x2-dy2=kn, with small k.
  35. sqroot: sqroot(d,n), n > 1, finds all solutions of x2 ≡ 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 logba 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=xp and p is the least prime with this property.
  40. primes: primes(m,n) prints the primes in the interval [m,n], if ma and n lie between 1 and 1010.
  41. nprime: nprime(m) finds the least Lucas-base2 strong pseudoprime p satisfying p ≥ m.
  42. 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.
  43. sturm: sturm(a[],n,b,e) prints a sturm polynomial sequence for the (squarefree) polynomial a[n]xn+···+a[0], evaluates the sequence at x=b and calculates its sign variation. e=0 suppresses printing.
  44. 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.
  45. reducepos: reduce(a,b,c) takes as input an indefinite binary quadratic form ax2+bxy+cy2 and uses the PQa algorithm to find a cycle of reduced forms. (See explanation.)
  46. 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 squarefree d in the range m ≤ d < n, where n<106.
  47. reduceneg: reduce(a,b,c) takes as input a positive definite binary quadratic form ax2+bxy+cy2 and uses an algorithm of Gauss to find the equivalent reduced form and unimodular transforming matrix.
  48. classnopos: class_number(d) (1 < d < 106 and squarefree) finds the class number of the real quadratic field Q() 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 < 106.
    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 x2-d*y2=-4 is determined.
    This is basically the same program as class_number(d), except that in the case of non-solubility of x2-d*y2=-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.
  49. unimodular: unimodular(p,q,r,s) expresses a unimodular matrix A ≠ I2 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 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.
    The number n of matrices in the product U0 Un-1 is returned.
  50. 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 .
  51. factorial: factorial_p(n,p) finds the power of a prime p dividing n!
    factorial(n) finds n!
  52. 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 b2=a (mod p) and 0 < b < p.
  53. 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 Ua=[a,1;1,0], note that Ua0...Ua2n=Ra0La1...Ra2nU0 and that Ua0...Ua2n+1=Ra0La1...La2n+1I2.
  54. davison: davison(l,m,n) performs the algorithm of J.L. Davison's paper An algorithm for the continued fraction of el/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 A0...Ak, for n* ≤ k ≤ n*+n.
    The number (count) of partial quotients of el/m found is returned. count becomes positive for all large n.
  55. squareroot: This is an improved version of sqroot and contains cornacchia(a,b,m). This finds all positive primitive solutions of ax2+by2=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 x2=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.
  56. 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.
  57. patz: patz(d,n,e) finds fundamental solutions for the diophantine equation x2-dy2=n, d > 0, d not a perfect square. We only find the fundamental solutions for classes P satisfying P2 ≡ 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.
  58. squareroot now contains quadratic(a,b,c,n,flag). This returns the number k of solutions of the congruence ax2+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.
  59. patz now contains binary(a,b,c,n). This solves ax2+bxy+cy2=n, where n is non-zero and b2-4ac is positive and not a perfect square. The method is from the paper The Diophantine equation ax2+bxy+cy2=N, D=b2-4ac > 0, Journal de Théorie des Nombres de Bordeaux, 14 (2002) 257-270 by K.R. Matthews.
  60. 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.
  61. sqrt_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.
  62. nipell. (a) Typing nipell(61,0) finds the smallest solution of Pell equation x2-61y2=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.
  63. nscf_pell. This program solves the Pell equations x2 - dy2 = ±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.

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

Last modified 11th March 2008