librarywork/xetal-libraries/libs/Trig/src/Trig.xtl
Trig: trigonometric and hyperbolic functions not in the base language, written in closed form from its builtins (s_in, c_os, a_tan, s_qrt, e_xp, e_xpm1, l_og, l_og1p). The hyperbolic tangent is the builtin t_anh; Trig keeps no copy of it. Import it with an alias of your choice: "tr:" u_se< "Trig". Put libs/Trig/src on XETAL_PATH ("just path"); the reference is libs/Trig/docs. Names with l: are exported; those under h: are private to this file.
Angles are in radians. Every function works on whole arrays, one element at a time, and its results are Floats.
Circular
ʰq̲uarter : Unit -> Float
Half of pi, the angle of a quarter turn.
ʰq̲uarter ← { u → 0.5 × p̲i u }
ˡt̲an : Num a => a -> Float
t_an x: the tangent, the sine over the cosine. Near an odd multiple of a quarter turn it is as large as the cosine there is small.
ˡt̲an ← { x → (s̲in x) ÷ c̲os x }
ˡa̲sin : Num a => a -> Float
a_sin x: the angle whose sine is x, from -pi/2 to pi/2, by the half-angle form 2 atan(x / (1 + sqrt(1 - x^2))), whose divisor is never below 1. Outside -1..1 the square root (s_qrt) refuses (error[domain]).
ˡa̲sin ← { x → y ← f̲loat x◆ 2.0 × a̲tan y ÷ 1.0 + s̲qrt 1.0 − y × y }
ˡa̲cos : Num a => a -> Float
a_cos x: the angle whose cosine is x, from 0 to pi, as 2 asin(sqrt((1 - x) / 2)), exact near both ends. Outside -1..1 the square root (s_qrt) refuses (error[domain]).
ˡa̲cos ← { x → 2.0 × ˡa̲sin s̲qrt 0.5 × 1.0 − f̲loat x }
ˡa̲tan2 : (Num a, Num b) => a -> b -> Float
y a_tan2 x: the angle of the point (x, y) from the positive x axis, from -pi to pi: a_tan (y / x) turned by pi into the left half plane (up for y >= 0, down for y < 0); on the y axis, pi/2 or -pi/2; at the origin, 0.
ˡa̲tan2 ← { y x → b ← f̲loat y a ← f̲loat x z ← f̲loat a = 0.0 t ← (1.0 − z) × a̲tan b ÷ a + z up ← (f̲loat b ≥ 0.0) − f̲loat b < 0.0 side ← (f̲loat b > 0.0) − f̲loat b < 0.0 t + ((p̲i @) × up × f̲loat a < 0.0) + (ʰq̲uarter @) × side × z }
Hyperbolic
ʰp̲ick : Float -> Float -> Float -> Float
((k p_ick lo)_ hi): lo where the mask k is 1.0, hi where it is 0.0. The branch not picked must be finite (0 times an infinity is NaN), so each caller clamps the argument it gives a branch.
ʰp̲ick ← { k lo hi → (k × lo) + (1.0 − k) × hi }
ʰs̲ign : Float -> Float
The sign of x as a Float, 1.0 for 0: the odd functions take it from their argument and work on its magnitude.
ʰs̲ign ← { x → 1.0 − 2.0 × f̲loat x < 0.0 }
ʰb̲ig : Any a => a -> Float
Where the large-argument forms take over: past it e^-a is below 2e-9 of e^a (and below 1e-16 of 1/a^2 for the inverses).
ʰb̲ig ← { u → 20.0 }
ʰh̲alfExp : Float -> Float
e^a / 2, as e^(a - ln 2) so it does not overflow until the result does (a up to about 710.47, where e^a alone stops at 709.78).
ʰh̲alfExp ← { a → e̲xp a − l̲og 2.0 }
ˡs̲inh : Num a => a -> Float
s_inh x: the hyperbolic sine, (e^x - e^-x) / 2. Up to |x| = 20 as (m + m / (m + 1)) / 2 with m = e^|x| - 1 (e_xpm1), so no cancellation near 0; past it e^|x| / 2. Infinite past about 710.47.
ˡs̲inh ← { x → y ← f̲loat x a ← a̲bs y m ← e̲xpm1 a m̲in ʰb̲ig @ k ← f̲loat a ≤ ʰb̲ig @ lo ← 0.5 × m + m ÷ m + 1.0 hi ← ʰh̲alfExp a m̲ax ʰb̲ig @ (ʰs̲ign y) × ((k ʰp̲ick lo)_ hi) }
ˡc̲osh : Num a => a -> Float
c_osh x: the hyperbolic cosine, (e^x + e^-x) / 2, never below 1. Up to |x| = 20 as 1 + m^2 / (2 (m + 1)) with m = e^|x| - 1, past it e^|x| / 2. Infinite past about 710.47.
ˡc̲osh ← { x → a ← a̲bs f̲loat x m ← e̲xpm1 a m̲in ʰb̲ig @ k ← f̲loat a ≤ ʰb̲ig @ lo ← 1.0 + (m × m) ÷ 2.0 × m + 1.0 hi ← ʰh̲alfExp a m̲ax ʰb̲ig @ ((k ʰp̲ick lo)_ hi) }
ʰh̲uge : Any a => a -> Float
Where the inverses' large-argument form takes over: past it sqrt(a^2 + 1) is a to the last bit (and a^2 cannot overflow below).
ʰh̲uge ← { u → 100000000.0 }
ʰl̲ogTwice : Num a => a -> Float
ln(2a), for a past h:h_uge: the inverses' large-argument form.
ʰl̲ogTwice ← { a → (l̲og a) + l̲og 2.0 }
ˡa̲sinh : Num a => a -> Float
a_sinh x: the inverse hyperbolic sine, ln(x + sqrt(x^2 + 1)), for every x. Up to |x| = 1e8 as log1p(a + a^2 / (1 + sqrt(1 + a^2))) (l_og1p), so no cancellation near 0; past it ln(2a).
ˡa̲sinh ← { x → y ← f̲loat x a ← a̲bs y b ← a m̲in ʰh̲uge @ k ← f̲loat a ≤ ʰh̲uge @ lo ← l̲og1p b + (b × b) ÷ 1.0 + s̲qrt 1.0 + b × b hi ← ʰl̲ogTwice a m̲ax ʰh̲uge @ (ʰs̲ign y) × ((k ʰp̲ick lo)_ hi) }
ˡa̲cosh : Num a => a -> Float
a_cosh x: the inverse hyperbolic cosine, ln(x + sqrt(x^2 - 1)), for x of 1 or more, from 0 up. Up to 1e8 as log1p(t + sqrt(t) sqrt(t + 2)) with t = x - 1, exact near 1; past it ln(2x). Below 1 the square root of t (s_qrt) refuses (error[domain]).
ˡa̲cosh ← { x → y ← f̲loat x t ← (y m̲in ʰh̲uge @) − 1.0 k ← f̲loat y ≤ ʰh̲uge @ lo ← l̲og1p t + (s̲qrt t) × s̲qrt t + 2.0 hi ← ʰl̲ogTwice y m̲ax ʰh̲uge @ ((k ʰp̲ick lo)_ hi) }
ˡa̲tanh : Num a => a -> Float
a_tanh x: the inverse hyperbolic tangent, ln((1 + x) / (1 - x)) / 2, for x strictly between -1 and 1, as log1p(2a / (1 - a)) / 2 on the magnitude a (no cancellation near 0). At -1 or 1 the division refuses (error[division-by-zero]); beyond them l_og1p does (error[domain]).
ˡa̲tanh ← { x → y ← f̲loat x a ← a̲bs y (ʰs̲ign y) × 0.5 × l̲og1p (2.0 × a) ÷ 1.0 − a }
Units and lengths
ˡd̲egrees : Num a => a -> Float
d_egrees x: the angle x (radians) in degrees, x * 180 / pi.
ˡd̲egrees ← { x → (f̲loat x) × 180.0 ÷ p̲i @ }
ˡr̲adians : Num a => a -> Float
r_adians x: the angle x (degrees) in radians, x * pi / 180.
ˡr̲adians ← { x → (f̲loat x) × (p̲i @) ÷ 180.0 }
ˡh̲ypot : (Num a, Num b) => a -> b -> Float
y h_ypot x: the length of (x, y), sqrt(x^2 + y^2), without overflow or underflow: both are divided by the larger magnitude m before squaring (so the squares are at most 1) and the root is multiplied by m again; 0 when both are 0.
ˡh̲ypot ← { y x → a ← a̲bs f̲loat x b ← a̲bs f̲loat y m ← a m̲ax b d ← m + f̲loat m = 0.0 m × s̲qrt ((a ÷ d) × a ÷ d) + (b ÷ d) × b ÷ d }