Module: math.basic

Source: ./math/basic.reef

Overview

math/basic module

Basic mathematical operations for integers and floating-point numbers. All floating-point functions are implemented in pure Reef using IEEE 754 bit manipulation - no libc/libm dependency.

Algorithms ported from A2 Oberon MathL.Mod (ETH Oberon System)

Functions

fn math_LN2(): float

ln(2) - natural log of 2

fn math_LN10(): float

ln(10) - natural log of 10

fn math_TWO_PI(): float

2*pi

fn math_PI(): float

pi

fn math_PI_2(): float

pi/2

fn math_PI_4(): float

pi/4

fn pow2_200(): float

Exact powers of two used by sqrt_f's exponent scaling. Built from bits so they are exact, and both are comfortably normal.

fn pow2_minus_200(): float

fn math_EPS(): float

Epsilon for convergence test (approximately 2.2e-16)

fn math_HUGE(): float

Large positive value (near max float)

fn math_TINY(): float

Small positive value for comparisons

fn float_equal(x: float, y: float): bool

fn int_max(): int

INT_MAX for a 32-bit int. INT_MIN is -(INT_MAX) - 1 and is spelled out where needed; it has no positive counterpart, which is the whole reason abs is partial.

fn abs(x: int): option.Option[int]

|x|, or None for INT_MIN.

0 - INT_MIN is INT_MIN again in two's complement, so the old if x < 0 then 0 - x returned a NEGATIVE number from a function documented to return a magnitude -- and gcd and lcm both called it.

fn checked_mul(a: int, b: int): option.Option[int]

a * b, or None when it would not fit. The guard divides rather than multiplies, so it cannot itself overflow.

fn min(a: int, b: int): int

Returns the smaller of a and b

fn max(a: int, b: int): int

Returns the larger of a and b

fn square(x: int): option.Option[int]

x * x, or None when that overflows (|x| > 46340).

fn cube(x: int): option.Option[int]

x * x * x, or None when any step overflows.

fn pow_int(base: int, exp: int): option.Option[int]

base ** exp by squaring, or None on overflow.

A negative exp is None rather than 1. The old body documented "exp must be non-negative" and then let the while e > 0 loop fall through, so pow_int(2, -1) returned 1 -- a wrong answer for an input the comment already called out as unsupported.

fn clamp(value: int, lo: int, hi: int): int

Clamps value to the range [lo, hi]

fn sign(x: int): int

Returns -1 if x < 0, 0 if x == 0, 1 if x > 0

fn gcd(a: int, b: int): option.Option[int]

Greatest common divisor, or None when the answer is |INT_MIN|.

Euclid runs on the SIGNED values and takes the magnitude only at the end. The old body called abs(a) and abs(b) up front, which handed INT_MIN to a function that could not represent its magnitude; the remainder sequence shrinks, so no intermediate here ever needs negating. None is therefore reachable only for gcd(INT_MIN, 0) and gcd(INT_MIN, INT_MIN), whose answer genuinely is 2147483648.

fn lcm(a: int, b: int): option.Option[int]

Least common multiple, or None on overflow.

reef-lang#65: this was abs(a * b) / gcd(a, b), which multiplies BEFORE dividing and so overflows on pairs whose LCM fits comfortably. Measured: lcm(196608, 65536) returned 0 -- because 196608 * 65536 is exactly 3 * 2^32, which wraps to 0 -- and 0 is also what lcm returns when an input is 0. The wrong answer was indistinguishable from a legitimate one.

Dividing first is the standard fix: a/g is exact because g divides a, and the product (a/g) * b overflows only when the true LCM is unrepresentable.

fn is_even(x: int): bool

Returns true if x is even

fn is_odd(x: int): bool

Returns true if x is odd

fn abs_f(x: float): float

Returns the absolute value of x

fn min_f(a: float, b: float): float

Returns the smaller of a and b

fn max_f(a: float, b: float): float

Returns the larger of a and b

fn trunc_f(x: float): float

Returns x truncated toward zero

fn floor_f(x: float): float

Returns the largest integer <= x (as float)

fn ceil_f(x: float): float

Returns the smallest integer >= x (as float)

fn round_f(x: float): float

Returns x rounded to the nearest integer (as float) Rounds half away from zero

fn exp_f(x: float): float

fn log_f(x: float): float

fn log10_f(x: float): float

fn next_up_pos(x: float): float

The next representable double above / below a POSITIVE finite value. For positive doubles the IEEE-754 bit pattern is monotone, so stepping the bits by one steps the value by exactly one ulp -- no magnitude constant and no special case per exponent, which matters because an ulp is not a fixed quantity. Only valid for x > 0 and finite; sqrt_correct guarantees both.

fn next_down_pos(x: float): float

fn square_err(a: float, p: float): float

Rounding error of aa, so that aa == (aa) + square_err(a, aa) EXACTLY.

Veltkamp splitting: multiplying by 2^27 + 1 and subtracting back splits a double into a high part with at most 26 significant bits and a low remainder, both exact. Products of the halves then fit a double with no rounding, so their combination recovers the bits the single rounding of a*a discarded. Needs only multiply and subtract -- no fma, no runtime help.

Safe against overflow here: sqrt_correct only calls it with c ~ sqrt(x), so c <= ~1.3e154 and 134217729 * c stays well under the double maximum.

fn sqrt_residual(x: float, c: float): float

x - cc, computed exactly enough to order candidates that differ by one ulp. cc is recovered exactly as p + err; x - p is itself exact whenever p is within a factor of two of x (Sterbenz), which holds because every candidate is within an ulp of the true root.

fn sqrt_correct(x: float, y: float): float

Pick the candidate whose square is nearest x. Comparing |cc - x| orders the candidates the same way as |c - sqrt(x)| does, because cc - x ~ 2c(c - sqrt(x)) and all three candidates share the same c to within an ulp.

fn sqrt_newton(x: float): float

Newton-Raphson to the floating-point fixed point.

reef-lang#57: this used to stop on float_equal(y, prev), an ABSOLUTE 2.2e-16. Below 1 that tolerance is wider than an ulp, so it exited while still 1-2 ulp out. The test is now exact equality, so the iteration runs until it genuinely stops moving.

fn sqrt_f(x: float): float

fn cbrt_f(x: float): float

fn is_inf_f(x: float): bool

True for +/-infinity only. Exponent 2047 with a zero mantissa is infinity; a non-zero mantissa there is NaN.

fn hypot_f(x: float, y: float): float

fn fmod_f(x: float, y: float): float

fn pow_f(base: float, exponent: float): float

fn sin_f(x: float): float

Sine function using Taylor series sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ...

fn cos_f(x: float): float

Cosine function using Taylor series cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ...

fn tan_f(x: float): float

Tangent function: tan(x) = sin(x) / cos(x)

fn atan_f(x: float): float

Arctangent: atan(x)

Reduction, in this order, so the Taylor series below is only ever evaluated where it converges fast: x < 0 : atan is odd, atan(x) = -atan(-x) x > 1 : atan(x) = pi/2 - atan(1/x) x > tan(pi/8): atan(x) = pi/4 + atan((x-1)/(x+1)) otherwise : Taylor series, |x| <= tan(pi/8)

reef-lang#56: the third step is the one that was missing, and without it this function was wrong by 5.0e-4 near |x| = 1. The series atan(x) = x - x^3/3 + x^5/5 - ... converges slower the closer |x| gets to 1, and AT x = 1 it IS the Leibniz series for pi/4, whose error after n terms is ~1/(2n). The loop below caps at 500 terms, so it exited on the cap rather than on convergence and returned ~1e-3 short. The function was not even monotonic across the seam: atan_f(1.0) returned a value below the true atan(0.999). Because x > 1 folds back into [0, 1] via the pi/2 identity, the whole neighbourhood of 1 was affected, not just the point.

It propagated structurally, which is why the fix belongs here and not in each caller: asin_f(x) = atan_f(x / sqrt(1 - x^2)) maps x = sqrt(2)/2 onto EXACTLY atan(1.0); acos_f is pi/2 - asin_f; and atan2_f(y, x) for x > 0 is atan_f(y/x), so y == x hits the same point.

Reducing at tan(pi/8) maps [tan(pi/8), 1] onto |t| <= tan(pi/8) ~ 0.4142, where successive terms shrink by t^2 ~ 0.17

The threshold literal below is the correctly rounded double for tan(pi/8). tan(pi/8) equals sqrt(2) - 1 in the reals, but do NOT substitute a computed sqrt_f(2.0) - 1.0: that lands 2 ulp away, and the threshold is compared against directly. Keep the literal. -- about 22 terms for 1e-17, comfortably inside the 500-term cap.

fn atan_series(x: float): float

Taylor series for atan, valid for any |x| <= 1 but only used by atan_f where |x| <= tan(pi/8) so that convergence is fast. Kept as its own function so the reduced-argument call above cannot accidentally re-enter the reduction and recurse.

fn atan2_f(y: float, x: float): float

Two-argument arctangent: atan2(y, x) = angle of point (x, y) Returns angle in radians in range (-pi, pi]

fn asin_f(x: float): float

Arcsine: asin(x) asin(x) = atan(x / sqrt(1 - x^2)) Domain: [-1, 1]

fn acos_f(x: float): float

Arccosine: acos(x) acos(x) = pi/2 - asin(x) Domain: [-1, 1]

fn sinh_f(x: float): float

Hyperbolic sine: sinh(x) = (e^x - e^(-x)) / 2

fn cosh_f(x: float): float

Hyperbolic cosine: cosh(x) = (e^x + e^(-x)) / 2

fn tanh_f(x: float): float

Hyperbolic tangent: tanh(x) = sinh(x) / cosh(x) = (e^x - e^(-x)) / (e^x + e^(-x))

fn log2_f(x: float): float

Generated by reefc doc