Bacovia: a symbolic math library

By
Named after the great symbolist poet, George Bacovia, I created this library to symbolically manipulate mathematical expressions in a simple and elegant way.


Before writing a symbolic math library, this was a somewhat mysterious subject to me, but in this post I would like to try to demystify and illustrate the beauty and satisfaction that comes from writing this simple, but powerful, symbolic math library.

The first thing in creating a new project, is the selection of the right programming language for the project. Just for fun, I decided to implement this library in the Sidef programming language, using recursion and multiple dispatch; two very powerful features that turned out to be just perfect for this task.

# Design

The library has a simple class hierarchy, with one base class, and several other classes inheriting from it.
class Bacovia              {}
class Power      < Bacovia {}
class Fraction   < Bacovia {}
class Difference < Bacovia {}
class Sum        < Bacovia {}
class Product    < Bacovia {}
class Log        < Bacovia {}
class Exp        < Bacovia {}
class Symbol     < Bacovia {}

# Classes

The classes provided by the Bacovia library, are described bellow:

Symbol(name, value=nil)

Represents a symbolic value. Optionally, it can have a numeric value.

Fraction(numerator, denominator)

Represents a symbolic fraction.

Difference(minuend, subtrahend)

Represents a symbolic difference.

Power(base, power)

Represents a symbolic exponentiation in a symbolic base.

Log(x)

Represents the natural logarithm of a symbolic value.

Exp(x)

Represents the natural exponentiation of a symbolic value.

Sum(a, b, c, ...)

Represents a summation of an arbitrary (finite) number of symbolic values.

Product(a, b, c, ...)

Represents a product of an arbitrary (finite) number of symbolic values.

# Special methods

To make the library more interesting, I included some special methods that do some interesting stuff with the self-expression. This methods are described bellow.

* alternatives()

The most exciting feature is the support for alternative representations, which uses common mathematical identities to create symbolically equivalent expressions from the self-expression, returning an array with the alternative representations, each as a distinct object.
Example:
  Exp(Log(3) * 2) -> alternatives.each { .say }
Output:
  Exp(Product(Log(3), 2))
  Power(3, 2)
  9

* pretty()

This method returns a human-readable representation for the self-expression, recursively calling the pretty() method on each object value.

Example: 
  Exp(4)**Log(Power(Fraction(1, 2), 3)) -> pretty.say
Output: 
  exp(log((1 / 2)^3) * 4)

* simple()

Based on the alternatives() method, the simple() method returns the shortest expression from the list of all the possible alternative expressions.

In selecting the shortest expressions, it also calls the pretty() method and computes the length of the human-readable text of each alternative expression. Then, based on this lengths, it selects the shortest alternative expression and returns it.
Example:
  Exp(Log(Log(Exp(Exp(Log(Symbol('x'))))))) -> simple.say
Output:
  Symbol('x')

* numeric()

The library also supports numerical evaluation, recursively evaluating each expression numerically.

Example:
  Fraction(Power(3, 42), 5) -> numeric.say
Output:
  21883797826302471841.8

# Bacovia

The Bacovia base class implements all the generic operations, which are the default operations for all the other classes that inherit from it.
class Bacovia {
    method +(Bacovia o) {
        Sum(self, o)
    }
    method -(Bacovia o) {
        Difference(self, o)
    }
    method *(Bacovia o) {
        Product(self, o)
    }
    method /(Bacovia o) {
        Fraction(self, o)        
    }
    method **(Bacovia o) {
        Power(self, o)
    }
    method neg {
        Difference(0, self)
    }
    method inv {
        Fraction(1, self)
    }
    method alternatives {
        [self]
    }
    method simple {
        self.alternatives.min_by {
            .pretty.len
        }
    }
}

# Power

The Power class takes two arguments: the base and the exponent, which can be implemented something like this:
class Power(v, n) {
    method **(Bacovia o) {
        Power(v, n * o)
    }
    method ==(Power o) {
        (v == o.v) &&
        (n == o.n)
    }
    method ==(_) {
        false
    }
    method inv {
        Power(v, n.neg)
    }
    method alternatives() is cached {
        gather {
            for a,b in (v.alternatives ~X n.alternatives) {
                take(Exp(Log(a) * b))
                take(Power(a, b))
                take(a ** b)

                # Identity: x^log(y) = y^log(x)
                if (b.kind_of(Log)) {
                    take(b.v**Log(a))
                }
            }
        }.uniq_by { .to_s }
    }
    method numeric {
        v.numeric ** n.numeric
    }
    method pretty() is cached {
        "(#{v.pretty})^(#{n.pretty})"
    }
    method to_s {
        "Power(#{v}, #{n})"
    }
}
Most of the other classes are implemented in a similar way, which for brevity, I will not include them in this post.

# Product

The Product and the Sum classes are slightly different that the other classes: they can take an arbitrary number of arguments which are internally stored inside an array. The interesting part comes in manipulating all of this values symbolically.
class Product(*values) {
    method *(Product o) {
        Product(values..., o.values...)
    }
    method /(Product o) {
        Product(values..., o.values.map{.inv}...)
    }
    method /(Bacovia o) {
        var copy = [values...]
        if (copy.remove_first(o)) {
            Product(copy...)
        }
        else {
            Product(copy..., o.inv)
        }
    }
    method *(Bacovia o) {
        var copy = [values...]
        if (copy.remove_first(o.inv)) {
            Product(copy...)
        }
        else {
            Product(copy..., o)
        }
    }
    method **(Bacovia o) {
        Product(values.map{ _ ** o }...)
    }
    method neg {
        Product(-1, values...)
    }
    method alternatives() is cached {
        gather {
            values.map{.alternatives}.cartesian { |*v|
                take(v.reduce('*'))
            }
        }.uniq_by { .to_s }
    }
    method ==(Product o) {
        values == o.values    
    }
    method ==(_) {
        false
    }
    method inv {
        Product(values.map { .inv }...)
    }
    method numeric {
        values.map { .numeric }.prod
    }
    method pretty() is cached {
        '(' + values.map{.pretty}.join(' * ') + ')'
    }
    method to_s {
        "Product(#{values.join(', ')})"
    }
}

Here, in the alternatives() method, we use the Cartesian product, which generates all the possible combinations of alternative representations. In combination with simple(), it can find the shortest possible alternative representation, based on the identities specified in each class.

# Examples

The library can be used for various things, such as finding alternative representations for a certain expression, or for symbolically simplifying expressions, making them shorter and easier to write them down.

Bellow I included a few examples to give the reader an idea for what can be done with this library. First, the library can be included using the include() keyword and the path to the main file, "bacovia.sf":

include('lib/bacovia.sf')

As a first example, we implement a closed-form to the Riemann zeta function for positive even integers, as defined here:
const τ = (Log(-1) * -2i)

func zeta_2n(n) {
    (-1)**(n+1) * bernoulli(2*n) * τ**(2*n) / ((2*n)! * 2)
}

for n in (1..5) {
    say zeta_2n(n).pretty
}

Output:
log(-1)^2  *    -4 * (1 /  6) * (1 / 4)
log(-1)^4  *    16 * (1 / 30) * (1 / 48)
log(-1)^6  *   -64 * (1 / 42) * (1 / 1440)
log(-1)^8  *   256 * (1 / 30) * (1 / 80640)
log(-1)^10 * -1024 * (5 / 66) * (1 / 7257600)

The second example illustrates the Fraction class, which includes some interesting identities for manipulating fractions in a non-reducible way.
var sum = Fraction(0, 1)

for n in (0..3) {
    sum += Fraction(1i**n, n+1)
    sum -> pretty.say
}
Output:
( 1 + 0i) / 1
( 2 + 1i) / 2
( 4 + 3i) / 6
(16 + 6i) / 24
(The real part of the numerator is A281964, and the imaginary part is A282132)

In this third example we take a look at some trigonometric functions applied on symbolic values. First, let's define a symbolic value with an associated numerical value.
var n = Symbol('n', 42)
Using the symbolic value defined above, we can start investigating two well-known trigonometric functions.
cos(Log(n)) -> simple.pretty.say
sin(Log(n)) -> simple.pretty.say
Output:
(n^-i + n^i) / 2
(n^i - n^-i) / (2i)
The library also has built-in support for hyperbolic trigonometric functions (including an inverse to each trigonometric function).
cosh(Log(n)) -> simple.pretty.say
sinh(Log(n)) -> simple.pretty.say
Output:
(1 + n^2) / (2 * n)
(n^2 - 1) / (2 * n)
When a symbolic value has a numeric value associated with it, we have the possibility of evaluating the expression numerically.
cos(n) -> numeric.say   #=> -0.3999853149883512939...
sin(n) -> numeric.say   #=> -0.9165215479156337858...
Inverses to hyperbolic trigonometric functions:
asinh(sinh(n)) -> numeric.say      #=> 42
acosh(cosh(n)) -> numeric.say      #=> 42

# Implementations

The code of the library is freely available on GitHub.

 Sidef: https://github.com/trizen/bacovia
Perl: https://github.com/trizen/Math-Bacovia

The first library requires the Sidef programming language.

Comments

Post a Comment