Rationalising Denominators 6: Rationalising Denominators

Posted on by Chris Warburton

Posts in this series:

Powers

Products

Sums

Ratios

Roots of unity

Rationalising denominators

Introduction

In an earlier installment we introduced Ratio values, but they had a problem: the same numerical quantity could be represented by multiple different Ratio values. For example (1 + √1∛2 + ∛2²)/3 is numerically equal to 1/(1 + ∛2), which we can demonstrate by multiplying the second expression by (1 + √1∛2 + ∛2²)/(1 + √1∛2 + ∛2²) (which is equal to 1, since it’s the Ratio of a non-zero number to itself):

1/(1 + ∛2) × (1 + √1∛2 + ∛2²)/(1 + √1∛2 + ∛2²)
= (1 × (1 + √1∛2 + ∛2²) / ((1 + ∛2) × (1 + √1∛2 + ∛2²))
= (1 + √1∛2 + ∛2²) / (1 + √1∛2 + ∛2² + ∛2 + √1∛2² + 2)
= (1 + √1∛2 + ∛2²) / 3

The result is the first Ratio, proving they are the same. This ambiguity means our Ratio representation does not yet have a normal form. In this post we’ll identify the normal form for Ratio as that with a rational denominator; and implement an algorithm to “rationalise the denominator” of any Ratio of Real numbers.

This algorithm uses the roots of unity introduced in the previous post, to construct the algebraic conjugates of a number: these are combined into a “rationalising factor” (like 1 + √1∛2 + ∛2² in the example above), such that they all cancel-out to give a Real number result.

Our approach is guaranteed to terminate with a normal form if the given Ratio is a Real number (i.e. containing no roots of unity other than 1 and √1). However, given a non-Real Ratio then it may fail to terminate at all! This comes from a reliance on unique factorisation, which breaks down when the Real numbers are extended with certain Complex roots of unity. Whilst interesting, this isn’t a big deal for us; since I only set out to normalise (a sub-set of) Real numbers in the first place!

Radicals don’t factorise uniquely

To normalise a Ratio we can try dividing the numerator and denominator by any factors they have in common. For Natural numbers this is easy: each factorises into a unique, finite set of primes; and those which don’t share a common prime factor do not share any common factor.

We lose this uniqueness if we allow factors to be arbitrary rationals, e.g. a prime like 2 can be factorised further into rationals like 6 × ⅓, or many other combinations. Thankfully we don’t need to do that, since it doesn’t allow us to represent any more numbers than our top-level Ratio already can. However, once we introduce radicals then fractional factors become unavoidable.

For example, consider √2 / ⁶√2⁵. Its numerator and denominator do not have any common prime factors, but they do have radical factors in common, like ⁶√2; since √2 = ⁶√2 × ∛2 and ⁶√2⁵ = ⁶√2 × ∛2². Dividing numerator and denominator by this common factor gives us ∛2 / ∛2²

These have a common factor of ∛2, which we can divide through to get 1 / ∛2

Alternatively, we could have divided-through by 1/∛2 to get ∛2² / 2

If we allow roots of unity, then we could keep going forever, dividing by more and more radical factors, without ever reaching a “base case” that no longer factors. We need a more systematic approach, to avoid such endless rearranging. For this, we’ll need to learn about algebraic conjugates.

Algebraic conjugates

All of the (Real) values we’ve been dealing with are Algebraic numbers, and those can be characterised as the zero-values of a univariate polynomial. However, that isn’t a unique definition, since polynomials can have multiple zero-values (depending on their degree). For example, the polynomial x² - 7 has zero-values of √7 and √1√7. These are called the algebraic conjugates of each other.

Conjugates of Powers

For a simple radical like ⁿ√b, the conjugates are found by multiplying by the n-th roots of unity. For example, the conjugates of √7 are √7 × 1 and √7 × √1 (where 1 and √1 are the square roots of unity).

When n = 1 we have a rational value. There is only a single “first root of unity”, which is the number 1; so these values have no other conjugates, just themselves.

The other special-case is √1 itself: this is a square root, so n = 2; but nevertheless it is still rational. In this case its conjugates consist of itself √1 × 1 = √1, but also its negation √1 × √1 = 1. Thankfully, this “extra” conjugate won’t cause us any issues, since we’ll be multiplying these conjugates together; so having an extra 1 in that Product makes no difference!

def get_power_conjugates(p: Power) -> List[Product]:
    """For a power with denominator n, multiply by distinct powers of nth
    roots of unity. For example, n = 3 uses the cube root of unity 1^1/3,
    whose powers are 1^0/3, 1^1/3, 1^2/3, 1^3/3, 1^4/3, etc. Since they're
    cyclic, 1^k/3 = 1^(k+3)/3, so distinct powers occur for k in [0,1,2].
    In general, 1^(k+n)/n = 1^k/n, so distinct nth roots occur for k in
    range(n). Note when n = 1, range(n) = [0], giving only 1^0/n = 1^0 = 1"""
    base, exp = p
    n = exp.denominator
    return [Product(p, Power(1, Fraction(k, n))) for k in range(n)]

Conjugates of Products

The conjugates of a Product are found by taking all possible combinations of the conjugates of its constituent Powers:

def get_product_conjugates(p: Product) -> List[Product]:
    """Get all conjugates of a product using roots of unity."""
    if product_is_rational(p):
        return [p]

    # For each power b^e in the product, we need all eth roots of unity
    return reduce(
        lambda conjugates, power: [
            multiply_products(conj, root)
            for conj in conjugates
            for root in get_power_conjugates(power)
        ],
        p.items(),
        [product_one],
    )

Conjugates of Sums

The conjugates of a Sum are found by taking all possible combinations of the conjugates of its constituent Products. This can get complicated, but the principle is the same:

def get_sum_conjugates(s: Sum) -> List[Sum]:
    """Get all conjugates of a sum by getting conjugates of each product
    and taking all possible combinations."""
    products = _sum_to_products(s)
    if len(products) == 0:
        return [s]  # Base case

    # Recursively combine conjugates of head and tail
    first, *rest = products
    return [
        add_sums(Sum(f), r)
        for f in get_product_conjugates(first)
        for r in get_sum_conjugates(Sum(*rest))
    ]

Rationalising denominators (finally!)

Our aim is to manipulate a Ratio, without changing its overall value, until its denominator becomes a Rational number. That will be our normal form. The key to achieving this is a theorem from Galois theory (which I admit is a little beyond my current understanding!), which tells us that multiplying an Algebraic number by all of its algebraic conjugates will always give a Rational result.

We can see this with our running example 1 / (1 + ∛2). Applying get_sum_conjugates to its denominator 1 + ∛2 gives us a set of three algebraic conjugates:

If we multiply together all of the conjugates except the denominator itself, we get a “rationalising factor” which (according to Galois theory) multiplies with the denominator to give a Rational number.

In our example, we have (1 + ∛2 × ∛1) × (1 + ∛2 × ∛1²)

Expanding this we get 1 + ∛2 × ∛1² + ∛2 × ∛1 + ∛2² × ∛1 × ∛1²

That simplifies to 1 + ∛2 × (∛1 + ∛1²) + ∛2² × ∛1³

We know that ∛1³ = 1, and on the roots of unity page we saw that 1 + ∛1 + ∛1² = 0 (since roots of unity form regular polygons centred on 0), which in this case implies ∛1 + ∛1² = √1.

Substituting these into the expression for our “rationalising factor”, we get 1 + ∛2 × √1 + ∛2² × 1, which is 1 + √1∛2 + ∛2² that we saw in the introduction.

Notice how the roots of unity have cancelled out to give us a (non-zero) Real number. If we multiply our numerator and denominator by this factor, the overall value of our Ratio will stay the same; yet we will have a Rational denominator, like we were aiming for!

Continuing our example of 1 / (1 + ∛2), the numerator is just 1, so multiplying it by the above factor leaves that factor unchanged. As we saw in the introduction, multiplying our denominator by this factor simplifies to 3 (supporting our assumption that the product of all conjugates is rational!).

Hence the final normal form for this example is (1 + √1∛2 + ∛2²) / 3

def rationalise_denominator(e: Ratio) -> Ratio:
    """Rationalise the denominator of an ratio by multiplying
    numerator and denominator by all conjugates of denominator."""
    num, den = e

    if sum_is_rational(den):
        return e

    # Get all conjugates of denominator except itself
    conjugates = get_sum_conjugates(den)[1:]

    # Multiply together all conjugates except the denominator itself
    multiplier = multiply_sums(*conjugates)

    # Multiply both numerator and denominator
    new_num = multiply_sums(num, multiplier)
    new_den = multiply_sums(den, multiplier)

    if new_den == sum_zero:
        error_message = 'Made a zero denominator'
        raise Exception(repr(locals()))

    return (new_num, new_den)

We can incorporate this step into our Ratio function as follows (making sure to recurse, in case this rationalisation algorithm introduces common rational factors we need to cancel out):

    if not sum_is_rational(denominator):
        return Ratio(
            *rationalise_denominator((numerator, denominator)),
            negating=negating,
        )

Conclusion

This is the final post in this series. We’ve finally arrived at a normal form for a large sub-set of the Real numbers (radical expressions with no nesting) which seems like it’s as large as we can get (or pretty close; since I don’t think a normal form is possible for nested radicals, or more general Algebraic numbers).

In the process, we’ve seen a general implementation of the classic “rationalise the denominator” algorithm; going well beyond the quadratic version often taught in school, which can only deal with square-roots.

We’ve taken a detour into Complex numbers, and roots of unity, in order to reach the idea of algebraic conjugates; which have ultimately cancelled themselves out to give us a relatively simple normal form for Ratio values, which is what we set out to do!

At some point I aim to re-implement this in Racket, and incorporate it into my Ivory Tower project; which can use it for the coefficients of Complex numbers, Dual numbers, and more!

Here is the full implementation of the Ratio function, which incorporates all of the normalisation steps we’ve discussed:

def Ratio(numerator: Sum, denominator: Sum, negating=False) -> Tuple[Sum, Sum]:
    assert denominator != sum_zero

    # 0/n normalises to 0/1
    if numerator == sum_zero:
        return (sum_zero, sum_one)

    # Find common factors in the numerator's terms
    num_factors, remainder = common_factors_of_sum(
        numerator,
        allow_remainder=False
    )
    assert remainder == 0

    # Find common factors in the denominator's terms, and keep those which were
    # also in the numerator's
    common_product = Product(*[
        (base, min(exp, num_factors[base]))
        for base, exp in common_factors_of_sum(
                denominator,
                allow_remainder=False
        )[0].items()
        if base in num_factors
    ]).copy()

    # Rationalising the denominator scales poorly, so we remove common factors
    # first. However, we want to avoid re-introducing irrationalities once the
    # denominator is rational!
    denominator_is_rational = sum_is_rational(denominator)
    if denominator_is_rational:
        # We keep the denominator rational by sticking to rational factors, i.e.
        # those with exponent denominator 1.
        for base, exp in list(common_product.items()):
            # Round-down fractional exp values to the nearest natural
            common_product[base] -= exp % 1
            # If it rounded down to zero, we can remove it (since base^0 == 1)
            if common_product[base] == 0:
                del(common_product[base])

    # Remove common factor, if it's not trivial
    if common_product not in [product_one, product_neg]:
        # Divide-through by common factor
        return Ratio(
            remove_common_product(numerator, common_product),
            remove_common_product(denominator, common_product),
            negating=negating,
        )

    # If we're here, our denominator is as small as possible; so now rationalise
    # it if needed
    if not denominator_is_rational:
        return Ratio(
            *rationalise_denominator((numerator, denominator)),
            negating=negating,
        )

    # Finally, try negating the whole thing (unless we've already done so) and
    # seeing if that's simpler.
    given = (numerator, denominator)
    negated = given if negating else Ratio(
        multiply_sums(numerator, sum_neg),
        multiply_sums(denominator, sum_neg),
        negating=True
    )
    return Ratio(*negated) if r_complexity(negated) < r_complexity(given) \
        else given

The full code for all the posts in this series, along with optimisations and tests, can be found in my conjugate repo.