> For the complete documentation index, see [llms.txt](https://docs.algenta.ai/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://docs.algenta.ai/numerical-and-mathematical-computing/polynomial-root-finding-and-number-theory-kernels.md).

# Polynomial, root-finding & number theory

Dense-coefficient polynomial arithmetic and calculus, five classical root-finding algorithms run against a fixed menu of test functions, and integer number theory -- primality, GCD, modular exponentia

Three small, self-contained modules that replace pieces of `numpy.polynomial`, `scipy.optimize`, and `sympy.ntheory` with pure Mojo. `polynomial` stores a polynomial as a flat `List[Float64]` of coefficients (`coeffs[i]` is the coefficient of `xⁱ`, constant term first) and does evaluation, arithmetic, calculus, and composition on that representation. `root_finding` implements five classical scalar root-finders -- bisection, Newton-Raphson, secant, regula falsi, and a Newton-flavored fixed-point variant -- run not against an arbitrary user-supplied function but against a closed menu of four hardcoded test functions selected by integer `func_id`, plus a bracket-search helper and a convergence-order estimator for comparing them. `number_theory` is integer-only: primality testing, GCD/LCM, modular exponentiation and inverse, Euler's totient, prime factorization, a sieve, Fibonacci, and integer perfect-square testing. None of the three share code with each other; they are grouped on this page because they are the engine's classical numerical-methods primitives, not because one calls another.

```python
from algenta import MojoRuntime

runtime = MojoRuntime(mode="local")

# roots of x^2 - 5x + 6 = (x-2)(x-3)
runtime.execute("polynomial", "quadratic_roots", [1.0, -5.0, 6.0])["result"]
# [2.0, 3.0]

# Newton-Raphson on func_id 0 (x^2 - 2), starting from x0 = 1.0
runtime.execute("root_finding", "newton_raphson", [1.0, 0])["result"]
# 1.4142135623746899 -- sqrt(2) to ~11 digits; see below for why not more

runtime.execute("number_theory", "is_prime", [561])["result"]
# False -- 561 = 3 x 11 x 17, a Carmichael number that fools a naive Fermat test
```

## The modules

| Module          | What it does                                                                                                                                                                                                                              |
| --------------- | ----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `polynomial`    | Dense-coefficient polynomial arithmetic: Horner evaluation, add, multiply, derivative, integral, degree, scaling, composition, building a polynomial from its roots, and the closed-form quadratic formula.                               |
| `root_finding`  | Bisection, Newton-Raphson, secant, regula falsi, and a Newton-form fixed-point iteration, all run against four hardcoded test functions selected by `func_id`; plus a bracket-search helper and an empirical convergence-order estimator. |
| `number_theory` | Primality (deterministic Miller-Rabin), GCD/LCM, modular exponentiation and inverse, extended GCD, Euler's totient, prime factorization, a sieve of Eratosthenes, coprimality, Fibonacci, and integer perfect-square testing.             |

## `polynomial`: dense-coefficient arithmetic and calculus

| Function                       | What it does                                                                                                                                                                                                                    |
| ------------------------------ | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `poly_eval(coeffs, x)`         | Evaluates `p(x)` by Horner's method, working from the highest-degree term down. `poly_eval([], x)` returns `0.0`.                                                                                                               |
| `poly_add(a, b)`               | Coefficient-wise sum; the shorter input is treated as zero-padded rather than requiring equal lengths.                                                                                                                          |
| `poly_multiply(a, b)`          | Convolution: `result[i+j] += a[i]*b[j]`. Output has `len(a)+len(b)-1` terms; either empty input yields an empty result.                                                                                                         |
| `poly_derivative(coeffs)`      | `d/dx`: term `i` of the input (`i >= 1`) becomes term `i-1` of the output, scaled by `i`. A constant or empty input returns `[0.0]`.                                                                                            |
| `poly_integral(coeffs, C=0.0)` | Indefinite integral: prepends the constant `C`, then divides term `i` of the input by `i+1`.                                                                                                                                    |
| `poly_degree(coeffs)`          | Index of the highest coefficient whose magnitude exceeds `1e-15`. Returns `-1` only for a genuinely empty list -- see the edge case below.                                                                                      |
| `quadratic_roots(a, b, c)`     | Roots of `a*x^2 + b*x + c` by the closed-form formula. Degrades to the single linear root when `a` is \~0, and returns an empty list both for a fully degenerate `a=b=0` input and for a negative discriminant (complex roots). |
| `poly_scale(coeffs, s)`        | Multiplies every coefficient by scalar `s`.                                                                                                                                                                                     |
| `poly_compose(f, g)`           | Substitutes `g(x)` into `f`, producing `f(g(x))`, via repeated Horner-style multiplication by `g`.                                                                                                                              |
| `poly_from_roots(roots)`       | Builds `∏(x - rᵢ)` from a list of real roots. An empty root list returns the constant polynomial `[1.0]` -- the empty product.                                                                                                  |

**Worked arithmetic and calculus, hand-verified.** Evaluating `p(x) = 1 - 2x + 3x^2` at `x = 2` by Horner's method gives `1 - 4 + 12 = 9.0`, matching `poly_eval([1.0, -2.0, 3.0], 2.0)`. Its derivative is `-2 + 6x`, i.e. `poly_derivative([1.0, -2.0, 3.0]) == [-2.0, 6.0]`, and integrating `[1.0, -2.0, 3.0]` with `C=0.0` gives `[0.0, 1.0, -1.0, 1.0]` -- the polynomial `x - x^2 + x^3`. Differentiating that result recovers the original `[1.0, -2.0, 3.0]` exactly, since `poly_derivative` multiplies term `i` of `[C, 1, -1, 1]` by `i` for `i = 1, 2, 3`, giving `[1*1, 2*-1, 3*1] = [1.0, -2.0, 3.0]` -- integral and derivative are exact inverses of each other up to the constant term, on this representation.

`poly_from_roots([2.0, 3.0])` builds `(x-2)(x-3) = x^2 - 5x + 6`, returning `[6.0, -5.0, 1.0]`, and `poly_compose` substitutes one polynomial into another: composing `f = 1 + x^2` (`[1.0, 0.0, 1.0]`) with `g = x + 1` (`[1.0, 1.0]`) gives `f(g(x)) = 1 + (x+1)^2 = x^2 + 2x + 2`, i.e. `poly_compose([1.0, 0.0, 1.0], [1.0, 1.0]) == [2.0, 2.0, 1.0]`.

**`quadratic_roots` cannot represent a complex pair -- it returns nothing instead.** For `x^2 + 1 = 0` (`a=1, b=0, c=1`), the discriminant `b^2 - 4ac = -4` is negative, and the function returns an empty list rather than a complex pair, because the return type is `List[Float64]`. A caller cannot distinguish "no real roots" from certain other empty-list cases (see below) without also knowing which branch produced it. The degenerate `a ≈ 0` branch is a real fallback, not an error: `quadratic_roots(0.0, 2.0, -4.0)` returns `[2.0]`, the single linear root `-c/b`, and `quadratic_roots(0.0, 0.0, 5.0)` -- both leading coefficients zero, an inconsistent equation -- returns an empty list, the same empty result a genuine complex-root case produces.

**`poly_degree` of an all-zero, non-empty list is `0`, not `-1`.** The loop that walks down from the top looking for a non-negligible coefficient is guarded by `deg > 0`, so it stops at index `0` even if `coeffs[0]` is also below `1e-15` in magnitude: `poly_degree([0.0, 0.0, 0.0])` returns `0`. Only a literally empty `List[Float64]` (`poly_degree([])`) produces the `-1` sentinel that means "no coefficients at all." The identically-zero polynomial and a nonzero constant polynomial both report degree `0`.

## `root_finding`: five algorithms sharing one closed menu of test functions

| Function                                                    | What it does                                                                                                                                                                                                       |
| ----------------------------------------------------------- | ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ |
| `eval_func(x, func_id)`                                     | Evaluates test function `func_id` at `x`. Unrecognized `func_id` values silently return `0.0` rather than raising.                                                                                                 |
| `eval_deriv(x, func_id)`                                    | The matching hand-derived derivative for the same four functions; same silent-`0.0` fallback.                                                                                                                      |
| `bisection(a, b, func_id, tol=1e-10, max_iter=100)`         | Standard bracketing bisection. If `f(a)` and `f(b)` don't have opposite signs, returns the midpoint immediately rather than searching.                                                                             |
| `newton_raphson(x0, func_id, tol=1e-10, max_iter=100)`      | `x_{n+1} = x_n - f(x_n)/f'(x_n)`, stopping when `\|f(x)\| < tol` or the derivative underflows below `1e-15`.                                                                                                       |
| `secant_method(x0, x1, func_id, tol=1e-10, max_iter=100)`   | Newton's method with a finite-difference secant in place of the analytic derivative.                                                                                                                               |
| `fixed_point(x0, func_id, tol=1e-10, max_iter=100)`         | See the callout below -- this is Newton's method again, not generic fixed-point iteration.                                                                                                                         |
| `regula_falsi(a, b, func_id, tol=1e-10, max_iter=100)`      | False position: like bisection, but the next probe point is a secant intersection instead of the midpoint. Refuses a degenerate or non-finite bracket by returning the plain midpoint.                             |
| `bracket_root(func_id, start=0.0, step=0.5, max_steps=100)` | Walks outward from `start` in both directions looking for a sign change. See the caveats below.                                                                                                                    |
| `convergence_order(errors)`                                 | Estimates the convergence order from the last usable triple of consecutive errors via `log(e_{n+1}/e_n) / log(e_n/e_{n-1})`. Returns `0.0` if fewer than 3 usable (positive, non-degenerate) errors are available. |

All five solvers and `eval_func`/`eval_deriv` take an `Int func_id` rather than a callback, because the four test functions are hardcoded in `eval_func` and `eval_deriv` themselves:

| `func_id` | Function      | Derivative    | Real root(s)                      |
| --------- | ------------- | ------------- | --------------------------------- |
| `0`       | `x^2 - 2`     | `2x`          | `√2 ≈ 1.41421356`                 |
| `1`       | `x^3 - x - 2` | `3x^2 - 1`    | `≈ 1.52137970`                    |
| `2`       | `cos(x) - x`  | `-sin(x) - 1` | `≈ 0.73908513`                    |
| `3`       | `exp(x) - 3x` | `exp(x) - 3`  | `≈ 0.61906129` and `≈ 1.51213455` |

This is a real, narrower scope than a general-purpose numerical solver library: there is no way to hand these functions an arbitrary user-defined `f`; every solver here is exercised against this fixed menu of four textbook test cases, which is what "test functions identified by `func_id`" in the module's own docstring means literally.

**`fixed_point` is not generic fixed-point iteration -- it is Newton's method under another name.** Classical fixed-point iteration repeatedly applies an arbitrary map `g` and converges when `x_{n+1} = g(x_n) ≈ x_n`, for any `g` the caller chooses. Reading the body of `fixed_point`, its update rule is hardcoded to `x - f(x)/f'(x)` -- calling `eval_deriv` exactly as `newton_raphson` does -- with the only difference being an added `|x_{n+1} - x_n| < tol` stopping check alongside the `|f(x)| < tol` one. For every `func_id` in the table above, `fixed_point(x0, func_id)` and `newton_raphson(x0, func_id)` converge to the same value by the same iteration; there is no generic `g(x) = x` map here to substitute a different one in.

**Newton, bisection, secant, and regula falsi all land on the same root, to comparable precision.** Solving `func_id = 0` (root `√2`) from reasonable starting brackets/points: `newton_raphson(1.0, 0)` returns `1.4142135623746899`, `bisection(0.0, 2.0, 0)` returns `1.4142135623842478`, `secant_method(1.0, 2.0, 0)` returns `1.4142135623730954`, and `regula_falsi(0.0, 2.0, 0)` returns `1.4142135623637997` -- all four agree with `√2 = 1.4142135623730951` to at least 9 significant digits, and all four were hand-traced against the real formulas in `mojo/engine/root_finding.mojo` and independently reproduced in Python. Note that the default `tol=1e-10` bounds `|f(x)|`, not `|x - root|`; near `√2` the derivative is `2x ≈ 2.83`, so an `f(x)` just under `1e-10` corresponds to an `x` error on the order of `1e-10 / 2.83`, which is why Newton's actual returned value differs from the true `√2` by only `1.6e-12` rather than landing exactly on `tol`.

**`convergence_order`, measured on Newton's own iterates, reports ≈2.0 -- textbook quadratic convergence.** Tracing `newton_raphson`'s five iterates starting from `x0 = 1.0` on `func_id = 0` and computing the true error `|x_i - √2|` at each step gives the sequence `[0.41421356, 0.08578644, 0.00245310, 2.1239e-06, 1.5947e-12]` -- each error is roughly the square of the previous one, the signature of quadratic convergence. Feeding that sequence into `convergence_order` returns `1.9997666710861988`, matching the theoretical order-2 convergence of Newton's method on a simple root to three decimal places; this number was computed by re-running the exact recurrence in `convergence_order`'s docstring formula in Python, not invented.

**`bracket_root`'s fallback is not validated, and a coarse step size can straddle a pair of roots without detecting either.** When the search loop exhausts `max_steps` without finding a sign change, the function falls back to returning `[start, start + step]` unconditionally -- it does not check `eval_func(start, func_id) * eval_func(start + step, func_id) < 0.0` before returning, so a caller must verify that product itself before trusting the result as an actual bracket. This is not a hypothetical: `func_id = 3` (`exp(x) - 3x`) has two real roots, `≈0.619` and `≈1.512`, both strictly between `x = 0.6` and `x = 1.6`. Calling `bracket_root(3, 0.6, 1.0, 60)` samples only `f(0.6) ≈ 0.0221` and `f(1.6) ≈ 0.153` -- both positive, because the step of `1.0` jumps clean over both roots without sampling the negative dip between them -- so no sign change is ever seen, and the function falls through every iteration to the unvalidated `[0.6, 1.6]` fallback, which is not a valid bracket for either root.

## `number_theory`: primality, modular arithmetic, and small integer sequences

| Function                       | What it does                                                                                                                                                                                                   |
| ------------------------------ | -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `is_prime(n)`                  | Deterministic Miller-Rabin using witnesses `{2, 3, 5, 7}`, valid for `n < 2^31` (the module docstring's claimed range fits inside the proven deterministic bound of `n < 3,215,031,751` for this witness set). |
| `gcd(a, b)`                    | Euclidean algorithm; negative inputs are converted to their absolute value first.                                                                                                                              |
| `lcm(a, b)`                    | `\|a*b\| / gcd(a, b)`; `0` if either input is `0`.                                                                                                                                                             |
| `mod_pow(base, exp, modulus)`  | Modular exponentiation by repeated squaring. See the negative-exponent caveat below.                                                                                                                           |
| `euler_totient(n)`             | Euler's φ(n) by trial division up to `√n`. `0` for `n <= 0`, `1` for `n == 1`.                                                                                                                                 |
| `prime_factors(n)`             | Prime factorization with repetition, smallest factor first, via trial division. Empty for `n <= 1`.                                                                                                            |
| `sieve_of_eratosthenes(limit)` | All primes up to and including `limit`. Empty for `limit < 2`.                                                                                                                                                 |
| `is_coprime(a, b)`             | `gcd(a, b) == 1`.                                                                                                                                                                                              |
| `extended_gcd(a, b)`           | Returns `[gcd, x, y]` satisfying Bezout's identity `a*x + b*y = gcd(a, b)`. Unlike `gcd`, does not normalize the sign of its inputs first.                                                                     |
| `mod_inverse(a, m)`            | Modular multiplicative inverse of `a` mod `m`, built on `extended_gcd(a, m)`; returns `-1` if `gcd(a, m) != 1`.                                                                                                |
| `fibonacci(n)`                 | Nth Fibonacci number, iterative, `fibonacci(0) = 0`, `fibonacci(1) = 1`.                                                                                                                                       |
| `is_perfect_square(n)`         | Integer square root via Newton's method, then checks `x*x == n`. `False` for negative `n`.                                                                                                                     |

**`is_prime` correctly rejects a Carmichael number, which a naive Fermat test would not.** `561 = 3 * 11 * 17` is the smallest Carmichael number -- it passes a plain Fermat primality test for every base coprime to it, which is exactly the failure mode Miller-Rabin's extra squaring step exists to catch. Tracing the real algorithm: `n - 1 = 560 = 2^4 * 35`, so `d = 35`, `r = 4`. The first witness, `a = 2`, gives `x = 2^35 mod 561 = 263`, which is neither `1` nor `560`; squaring `r - 1 = 3` more times never lands on `560` either, so the witness proves `561` composite and `is_prime(561)` correctly returns `False` -- confirmed independently in Python against the same witness set and against brute-force trial division.

**`mod_pow` silently returns `1` for a negative exponent -- it does not compute an inverse.** The squaring loop is `while e > 0`, so a negative `exp` never enters it and `result` is returned unchanged at its initial value of `1`, identical to what `exp = 0` produces: `mod_pow(3, -1, 7)` returns `1`, not the true inverse `5` (since `3 * 5 = 15 ≡ 1 mod 7`). A caller wanting a modular inverse needs the dedicated `mod_inverse` function, not a negative exponent passed to `mod_pow`.

**`gcd`/`lcm` normalize sign; `extended_gcd`/`mod_inverse` do not, except at the very end.** `gcd` and `lcm` explicitly take the absolute value of both inputs before running the Euclidean algorithm, so `gcd(-48, 18) == gcd(48, 18) == 6`. `extended_gcd` runs the extended Euclidean algorithm directly on the signed inputs it is given with no such normalization step, so its Bezout coefficients (and occasionally the sign of the reported gcd itself) depend on the actual signs of `a` and `b`, not their magnitudes. `mod_inverse` builds on `extended_gcd(a, m)` and only normalizes at the very end, folding a negative Bezout coefficient into `[0, m)` with a single `+= m` if needed -- e.g. `mod_inverse(3, 11)` returns `4` (since `3*4 = 12 ≡ 1 mod 11`), and `mod_inverse(4, 8)` returns `-1` (since `gcd(4, 8) = 4 != 1`, no inverse exists), both confirmed by re-running the same recurrence independently.

**Worked chain: `euler_totient`, `prime_factors`, and `sieve_of_eratosthenes` on the same numbers.** `36 = 2^2 * 3^2`, so `prime_factors(36)` gives `[2, 2, 3, 3]` and `euler_totient(36) = 36 * (1 - 1/2) * (1 - 1/3) = 12` -- both independently re-derived and matching the trial-division formula in the source. `sieve_of_eratosthenes(30)` returns `[2, 3, 5, 7, 11, 13, 17, 19, 23, 29]`, the ten primes up to 30. `fibonacci(10) = 55` and `is_perfect_square(144)` is `True` (`12^2 = 144`) while `is_perfect_square(145)` is `False`, both traced through the integer-Newton loop by hand.

## Verification

All three `mojo/tests/libs/test_*.mojo` files behind this page (`test_polynomial.mojo`, `test_root_finding.mojo`, `test_number_theory.mojo`) are auto-generated coverage smoke tests: each calls every function once with a literal fixture (mostly `0.5`-valued floats or the integer `5`) and checks only that a numeric result is finite or a list result has non-negative length -- none of the three pin any specific expected output. Every numeric claim on this page -- the Horner evaluation, the derivative/integral round trip, the quadratic and composed-polynomial coefficients, the four solvers' agreement on `√2`, the Newton convergence-order measurement, the Carmichael-number primality trace, the negative-exponent `mod_pow` behavior, and the totient/factorization/sieve chain -- was instead hand-traced against the real formulas in `mojo/engine/polynomial.mojo`, `root_finding.mojo`, and `number_theory.mojo`, and independently re-run in Python to confirm the value before being presented here.

## Relationship to other pages

Grepping `gitbook/` for `polynomial`, `root_finding`, and `number_theory` as module names turns up nothing besides this page and the auto-generated [Runtime library catalog](/engine-runtime/runtime-library-catalog.md), which lists them with 10, 9, and 12 functions respectively -- matching the counts above exactly. Two near-miss names are worth ruling out explicitly: [`ml-classical-algorithms.md`](/machine-learning/ml-classical-algorithms.md)'s `ml.preprocessing.polynomial_features` shares the word "polynomial" but is a feature-engineering transform that expands a numeric vector into monomial cross-terms for a downstream model -- it does not touch this page's `List[Float64]` coefficient representation or arithmetic at all. And [`optimization.md`](/optimization-and-simulation/optimization.md)'s six modules (`optimization.model`, `.linear`, `.integer`, `.assignment`, `.flow`, `.explain`) solve linear and integer programs, assignment, and network flow -- none of them do scalar root-finding or expose anything resembling `newton_raphson` or `bisection`; the two families do not share code or call each other.

**Next:** see the [Runtime library catalog](/engine-runtime/runtime-library-catalog.md) for every module and function, or [Optimization & decision programs](/optimization-and-simulation/optimization.md) for the engine's other family of classical numerical methods.


---

# Agent Instructions
This documentation is published with GitBook. GitBook is the documentation platform designed so that both humans and AI agents can read, navigate, and reason over technical content effectively. Learn more at gitbook.com.

## Querying This Documentation
If you need additional information that is not directly available in this page, you can query the documentation by asking a question.

Perform an HTTP GET request on the following URL with the `ask` and `goal` query parameters:

```
GET https://docs.algenta.ai/numerical-and-mathematical-computing/polynomial-root-finding-and-number-theory-kernels.md?ask=<question>&goal=<user_goal>
```

`ask` is the immediate question: it should be specific, self-contained, and written in natural language.
`goal` is what the user is ultimately trying to achieve, the reason they need the answer. Sharing it helps GitBook give you a better, more relevant answer. A goal is most helpful when it describes the outcome the user wants rather than restating the question. For example, with `ask=how do I create an API token`, a goal like `automate deployments from our CI pipeline` lets GitBook tailor the answer to that use case.

The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.