Skip to content

Common Number Theory Algorithms in Python

My number theory library provides the following common number theory functions and algorithms:

     gcd greatest common divisor
lcm least common multiple
xgcd     extended greatest common divisor
cong solve a.x = p (mod n) for x
mod_inverse inverse of a (mod n)
crt chinese remainder theorem
isqrt integer square root
iroot integer k’th root 
miller_rabin   Miller Rabin probable prime test 
Sieve a class for the sieve of Eratosthenes
Primes a class for working with primes
is_prime test if a positive integer is prime
factor prime factors and their multiplicity 
factors prime factors (with repeated primes)
divisors divisors of an integer
divisor_pairs divisor pairs of an integer
tau number of divisors of an integer
omega sum of the divisors of an integer
phi the Euler totient function
mu the Mobius function
jacobi the Jacobi symbol for an odd integer
pythag_primitive generate primitive Pythagorean triangles in order up to maximum side
pythag generate all Pythagorean triangles in order up to maximum side
hensel_lift_2 Hensel root lifting for powers of 2
hensel_lift Hensel root lifting for powers of odd primes
sqrt_mod_p the modular square root for a prime modulus
sqrt_mod_pn   the modular square root for prime power modulus
sqrt_mod_m the modular square root for a composite modulus
PQa continued fraction expansion of a quadratic irrational
Qde class for solving x^2 – d.y^2 = n
frobenius_number   the Frobenius number for an integer sequence
frobenius_solve Frobenius solutions for an integer sequence

Please note that these routines are still under development so they may contain bugs

My polynomial library provides facilities for manipulating polynomials with rational coefficients (this is used in the Funny Dice Sunday Times Teaser).

9 Comments Leave one →
  1. This is cool code. I finally got around to trying the Frobenius function and it works great!… as long as at least two of the items are coprime. When they all have a common factor > 1, the funcntion returns an answer that doesn’t make sense to me. As far as I can tell this is not the result of anything I did but I wouldn’t rule it out.

    Here are a few examples of my program using your function:

    d:dev>rpn [ 6 9 20 ] frobenius
    43

    d:dev>rpn [ 6 9 24 ] frobenius
    3

    d:dev>rpn [ 6 8 13 ] frobenius
    23

    d:dev>rpn [ 6 8 12 ] frobenius
    10

    I haven’t had a chance to go through the function to see if I can figure out why it does this. It will probably take me some time to understand it well enough, but I thought I would drop you a note.

    p.s. I’m looking forward to implementing operators for some of the other functions as well.
    Now I am off to McDonald’s to order 43 Chicken McNuggets!

    • Brian Gladman permalink

      Frobenius numbers are only defined for integer sequenes for which:

      \[{gcd}(a_1, a_2, …, a_n) = 1\]

      You can add a check by changing the last line in the frobenius_number function with the following:

      Alternatively add this before the last line:

      I have been updating other parts of the library today so I have added this check as well.

  2. I figured if input with a gcd > 1 was disallowed for some reason, I would make the check myself and throw an error.

    Thanks again. I still want to dig through and understand how it works. I’m still a number theory newbie.

    Rick

  3. Lucas A. Brown permalink

    The sqrtmod code is incorrect — specifically, for each m divisible by 3 or more distinct primes, sqrt_mod_m(a, m) returns incorrect results for a handful of values of a.

    • Brian Gladman permalink

      Thanks for reporting the bug, which I have now fixed.

  4. Lucas A. Brown permalink

    Also, could you explain the output of PQa()?

  5. Jurjen permalink

    For an implementation of Martin Seysen’s EQFT test (probable prime test with very high reliability), see eqft.py. This implementation includes the full SQFT3 test.

    I Hope this is useful.

  6. Brian Gladman permalink

    Thank you for the information on your implementation of a high reliability Miller-Rabin alternative.

Leave a Reply

Note: HTML is allowed. Your email address will not be published.

Subscribe to this comment feed via RSS