# Myrddin: libmath

## Summary

Math is hard, let's go shopping.

``````pkg math =
trait fpmath @f =
exp : (x : @f -> @f)
expm1 : (x : @f -> @f)

fma : (x : @f, y : @f, z : @f -> @f)

log : (x : @f -> @f)
log1p : (x : @f -> @f)

horner_poly : (x : @f, a : @f[:] -> @f)
horner_polyu : (x : @f, a : @u[:] -> @f)

powr : (x : @f, y : @f -> @f)

scale2 : (x : @f, m : @i -> @f)

sin : (x : @f -> @f)
cos : (x : @f -> @f)
sincos : (x : @f -> (@f, @f))

sqrt : (x : @f -> @f)

kahan_sum : (a : @f[:] -> @f)
priest_sum : (a : @f[:] -> @f)

tan : (x : @f -> @f)
cot : (x : @f -> @f)

trunc : (f : @f -> @f)
ceil  : (f : @f -> @f)
floor : (f : @f -> @f)
;;

trait roundable @f -> @i =
rn : (f : @f -> @i)
;;

impl std.equatable flt32
impl std.equatable flt64
impl roundable flt64 -> int64
impl roundable flt32 -> int32
impl fpmath flt32
impl fpmath flt64
;;
``````

Libmath currently does not raise floating point exceptions. Only the rounding mode round-to-nearest (RN) is supported, and the only base is 2, in binary32 and binary64 formats.

## Exponentiation

``````const exp : (x : @f -> @f)
const expm1 : (x : @f -> @f)
``````

The `exp` function computes the base e exponential e^x. `expm1` computes e^x-1. `expm1` is accurate even for small values of `x`. Results should be correct to within 0.77 ulps.

``````const fma : (x : @f, y : @f, z : @f -> @f)
``````

The `fma` function computes `x*y+z` as a single operation, avoiding a rounding error and subtractive cancellation on the intermediate value. When supported, a specialized CPU instruction is used.

## Logarithms

``````log : (x : @f -> @f)
log1p : (x : @f -> @f)
``````

The `log` function computes the natural logarithm log(x). `log1p` computes log(1 + x). `log1p` is accurate even for small values of `x`. Results should be correct to within 0.6 ulps.

## Horner Polynomials

``````const horner_poly : (x : @f, a : @f[:] -> @f)
const horner_polyu : (x : @f, a : @u[:] -> @f)
``````

Evaluates a polynomial represented by the coefficients `a` at `x`. The `polyu` version accepts coefficients in bit form (as output from `flt64bits`).

Powr

``````powr : (x : @f, y : @f -> @f)
``````

The `powr` function computes `x^y`, with exceptions arising from considering the formula e^(log(y) * x). This is not the `pow` function of the C standard; negative integer values are not treated specially.

This function is relatively accurate for real-world values of `x` and `y`, but exhibits poor behavior when `x` is extremely high and `y` is extremely small, with errors up to 16 ulps observed.

## Scaling

``````scale2 : (x : @f, m : @i -> @f)
``````

Computes `2^m * x`, including cases when `2^m` would overflow or underflow floating-point representation, but the product is representable.

## Sine and Cosine

``````sin : (x : @f -> @f)
cos : (x : @f -> @f)
sincos : (x : @f -> (@f, @f))
``````

Compute sine and cosine of `x`. The calculation is exhaustively correct for `flt32` arguments, and is usually within 1 ulp for `flt64` arguments. All calculation, however, is performed in 64-bit precision, so `flt32` calculations are relatively slow.

## Sqrt

``````sqrt : (x : @f -> @f)
``````

Computes the square root of `x`. The calculation should be correctly rounded.

## Summation

``````kahan_sum : (a : @f[:] -> @f)
priest_sum : (a : @f[:] -> @f)
``````

Implements Kahan and Priest summation, respectively. These functions provide more accurate summation than simply adding each number in the sequence naively.

Each function sums a sequence of floats, returning the computed sum.

## Tangent and Cotangent

``````tan : (x : @f -> @f)
cot : (x : @f -> @f)
``````

Compute tangent and cotangent of `x`. The calculation is exhaustively correct for `flt32` arguments, and is usually within 1 ulp for `flt64` arguments, although errors of up to 4 ulps have been observed. All calculation is performed in 64-bit precision, so `flt32` calculations are relatively slow.

## Rounding

``````trunc : (x : @f -> @f)
ceil  : (x : @f -> @f)
floor : (x : @f -> @f)
rn : (x : @f -> @i)
``````

Trunc truncates its argument, or equivalently, rounds its argument downwards towards zero.

Floor and ceil compute the floor and ceiling of their argument, rounding towards negative infinity and infinity, respectively.

Rn rounds its argument to the nearest integer, rounding towards even numbers in order to break ties.