Rationalising Denominators 2: Products

Posts in this series:

Powers

Products

Sums

Ratios

Roots of unity

Rationalising denominators

Introduction

Last time we used Python’s int and Fraction types to implement a new type called Power, to represent nth roots.

We implemented a few normalisation steps, to avoid redundant values like 1² (which is just 1) and 0³ (which is just 0).

We were still left with some redundancy though, like 4¹ and 2² both being the number 4; or, conversely, 2¹ and √4 both being the number 2.

This time we’ll use Power to implement another type, called Product, which completely avoids these redundancies.

Products

A product is a collection of numbers multiplied together. We can define a Product type in Python as a container of Power values:

def Product(*powers: Power) -> List[Power]:
    return list(powers)

This can represent all the same values as Power (by using a list with a single value), but it lets us represent some values in a different way. For example, Power can only represent some numbers by taking a large root of a large base, like ⁶√72; whilst Product can break that down into ∛3×√2.

Here are some useful constants for this type:

product_zero = Product(power_zero)
product_one = Product(power_one)

Representing negatives

Another advantage of Product is that it can represent negative numbers, if we’re somehow able to represent -1. To do this, we’ll use a neat insight: that square roots have two solutions, e.g. 4 = 2² and 4 = (-2)². When we interpret a Power representing a square root, we’ve been assuming that it refers to the positive root; and that’s a good convention to stick with. However, we’ll make a special exception for the value √1: in that case, the positive root is 1, which is perfectly represented already, without the need to take a square root. That redundancy allows us to instead use √1 to represent -1, which has no other encoding in our representation.

Whilst this interpretation only expands the range of Power by a single extra value, the inclusion of -1 in a Product allows us to represent the negative of all our values!

This is why, in the previous post, we took the seemingly-strange decision to use the modulo operation when normalising Power values with base 1. That was overkill for normalising whole exponents of 1, but is exactly what we need to handle this fractional representation of -1. For example, the improper fraction ³⁄₂, as an exponent of 1, would mean -1³, or -1 × -1 × -1, which normalises to -1; with exponent ½.

We’ll define a constant for this special value:

power_neg = Power(1, half)

Since -1 is a rational number, we need to add it as a special-case to the power_is_rational function we defined last time (since its exp.denominator == 1 check isn’t compatible with this -1 = √1 trick):

def power_is_rational(p: Power) -> bool:
    base, exp = p
    return p == power_neg or exp.denominator == 1

We also need a special-case in eval_power_int, to interpret power_neg as the integer -1:

def eval_power_int(p: Power) -> int:
    assert power_is_rational(p), f"Can't eval {p} as int"
    if p == power_neg:
        return -1
    base, exp = p
    return base**exp.numerator

Improving our representation

Since we’re dealing with commutative multiplication (for the moment), changing the order of arguments given to Product should not affect its result; yet the order does matter if we represent them with a list. Instead, we could use a bag (AKA a multiset); which can be implemented as a dictionary which counts the number of times a particular Power appears:

def Product(*powers: Power) -> dict[Power, int]:
    counts = {}
    for power in powers:
        counts[power] = counts.get(power, 0) + 1
    return counts

However, we can exploit the semantics of Power to get a simpler representation: multiplying two Power values with the same base results in that base raised to the sum of their exponents. This means our dictionary can directly associate bases with (sums of) exponents:

def Product(*powers: Power) -> dict[Base, Exponent]:
    result = {}
    for base, exp in powers:
        result[base] = result.get(base, zero) + exp
    return result

To turn such a Product back into a list of Power values, we can call its .items() method.

Note that this will automatically handle double-negatives, since √1 has an exponent of ½, so a pair of them will have their exponents added to give a positive factor of 1¹, and this will keep cycling for each √1 in the Product.

Multiplication

Now we’ve chosen the final structure of Product, we can write a multiplication function for it. Since Product already represents a multiplication of its contents, we can multiply Product values together by collecting up all their elements into one big Product:

def multiply_products(*ps: Product) -> Product:
    return Product(*sum([list(p.items()) for p in ps], []))

Normalisation

There is plenty of redundancy in this simple definition, so we’ll build up a normalisation algorithm which consolidates such values to a single representation.

Absorbing Zero

All products which contain 0 are equivalent (since 0 is the absorbing value for multiplication). The simplest of these is the product which only contains 0, so we’ll normalise the others to that:

        if base == 0:
            # If any element is zero, the whole Product is zero
            return dict([power_zero])

Discarding Ones

The number 1 is the identity value for multiplication, so its presence in a product does not alter the result. Hence the normalisation of a Product should remove any such Power. We can check for this in our loop, by ignoring all exponents that are 0:

        if exp != zero:
            result[base] = result.get(base, zero) + exp

Note that, if every Power given to Product is 1, then the result will be the empty product; and this is the normal form for the Product representing the number 1.

Fixed Points

Performing multiple checks in sequence may cause us to miss redundancies, where earlier checks could normalise the results of later ones. There are many ways we could adjust our code to fix this, but one particularly simple approach is to call ourselves recursively, until our result matches our input (i.e. until we reach a fixed point):

    return result if len(powers) == len(result) and result == dict(powers) \
        else Product(*result.items())

Factorising

Now we confront the main source of redundancy in both Product and Power: that different sets of powers can multiply to the same result; and indeed that raising different base values to appropriate exponents can also produce the same result.

We will avoid this redundancy by only allowing elements of a Product to have base values which are prime, or the special values 0 and 1.

To implement this normalisation, we need some helper functions to factorise base values. There are many algorithms we could use, and the details aren’t too important for our purposes; you can expand the following box to see the implementation I’ve used.

from functools import lru_cache

def factorise() -> Callable[[Base], Product]:
    # Precompute primes up to 1000 using Sieve of Eratosthenes
    sieve = [True] * 1001
    sieve[0] = sieve[1] = False
    for i in range(2, 32):  # sqrt(1000) ≈ 31.6
        if sieve[i]:
            for j in range(i * i, 1001, i):
                sieve[j] = False
    SMALL_PRIMES = [i for i in range(1001) if sieve[i]]
    del(sieve)

    # Slow, but works for arbitrary n
    def wheel_factorise(n: int) -> Product:
        factors = {}

        # Handle 2, 3 and 5 separately for efficiency
        for p in (2, 3, 5):
            if n % p == 0:
                factors[p] = Fraction(0, 1)
                while n % p == 0:
                    factors[p] += 1
                    n //= p

        # Wheel increments for numbers coprime to 2, 3 and 5
        increments = [4, 2, 4, 2, 4, 6, 2, 6]
        i = 0
        f = 7

        while f * f <= n:
            if n % f == 0:
                factors[f] = Fraction(0, 1)
                while n % f == 0:
                    factors[f] += 1
                    n //= f
            f += increments[i]
            i = (i + 1) % 8

        if n > 1:
            factors[n] = Fraction(1, 1)

        return factors

    @lru_cache(maxsize=1024*1024)
    def factorise(n: int) -> Product:
        if n == 0:
            return dict([power_zero])
        if n == 1:
            return dict([power_one])
        if n < 0:
            # Negatives have negative one as a factor
            factors = factorise(-n)
            factors[1] = factors.get(1, Fraction(0, 1)) + Fraction(1, 2)
            return factors
        if n > 1000:
            # Fallback method for large numbers
            return wheel_factorise(n)
        else:
            # Small numbers can use precomputed primes
            factors = {}
            for p in SMALL_PRIMES:
                if p * p > n:
                    if n > 1:
                        factors[n] = one
                    break

                if n % p == 0:
                    factors[p] = zero
                    while n % p == 0:
                        factors[p] += 1
                        n //= p
            return factors
    return factorise
factorise = factorise()

We use that factorise function to implement the following split_power function, which adjusts factors of a base to account for the given exponent:

def split_power(power: Power) -> List[Power]:
    base, exp = power
    return [power] if base == 1 else [
        Power(b, exp * e) for b, e in factorise(base).items()
    ]

Then, rather than Product accumulating its given Power arguments as-is, it can instead accumulate the results of split_power called on those arguments:

def Product(*powers: Power) -> dict[Base, Exponent]:
    result = {}
    for power in powers:
        for base, exp in split_power(Power(*power)):
            if base == 0:
                # If any element is zero, the whole Product is zero
                return dict([power_zero])
            b, e = Power(base, result.pop(base, zero) + exp)
            if e != zero:
                result[b] = e
    return result if len(powers) == len(result) and result == dict(powers) \
        else Product(*result.items())

The redundancy in our earlier examples is now avoided: since 4 is not prime, any Power with base 4 will be replaced by an equivalent involving base 2:

• 4¹ will normalise to 2².
• √4 will normalise to 2¹.

Note that we don’t perform this factorising in Power, since it may require a Product of multiple elements. For example, normalising 18² produces 2²×3⁴.

Also note that the use of b, e = Power(…) will normalise any powers of 1, which may have accumulated due to repeated negation, into having an exponent that’s strictly less than one. In particular, factors of 1¹ (caused by double-negatives) will get normalised to 1⁰, and hence get skipped by the if e != zero: check.

Extra Helpers

We’ll be needing the following helper functions, so might as well define them here alongside Product:

def eval_product_int(p: Product) -> int:
    return reduce(
        lambda x, y: x * eval_power_int(y),
        p.items(),
        1
    )

def product_is_rational(p: Product) -> bool:
    return all(power_is_rational(x) for x in p.items())

Conclusion

We’ve extended our representation of numbers from Power, which suffered from redundancy; to Product, which does not have any redundancies. Normalising this Product representation has exposed some structure that will be useful in future installments: namely the use of prime factors.

In the next post we’ll extend our representation again: to sums of products of powers!

The full code for this post is available in the fraction_products module of my conjugate repo.