math.oak
// libmath implements basic arithmetic and algebraic functions
//
// For functions dealing with coordinate pairs and angles, the coordinate plane
// is assumed to be a Cartesian plane with +x to the east and +y to the north,
// where the angle is measured in radians from the +x axis, counterclockwise.
{
default: default
map: map
reduce: reduce
} := import('std')
{
sort: sort
} := import('sort')
// Pi, the circle constant
Pi := 3.14159265358979323846264338327950288419716939937510
// E, the base of the natural logarithm
E := 2.71828182845904523536028747135266249775724709369995
// sign returns -1 for all negative numbers, and 1 otherwise
fn sign(n) if n >= 0 {
true -> 1
_ -> -1
}
// abs returns the absolute value of a real number
fn abs(n) if n >= 0 {
true -> n
_ -> -n
}
// sqrt returns the principal square root of a real number, or ? if the number
// is negative.
fn sqrt(n) if n >= 0 -> pow(n, 0.5)
// hypot returns the Euclidean distance between two points, equivalent to the
// hypotenuse of a right triangle with the given two points as vertices.
fn hypot(x0, y0, x1, y1) {
if x1 = ? & y1 = ? -> x1 <- y1 <- 0
sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1))
}
// scale maps the value x in the range [a, b] to the range [c, d]. If [c, d]
// are not provided, they are assumed to be [0, 1]. x may be outside the range
// [a, b], in which case the value is scaled to be an equal amount outside of
// the range [c, d].
fn scale(x, a, b, c, d) {
normed := (x - a) / (b - a)
if {
c = ? & d = ? -> normed
_ -> (d - c) * normed + c
}
}
// bearing returns the point [x', y'] at the other end of a line segment
// starting at (x, y) and extending by distance d at angle t.
fn bearing(x, y, d, t) [
x + d * cos(t)
y + d * sin(t)
]
// orient returns the angle of the line extending from (x0, y0) to (x1, y1). If
// (x1, y1) is not provided, the given coordinate point is assumed to be (x1,
// y1) and (x0, y0) is assumed to be the origin (0, 0). Return values are in
// the range (-Pi, Pi]. This function is more commonly known in the form
// `atan2(y, x)` (note the reversed argument order).
fn orient(x0, y0, x1, y1) {
[x, y] := if x1 = ? & y1 = ? {
true -> [x0, y0]
_ -> [x1 - x0, y1 - y0]
}
if {
x > 0 -> 2 * atan(y / (hypot(x, y) + x))
x <= 0 & y != 0 -> 2 * atan((hypot(x, y) - x) / y)
x < 0 & y = 0 -> Pi
}
}
// sum takes a sequence of values and returns their sum
fn sum(xs...) xs |> reduce(0, fn(a, b) a + b)
// prod takes a sequence of values and returns their product
fn prod(xs...) xs |> reduce(1, fn(a, b) a * b)
// min returns the minimum value of all given values
fn min(xs...) xs |> reduce(xs.0, fn(acc, n) if n < acc {
true -> n
_ -> acc
})
// max returns the maximum value of all given values
fn max(xs...) xs |> reduce(xs.0, fn(acc, n) if n > acc {
true -> n
_ -> acc
})
// clamp returns a value bounded by some upper and lower bounds a and b. If the
// given x is between a and b, it is returned as-is; if it is outside the
// bounds, the closer of the two bounds is returned.
fn clamp(x, a, b) if {
x < a -> a
x > b -> b
_ -> x
}
// mean returns the arithmetic mean, or average, of all given values. If the
// list is empty, mean returns ?.
fn mean(xs) if len(xs) {
0 -> ?
_ -> sum(xs...) / len(xs)
}
// median returns the median, or "middle value", of all given values. If there
// is an even number of values given, median computes the mean of the middle
// two values in the list. If the list is empty, median returns ?.
fn median(xs) {
xs := sort(xs)
count := len(xs)
half := int(count / 2)
if 0 {
count -> ?
count % 2 -> (xs.(half - 1) + xs.(half)) / 2
_ -> xs.(half)
}
}
// stddev returns the population standard deviation computed from the given
// list of values. If the list is empty, stddev returns ?.
fn stddev(xs) if ? != xmean := mean(xs) -> {
xs |> map(fn(x) pow(xmean - x, 2)) |> mean() |> sqrt()
}
// round takes a number `n` and returns a floating-point number that represents
// `n` round to the nearest `decimals`-th decimal place. For negative values of
// `decimals`, no rounding occurs and `n` is returned exactly.
fn round(n, decimals) {
decimals := int(decimals) |> default(0)
if decimals < 0 {
true -> n
_ -> {
order := pow(10, decimals)
if n >= 0 {
true -> int(n * order + 0.5) / order
_ -> -int(-n * order + 0.5) / order
}
}
}
}