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