Libmath updates merged.

[Ori Bernstein <ori@xxxxxxxxxxxxxx>]

Hey,

For those watching the repo, you may have noticed that I merged npnth's work
on libmath, which contains the bulk of the trig functions. I think that it's
reached the point where we've actually got a usable libmath!

By way of summarizing changes, I'll just post the update to the libmath
documentation.


-- 
    Ori BernsteinFrom f051a05d9583d0f4327e403538e7addeaf39b60e Mon Sep 17 00:00:00 2001
From: "S. Gilles" <sgilles@xxxxxxxxxxxx>
Date: Mon, 23 Jul 2018 00:37:33 -0400
Subject: [PATCH] Document current libmath progress.

---
 doc/libmath/index.txt | 116 +++++++++++++++++++++++++++++++++---------
 1 file changed, 92 insertions(+), 24 deletions(-)

diff --git a/doc/libmath/index.txt b/doc/libmath/index.txt
index 5214068..61d5f1a 100644
--- a/doc/libmath/index.txt
+++ b/doc/libmath/index.txt
@@ -10,28 +10,39 @@ Math is hard, let's go shopping.
 
 	pkg math =
 		trait fpmath @f =
-			exp : (f : @f -> @f)
-			expm1 : (f : @f -> @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)
 
-			scale2 : (f : @f, m : @i -> @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 : (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 =
-			/* round-impl */
 			rn : (f : @f -> @i)
 		;;
 
@@ -43,25 +54,38 @@ Math is hard, let's go shopping.
 		impl fpmath flt64
 	;;
 
-Libmath currently does not raise floating point exceptions.
+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 : (f : @f -> @f)
-	const expm1 : (f : @f -> @f)
+	const exp : (x : @f -> @f)
+	const expm1 : (x : @f -> @f)
 
-The `exp` function computes the base *e* exponential
-e^f. `expm1` computes e^f-1. `expm1` is accurate even
-for small values of `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.
 
 Fused-mul-add
 -------------
 
 	const fma : (x : @f, y : @f, z : @f -> @f)
 
-The `fma` function computes `x*y+z` as a single operation,
-avoiding a rounding error on the intermediate value.
+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
 ------------------
@@ -70,19 +94,51 @@ Horner Polynomials
 	const horner_polyu : (x : @f, a : @u[:] -> @f)
 
 Evaluates a polynomial represented by the coefficients `a` at
-`x`.
+`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 : (f : @f, m : @i -> @f)
-	
+	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 : (f : @f -> @f)
+	sqrt : (x : @f -> @f)
 
-Computes the square root of `f`.
+Computes the square root of `x`. The calculation should be correctly
+rounded.
 
 Summation
 ---------
@@ -97,19 +153,31 @@ 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 : (f : @f -> @f)
-	ceil  : (f : @f -> @f)
-	floor : (f : @f -> @f)
-	rn : (f : @f -> @i)
+	trunc : (x : @f -> @f)
+	ceil  : (x : @f -> @f)
+	floor : (x : @f -> @f)
+	rn : (x : @f -> @i)
 
-Trunc truncates its argument, or equivalently, rounds it's argument downwards
+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.
+in order to break ties.

-- 
2.18.0