Math

Package Objects

This page describes functions in axle.logic and axle.math package objects.

Imports

import cats.implicits._

import spire.algebra._

import axle.logic._
import axle.math._

implicit val rngInt: Rng[Int] = spire.implicits.IntAlgebra
implicit val ringLong: Ring[Long] = spire.implicits.LongAlgebra
implicit val boolBoolean: Bool[Boolean] = spire.implicits.BooleanStructure

Logic aggregators ∃ and ∀:

∃(List(1, 2, 3)) { i: Int => i % 2 == 0 }
// res1: Boolean = true

∀(List(1, 2, 3)) { i: Int => i % 2 == 0 }
// res2: Boolean = false

Sum and multiply aggregators Σ and Π. Note that Σ and Π are also available in spire.optional.unicode._.

Σ((1 to 10) map { _ * 2 })
// res3: Int = 110

Π((1L to 10L) map { _ * 2 })
// res4: Long = 3715891200L

Doubles, triples, and cross-products

doubles(Set(1, 2, 3))
// res5: Seq[(Int, Int)] = List((1, 2), (1, 3), (2, 1), (2, 3), (3, 1), (3, 2))

triples(Set(1, 2, 3))
// res6: Seq[(Int, Int, Int)] = List(
//   (1, 2, 3),
//   (1, 3, 2),
//   (2, 1, 3),
//   (2, 3, 1),
//   (3, 1, 2),
//   (3, 2, 1)
// )

⨯(List(1, 2, 3))(List(4, 5, 6)).toList
// res7: List[(Int, Int)] = List(
//   (1, 4),
//   (1, 5),
//   (1, 6),
//   (2, 4),
//   (2, 5),
//   (2, 6),
//   (3, 4),
//   (3, 5),
//   (3, 6)
// )

Powerset

℘(0 until 4)
// res8: IndexedPowerSet[Int] = Iterable(
//   Set(),
//   Set(0),
//   Set(1),
//   Set(0, 1),
//   Set(2),
//   Set(0, 2),
//   Set(1, 2),
//   Set(0, 1, 2),
//   Set(3),
//   Set(0, 3),
//   Set(1, 3),
//   Set(0, 1, 3),
//   Set(2, 3),
//   Set(0, 2, 3),
//   Set(1, 2, 3),
//   Set(0, 1, 2, 3)
// )

val ps = ℘(Vector("a", "b", "c"))
// ps: IndexedPowerSet[String] = Iterable(
//   Set(),
//   Set("a"),
//   Set("b"),
//   Set("a", "b"),
//   Set("c"),
//   Set("a", "c"),
//   Set("b", "c"),
//   Set("a", "b", "c")
// )

ps.size
// res9: Int = 8

ps(7)
// res10: Set[String] = Set("a", "b", "c")

Permutations

permutations(0 until 4)(2).toList
// res11: List[IndexedSeq[Int]] = List(
//   Vector(0, 1),
//   Vector(0, 2),
//   Vector(0, 3),
//   Vector(1, 0),
//   Vector(1, 2),
//   Vector(1, 3),
//   Vector(2, 0),
//   Vector(2, 1),
//   Vector(2, 3),
//   Vector(3, 0),
//   Vector(3, 1),
//   Vector(3, 2)
// )

Combinations

combinations(0 until 4)(2).toList
// res12: List[IndexedSeq[Int]] = List(
//   Vector(0, 1),
//   Vector(0, 2),
//   Vector(0, 3),
//   Vector(1, 2),
//   Vector(1, 3),
//   Vector(2, 3)
// )

Indexed Cross Product

val icp = IndexedCrossProduct(Vector(
  Vector("a", "b", "c"),
  Vector("d", "e"),
  Vector("f", "g", "h")))
// icp: IndexedCrossProduct[String] = Iterable(
//   List("a", "d", "f"),
//   List("a", "d", "g"),
//   List("a", "d", "h"),
//   List("a", "e", "f"),
//   List("a", "e", "g"),
//   List("a", "e", "h"),
//   List("b", "d", "f"),
//   List("b", "d", "g"),
//   List("b", "d", "h"),
//   List("b", "e", "f"),
//   List("b", "e", "g"),
//   List("b", "e", "h"),
//   List("c", "d", "f"),
//   List("c", "d", "g"),
//   List("c", "d", "h"),
//   List("c", "e", "f"),
//   List("c", "e", "g"),
//   List("c", "e", "h")
// )

icp.size
// res13: Int = 18

icp(4)
// res14: Seq[String] = List("a", "e", "g")

Pi

Two estimators for π

import axle.math._

Wallis

The first is attributed to Englishman John Wallis (1616 - 1703) who published this function in 1655. It is quite slow.

wallisΠ(100).toDouble
// res16: Double = 3.1337874906281624

wallisΠ(200).toDouble
// res17: Double = 3.137677900950936

wallisΠ(400).toDouble
// res18: Double = 3.1396322219293964

wallisΠ(800).toDouble
// res19: Double = 3.1406116723489452

wallisΠ(1600).toDouble
// res20: Double = 3.1411019714193746

wallisΠ(3200).toDouble
// res21: Double = 3.141347264592393

Monte Carlo

import cats.implicits._
import spire.algebra.Field

implicit val fieldDouble: Field[Double] = spire.implicits.DoubleAlgebra

See the Wikipedia page on Monte Carlo Methods

This particular implementation requires that the number of trials be passed as a type F such that witnesses for typeclasses Aggregatable, Finite, and Functor are available in implicit scope.

While this may may seem initially over-engineered, it allows F as varied as List and Spark's RDD to be used to represent the number of trials and support the Monte Carlo simulation and resulting aggregation.

monteCarloPiEstimate((1 to 10000).toList, (n: Int) => n.toDouble)
// res22: Double = 3.1496

Fibonacci

import axle.math._

Linear using foldLeft

fibonacciByFold(10)
// res24: Int = 89

Recursive

fibonacciRecursively(10)
// res25: Int = 89

Some alternatives that are not in Axle include

Recursive with memoization

val memo = collection.mutable.Map(0 -> 0L, 1 -> 1L)

def fibonacciRecursivelyWithMemo(n: Int): Long = {
  if (memo.contains(n)) {
    memo(n)
  } else {
    val result = fibonacciRecursivelyWithMemo(n - 2) + fibonacciRecursivelyWithMemo(n - 1)
    memo += n -> result
    result
  }
}
fibonacciRecursivelyWithMemo(10)
// res26: Long = 55L

Recursive squaring

Imports

import org.jblas.DoubleMatrix

import cats.implicits._

import spire.algebra.EuclideanRing
import spire.algebra.NRoot
import spire.algebra.Rng

import axle._
import axle.jblas._

implicit val eucRingInt: EuclideanRing[Int] = spire.implicits.IntAlgebra
implicit val rngDouble: Rng[Double] = spire.implicits.DoubleAlgebra
implicit val nrootDouble: NRoot[Double] = spire.implicits.DoubleAlgebra
implicit val laJblasDouble = axle.jblas.linearAlgebraDoubleMatrix[Double]
import laJblasDouble._

The fibonacci sequence at N can be generated by taking the Nth power of a special 2x2 matrix. By employing the general-purpose strategy for exponentiation called "recursive squaring", we can achieve sub-linear time.

val base = fromColumnMajorArray(2, 2, List(1d, 1d, 1d, 0d).toArray)

def fibonacciSubLinear(n: Int): Long = n match {
  case 0 => 0L
  case _ => exponentiateByRecursiveSquaring(base, n).get(0, 1).toLong
}

Demo:

fibonacciSubLinear(78)
// res27: Long = 8944394323791464L

Note: Beyond 78 inaccuracies creep in due to the limitations of the Double number type.

Ackermann

See the Wikipedia page on the Ackermann function

import axle.math._

The computational complexity is enormous. Only for very small m and n can the function complete:

ackermann(1, 1)
// res29: Long = 3L

ackermann(3, 3)
// res30: Long = 61L

Future Work