[libre-riscv-dev] algebraic numbers library

[Jacob Lifshay]

On Tue, Sep 17, 2019, 06:58 Hendrik Boom <hendrik at topoi.pooq.com> wrote:

> On Mon, Sep 16, 2019 at 04:11:12AM -0700, Jacob Lifshay wrote:
> > I've been writing an algebraic numbers library for a while with the
> > intention that it could be used to create a reference by which we
> > could test our IEEE 754 operations, since the numbers it represents
> > are exact instead of approximated.
> >
> > The main API is using RealAlgebraicNumber in
> >
> https://salsa.debian.org/Kazan-team/algebraics/blob/master/src/algebraic_numbers.rs
> >
> > the RealAlgebraicNumber type supports add, sub, mul, div (but not by
> > zero), comparisons, raising to rational powers, floor, ceil, trunc,
> > fract, conversion to rational numbers and integers (returning None if
> > the number is not rational or integer), and conversions from rational
> > numbers and integers.
> >
> > It is based on polynomials over bigints, so the only restrictions are
> > execution time and memory usage.
>
> Is there any documentation about the algorithms you used and the
> mathematical theory behind them?
>

not much, but they're well-known algorithms.

I based it on storing the content-free minimal polynomial for the
particular algebraic number as well as an interval (DyadicFractionInterval
-- basically intervals of the form a/2^n <= x <= b/2^n where x is a real
number, a and b are integers with a <= b, and n is a unsigned integer) in
which the number is located to distinguish which root it is.

having it based on the content-free minimal polynomial means that the
polynomial is irreducible (like a prime number, but for polynomials), so it
doesn't ever have any repeated roots. the other benefit of content-free
minimal polynomials is that for all algebraic numbers, two equal algebraic
numbers always have the same minimal polynomial (but not the converse,
polynomials can have multiple roots).

the algorithms are based on resultants, polynomial factorization, sturm
sequences, and interval arithmetic.

for polynomial factorization, I based it on dividing out the content to get
the content-free part (generally called primitive part) then using yun's
algorithm for square-free factorization then using a variation of
Zassenhaus's algorithm (
https://en.wikipedia.org/wiki/Factorization_of_polynomials#Factoring_univariate_polynomials_over_the_integers
-- not to be confused with
https://en.wikipedia.org/wiki/Zassenhaus_algorithm). I'm planning on
eventually switching the exponential-time last step (finding which modular
factors combine to form the integer factors by trying all subsets) to using
the polynomial-time LLL algorithm (which I implemented), but have not done
that yet.

for square-free modular polynomial factorization, I used the combination of
https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Distinct-degree_factorization
and
https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Cantor%E2%80%93Zassenhaus_algorithm
(note that Fq is the quotient ring Fp[x]/(e) where e is an element of Fp[x]
-- basically, polynomial remainders after dividing by e. I think q is e,
but not sure. Fp[x] is the ring of polynomials in the variable x over a
finite field Fp where p is a prime number. Fp is the ring of integers mod
p. The algorithms require p to be an odd prime. It took me a while to
figure it all out.)

I typed all the following out, then realised i'm just paraphrasing the
source -- keeping anyway.

definitions:
content-free sturm sequence -- a sturm sequence where each polynomial is
scaled by a potentially-different positive fraction such that each
polynomial is a content-free polynomial with a denominator of 1 -- has same
sign changing behavior as a standard sturm sequence, but the numbers are
much smaller and integers rather than fractions.

the operations are:
shrink RealAlgebraicNumber's interval:
-- calculates a smaller interval containing the algebraic number
1. if the interval's lower and upper bounds are equal, return
2. if the interval's lower and upper bounds are close enough, increase the
interval's denominator and numerators proportionally
3. let middle be the mid-point of the interval
4. let cfss be the content-free sturm sequence for the minimal polynomial.
5. compute the number of sign changes for the lower and upper bounds of the
interval and for middle using cfss
6. if any of the lower bound, upper bound, or middle are a root, set
interval to that number and return
7. if the number of sign changes for the lower bound equals the number of
sign changes for middle, then:
    a. the root is not in the interval between the lower bound and middle
    b. set interval's lower bound to middle
8. else:
    a. assert that the number of sign changes for the upper bound equals
the number
        of sign changes for zero -- the root is not in the interval between
the upper
        bound and zero
    b. set interval's upper bound to middle
9. return

RootSelector is a trait with the following operations:
    get_interval(&self): returns an interval containing the root to be
selected
    shrink_interval(&mut self): shrinks the intervals, causing get_interval
to eventually
        return a more accurate interval
    select_root(self, p(x)): selects the correct root of the polynomial
p(x), takes
        ownership of self

RootSelector::select_root(p(x)):
1. factor p(x) into a integer part and a list of irreducible polynomials
with the polynomial's power -- so the polynomial 2 * (x)^4 * (2*x - 1)^3
has 2 as the integer part and [(x, 4), (2*x - 1, 3)] as the list
2. let factors2 be the list of just the irreducible polynomials, discarding
the powers and integer part
3. let factors be the list of content-free sturm sequences for each
polynomial in factors2
4. let factors_temp be an empty list
3. loop:
    a. let interval be self.get_interval()
    b. swap factors and factors_temp
    c. clear factors
    d. let roots_left = 0
    e. for each element `factor` of factors_temp:
        1. let lower_bound_sign_changes be the number of sign changes for
interval's
            lower bound using factor
        2. if interval's lower bound is a root, return lower bound
converted to
            RealAlgebraicNumber
        3. let upper_bound_sign_changes be the number of sign changes for
interval's
            upper bound using factor
        4. if interval's upper bound is a root, return upper bound
converted to
            RealAlgebraicNumber
        5. if lower_bound_sign_changes != upper_bound_sign_changes:
            a. roots_left += abs(lower_bound_sign_changes -
upper_bound_sign_changes)
            b. append factor to factors
    f. if roots_left <= 1:
        1. assert factors is not empty (implicit in .remove(0) in source
code)
        2. let factor be the corresponding minimal polynomial for the only
element of
            factors
        3. return the RealAlgebraicNumber for factor and interval
    g. call self.shrink_interval()


convert from rational number (or integer):
1. calculate an interval containing the rational number
2. the minimal polynomial is d*x - n where d is positive and n/d is the
rational number
3. return the corresponding RealAlgebraicNumber

negate:
based on calculating the polynomial p(-x) where the input polynomial is
p(x) and computing the negative of the interval

compare with zero:
1. if the number is rational (minimal polynomial is degree 1), compare
rational with zero directly (this is the only case that needs to handle
equality since zero is a rational number)
2. if zero is not inside the interval, return less or greater based on
comparing the interval with zero
3. let cfss be the content-free sturm sequence for the minimal polynomial
4. compute the number of sign changes for the lower and upper bounds of the
interval and for zero using cfss
5. if the number of sign changes for the lower bound equals the number of
sign changes for zero, return greater -- the root is not in the interval
between the lower bound and zero
6. assert that the number of sign changes for the upper bound equals the
number of sign changes for zero -- the root is not in the interval between
the upper bound and zero
7. return less

add lhs + rhs:
1. if lhs and rhs are rational, add them as rational numbers and return the
result converted to RealAlgebraicNumber
2. let lp(x) be the lhs's minimal polynomial and rp(x) be rhs's.
3. let p(x) be the resultant of lp(y) and rp(x-y) with respect to y (so, in
Maxima terms, resultant(lp(y), rp(x-y), y))
4. AddRootSelector is a RootSelector where:
    RootSelector::get_interval returns lhs.interval + rhs.interval
    RootSelector::shrink_interval runs shrink on lhs and runs shrink on rhs
5. return AddRootSelector.select_root(p(x))

subtract lhs - rhs:
1. return lhs + -rhs

compare lhs with rhs:
1. return (lhs - rhs).compare_with_zero()

Jacob Lifshay