Improved documentation

book/build.py

@@ -69,6 +69,9 @@ def list_of_functions(file_name, document):
     with open(path, "w") as f:
         print(f"# {document['title']}\n", file=f, flush=True)
 
+        if introduction := document.get("introduction"):
+            print(introduction + "\n", file=f, flush=True)
+
         sections = document["sections"]
 
         if len(sections) >= 3:
@@ -115,22 +118,42 @@ list_of_functions(
                 "title": "Basics",
                 "modules": ["core::functions"],
             },
+            {
+                "title": "Transcendental functions",
+                "modules": ["math::transcendental"],
+            },
             {
                 "title": "Trigonometry",
-                "modules": ["math::functions", "math::trigonometry_extra"],
+                "modules": ["math::trigonometry"],
+            },
+            {
+                "title": "Statistics",
+                "modules": ["math::statistics"],
+            },
+            {
+                "title": "Random sampling, distributions",
+                "modules": ["core::random", "math::distributions"],
             },
             {
-                "title": "Random numbers",
-                "modules": ["core::random", "math::statistics"],
+                "title": "Number theory",
+                "modules": ["core::number_theory"],
             },
             {
                 "title": "Numerical methods",
                 "modules": ["numerics::diff", "numerics::solve"],
             },
+            {
+                "title": "Geometry",
+                "modules": ["math::geometry"],
+            },
             {
                 "title": "Algebra",
                 "modules": ["extra::algebra"],
             },
+            {
+                "title": "Trigonometry (extra)",
+                "modules": ["math::trigonometry_extra"],
+            },
         ],
     },
 )
@@ -163,14 +186,10 @@ list_of_functions(
     "datetime",
     {
         "title": "Date and time",
+        "introduction": "See [this page](./date-and-time.md) for a general introduction to date and time handling in Numbat.",
         "sections": [
             {
-                "title": "Basics",
-                "modules": ["datetime::functions"],
-            },
-            {
-                "title": "Human-readable",
-                "modules": ["datetime::human"],
+                "modules": ["datetime::functions", "datetime::human"],
             },
         ],
     },
@@ -186,6 +205,10 @@ list_of_functions(
                 "title": "Error handling",
                 "modules": ["core::error"],
             },
+            {
+                "title": "Floating point",
+                "modules": ["core::numbers"],
+            },
             {
                 "title": "Quantities",
                 "modules": ["core::quantities"],

book/src/list-functions-datetime.md

@@ -1,28 +1,25 @@
 # Date and time
 
-## Basics
+See [this page](./date-and-time.md) for a general introduction to date and time handling in Numbat.
 
-Defined in: `datetime::functions`
+Defined in: `datetime::functions`, `datetime::human`
 
 ### `now`
 Returns the current date and time.
-More information [here](https://numbat.dev/doc/date-and-time.html).
 
 ```nbt
 fn now() -> DateTime
 ```
 
 ### `datetime`
-Parses a string (date and time) into a DateTime object. See [this page](/doc/date-and-time.html#date-time-formats) for an overview of the supported formats.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Parses a string (date and time) into a `DateTime` object. See [here](./date-and-time.md#date-time-formats) for an overview of the supported formats.
 
 ```nbt
 fn datetime(input: String) -> DateTime
 ```
 
 ### `format_datetime`
-Formats a DateTime object as a string.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Formats a `DateTime` object as a string.
 
 ```nbt
 fn format_datetime(format: String, input: DateTime) -> String
@@ -30,7 +27,6 @@ fn format_datetime(format: String, input: DateTime) -> String
 
 ### `get_local_timezone`
 Returns the users local timezone.
-More information [here](https://numbat.dev/doc/date-and-time.html).
 
 ```nbt
 fn get_local_timezone() -> String
@@ -38,23 +34,20 @@ fn get_local_timezone() -> String
 
 ### `tz`
 Returns a timezone conversion function, typically used with the conversion operator.
-More information [here](https://numbat.dev/doc/date-and-time.html).
 
 ```nbt
 fn tz(tz: String) -> Fn[(DateTime) -> DateTime]
 ```
 
 ### `unixtime`
-Converts a DateTime to a UNIX timestamp.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Converts a `DateTime` to a UNIX timestamp. Can be used on the right hand side of a conversion operator: `now() -> unixtime`.
 
 ```nbt
 fn unixtime(input: DateTime) -> Scalar
 ```
 
 ### `from_unixtime`
-Converts a UNIX timestamp to a DateTime object.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Converts a UNIX timestamp to a `DateTime` object.
 
 ```nbt
 fn from_unixtime(input: Scalar) -> DateTime
@@ -62,32 +55,25 @@ fn from_unixtime(input: Scalar) -> DateTime
 
 ### `today`
 Returns the current date at midnight (in the local time).
-More information [here](https://numbat.dev/doc/date-and-time.html).
 
 ```nbt
 fn today() -> DateTime
 ```
 
 ### `date`
-Parses a string (only date) into a DateTime object.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Parses a string (only date) into a `DateTime` object.
 
 ```nbt
 fn date(input: String) -> DateTime
 ```
 
 ### `time`
-Parses a string (only time) into a DateTime object.
-More information [here](https://numbat.dev/doc/date-and-time.html).
+Parses a string (time only) into a `DateTime` object.
 
 ```nbt
 fn time(input: String) -> DateTime
 ```
 
-## Human-readable
-
-Defined in: `datetime::human`
-
 ### `human` (Human-readable time duration)
 Converts a time duration to a human-readable string in days, hours, minutes and seconds.
 More information [here](https://numbat.dev/doc/date-and-time.html).
@@ -95,3 +81,4 @@ More information [here](https://numbat.dev/doc/date-and-time.html).
 ```nbt
 fn human(time: Time) -> String
 ```
+

book/src/list-functions-math.md

@@ -1,83 +1,87 @@
 # Mathematical functions
 
-Jump to: [Basics](#basics) · [Trigonometry](#trigonometry) · [Random numbers](#random-numbers) · [Numerical methods](#numerical-methods) · [Algebra](#algebra)
+[Basics](#basics) · [Transcendental functions](#transcendental-functions) · [Trigonometry](#trigonometry) · [Statistics](#statistics) · [Random sampling, distributions](#random-sampling,-distributions) · [Number theory](#number-theory) · [Numerical methods](#numerical-methods) · [Geometry](#geometry) · [Algebra](#algebra) · [Trigonometry (extra)](#trigonometry-(extra))
 
 ## Basics
 
 Defined in: `core::functions`
 
 ### `id` (Identity function)
+Return the input value.
 
 ```nbt
 fn id<A>(x: A) -> A
 ```
 
 ### `abs` (Absolute value)
+Return the absolute value of the input, \\( |x| \\). This works for quantities, too: `abs(-5 m) = 5 m`.
 More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.abs).
 
 ```nbt
 fn abs<T: Dim>(x: T) -> T
 ```
 
-### `round` (Round)
-Round to the nearest integer.
-More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.round).
+### `sqrt` (Square root)
+Return the square root of the input, \\( \sqrt{x} \\): `sqrt(121 m^2) = 11 m`.
+More information [here](https://en.wikipedia.org/wiki/Square_root).
 
 ```nbt
-fn round<T: Dim>(x: T) -> T
+fn sqrt<D: Dim>(x: D^2) -> D
 ```
 
-### `floor` (Floor function)
-More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.floor).
+### `cbrt` (Cube root)
+Return the cube root of the input, \\( \sqrt[3]{x} \\): `cbrt(8 m^3) = 2 m`.
+More information [here](https://en.wikipedia.org/wiki/Cube_root).
 
 ```nbt
-fn floor<T: Dim>(x: T) -> T
+fn cbrt<D: Dim>(x: D^3) -> D
 ```
 
-### `ceil` (Ceil function)
-More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.ceil).
+### `sqr` (Square function)
+Return the square of the input, \\( x^2 \\): `sqr(5 m) = 25 m^2`.
 
 ```nbt
-fn ceil<T: Dim>(x: T) -> T
+fn sqr<D: Dim>(x: D) -> D^2
 ```
 
-### `mod` (Modulo)
-More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.rem_euclid).
+### `round` (Rounding)
+Round to the nearest integer. If the value is half-way between two integers, round away from 0.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.round).
 
 ```nbt
-fn mod<T: Dim>(a: T, b: T) -> T
+fn round<T: Dim>(x: T) -> T
 ```
 
-## Trigonometry
-
-Defined in: `math::functions`, `math::trigonometry_extra`
-
-### `is_nan`
+### `floor` (Floor function)
+Returns the largest integer less than or equal to `x`.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.floor).
 
 ```nbt
-fn is_nan<T: Dim>(n: T) -> Bool
+fn floor<T: Dim>(x: T) -> T
 ```
 
-### `is_infinite`
+### `ceil` (Ceil function)
+Returns the smallest integer greater than or equal to `x`.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.ceil).
 
 ```nbt
-fn is_infinite<T: Dim>(n: T) -> Bool
+fn ceil<T: Dim>(x: T) -> T
 ```
 
-### `sqrt` (Square root)
-More information [here](https://en.wikipedia.org/wiki/Square_root).
+### `mod` (Modulo)
+Calculates the least nonnegative remainder of `a (mod b)`.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.rem_euclid).
 
 ```nbt
-fn sqrt<D: Dim>(x: D^2) -> D
+fn mod<T: Dim>(a: T, b: T) -> T
 ```
 
-### `sqr` (Square function)
+## Transcendental functions
 
-```nbt
-fn sqr<D: Dim>(x: D) -> D^2
-```
+Defined in: `math::transcendental`
 
 ### `exp` (Exponential function)
+The exponential function, \\( e^x \\).
 More information [here](https://en.wikipedia.org/wiki/Exponential_function).
 
 ```nbt
@@ -85,6 +89,7 @@ fn exp(x: Scalar) -> Scalar
 ```
 
 ### `ln` (Natural logarithm)
+The natural logarithm with base \\( e \\).
 More information [here](https://en.wikipedia.org/wiki/Natural_logarithm).
 
 ```nbt
@@ -92,6 +97,7 @@ fn ln(x: Scalar) -> Scalar
 ```
 
 ### `log` (Natural logarithm)
+The natural logarithm with base \\( e \\).
 More information [here](https://en.wikipedia.org/wiki/Natural_logarithm).
 
 ```nbt
@@ -99,6 +105,7 @@ fn log(x: Scalar) -> Scalar
 ```
 
 ### `log10` (Common logarithm)
+The common logarithm with base \\( 10 \\).
 More information [here](https://en.wikipedia.org/wiki/Common_logarithm).
 
 ```nbt
@@ -106,12 +113,25 @@ fn log10(x: Scalar) -> Scalar
 ```
 
 ### `log2` (Binary logarithm)
+The binary logarithm with base \\( 2 \\).
 More information [here](https://en.wikipedia.org/wiki/Binary_logarithm).
 
 ```nbt
 fn log2(x: Scalar) -> Scalar
 ```
 
+### `gamma` (Gamma function)
+The gamma function, \\( \Gamma(x) \\).
+More information [here](https://en.wikipedia.org/wiki/Gamma_function).
+
+```nbt
+fn gamma(x: Scalar) -> Scalar
+```
+
+## Trigonometry
+
+Defined in: `math::trigonometry`
+
 ### `sin` (Sine)
 More information [here](https://en.wikipedia.org/wiki/Trigonometric_functions).
 
@@ -203,28 +223,26 @@ More information [here](https://en.wikipedia.org/wiki/Hyperbolic_functions).
 fn atanh(x: Scalar) -> Scalar
 ```
 
-### `gamma` (Gamma function)
-More information [here](https://en.wikipedia.org/wiki/Gamma_function).
+## Statistics
 
-```nbt
-fn gamma(x: Scalar) -> Scalar
-```
+Defined in: `math::statistics`
 
 ### `maximum` (Maxmimum)
-Get the largest element of a list.
+Get the largest element of a list: `maximum([30 cm, 2 m]) = 2 m`.
 
 ```nbt
 fn maximum<D: Dim>(xs: List<D>) -> D
 ```
 
 ### `minimum` (Minimum)
-Get the smallest element of a list.
+Get the smallest element of a list: `minimum([30 cm, 2 m]) = 30 cm`.
 
 ```nbt
 fn minimum<D: Dim>(xs: List<D>) -> D
 ```
 
 ### `mean` (Arithmetic mean)
+Calculate the arithmetic mean of a list of quantities: `mean([1 m, 2 m, 300 cm]) = 2 m`.
 More information [here](https://en.wikipedia.org/wiki/Arithmetic_mean).
 
 ```nbt
@@ -255,268 +273,261 @@ More information [here](https://en.wikipedia.org/wiki/Median).
 fn median<D: Dim>(xs: List<D>) -> D
 ```
 
-### `gcd` (Greatest common divisor)
-More information [here](https://en.wikipedia.org/wiki/Greatest_common_divisor).
+## Random sampling, distributions
 
-```nbt
-fn gcd(a: Scalar, b: Scalar) -> Scalar
-```
+Defined in: `core::random`, `math::distributions`
 
-### `lcm` (Least common multiple)
-More information [here](https://en.wikipedia.org/wiki/Least_common_multiple).
+### `random` (Standard uniform distribution sampling)
+Uniformly samples the interval \\( [0,1) \\).
 
 ```nbt
-fn lcm(a: Scalar, b: Scalar) -> Scalar
+fn random() -> Scalar
 ```
 
-### `hypot2`
+### `rand_uniform` (Continuous uniform distribution sampling)
+Uniformly samples the interval \\( [a,b) \\) if \\( a \le b \\) or \\( [b,a) \\) if \\( b<a \\) using inversion sampling.
+More information [here](https://en.wikipedia.org/wiki/Continuous_uniform_distribution).
 
 ```nbt
-fn hypot2<T: Dim>(x: T, y: T) -> T
+fn rand_uniform<T: Dim>(a: T, b: T) -> T
 ```
 
-### `hypot3`
+### `rand_int` (Discrete uniform distribution sampling)
+Uniformly samples the integers in the interval \\( [a, b] \\).
+More information [here](https://en.wikipedia.org/wiki/Discrete_uniform_distribution).
 
 ```nbt
-fn hypot3<T: Dim>(x: T, y: T, z: T) -> T
+fn rand_int(a: Scalar, b: Scalar) -> Scalar
 ```
 
-### `circle_area`
+### `rand_bernoulli` (Bernoulli distribution sampling)
+Samples a Bernoulli random variable, that is, \\( 1 \\) with probability \\( p \\), \\( 0 \\) with probability \\( 1-p \\). The parameter \\( p \\) must be a probability (\\( 0 \le p \le 1 \\)).
+More information [here](https://en.wikipedia.org/wiki/Bernoulli_distribution).
 
 ```nbt
-fn circle_area<L: Dim>(radius: L) -> L^2
+fn rand_bernoulli(p: Scalar) -> Scalar
 ```
 
-### `circle_circumference`
+### `rand_binom` (Binomial distribution sampling)
+Samples a binomial distribution by doing \\( n \\) Bernoulli trials with probability \\( p \\).
+              The parameter \\( n \\) must be a positive integer, the parameter \\( p \\) must be a probability (\\( 0 \le p \le 1 \\)).
+More information [here](https://en.wikipedia.org/wiki/Binomial_distribution).
 
 ```nbt
-fn circle_circumference<L: Dim>(radius: L) -> L
+fn rand_binom(n: Scalar, p: Scalar) -> Scalar
 ```
 
-### `sphere_area`
+### `rand_norm` (Normal distribution sampling)
+Samples a normal distribution with mean \\( \mu \\) and standard deviation \\( \sigma \\) using the Box-Muller transform.
+More information [here](https://en.wikipedia.org/wiki/Normal_distribution).
 
 ```nbt
-fn sphere_area<L: Dim>(radius: L) -> L^2
+fn rand_norm<T: Dim>(μ: T, σ: T) -> T
 ```
 
-### `sphere_volume`
+### `rand_geom` (Geometric distribution sampling)
+Samples a geometric distribution (the distribution of the number of Bernoulli trials with probability \\( p \\) needed to get one success) by inversion sampling. The parameter \\( p \\) must be a probability (\\( 0 \le p \le 1 \\)).
+More information [here](https://en.wikipedia.org/wiki/Geometric_distribution).
 
 ```nbt
-fn sphere_volume<L: Dim>(radius: L) -> L^3
+fn rand_geom(p: Scalar) -> Scalar
 ```
 
-### `cot`
+### `rand_poisson` (Poisson distribution sampling)
+Sampling a poisson distribution with rate \\( \lambda \\), that is, the distribution of the number of events occurring in a fixed interval if these events occur with mean rate \\( \lambda \\). The rate parameter \\( \lambda \\) must be non-negative.
+More information [here](https://en.wikipedia.org/wiki/Poisson_distribution).
 
 ```nbt
-fn cot(x: Scalar) -> Scalar
+fn rand_poisson(λ: Scalar) -> Scalar
 ```
 
-### `acot`
+### `rand_expon` (Exponential distribution sampling)
+Sampling an exponential distribution (the distribution of the distance between events in a Poisson process with rate \\( \lambda \\)) using inversion sampling. The rate parameter \\( \lambda \\) must be positive.
+More information [here](https://en.wikipedia.org/wiki/Exponential_distribution).
 
 ```nbt
-fn acot(x: Scalar) -> Scalar
+fn rand_expon<T: Dim>(λ: T) -> 1 / T
 ```
 
-### `coth`
+### `rand_lognorm` (Log-normal distribution sampling)
+Sampling a log-normal distribution, that is, a distribution whose logarithm is a normal distribution with mean \\( \mu \\) and standard deviation \\( \sigma \\).
+More information [here](https://en.wikipedia.org/wiki/Log-normal_distribution).
 
 ```nbt
-fn coth(x: Scalar) -> Scalar
+fn rand_lognorm(μ: Scalar, σ: Scalar) -> Scalar
 ```
 
-### `acoth`
+### `rand_pareto` (Pareto distribution sampling)
+Sampling a Pareto distribution with minimum value `min` and shape parameter \\( \alpha \\) using inversion sampling. Both parameters must be positive.
+More information [here](https://en.wikipedia.org/wiki/Pareto_distribution).
 
 ```nbt
-fn acoth(x: Scalar) -> Scalar
+fn rand_pareto<T: Dim>(α: Scalar, min: T) -> T
 ```
 
-### `secant`
+## Number theory
+
+Defined in: `core::number_theory`
+
+## Numerical methods
+
+Defined in: `numerics::diff`, `numerics::solve`
+
+### `diff` (Numerical differentiation)
+Compute the numerical derivative of the function \\( f \\) at point \\( x \\) using the central difference method.
+More information [here](https://en.wikipedia.org/wiki/Numerical_differentiation).
 
 ```nbt
-fn secant(x: Scalar) -> Scalar
+fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X) -> Y / X
 ```
 
-### `arcsecant`
+### `root_bisect` (Bisection method)
+Find the root of the function \\( f \\) in the interval \\( [x_1, x_2] \\) using the bisection method. The function \\( f \\) must be continuous and \\( f(x_1) \cdot f(x_2) < 0 \\).
+More information [here](https://en.wikipedia.org/wiki/Bisection_method).
 
 ```nbt
-fn arcsecant(x: Scalar) -> Scalar
+fn root_bisect<A: Dim, B: Dim>(f: Fn[(A) -> B], x1: A, x2: A, x_tol: A, y_tol: B) -> A
 ```
 
-### `cosecant`
+### `root_newton` (Newton's method)
+Find the root of the function \\( f(x) \\) and its derivative \\( f'(x) \\) using Newton's method.
+More information [here](https://en.wikipedia.org/wiki/Newton%27s_method).
 
 ```nbt
-fn cosecant(x: Scalar) -> Scalar
+fn root_newton<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B / A], x0: A, y_tol: B) -> A
 ```
 
-### `csc`
+## Geometry
+
+Defined in: `math::geometry`
+
+### `hypot2`
 
 ```nbt
-fn csc(x: Scalar) -> Scalar
+fn hypot2<T: Dim>(x: T, y: T) -> T
 ```
 
-### `acsc`
+### `hypot3`
 
 ```nbt
-fn acsc(x: Scalar) -> Scalar
+fn hypot3<T: Dim>(x: T, y: T, z: T) -> T
 ```
 
-### `sech`
+### `circle_area`
 
 ```nbt
-fn sech(x: Scalar) -> Scalar
+fn circle_area<L: Dim>(radius: L) -> L^2
 ```
 
-### `asech`
+### `circle_circumference`
 
 ```nbt
-fn asech(x: Scalar) -> Scalar
+fn circle_circumference<L: Dim>(radius: L) -> L
 ```
 
-### `csch`
+### `sphere_area`
 
 ```nbt
-fn csch(x: Scalar) -> Scalar
+fn sphere_area<L: Dim>(radius: L) -> L^2
 ```
 
-### `acsch`
+### `sphere_volume`
 
 ```nbt
-fn acsch(x: Scalar) -> Scalar
+fn sphere_volume<L: Dim>(radius: L) -> L^3
 ```
 
-## Random numbers
+## Algebra
 
-Defined in: `core::random`, `math::statistics`
+Defined in: `extra::algebra`
 
-### `random` (Standard uniform distribution sampling)
-Uniformly samples the interval [0,1).
+### `quadratic_equation` (Solve quadratic equations)
+Returns the solutions of the equation a x² + b x + c = 0.
+More information [here](https://en.wikipedia.org/wiki/Quadratic_equation).
 
 ```nbt
-fn random() -> Scalar
+fn quadratic_equation<A: Dim, B: Dim>(a: A, b: B, c: B^2 / A) -> List<B / A>
 ```
 
-### `rand_uniform` (Continuous uniform distribution sampling)
-Uniformly samples the interval [a,b) if a<=b or [b,a) if b<a using inversion sampling.
-More information [here](https://en.wikipedia.org/wiki/Continuous_uniform_distribution).
+## Trigonometry (extra)
 
-```nbt
-fn rand_uniform<T: Dim>(a: T, b: T) -> T
-```
+Defined in: `math::trigonometry_extra`
 
-### `rand_int` (Discrete uniform distribution sampling)
-Uniformly samples the integers in the interval [a, b].
-More information [here](https://en.wikipedia.org/wiki/Discrete_uniform_distribution).
+### `cot`
 
 ```nbt
-fn rand_int(a: Scalar, b: Scalar) -> Scalar
+fn cot(x: Scalar) -> Scalar
 ```
 
-### `rand_bernoulli` (Bernoulli distribution sampling)
-Samples a Bernoulli random variable, that is, 1 with probability p, 0 with probability 1-p.
-              Parameter p must be a probability (0 <= p <= 1).
-More information [here](https://en.wikipedia.org/wiki/Bernoulli_distribution).
+### `acot`
 
 ```nbt
-fn rand_bernoulli(p: Scalar) -> Scalar
+fn acot(x: Scalar) -> Scalar
 ```
 
-### `rand_binom` (Binomial distribution sampling)
-Samples a binomial distribution by doing n Bernoulli trials with probability p.
-              Parameter n must be a positive integer, parameter p must be a probability (0 <= p <= 1).
-More information [here](https://en.wikipedia.org/wiki/Binomial_distribution).
+### `coth`
 
 ```nbt
-fn rand_binom(n: Scalar, p: Scalar) -> Scalar
+fn coth(x: Scalar) -> Scalar
 ```
 
-### `rand_norm` (Normal distribution sampling)
-Samples a normal distribution with mean μ and standard deviation σ using the Box-Muller transform.
-More information [here](https://en.wikipedia.org/wiki/Normal_distribution).
+### `acoth`
 
 ```nbt
-fn rand_norm<T: Dim>(μ: T, σ: T) -> T
+fn acoth(x: Scalar) -> Scalar
 ```
 
-### `rand_geom` (Geometric distribution sampling)
-Samples a geometric distribution (the distribution of the number of Bernoulli trials with probability p needed to get one success) by inversion sampling.
-              Parameter p must be a probability (0 <= p <= 1).
-More information [here](https://en.wikipedia.org/wiki/Geometric_distribution).
+### `secant`
 
 ```nbt
-fn rand_geom(p: Scalar) -> Scalar
+fn secant(x: Scalar) -> Scalar
 ```
 
-### `rand_poisson` (Poisson distribution sampling)
-Sampling a poisson distribution with rate λ, that is, the distribution of the number of events occurring in a fixed interval if these events occur with mean rate λ.
-              The rate parameter λ must not be negative.
-More information [here](https://en.wikipedia.org/wiki/Poisson_distribution).
+### `arcsecant`
 
 ```nbt
-fn rand_poisson(λ: Scalar) -> Scalar
+fn arcsecant(x: Scalar) -> Scalar
 ```
 
-### `rand_expon` (Exponential distribution sampling)
-Sampling an exponential distribution (the distribution of the distance between events in a Poisson process with rate λ) using inversion sampling.
-              The rate parameter λ must be positive.
-More information [here](https://en.wikipedia.org/wiki/Exponential_distribution).
+### `cosecant`
 
 ```nbt
-fn rand_expon<T: Dim>(λ: T) -> 1 / T
+fn cosecant(x: Scalar) -> Scalar
 ```
 
-### `rand_lognorm` (Log-normal distribution sampling)
-Sampling a log-normal distribution, that is, a distribution whose log is a normal distribution with mean μ and standard deviation σ.
-More information [here](https://en.wikipedia.org/wiki/Log-normal_distribution).
+### `csc`
 
 ```nbt
-fn rand_lognorm(μ: Scalar, σ: Scalar) -> Scalar
+fn csc(x: Scalar) -> Scalar
 ```
 
-### `rand_pareto` (Pareto distribution sampling)
-Sampling a Pareto distribution with minimum value min and shape parameter α using inversion sampling.
-              Both parameters α and min must be positive.
-More information [here](https://en.wikipedia.org/wiki/Pareto_distribution).
+### `acsc`
 
 ```nbt
-fn rand_pareto<T: Dim>(α: Scalar, min: T) -> T
+fn acsc(x: Scalar) -> Scalar
 ```
 
-## Numerical methods
-
-Defined in: `numerics::diff`, `numerics::solve`
-
-### `diff` (Numerical differentiation)
-Compute the numerical derivative of a function at a point using the central difference method.
-More information [here](https://en.wikipedia.org/wiki/Numerical_differentiation).
+### `sech`
 
 ```nbt
-fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X) -> Y / X
+fn sech(x: Scalar) -> Scalar
 ```
 
-### `root_bisect` (Bisection method)
-Find the root of the function f in the interval [x1, x2] using the bisection method. The function f must be continuous and f(x1) × f(x2) < 0.
-More information [here](https://en.wikipedia.org/wiki/Bisection_method).
+### `asech`
 
 ```nbt
-fn root_bisect<A: Dim, B: Dim>(f: Fn[(A) -> B], x1: A, x2: A, x_tolerance: A, y_tolerance: B) -> A
+fn asech(x: Scalar) -> Scalar
 ```
 
-### `root_newton` (Newton's method)
-Find the root of the function f(x) and its derivative f'(x) using Newton's method.
-More information [here](https://en.wikipedia.org/wiki/Newton%27s_method).
+### `csch`
 
 ```nbt
-fn root_newton<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B / A], x0: A, y_tolerance: B) -> A
+fn csch(x: Scalar) -> Scalar
 ```
 
-## Algebra
-
-Defined in: `extra::algebra`
-
-### `quadratic_equation` (Solve quadratic equations)
-Returns the solutions of the equation a x² + b x + c = 0.
-More information [here](https://en.wikipedia.org/wiki/Quadratic_equation).
+### `acsch`
 
 ```nbt
-fn quadratic_equation<A: Dim, B: Dim>(a: A, b: B, c: B^2 / A) -> List<B / A>
+fn acsch(x: Scalar) -> Scalar
 ```
 

book/src/list-functions-other.md

@@ -1,6 +1,6 @@
 # Other functions
 
-Jump to: [Error handling](#error-handling) · [Quantities](#quantities) · [Chemical elements](#chemical-elements) · [Temperature conversion](#temperature-conversion)
+[Error handling](#error-handling) · [Floating point](#floating-point) · [Quantities](#quantities) · [Chemical elements](#chemical-elements) · [Temperature conversion](#temperature-conversion)
 
 ## Error handling
 
@@ -13,6 +13,26 @@ Throw an error with the specified message. Stops the execution of the program.
 fn error<T>(message: String) -> T
 ```
 
+## Floating point
+
+Defined in: `core::numbers`
+
+### `is_nan`
+Returns true if the input is `NaN`.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.is_nan).
+
+```nbt
+fn is_nan<T: Dim>(n: T) -> Bool
+```
+
+### `is_infinite`
+Returns true if the input is positive infinity or negative infinity.
+More information [here](https://doc.rust-lang.org/std/primitive.f64.html#method.is_infinite).
+
+```nbt
+fn is_infinite<T: Dim>(n: T) -> Bool
+```
+
 ## Quantities
 
 Defined in: `core::quantities`

numbat/examples/inspect.rs

@@ -56,8 +56,21 @@ fn inspect_functions_in_module(ctx: &Context, module: String) {
             println!("### `{fn_name}`");
         }
 
-        if let Some(ref description) = description {
-            let description = description.trim();
+        if let Some(ref description_raw) = description {
+            let description_raw = description_raw.trim().to_string();
+
+            // Replace $..$ with \\( .. \\) for mdbook.
+            let mut description = String::new();
+            for (i, part) in description_raw.split('$').enumerate() {
+                if i % 2 == 0 {
+                    description.push_str(part);
+                } else {
+                    description.push_str("\\\\( ");
+                    description.push_str(part);
+                    description.push_str(" \\\\)");
+                }
+            }
+
             if description.ends_with('.') {
                 println!("{description}");
             } else {

numbat/modules/all.nbt

@@ -6,7 +6,5 @@ use units::hartree
 
 use extra::algebra
 
-use math::trigonometry_extra
-
 use numerics::diff
 use numerics::solve

numbat/modules/core/functions.nbt

@@ -1,23 +1,42 @@
 @name("Identity function")
+@description("Return the input value.")
 fn id<A>(x: A) -> A = x
 
 @name("Absolute value")
+@description("Return the absolute value of the input, $|x|$. This works for quantities, too: `abs(-5 m) = 5 m`.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.abs")
 fn abs<T: Dim>(x: T) -> T
 
-@name("Round")
+@name("Square root")
+@description("Return the square root of the input, $\\sqrt\{x\}$: `sqrt(121 m^2) = 11 m`.")
+@url("https://en.wikipedia.org/wiki/Square_root")
+fn sqrt<D: Dim>(x: D^2) -> D = x^(1/2)
+
+@name("Cube root")
+@description("Return the cube root of the input, $\\sqrt[3]\{x\}$: `cbrt(8 m^3) = 2 m`.")
+@url("https://en.wikipedia.org/wiki/Cube_root")
+fn cbrt<D: Dim>(x: D^3) -> D = x^(1/3)
+
+@name("Square function")
+@description("Return the square of the input, $x^2$: `sqr(5 m) = 25 m^2`.")
+fn sqr<D: Dim>(x: D) -> D^2 = x^2
+
+@name("Rounding")
+@description("Round to the nearest integer. If the value is half-way between two integers, round away from 0.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.round")
-@description("Round to the nearest integer.")
 fn round<T: Dim>(x: T) -> T
 
 @name("Floor function")
+@description("Returns the largest integer less than or equal to `x`.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.floor")
 fn floor<T: Dim>(x: T) -> T
 
 @name("Ceil function")
+@description("Returns the smallest integer greater than or equal to `x`.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.ceil")
 fn ceil<T: Dim>(x: T) -> T
 
 @name("Modulo")
+@description("Calculates the least nonnegative remainder of `a (mod b)`.")
 @url("https://doc.rust-lang.org/std/primitive.f64.html#method.rem_euclid")
 fn mod<T: Dim>(a: T, b: T) -> T

numbat/modules/core/numbers.nbt

@@ -0,0 +1,7 @@
+@description("Returns true if the input is `NaN`.")
+@url("https://doc.rust-lang.org/std/primitive.f64.html#method.is_nan")
+fn is_nan<T: Dim>(n: T) -> Bool
+
+@description("Returns true if the input is positive infinity or negative infinity.")
+@url("https://doc.rust-lang.org/std/primitive.f64.html#method.is_infinite")
+fn is_infinite<T: Dim>(n: T) -> Bool

numbat/modules/core/random.nbt

@@ -1,5 +1,5 @@
 use core::scalar
 
 @name("Standard uniform distribution sampling")
-@description("Uniformly samples the interval [0,1).")
+@description("Uniformly samples the interval $[0,1)$.")
 fn random() -> Scalar

numbat/modules/datetime/functions.nbt

@@ -1,56 +1,44 @@
 use core::strings
 use units::si
 
-@url("https://numbat.dev/doc/date-and-time.html")
 @description("Returns the current date and time.")
 fn now() -> DateTime
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Parses a string (date and time) into a DateTime object. See [this page](/doc/date-and-time.html#date-time-formats) for an overview of the supported formats.")
+@description("Parses a string (date and time) into a `DateTime` object. See [here](./date-and-time.md#date-time-formats) for an overview of the supported formats.")
 fn datetime(input: String) -> DateTime
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Formats a DateTime object as a string.")
+@description("Formats a `DateTime` object as a string.")
 fn format_datetime(format: String, input: DateTime) -> String
 
-@url("https://numbat.dev/doc/date-and-time.html")
 @description("Returns the users local timezone.")
 fn get_local_timezone() -> String
 
-@url("https://numbat.dev/doc/date-and-time.html")
 @description("Returns a timezone conversion function, typically used with the conversion operator.")
 fn tz(tz: String) -> Fn[(DateTime) -> DateTime]
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Timezone conversion function targeting the users local timezone (datetime -> local).")
+@description("Timezone conversion function targeting the users local timezone (`datetime -> local`).")
 let local: Fn[(DateTime) -> DateTime] = tz(get_local_timezone())
 
-@url("https://numbat.dev/doc/date-and-time.html")
 @description("Timezone conversion function to UTC.")
 let UTC: Fn[(DateTime) -> DateTime] = tz("UTC")
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Converts a DateTime to a UNIX timestamp.")
+@description("Converts a `DateTime` to a UNIX timestamp. Can be used on the right hand side of a conversion operator: `now() -> unixtime`.")
 fn unixtime(input: DateTime) -> Scalar
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Converts a UNIX timestamp to a DateTime object.")
+@description("Converts a UNIX timestamp to a `DateTime` object.")
 fn from_unixtime(input: Scalar) -> DateTime
 
 fn _today_str() = format_datetime("%Y-%m-%d", now())
 
-@url("https://numbat.dev/doc/date-and-time.html")
 @description("Returns the current date at midnight (in the local time).")
 fn today() -> DateTime = datetime("{_today_str()} 00:00:00")
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description("Parses a string (only date) into a DateTime object.")
+@description("Parses a string (only date) into a `DateTime` object.")
 fn date(input: String) -> DateTime =
   if str_contains(input, " ")
     then datetime(str_replace(input, " ", " 00:00:00 "))
     else datetime("{input} 00:00:00")
 
-@url("https://numbat.dev/doc/date-and-time.html")
-@description(" Parses a string (only time) into a DateTime object.")
+@description("Parses a string (time only) into a `DateTime` object.")
 fn time(input: String) -> DateTime =
   datetime("{_today_str()} {input}")

numbat/modules/extra/algebra.nbt

@@ -1,5 +1,5 @@
 use core::error
-use math::functions
+use core::functions
 
 fn _qe_solution<A: Dim, B: Dim>(a: A, b: B, c: B² / A, sign: Scalar) -> B / A =
   (-b + sign × sqrt(b² - 4 a c)) / 2 a

numbat/modules/math/distributions.nbt

@@ -0,0 +1,96 @@
+use core::scalar
+use core::random
+use core::quantities
+use core::error
+use core::functions
+use math::constants
+use math::transcendental
+use math::trigonometry
+
+@name("Continuous uniform distribution sampling")
+@url("https://en.wikipedia.org/wiki/Continuous_uniform_distribution")
+@description("Uniformly samples the interval $[a,b)$ if $a \\le b$ or $[b,a)$ if $b<a$ using inversion sampling.")
+fn rand_uniform<T: Dim>(a: T, b: T) -> T =
+    if a <= b
+    then random() * (b - a) + a
+    else random() * (a - b) + b
+
+@name("Discrete uniform distribution sampling")
+@url("https://en.wikipedia.org/wiki/Discrete_uniform_distribution")
+@description("Uniformly samples the integers in the interval $[a, b]$.")
+fn rand_int(a: Scalar, b: Scalar) -> Scalar =
+    if a <= b
+    then floor( random() * (floor(b) - ceil(a) + 1) ) + ceil(a)
+    else floor( random() * (floor(a) - ceil(b) + 1) ) + ceil(b)
+
+@name("Bernoulli distribution sampling")
+@url("https://en.wikipedia.org/wiki/Bernoulli_distribution")
+@description("Samples a Bernoulli random variable, that is, $1$ with probability $p$, $0$ with probability $1-p$. The parameter $p$ must be a probability ($0 \le p \le 1$).")
+fn rand_bernoulli(p: Scalar) -> Scalar =
+    if p>=0 && p<=1
+    then (if random() < p
+        then 1
+        else 0)
+    else error("Argument p must be a probability (0 <= p <= 1).")
+
+@name("Binomial distribution sampling")
+@url("https://en.wikipedia.org/wiki/Binomial_distribution")
+@description("Samples a binomial distribution by doing $n$ Bernoulli trials with probability $p$.
+              The parameter $n$ must be a positive integer, the parameter $p$ must be a probability ($0 \le p \le 1$).")
+fn rand_binom(n: Scalar, p: Scalar) -> Scalar =
+    if n >= 1
+    then rand_binom(n-1, p) + rand_bernoulli(p)
+    else if n == 0
+    then 0
+    else error("Argument n must be a positive integer.")
+
+@name("Normal distribution sampling")
+@url("https://en.wikipedia.org/wiki/Normal_distribution")
+@description("Samples a normal distribution with mean $\\mu$ and standard deviation $\\sigma$ using the Box-Muller transform.")
+fn rand_norm<T: Dim>(μ: T, σ: T) -> T =
+    μ + sqrt(-2 σ² × ln(random())) × sin(2π × random())
+
+@name("Geometric distribution sampling")
+@url("https://en.wikipedia.org/wiki/Geometric_distribution")
+@description("Samples a geometric distribution (the distribution of the number of Bernoulli trials with probability $p$ needed to get one success) by inversion sampling. The parameter $p$ must be a probability ($0 \le p \le 1$).")
+fn rand_geom(p: Scalar) -> Scalar =
+    if p>=0 && p<=1
+    then ceil( ln(1-random()) / ln(1-p) )
+    else error("Argument p must be a probability (0 <= p <= 1).")
+
+# A helper function for rand_poisson, counts how many samples of the standard uniform distribution need to be multiplied to fall below lim.
+fn _poisson(lim: Scalar, prod: Scalar) -> Scalar =
+    if prod > lim
+    then _poisson(lim,  prod × random()) + 1
+    else -1
+
+@name("Poisson distribution sampling")
+@url("https://en.wikipedia.org/wiki/Poisson_distribution")
+@description("Sampling a poisson distribution with rate $\\lambda$, that is, the distribution of the number of events occurring in a fixed interval if these events occur with mean rate $\\lambda$. The rate parameter $\\lambda$ must be non-negative.")
+# This implementation is based on the exponential distribution of inter-arrival times. For details see L. Devroye, Non-Uniform Random Variate Generation, p. 504, Lemma 3.3.
+fn rand_poisson(λ: Scalar) -> Scalar =
+    if λ >= 0
+    then _poisson(exp(-λ), 1)
+    else error("Argument λ must not be negative.")
+
+@name("Exponential distribution sampling")
+@url("https://en.wikipedia.org/wiki/Exponential_distribution")
+@description("Sampling an exponential distribution (the distribution of the distance between events in a Poisson process with rate $\\lambda$) using inversion sampling. The rate parameter $\\lambda$ must be positive.")
+fn rand_expon<T: Dim>(λ: T) -> 1/T =
+    if value_of(λ) > 0
+    then - ln(1-random()) / λ
+    else error("Argument λ must be positive.")
+
+@name("Log-normal distribution sampling")
+@url("https://en.wikipedia.org/wiki/Log-normal_distribution")
+@description("Sampling a log-normal distribution, that is, a distribution whose logarithm is a normal distribution with mean $\\mu$ and standard deviation $\\sigma$.")
+fn rand_lognorm(μ: Scalar, σ: Scalar) -> Scalar =
+    exp( μ + σ × rand_norm(0, 1) )
+
+@name("Pareto distribution sampling")
+@url("https://en.wikipedia.org/wiki/Pareto_distribution")
+@description("Sampling a Pareto distribution with minimum value `min` and shape parameter $\\alpha$ using inversion sampling. Both parameters must be positive.")
+fn rand_pareto<T: Dim>(α: Scalar, min: T) -> T =
+    if value_of(min) > 0 && α > 0
+    then min / ((1-random())^(1/α))
+    else error("Both arguments α and min must be positive.")

numbat/modules/math/functions.nbt

@@ -1,173 +0,0 @@
-use core::scalar
-use core::lists
-use math::constants
-
-## Numerical properties
-
-fn is_nan<T: Dim>(n: T) -> Bool
-
-fn is_infinite<T: Dim>(n: T) -> Bool
-
-## Basics
-
-@name("Square root")
-@url("https://en.wikipedia.org/wiki/Square_root")
-fn sqrt<D: Dim>(x: D^2) -> D = x^(1/2)
-
-@name("Square function")
-fn sqr<D: Dim>(x: D) -> D^2 = x^2
-
-## Exponential and logarithm
-
-@name("Exponential function")
-@url("https://en.wikipedia.org/wiki/Exponential_function")
-fn exp(x: Scalar) -> Scalar
-
-@name("Natural logarithm")
-@url("https://en.wikipedia.org/wiki/Natural_logarithm")
-fn ln(x: Scalar) -> Scalar
-
-@name("Natural logarithm")
-@url("https://en.wikipedia.org/wiki/Natural_logarithm")
-fn log(x: Scalar) -> Scalar = ln(x)
-
-@name("Common logarithm")
-@url("https://en.wikipedia.org/wiki/Common_logarithm")
-fn log10(x: Scalar) -> Scalar
-
-@name("Binary logarithm")
-@url("https://en.wikipedia.org/wiki/Binary_logarithm")
-fn log2(x: Scalar) -> Scalar
-
-## Trigonometry
-
-@name("Sine")
-@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
-fn sin(x: Scalar) -> Scalar
-
-@name("Cosine")
-@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
-fn cos(x: Scalar) -> Scalar
-
-@name("Tangent")
-@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
-fn tan(x: Scalar) -> Scalar
-
-@name("Arc sine")
-@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
-fn asin(x: Scalar) -> Scalar
-
-@name("Arc cosine")
-@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
-fn acos(x: Scalar) -> Scalar
-
-@name("Arc tangent")
-@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
-fn atan(x: Scalar) -> Scalar
-
-@url("https://en.wikipedia.org/wiki/Atan2")
-fn atan2<T: Dim>(y: T, x: T) -> Scalar
-
-@name("Hyperbolic sine")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn sinh(x: Scalar) -> Scalar
-
-@name("Hyperbolic cosine")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn cosh(x: Scalar) -> Scalar
-
-@name("Hyperbolic tangent")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn tanh(x: Scalar) -> Scalar
-
-@name("Area hyperbolic sine")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn asinh(x: Scalar) -> Scalar
-
-@name("Area hyperbolic cosine")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn acosh(x: Scalar) -> Scalar
-
-@name("Area hyperbolic tangent ")
-@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
-fn atanh(x: Scalar) -> Scalar
-
-# Note: there are even more functions in `math::trigonmetry_extra`.
-
-## Others
-
-
-@name("Gamma function")
-@url("https://en.wikipedia.org/wiki/Gamma_function")
-fn gamma(x: Scalar) -> Scalar
-
-### Statistics
-
-# TODO: remove these helpers once we support local definitions
-fn _max<D: Dim>(x: D, y: D) -> D = if x > y then x else y
-fn _min<D: Dim>(x: D, y: D) -> D = if x < y then x else y
-
-@name("Maxmimum")
-@description("Get the largest element of a list")
-fn maximum<D: Dim>(xs: List<D>) -> D =
-  if len(xs) == 1
-    then head(xs)
-    else _max(head(xs), maximum(tail(xs)))
-
-@name("Minimum")
-@description("Get the smallest element of a list")
-fn minimum<D: Dim>(xs: List<D>) -> D =
-  if len(xs) == 1
-    then head(xs)
-    else _min(head(xs), minimum(tail(xs)))
-
-@name("Arithmetic mean")
-@url("https://en.wikipedia.org/wiki/Arithmetic_mean")
-fn mean<D: Dim>(xs: List<D>) -> D = if is_empty(xs) then 0 else sum(xs) / len(xs)
-
-@name("Variance")
-@url("https://en.wikipedia.org/wiki/Variance")
-@description("Calculate the population variance of a list of quantities")
-fn variance<D: Dim>(xs: List<D>) -> D^2 =
-  mean(map(sqr, xs)) - sqr(mean(xs))
-
-@name("Standard deviation")
-@url("https://en.wikipedia.org/wiki/Standard_deviation")
-@description("Calculate the population standard deviation of a list of quantities")
-fn stdev<D: Dim>(xs: List<D>) -> D = sqrt(variance(xs))
-
-@name("Median")
-@url("https://en.wikipedia.org/wiki/Median")
-@description("Calculate the median of a list of quantities")
-fn median<D: Dim>(xs: List<D>) -> D =  # TODO: this is extremely inefficient
-  if mod(len(xs), 2) == 1
-    then element_at((len(xs) - 1) / 2, sort(xs))
-    else mean([element_at(len(xs) / 2 - 1, sort(xs)), element_at(len(xs) / 2, sort(xs))])
-
-
-### Number theory
-
-@name("Greatest common divisor")
-@url("https://en.wikipedia.org/wiki/Greatest_common_divisor")
-fn gcd(a: Scalar, b: Scalar) -> Scalar =
-  if b == 0
-    then abs(a)
-    else gcd(b, mod(a, b))
-
-@name("Least common multiple")
-@url("https://en.wikipedia.org/wiki/Least_common_multiple")
-fn lcm(a: Scalar, b: Scalar) -> Scalar = abs(a * b) / gcd(a, b)
-
-### Geometry
-
-fn hypot2<T: Dim>(x: T, y: T) -> T = sqrt(x^2 + y^2)
-fn hypot3<T: Dim>(x: T, y: T, z: T) -> T = sqrt(x^2 + y^2 + z^2)
-
-# The following functions use a generic dimension instead of
-# 'Length' in order to allow for computations in pixels, for
-# example
-
-fn circle_area<L: Dim>(radius: L) -> L^2 = π × radius^2
-fn circle_circumference<L: Dim>(radius: L) -> L = 2 π × radius
-fn sphere_area<L: Dim>(radius: L) -> L^2 = 4 π × radius^2
-fn sphere_volume<L: Dim>(radius: L) -> L^3 = 4/3 × π × radius^3

numbat/modules/math/geometry.nbt

@@ -0,0 +1,14 @@
+use core::functions
+use math::constants
+
+fn hypot2<T: Dim>(x: T, y: T) -> T = sqrt(x^2 + y^2)
+fn hypot3<T: Dim>(x: T, y: T, z: T) -> T = sqrt(x^2 + y^2 + z^2)
+
+# The following functions use a generic dimension instead of
+# 'Length' in order to allow for computations in pixels, for
+# example
+
+fn circle_area<L: Dim>(radius: L) -> L^2 = π × radius^2
+fn circle_circumference<L: Dim>(radius: L) -> L = 2 π × radius
+fn sphere_area<L: Dim>(radius: L) -> L^2 = 4 π × radius^2
+fn sphere_volume<L: Dim>(radius: L) -> L^3 = 4/3 × π × radius^3

numbat/modules/math/number_theory.nbt

@@ -0,0 +1,13 @@
+use core::scalar
+use core::functions
+
+@name("Greatest common divisor")
+@url("https://en.wikipedia.org/wiki/Greatest_common_divisor")
+fn gcd(a: Scalar, b: Scalar) -> Scalar =
+  if b == 0
+    then abs(a)
+    else gcd(b, mod(a, b))
+
+@name("Least common multiple")
+@url("https://en.wikipedia.org/wiki/Least_common_multiple")
+fn lcm(a: Scalar, b: Scalar) -> Scalar = abs(a * b) / gcd(a, b)

numbat/modules/math/statistics.nbt

@@ -1,99 +1,43 @@
-use core::scalar
-use core::random
-use core::quantities
-use core::error
-use core::functions
-use math::functions
-
-@name("Continuous uniform distribution sampling")
-@url("https://en.wikipedia.org/wiki/Continuous_uniform_distribution")
-@description("Uniformly samples the interval [a,b) if a<=b or [b,a) if b<a using inversion sampling.")
-fn rand_uniform<T: Dim>(a: T, b: T) -> T =
-    if a <= b
-    then random() * (b - a) + a
-    else random() * (a - b) + b
-
-@name("Discrete uniform distribution sampling")
-@url("https://en.wikipedia.org/wiki/Discrete_uniform_distribution")
-@description("Uniformly samples the integers in the interval [a, b].")
-fn rand_int(a: Scalar, b: Scalar) -> Scalar =
-    if a <= b
-    then floor( random() * (floor(b) - ceil(a) + 1) ) + ceil(a)
-    else floor( random() * (floor(a) - ceil(b) + 1) ) + ceil(b)
-
-@name("Bernoulli distribution sampling")
-@url("https://en.wikipedia.org/wiki/Bernoulli_distribution")
-@description("Samples a Bernoulli random variable, that is, 1 with probability p, 0 with probability 1-p.
-              Parameter p must be a probability (0 <= p <= 1).")
-fn rand_bernoulli(p: Scalar) -> Scalar =
-    if p>=0 && p<=1
-    then (if random() < p
-        then 1
-        else 0)
-    else error("Argument p must be a probability (0 <= p <= 1).")
-
-@name("Binomial distribution sampling")
-@url("https://en.wikipedia.org/wiki/Binomial_distribution")
-@description("Samples a binomial distribution by doing n Bernoulli trials with probability p.
-              Parameter n must be a positive integer, parameter p must be a probability (0 <= p <= 1).")
-fn rand_binom(n: Scalar, p: Scalar) -> Scalar =
-    if n >= 1
-    then rand_binom(n-1, p) + rand_bernoulli(p)
-    else if n == 0
-    then 0
-    else error("Argument n must be a positive integer.")
-
-@name("Normal distribution sampling")
-@url("https://en.wikipedia.org/wiki/Normal_distribution")
-@description("Samples a normal distribution with mean μ and standard deviation σ using the Box-Muller transform.")
-fn rand_norm<T: Dim>(μ: T, σ: T) -> T =
-    μ + sqrt(-2 σ² × ln(random())) × sin(2π × random())
-
-@name("Geometric distribution sampling")
-@url("https://en.wikipedia.org/wiki/Geometric_distribution")
-@description("Samples a geometric distribution (the distribution of the number of Bernoulli trials with probability p needed to get one success) by inversion sampling.
-              Parameter p must be a probability (0 <= p <= 1).")
-fn rand_geom(p: Scalar) -> Scalar =
-    if p>=0 && p<=1
-    then ceil( ln(1-random()) / ln(1-p) )
-    else error("Argument p must be a probability (0 <= p <= 1).")
-
-# A helper function for rand_poisson, counts how many samples of the standard uniform distribution need to be multiplied to fall below lim.
-fn _poisson(lim: Scalar, prod: Scalar) -> Scalar =
-    if prod > lim
-    then _poisson(lim,  prod × random()) + 1
-    else -1
-
-@name("Poisson distribution sampling")
-@url("https://en.wikipedia.org/wiki/Poisson_distribution")
-@description("Sampling a poisson distribution with rate λ, that is, the distribution of the number of events occurring in a fixed interval if these events occur with mean rate λ.
-              The rate parameter λ must not be negative.")
-# This implementation is based on the exponential distribution of inter-arrival times. For details see L. Devroye, Non-Uniform Random Variate Generation, p. 504, Lemma 3.3.
-fn rand_poisson(λ: Scalar) -> Scalar =
-    if λ >= 0
-    then _poisson(exp(-λ), 1)
-    else error("Argument λ must not be negative.")
-
-@name("Exponential distribution sampling")
-@url("https://en.wikipedia.org/wiki/Exponential_distribution")
-@description("Sampling an exponential distribution (the distribution of the distance between events in a Poisson process with rate λ) using inversion sampling.
-              The rate parameter λ must be positive.")
-fn rand_expon<T: Dim>(λ: T) -> 1/T =
-    if value_of(λ) > 0
-    then - ln(1-random()) / λ
-    else error("Argument λ must be positive.")
-
-@name("Log-normal distribution sampling")
-@url("https://en.wikipedia.org/wiki/Log-normal_distribution")
-@description("Sampling a log-normal distribution, that is, a distribution whose log is a normal distribution with mean μ and standard deviation σ.")
-fn rand_lognorm(μ: Scalar, σ: Scalar) -> Scalar =
-    exp( μ + σ × rand_norm(0, 1) )
-
-@name("Pareto distribution sampling")
-@url("https://en.wikipedia.org/wiki/Pareto_distribution")
-@description("Sampling a Pareto distribution with minimum value min and shape parameter α using inversion sampling.
-              Both parameters α and min must be positive.")
-fn rand_pareto<T: Dim>(α: Scalar, min: T) -> T =
-    if value_of(min) > 0 && α > 0
-    then min / ((1-random())^(1/α))
-    else error("Both arguments α and min must be positive.")
+use core::lists
+
+# TODO: remove these helpers once we support local definitions
+fn _max<D: Dim>(x: D, y: D) -> D = if x > y then x else y
+fn _min<D: Dim>(x: D, y: D) -> D = if x < y then x else y
+
+@name("Maxmimum")
+@description("Get the largest element of a list: `maximum([30 cm, 2 m]) = 2 m`.")
+fn maximum<D: Dim>(xs: List<D>) -> D =
+  if len(xs) == 1
+    then head(xs)
+    else _max(head(xs), maximum(tail(xs)))
+
+@name("Minimum")
+@description("Get the smallest element of a list: `minimum([30 cm, 2 m]) = 30 cm`.")
+fn minimum<D: Dim>(xs: List<D>) -> D =
+  if len(xs) == 1
+    then head(xs)
+    else _min(head(xs), minimum(tail(xs)))
+
+@name("Arithmetic mean")
+@description("Calculate the arithmetic mean of a list of quantities: `mean([1 m, 2 m, 300 cm]) = 2 m`.")
+@url("https://en.wikipedia.org/wiki/Arithmetic_mean")
+fn mean<D: Dim>(xs: List<D>) -> D = if is_empty(xs) then 0 else sum(xs) / len(xs)
+
+@name("Variance")
+@url("https://en.wikipedia.org/wiki/Variance")
+@description("Calculate the population variance of a list of quantities")
+fn variance<D: Dim>(xs: List<D>) -> D^2 =
+  mean(map(sqr, xs)) - sqr(mean(xs))
+
+@name("Standard deviation")
+@url("https://en.wikipedia.org/wiki/Standard_deviation")
+@description("Calculate the population standard deviation of a list of quantities")
+fn stdev<D: Dim>(xs: List<D>) -> D = sqrt(variance(xs))
+
+@name("Median")
+@url("https://en.wikipedia.org/wiki/Median")
+@description("Calculate the median of a list of quantities")
+fn median<D: Dim>(xs: List<D>) -> D =  # TODO: this is extremely inefficient
+  if mod(len(xs), 2) == 1
+    then element_at((len(xs) - 1) / 2, sort(xs))
+    else mean([element_at(len(xs) / 2 - 1, sort(xs)), element_at(len(xs) / 2, sort(xs))])

numbat/modules/math/transcendental.nbt

@@ -0,0 +1,31 @@
+use core::scalar
+
+@name("Exponential function")
+@description("The exponential function, $e^x$.")
+@url("https://en.wikipedia.org/wiki/Exponential_function")
+fn exp(x: Scalar) -> Scalar
+
+@name("Natural logarithm")
+@description("The natural logarithm with base $e$.")
+@url("https://en.wikipedia.org/wiki/Natural_logarithm")
+fn ln(x: Scalar) -> Scalar
+
+@name("Natural logarithm")
+@description("The natural logarithm with base $e$.")
+@url("https://en.wikipedia.org/wiki/Natural_logarithm")
+fn log(x: Scalar) -> Scalar = ln(x)
+
+@name("Common logarithm")
+@description("The common logarithm with base $10$.")
+@url("https://en.wikipedia.org/wiki/Common_logarithm")
+fn log10(x: Scalar) -> Scalar
+
+@name("Binary logarithm")
+@description("The binary logarithm with base $2$.")
+@url("https://en.wikipedia.org/wiki/Binary_logarithm")
+fn log2(x: Scalar) -> Scalar
+
+@name("Gamma function")
+@description("The gamma function, $\\Gamma(x)$.")
+@url("https://en.wikipedia.org/wiki/Gamma_function")
+fn gamma(x: Scalar) -> Scalar

numbat/modules/math/trigonometry.nbt

@@ -0,0 +1,54 @@
+use core::scalar
+
+@name("Sine")
+@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
+fn sin(x: Scalar) -> Scalar
+
+@name("Cosine")
+@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
+fn cos(x: Scalar) -> Scalar
+
+@name("Tangent")
+@url("https://en.wikipedia.org/wiki/Trigonometric_functions")
+fn tan(x: Scalar) -> Scalar
+
+@name("Arc sine")
+@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
+fn asin(x: Scalar) -> Scalar
+
+@name("Arc cosine")
+@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
+fn acos(x: Scalar) -> Scalar
+
+@name("Arc tangent")
+@url("https://en.wikipedia.org/wiki/Inverse_trigonometric_functions")
+fn atan(x: Scalar) -> Scalar
+
+@url("https://en.wikipedia.org/wiki/Atan2")
+fn atan2<T: Dim>(y: T, x: T) -> Scalar
+
+@name("Hyperbolic sine")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn sinh(x: Scalar) -> Scalar
+
+@name("Hyperbolic cosine")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn cosh(x: Scalar) -> Scalar
+
+@name("Hyperbolic tangent")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn tanh(x: Scalar) -> Scalar
+
+@name("Area hyperbolic sine")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn asinh(x: Scalar) -> Scalar
+
+@name("Area hyperbolic cosine")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn acosh(x: Scalar) -> Scalar
+
+@name("Area hyperbolic tangent ")
+@url("https://en.wikipedia.org/wiki/Hyperbolic_functions")
+fn atanh(x: Scalar) -> Scalar
+
+# Note: there are even more functions in `math::trigonometry_extra`.

numbat/modules/math/trigonometry_extra.nbt

@@ -1,4 +1,8 @@
-use math::functions
+use core::scalar
+use core::functions
+use math::constants
+use math::trigonometry
+use math::transcendental
 
 fn cot(x: Scalar) -> Scalar = 1 / tan(x)
 fn acot(x: Scalar) -> Scalar = atan(1 / x)

numbat/modules/numerics/diff.nbt

@@ -5,6 +5,6 @@ fn _delta<X: Dim>(x: X) -> X = 1e-10 × unit_of(x)
 
 @name("Numerical differentiation")
 @url("https://en.wikipedia.org/wiki/Numerical_differentiation")
-@description("Compute the numerical derivative of a function at a point using the central difference method")
+@description("Compute the numerical derivative of the function $f$ at point $x$ using the central difference method.")
 fn diff<X: Dim, Y: Dim>(f: Fn[(X) -> Y], x: X) -> Y / X =
   (f(x + _delta(x)) - f(x - _delta(x))) / (2 _delta(x))

numbat/modules/numerics/solve.nbt

@@ -3,26 +3,26 @@ use core::error
 
 @name("Bisection method")
 @url("https://en.wikipedia.org/wiki/Bisection_method")
-@description("Find the root of the function f in the interval [x1, x2] using the bisection method. The function f must be continuous and f(x1) × f(x2) < 0.")
-fn root_bisect<A: Dim, B: Dim>(f: Fn[(A) -> B], x1: A, x2: A, x_tolerance: A, y_tolerance: B) -> A =
-  if abs(x1 - x2) < x_tolerance
+@description("Find the root of the function $f$ in the interval $[x_1, x_2]$ using the bisection method. The function $f$ must be continuous and $f(x_1) \cdot f(x_2) < 0$.")
+fn root_bisect<A: Dim, B: Dim>(f: Fn[(A) -> B], x1: A, x2: A, x_tol: A, y_tol: B) -> A =
+  if abs(x1 - x2) < x_tol
     then (x1 + x2) / 2
-    else if abs(f((x1 + x2) / 2)) < y_tolerance
+    else if abs(f((x1 + x2) / 2)) < y_tol
       then (x1 + x2) / 2
       else if f((x1 + x2) / 2) * f(x1) < 0
-        then root_bisect(f, x1, (x1 + x2) / 2, x_tolerance, y_tolerance)
-        else root_bisect(f, (x1 + x2) / 2, x2, x_tolerance, y_tolerance)
+        then root_bisect(f, x1, (x1 + x2) / 2, x_tol, y_tol)
+        else root_bisect(f, (x1 + x2) / 2, x2, x_tol, y_tol)
   # TODO: move (x1 + x2) / 2 to a local variable once we support them
 
-fn _root_newton_helper<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B/A], x0: A, y_tolerance: B, max_iterations: Scalar) -> A =
+fn _root_newton_helper<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B/A], x0: A, y_tol: B, max_iterations: Scalar) -> A =
   if max_iterations <= 0
     then error("root_newton: Maximum number of iterations reached. Try another initial guess?")
-    else if abs(f(x0)) < y_tolerance
+    else if abs(f(x0)) < y_tol
       then x0
-      else _root_newton_helper(f, f_prime, x0 - f(x0) / f_prime(x0), y_tolerance, max_iterations - 1)
+      else _root_newton_helper(f, f_prime, x0 - f(x0) / f_prime(x0), y_tol, max_iterations - 1)
 
 @name("Newton's method")
 @url("https://en.wikipedia.org/wiki/Newton%27s_method") 
-@description("Find the root of the function f(x) and its derivative f'(x) using Newton's method.")
-fn root_newton<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B/A], x0: A, y_tolerance: B) -> A =
-  _root_newton_helper(f, f_prime, x0, y_tolerance, 10_000)
+@description("Find the root of the function $f(x)$ and its derivative $f'(x)$ using Newton's method.")
+fn root_newton<A: Dim, B: Dim>(f: Fn[(A) -> B], f_prime: Fn[(A) -> B/A], x0: A, y_tol: B) -> A =
+  _root_newton_helper(f, f_prime, x0, y_tol, 10_000)

numbat/modules/physics/constants.nbt

@@ -1,4 +1,3 @@
-use math::functions
 use units::si
 
 @name("Speed of light in vacuum")

numbat/modules/prelude.nbt

@@ -6,11 +6,16 @@ use core::lists
 use core::strings
 use core::error
 use core::random
+use core::numbers
 
 use math::constants
-use math::functions
+use math::transcendental
+use math::trigonometry
 use math::trigonometry_extra
 use math::statistics
+use math::number_theory
+use math::distributions
+use math::geometry
 
 use units::si
 use units::time

numbat/modules/units/planck.nbt

@@ -1,5 +1,6 @@
 # https://en.wikipedia.org/wiki/Planck_units
 
+use core::functions
 use physics::constants
 
 @name("Planck length")

numbat/modules/units/stoney.nbt

@@ -1,3 +1,5 @@
+use core::functions
+use math::constants
 use physics::constants
 
 @name("Stoney length")