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