mathematicskit.number_theory#

mathematicskit.number_theory hub · History · Examples · API reference

The extended Euclidean algorithm and modular inverses; fast modular exponentiation; primality testing (trial division, Miller-Rabin) and prime generation (sieve of Eratosthenes); the Chinese Remainder Theorem; continued-fraction expansion and best rational approximations; Euler’s totient function and other multiplicative functions (Mobius, divisor-sum); and linear and Pell Diophantine equation solvers. Every algorithm is hand-rolled – exact-integer number theory has no numpy/scipy equivalent.

mathematicskit.number_theory: elementary and computational number theory.

Every algorithm here is hand-rolled from its textbook definition – exact-integer arithmetic has no numpy/scipy equivalent. The extended Euclidean algorithm and modular inverses; fast modular exponentiation; primality testing (trial division, Miller-Rabin) and prime generation (sieve of Eratosthenes); the Chinese Remainder Theorem; continued-fraction expansion and best rational approximations; Euler’s totient function and other multiplicative number-theoretic functions (Mobius, divisor-sum); and linear and Pell Diophantine equation solvers; the Lucas-Lehmer test for Mersenne primes; Pollard’s rho factorization; Legendre and Jacobi symbols and Tonelli-Shanks modular square roots; sums of two and four squares; Euler’s product for the zeta function (checked against scipy.special.zeta); and the prime-counting function, logarithmic integral (scipy.special.expi), and primes in arithmetic progressions.

class mathematicskit.number_theory.BezoutResult(gcd, x, y)[source]#

Bases: object

Container for the extended Euclidean algorithm’s output: Bezout’s identity.

Parameters:
gcd: int#

\(\gcd(a, b)\), always non-negative.

Type:

int

x: int#

The Bezout coefficient of a in \(ax + by = \gcd(a,b)\).

Type:

int

y: int#

The Bezout coefficient of b in \(ax + by = \gcd(a,b)\).

Type:

int

class mathematicskit.number_theory.CRTResult(residue, modulus)[source]#

Bases: object

Container for the Chinese Remainder Theorem’s combined solution.

Parameters:
modulus: int#

The combined modulus (the product of the pairwise-coprime input moduli).

Type:

int

residue: int#

The unique solution modulo modulus.

Type:

int

class mathematicskit.number_theory.ContinuedFractionResult(terms, convergents)[source]#

Bases: object

Container for a continued-fraction expansion and its convergents.

Parameters:
convergents: list[tuple[int, int]]#

Successive \((p_k, q_k)\) best rational approximations \(p_k/q_k\), in the same order as terms.

Type:

list of (int, int)

terms: list[int]#

\([a_0; a_1, a_2, \dots]\).

Type:

list of int

class mathematicskit.number_theory.LinearDiophantineResult(has_solution, x0=None, y0=None, gcd=0, x_step=0, y_step=0)[source]#

Bases: object

Container for a linear Diophantine equation’s solution family.

Parameters:
gcd: int = 0#

\(\gcd(a, b)\), or 0 in the degenerate a == b == 0 case, where there is no gcd to speak of.

Type:

int

has_solution: bool#

Whether \(ax + by = c\) has an integer solution (equivalently, whether \(\gcd(a,b)\) divides c).

Type:

bool

x0: int | None = None#

The x of one particular solution, or None when has_solution is False.

Type:

int, optional

x_step: int = 0#

\(b/\gcd(a,b)\), the amount x advances per step along the solution family.

Type:

int

y0: int | None = None#

The y of that same particular solution.

Type:

int, optional

y_step: int = 0#

\(a/\gcd(a,b)\). Together with x_step, the general solution is (x0 + k*x_step, y0 - k*y_step) for any integer k; note the minus sign on the y component.

Type:

int

class mathematicskit.number_theory.PellResult(x, y, d, extra=<factory>)[source]#

Bases: object

Container for the fundamental solution of Pell’s equation \(x^2 - Dy^2 = 1\).

Parameters:
d: int#

The D in \(x^2 - Dy^2 = 1\).

Type:

int

extra: dict#

Free-form diagnostics slot, unused by the current solver.

Type:

dict

x: int#

The x of the fundamental (smallest positive) solution.

Type:

int

y: int#

The y of the fundamental solution.

Type:

int

class mathematicskit.number_theory.PollardRhoResult(factor, cofactor, iterations, c)[source]#

Bases: object

Container for a factor found by Pollard’s rho method.

Parameters:
c: int#

The additive constant of \(x \mapsto x^2 + c\) that succeeded.

Type:

int

cofactor: int#

n // factor.

Type:

int

factor: int#

A nontrivial factor of n (not necessarily prime).

Type:

int

iterations: int#

Total tortoise-and-hare steps taken, across all restarts.

Type:

int

mathematicskit.number_theory.best_rational_approximation(x, max_denominator)[source]#

Best rational approximation to x with denominator \(\leq\) max_denominator.

Expands x’s continued fraction and takes the last convergent \(p_k/q_k\) with \(q_k \leq\) max_denominator, then compares it against the best semiconvergent (or “intermediate fraction”) \(\dfrac{p_{k-1} + t\,p_k}{q_{k-1} + t\,q_k}\) that still fits under the denominator bound, returning whichever is closer to x.

Both candidates are needed: a convergent is only guaranteed to be the best approximation among fractions with denominator at most its own \(q_k\), and the next convergent’s denominator can jump far past max_denominator, leaving room for a semiconvergent in between. For \(x=\pi\) with max_denominator=57, for instance, the last convergent is \(22/7\) but the semiconvergent \(179/57\) (between \(22/7\) and \(333/106\)) is genuinely closer. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Theorem 7.13 and Sec. 7.4.

Parameters:

  • x (float)

  • max_denominator (int) – >= 1.

Return type:

tuple

Returns:

(int, int) – (p, q) with q <= max_denominator.

Examples

>>> # The classic approximation pi ~ 355/113 (denominator <= 200).
>>> best_rational_approximation(3.14159265358979, max_denominator=200)
(355, 113)
>>> # Under 57, the best fraction is the semiconvergent 179/57, not 22/7.
>>> best_rational_approximation(3.14159265358979, max_denominator=57)
(179, 57)
>>> abs(179 / 57 - 3.14159265358979) < abs(22 / 7 - 3.14159265358979)
True
mathematicskit.number_theory.chinese_remainder_theorem(remainders, moduli)[source]#

Solve the system \(x \equiv r_i \pmod{m_i}\) for pairwise-coprime \(m_i\).

Combines the congruences two at a time: given a solution \(x \equiv r_1 \pmod{m_1}\) and a new congruence \(x \equiv r_2 \pmod{m_2}\), the combined solution modulo \(m_1 m_2\) is found via the extended Euclidean algorithm’s Bezout coefficients for \((m_1, m_2)\). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Theorem 2.18.

Parameters:

  • remainders (Sequence[int])

  • moduli (Sequence[int]) – Pairwise coprime, i.e. gcd(moduli[i], moduli[j]) == 1 for every i != j.

Return type:

CRTResult

Returns:

CRTResult

Examples

>>> # x = 2 mod 3, x = 3 mod 5, x = 2 mod 7 -> x = 23 mod 105 (a classic example).
>>> result = chinese_remainder_theorem([2, 3, 2], [3, 5, 7])
>>> result.residue
23
>>> result.modulus
105
mathematicskit.number_theory.continued_fraction_expansion(x, max_terms=20, tol=1e-10)[source]#

Expand x as a (simple) continued fraction \([a_0; a_1, a_2, \dots]\).

\(a_0 = \lfloor x\rfloor\), then repeatedly takes the reciprocal of the remaining fractional part and floors it again, stopping early once the fractional part is within tol of zero (an exact rational input) or after max_terms. Each successive convergent \(p_k/q_k\) (via the standard recurrence \(p_k = a_k p_{k-1} + p_{k-2}\), \(q_k = a_k q_{k-1} + q_{k-2}\)) is the best rational approximation to x among all fractions with denominator \(\leq q_k\). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 7.1-7.4.

Parameters:

  • x (float)

  • max_terms (int)

  • tol (float) – Stop expanding once the remaining fractional part is smaller than this (exact termination for rational x).

Return type:

ContinuedFractionResult

Returns:

ContinuedFractionResult

Examples

>>> result = continued_fraction_expansion(3.245, max_terms=10)
>>> result.terms[:3]
[3, 4, 12]
>>> p, q = result.convergents[-1]
>>> abs(p / q - 3.245) < 1e-9
True
>>> # The golden ratio's continued fraction is famously all 1s.
>>> phi = (1 + 5**0.5) / 2
>>> continued_fraction_expansion(phi, max_terms=8).terms
[1, 1, 1, 1, 1, 1, 1, 1]
mathematicskit.number_theory.divisor_sum(n, power=1)[source]#

The divisor-power-sum function \(\sigma_k(n) = \sum_{d \mid n} d^k\).

power=0 gives the number of divisors \(d(n)\); power=1 (the default) gives the ordinary sum of divisors \(\sigma(n)\) (a perfect number satisfies \(\sigma(n) = 2n\)). Multiplicative, computed from prime_factorization() via \(\sigma_k(p^e) = \sum_{j=0}^{e} p^{jk}\). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 4.2.

Parameters:

  • n (int) – n >= 1.

  • power (int)

Return type:

int

Returns:

int

Examples

>>> divisor_sum(6)  # 1 + 2 + 3 + 6 = 12 = 2*6: a perfect number
12
>>> divisor_sum(28)  # also perfect
56
>>> divisor_sum(12, power=0)  # 1, 2, 3, 4, 6, 12: 6 divisors
6
mathematicskit.number_theory.euler_product(s, limit)[source]#

Truncated Euler product \(\prod_{p \le N} (1 - p^{-s})^{-1}\).

Euler (1737) showed that for \(s > 1\)

\[\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s} = \prod_{p\ \mathrm{prime}} \frac{1}{1 - p^{-s}},\]

which is unique factorization written analytically. Truncating the product at primes \(p \le N\) converges to \(\zeta(s)\) as \(N \to \infty\); at \(s = 1\) it diverges like \(e^{\gamma}\ln N\) (Mertens’ third theorem), reproving that there are infinitely many primes.

Parameters:

  • s (float) – Exponent, s > 0.

  • limit (int) – Include primes up to and including limit.

Return type:

float

Returns:

float

Examples

>>> import math
>>> from scipy.special import zeta
>>> abs(euler_product(2.0, 10**5) - math.pi**2 / 6) < 1e-5  # the Basel problem
True
>>> bool(abs(euler_product(3.0, 1000) - zeta(3.0)) < 1e-6)
True
mathematicskit.number_theory.euler_totient(n)[source]#

Euler’s totient \(\varphi(n)\): count of integers in [1, n] coprime to n.

Multiplicative: \(\varphi(n) = n\prod_{p \mid n}(1 - 1/p)\) over n’s distinct prime factors (from prime_factorization()). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 4.2, Theorem 4.11.

Parameters:

n (int) – n >= 1.

Return type:

int

Returns:

int

Examples

>>> euler_totient(1)
1
>>> euler_totient(9)
6
>>> euler_totient(17)  # prime: phi(p) = p - 1
16
mathematicskit.number_theory.extended_gcd(a, b)[source]#

Extended Euclidean algorithm: find \(\gcd(a,b)\) and Bezout coefficients.

Returns integers \(g, x, y\) with \(ax + by = g = \gcd(a,b)\), computed by unwinding the ordinary Euclidean algorithm’s recursion. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 1.2, Theorem 1.3.

Parameters:
Return type:

BezoutResult

Returns:

BezoutResult

Examples

>>> result = extended_gcd(240, 46)
>>> result.gcd
2
>>> 240 * result.x + 46 * result.y == result.gcd
True
mathematicskit.number_theory.fast_mod_pow(base, exponent, modulus)[source]#

Fast modular exponentiation via repeated squaring: \(\text{base}^{\text{exponent}} \bmod \text{modulus}\).

\(O(\log(\text{exponent}))\) multiplications, by writing the exponent in binary and squaring the running result once per bit (multiplying it in whenever that bit is 1) – versus the naive \(O(\text{exponent})\) repeated multiplication. Equivalent to Python’s built-in three-argument pow(base, exponent, modulus) (used as a cross-check in this module’s tests), reimplemented here since exposing the square-and-multiply algorithm itself is the point. See Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 31.6.

Parameters:

  • base (int) – exponent >= 0.

  • exponent (int) – exponent >= 0.

  • modulus (int) – exponent >= 0.

Return type:

int

Returns:

int

Examples

>>> fast_mod_pow(7, 128, 13)
3
>>> fast_mod_pow(2, 10, 1000)
24
mathematicskit.number_theory.gcd(a, b)[source]#

Greatest common divisor via the Euclidean algorithm.

\(\gcd(a,b) = \gcd(b, a \bmod b)\), terminating when the second argument reaches 0. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 1.2.

Parameters:
Return type:

int

Returns:

int – Non-negative.

Examples

>>> gcd(48, 18)
6
>>> gcd(17, 5)
1
mathematicskit.number_theory.is_prime_miller_rabin(n, k=40, seed=0)[source]#

Probabilistic primality test (Miller-Rabin).

Writes \(n-1 = 2^r d\) with d odd, and for each of k random witnesses a checks whether \(a^d \equiv 1\) or \(a^{2^i d} \equiv -1 \pmod n\) for some \(0\leq i<r\); if neither holds, n is certainly composite (”a is a witness to compositeness”), otherwise n is declared probably prime. Each round has failure probability \(\leq 1/4\) for composite n, so \(k=40\) rounds give false-positive probability \(\leq 4^{-40}\) – negligible in practice, but this remains a probabilistic (not certificate) test, unlike is_prime_trial_division(). See Cormen et al., Introduction to Algorithms, 3rd ed., Ch. 31.8.

Parameters:

  • n (int)

  • k (int) – Number of random witnesses tested.

  • seed (int) – Random seed for witness selection.

Return type:

bool

Returns:

bool

Examples

>>> is_prime_miller_rabin(97)
True
>>> is_prime_miller_rabin(91)  # 91 = 7 * 13
False
>>> is_prime_miller_rabin(2**61 - 1)  # a known Mersenne prime
True
>>> is_prime_miller_rabin(2**89 - 1)  # arbitrary precision: no 64-bit ceiling
True
mathematicskit.number_theory.is_prime_trial_division(n)[source]#

Primality test by trial division up to \(\sqrt n\).

\(O(\sqrt n)\): a composite n must have a factor \(\leq\sqrt n\), so it suffices to check candidate divisors up to there (only 2 and odd numbers, after handling 2 separately). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 1.1.

Parameters:

n (int)

Return type:

bool

Returns:

bool

Examples

>>> [k for k in range(2, 30) if is_prime_trial_division(k)]
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
mathematicskit.number_theory.jacobi_symbol(a, n)[source]#

The Jacobi symbol \(\left(\frac{a}{n}\right)\) for odd n > 0.

The multiplicative extension of the Legendre symbol to odd composite moduli. Computed without factoring n by the Euclid-like algorithm that quadratic reciprocity makes possible: pull out factors of 2 with \(\left(\frac{2}{n}\right) = (-1)^{(n^2-1)/8}\) and flip the symbol with \(\left(\frac{a}{n}\right) \left(\frac{n}{a}\right) = (-1)^{\frac{a-1}{2}\frac{n-1}{2}}\). See Cohen, A Course in Computational Algebraic Number Theory, Algorithm 1.4.10.

Parameters:

  • a (int)

  • n (int) – Odd, n > 0.

Return type:

int

Returns:

int – One of -1, 0, 1.

Examples

>>> jacobi_symbol(1001, 9907)
-1
>>> jacobi_symbol(2, 15)  # = (2/3)(2/5) = (-1)(-1), yet 2 is not a square mod 15
1
mathematicskit.number_theory.lcm(a, b)[source]#

Least common multiple: \(\mathrm{lcm}(a,b) = |ab|/\gcd(a,b)\).

Parameters:
Return type:

int

Returns:

int

Examples

>>> lcm(4, 6)
12
>>> lcm(21, 6)
42
mathematicskit.number_theory.legendre_symbol(a, p)[source]#

The Legendre symbol \(\left(\frac{a}{p}\right)\) via Euler’s criterion.

For an odd prime p, \(\left(\frac{a}{p}\right) \equiv a^{(p-1)/2} \pmod p\) is 1 when a is a nonzero square modulo p, -1 when it is not, and 0 when p divides a. See Niven, Zuckerman & Montgomery, 5th ed., Sec. 3.1, Theorem 3.1.

Parameters:

  • a (int)

  • p (int) – An odd prime (not checked).

Return type:

int

Returns:

int – One of -1, 0, 1.

Examples

>>> [legendre_symbol(a, 7) for a in range(7)]  # squares mod 7: 1, 2, 4
[0, 1, 1, -1, 1, -1, -1]
mathematicskit.number_theory.logarithmic_integral(x)[source]#

The logarithmic integral \(\operatorname{li}(x) = \int_0^x \frac{dt}{\ln t}\).

The (Cauchy principal value) integral equals \(\operatorname{Ei}(\ln x)\), evaluated with scipy.special.expi(). Gauss’s conjectured approximation to \(\pi(x)\): the prime number theorem (Hadamard and de la Vallée Poussin, 1896) states \(\pi(x) \sim \operatorname{li}(x) \sim x/\ln x\).

Parameters:

x (float or array_like) – x > 0.

Returns:

float or ndarray

Examples

>>> round(float(logarithmic_integral(2.0)), 6)  # li(2), the Ramanujan-Soldner offset
1.045164
>>> round(float(logarithmic_integral(1e6)))  # vs. pi(10**6) = 78498
78628
mathematicskit.number_theory.lucas_lehmer(p)[source]#

Lucas-Lehmer test: is the Mersenne number \(M_p = 2^p - 1\) prime?

For an odd prime p, set \(s_0 = 4\) and \(s_{k+1} = s_k^2 - 2 \bmod M_p\); then \(M_p\) is prime if and only if \(s_{p-2} \equiv 0 \pmod{M_p}\) (Lucas, 1878; Lehmer, 1930). Only p - 2 modular squarings are needed, so the test is a deterministic proof of primality far faster than any general method – it is why the largest known primes are almost all Mersenne primes. See Crandall & Pomerance, Prime Numbers: A Computational Perspective, 2nd ed., Theorem 4.2.6.

Parameters:

p (int) – Exponent, p >= 2. A composite p gives a composite \(M_p\) and returns False immediately; p = 2 (where \(M_2 = 3\)) is handled as a special case.

Return type:

bool

Returns:

bool

Examples

>>> [p for p in range(2, 130) if lucas_lehmer(p)]  # exponents of Mersenne primes
[2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127]
mathematicskit.number_theory.mobius(n)[source]#

The Mobius function \(\mu(n)\).

\(\mu(1)=1\); \(\mu(n)=0\) if n has any squared prime factor; otherwise \(\mu(n)=(-1)^k\) for k distinct prime factors. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 4.2.

Parameters:

n (int) – n >= 1.

Return type:

int

Returns:

int – One of -1, 0, 1.

Examples

>>> [mobius(k) for k in range(1, 11)]
[1, -1, -1, 0, -1, 1, -1, 0, 0, 1]
mathematicskit.number_theory.mod_inverse(a, m)[source]#

Modular inverse of a modulo m: the x with \(ax \equiv 1 \pmod m\).

Exists iff \(\gcd(a,m)=1\), found via extended_gcd(). See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 2.1.

Parameters:
Return type:

int

Returns:

int – In the range [0, m).

Raises:

ValueError – If gcd(a, m) != 1 (no inverse exists).

Examples

>>> mod_inverse(3, 11)
4
>>> (3 * mod_inverse(3, 11)) % 11
1
mathematicskit.number_theory.pollard_rho(n, c=1, x0=2, max_iter=10000000)[source]#

Find a nontrivial factor of composite n with Pollard’s rho method.

Iterates the pseudo-random map \(x \mapsto x^2 + c \bmod n\). Reduced modulo an unknown prime factor \(p\) of n, the sequence must repeat within about \(\sqrt{\pi p/2}\) steps (the birthday paradox), tracing out the Greek letter \(\rho\). Floyd’s tortoise-and-hare cycle detection finds that repeat by comparing \(x_i\) with \(x_{2i}\) and testing \(\gcd(|x_i - x_{2i}|, n)\). The expected cost is \(O(\sqrt p) \le O(n^{1/4})\) steps, compared with \(O(\sqrt n)\) for trial division. If the cycle closes modulo every factor at once (the gcd is n itself), the constant c is incremented and the search restarted.

Parameters:

  • n (int) – Composite, n >= 4. An even n returns the factor 2 immediately.

  • c (int) – Initial additive constant of the iteration map.

  • x0 (int) – Starting value.

  • max_iter (int) – Total iteration budget across all restarts.

Return type:

PollardRhoResult

Returns:

PollardRhoResult

Examples

>>> result = pollard_rho(8051)  # Pollard's own worked example: 83 * 97
>>> sorted([result.factor, result.cofactor])
[83, 97]
>>> result = pollard_rho(2**64 + 1)  # Landry, 1880: 274177 * 67280421310721
>>> sorted([result.factor, result.cofactor])
[274177, 67280421310721]
mathematicskit.number_theory.prime_counting(x)[source]#

The prime-counting function \(\pi(x)\): the number of primes \(\le x\).

Sieves once up to max(x) and counts with numpy.searchsorted(), so a whole array of x costs a single sieve.

Parameters:

x (int or array_like of int) – Non-negative arguments.

Return type:

ndarray

Returns:

int or ndarray of int – Same shape as x.

Examples

>>> int(prime_counting(100))
25
>>> prime_counting([10, 100, 1000, 10000])
array([   4,   25,  168, 1229])
mathematicskit.number_theory.prime_factorization(n)[source]#

Prime factorization of n via trial division: \(n = \prod_i p_i^{e_i}\).

Parameters:

n (int) – n >= 1.

Return type:

dict

Returns:

dict – {prime: exponent}, empty for n == 1.

Examples

>>> prime_factorization(360)
{2: 3, 3: 2, 5: 1}
mathematicskit.number_theory.primes_in_progression(a, q, limit)[source]#

Primes \(p \le\) limit with \(p \equiv a \pmod q\).

Dirichlet’s theorem (1837) guarantees infinitely many such primes whenever \(\gcd(a, q) = 1\); the prime number theorem for arithmetic progressions sharpens this to each of the \(\varphi(q)\) coprime residue classes receiving an asymptotic share \(1/\varphi(q)\) of all primes.

Parameters:

  • a (int) – Residue class.

  • q (int) – Modulus, q >= 1.

  • limit (int)

Return type:

ndarray

Returns:

ndarray of int – The primes in [2, limit] congruent to a mod q, ascending.

Examples

>>> primes_in_progression(3, 4, 50)
array([ 3,  7, 11, 19, 23, 31, 43, 47])
>>> primes_in_progression(2, 4, 1000)  # gcd(2, 4) = 2: only the prime 2
array([2])
mathematicskit.number_theory.sieve_of_eratosthenes(limit)[source]#

All primes up to and including limit, via the sieve of Eratosthenes.

Marks composites by striking out multiples of each prime found, from 2 upward – \(O(n\log\log n)\), far faster than testing each number individually. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 1.1.

Parameters:

limit (int)

Return type:

ndarray

Returns:

ndarray, int – Primes in [2, limit], ascending.

Examples

>>> sieve_of_eratosthenes(30)
array([ 2,  3,  5,  7, 11, 13, 17, 19, 23, 29])
mathematicskit.number_theory.solve_linear_diophantine(a, b, c)[source]#

Solve \(ax + by = c\) in integers.

A solution exists iff \(\gcd(a,b) \mid c\); given a particular solution \((x_0, y_0)\) (scaled up from extended_gcd()’s Bezout coefficients), the general solution is \((x_0 + k\,b/g,\ y_0 - k\,a/g)\) for any integer k. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Theorem 2.9.

Parameters:
Return type:

LinearDiophantineResult

Returns:

LinearDiophantineResult

Examples

>>> result = solve_linear_diophantine(3, 5, 1)
>>> 3 * result.x0 + 5 * result.y0 == 1
True
>>> # Every (x0 + k*x_step, y0 - k*y_step) is also a solution.
>>> k = 7
>>> x, y = result.x0 + k * result.x_step, result.y0 - k * result.y_step
>>> 3 * x + 5 * y == 1
True
>>> solve_linear_diophantine(2, 4, 3).has_solution  # gcd(2,4)=2 does not divide 3
False
>>> solve_linear_diophantine(0, 0, 5).has_solution  # 0x + 0y = 5 is unsolvable
False
>>> solve_linear_diophantine(0, 0, 0).has_solution  # ... but 0x + 0y = 0 holds for every (x, y)
True
mathematicskit.number_theory.solve_pell_equation(d)[source]#

Fundamental solution of Pell’s equation \(x^2 - Dy^2 = 1\), D not a perfect square.

The continued fraction of \(\sqrt D\) is eventually periodic; the fundamental (smallest positive) solution is given by one of its convergents \(p/q\). Unlike continued_fraction_expansion() (which works from a floating-point value), this function tracks the continued fraction of \(\sqrt D\) with the standard exact-integer recurrence for quadratic irrationals, since Pell periods can run to dozens of terms (e.g. \(D=61\) has period 11) – far beyond float64’s precision. See Niven, Zuckerman & Montgomery, An Introduction to the Theory of Numbers, 5th ed., Sec. 7.8, Theorem 7.26.

Parameters:

d (int) – Not a perfect square.

Return type:

PellResult

Returns:

PellResult

Examples

>>> result = solve_pell_equation(2)
>>> (result.x, result.y)
(3, 2)
>>> result.x**2 - 2 * result.y**2
1
>>> result2 = solve_pell_equation(61)  # a famously large fundamental solution
>>> result2.x**2 - 61 * result2.y**2
1
mathematicskit.number_theory.sqrt_mod(a, p)[source]#

A square root of a modulo an odd prime p (Tonelli-Shanks).

Writes \(p - 1 = 2^s q\) with q odd, starts from the candidate \(a^{(q+1)/2}\), and repeatedly corrects it with powers of a quadratic non-residue until the error term has order 1 (Tonelli, 1891; Shanks, 1973). Runs in \(O(\log^2 p)\) modular multiplications. See Cohen, Algorithm 1.5.1.

Parameters:

  • a (int) – A quadratic residue modulo p (legendre_symbol(a, p) != -1).

  • p (int) – An odd prime.

Return type:

int

Returns:

int – The root x in [0, p // 2]; the other root is p - x.

Examples

>>> sqrt_mod(10, 13)  # 6*6 = 36 = 10 (mod 13)
6
>>> x = sqrt_mod(2, 2**61 - 1)
>>> x * x % (2**61 - 1)
2
mathematicskit.number_theory.sum_of_four_squares(n)[source]#

Write n as \(a^2 + b^2 + c^2 + d^2\) (Lagrange’s four-square theorem).

Lagrange proved in 1770 that every non-negative integer is a sum of four squares, so this function always succeeds. The search tries \(a \ge b \ge c\) in decreasing order and solves for d; for n up to a few million it returns almost immediately. The number of ordered, signed representations is given by Jacobi’s (1829) formula \(r_4(n) = 8\sum_{d \mid n,\ 4 \nmid d} d\).

Parameters:

n (int) – n >= 0.

Return type:

tuple[int, int, int, int]

Returns:

tuple of (int, int, int, int) – (a, b, c, d) with a >= b >= c >= d >= 0 and squares summing to n.

Examples

>>> sum_of_four_squares(7)  # 7 = 4 + 1 + 1 + 1 needs all four squares
(2, 1, 1, 1)
>>> sum_of_four_squares(310)
(17, 4, 2, 1)
mathematicskit.number_theory.sum_of_two_squares(n)[source]#

Write n as \(a^2 + b^2\) with \(0 \le a \le b\), if possible.

By Fermat’s two-squares theorem (stated 1640, first proved by Euler in 1749) an odd prime \(p\) is a sum of two squares exactly when \(p \equiv 1 \pmod 4\); more generally n is a sum of two squares exactly when every prime \(q \equiv 3 \pmod 4\) divides n to an even power (Hardy & Wright, 6th ed., Theorem 366). The search runs over \(a \le \sqrt{n/2}\) – \(O(\sqrt n)\) – and returns the representation with the smallest a.

Parameters:

n (int) – n >= 0.

Return type:

tuple[int, int] | None

Returns:

tuple of (int, int) or None – (a, b) with a*a + b*b == n and a <= b, or None when no representation exists.

Examples

>>> sum_of_two_squares(13)  # 13 = 1 (mod 4): 4 + 9
(2, 3)
>>> sum_of_two_squares(7) is None  # 7 = 3 (mod 4)
True
>>> sum_of_two_squares(2**31 - 1) is None  # a Mersenne prime is always 3 (mod 4)
True