# Common Number Theory Algorithms in Python

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

 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 iroot integer k’th root is_square check for an integer perfect square 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 lambdaf Carmichael function, min m for a^m == 1 (mod n) , all a co-prime to n mu the Mobius function jacobi the Jacobi symbol for an odd integer is_square_free check that an integer is not divisible by any square greater than 1 coprime_pairs Generate co-prime pairs pythag_primitive generate primitive Pythagorean triangles in order up to maximum side pythag generate all Pythagorean triangles in order up to maximum side rseq generate recursive sequences 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.

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!

• 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. 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.

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

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

• There is a nice description here.

5. 6. 