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.
9 Comments
Leave one →
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.
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
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.
Also, could you explain the output of PQa()?
There is a nice description here.
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.
Thank you for the information on your implementation of a high reliability Miller-Rabin alternative.