Trig.xtl

1⍝│ Trig: trigonometric and hyperbolic functions not in the base
2⍝│ language, written in closed form from its builtins (s_in, c_os,
3⍝│ a_tan, s_qrt, e_xp, e_xpm1, l_og, l_og1p). The hyperbolic tangent
4⍝│ is the builtin t_anh; Trig keeps no copy of it.
5⍝│ Import it with an alias of your choice: "tr:" u_se< "Trig".
6⍝│ Put libs/Trig/src on XETAL_PATH ("just path"); the reference is libs/Trig/docs.
7⍝│ Names with l: are exported; those under h: are private to this file.
8⍝│
9⍝│ Angles are in radians. Every function works on whole arrays, one
10⍝│ element at a time, and its results are Floats.
11
12⍝━━ Circular
13
14⍝│ Half of pi, the angle of a quarter turn.
15ʰq̲uarter ← { u → 0.5 × p̲i u }
16
17⍝│ t_an x: the tangent, the sine over the cosine. Near an odd multiple
18⍝│ of a quarter turn it is as large as the cosine there is small.
19ˡt̲an ← { x → (s̲in x) ÷ c̲os x }
20
21⍝│ a_sin x: the angle whose sine is x, from -pi/2 to pi/2, by the
22⍝│ half-angle form 2 atan(x / (1 + sqrt(1 - x^2))), whose divisor is
23⍝│ never below 1. Outside -1..1 the square root (s_qrt) refuses
24⍝│ (error[domain]).
25ˡa̲sin ← { x → y ← f̲loat x◆ 2.0 × a̲tan y ÷ 1.0 + s̲qrt 1.0 − y × y }
26
27⍝│ a_cos x: the angle whose cosine is x, from 0 to pi, as
28⍝│ 2 asin(sqrt((1 - x) / 2)), exact near both ends. Outside -1..1 the
29⍝│ square root (s_qrt) refuses (error[domain]).
30ˡa̲cos ← { x → 2.0 × ˡa̲sin s̲qrt 0.5 × 1.0 − f̲loat x }
31
32⍝│ y a_tan2 x: the angle of the point (x, y) from the positive x axis,
33⍝│ from -pi to pi: a_tan (y / x) turned by pi into the left half plane
34⍝│ (up for y >= 0, down for y < 0); on the y axis, pi/2 or -pi/2; at
35⍝│ the origin, 0.
36ˡa̲tan2 ← { y x →
37  b ← f̲loat y
38  a ← f̲loat x
39  z ← f̲loat a = 0.0
40  t ← (1.0 − z) × a̲tan b ÷ a + z
41  up ← (f̲loat b ≥ 0.0) − f̲loat b < 0.0
42  side ← (f̲loat b > 0.0) − f̲loat b < 0.0
43  t + ((p̲i @) × up × f̲loat a < 0.0) + (ʰq̲uarter @) × side × z
44}
45
46⍝━━ Hyperbolic
47
48⍝│ ((k p_ick lo)_ hi): lo where the mask k is 1.0, hi where it is 0.0.
49⍝│ The branch not picked must be finite (0 times an infinity is NaN),
50⍝│ so each caller clamps the argument it gives a branch.
51ʰp̲ick ← { k lo hi → (k × lo) + (1.0 − k) × hi }
52
53⍝│ The sign of x as a Float, 1.0 for 0: the odd functions take it from
54⍝│ their argument and work on its magnitude.
55ʰs̲ign ← { x → 1.0 − 2.0 × f̲loat x < 0.0 }
56
57⍝│ Where the large-argument forms take over: past it e^-a is below
58⍝│ 2e-9 of e^a (and below 1e-16 of 1/a^2 for the inverses).
59ʰb̲ig ← { u → 20.0 }
60
61⍝│ e^a / 2, as e^(a - ln 2) so it does not overflow until the result
62⍝│ does (a up to about 710.47, where e^a alone stops at 709.78).
63ʰh̲alfExp ← { a → e̲xp a − l̲og 2.0 }
64
65⍝│ s_inh x: the hyperbolic sine, (e^x - e^-x) / 2. Up to |x| = 20 as
66⍝│ (m + m / (m + 1)) / 2 with m = e^|x| - 1 (e_xpm1), so no
67⍝│ cancellation near 0; past it e^|x| / 2. Infinite past about 710.47.
68ˡs̲inh ← { x →
69  y ← f̲loat x
70  a ← a̲bs y
71  m ← e̲xpm1 a m̲in ʰb̲ig @
72  k ← f̲loat a ≤ ʰb̲ig @
73  lo ← 0.5 × m + m ÷ m + 1.0
74  hi ← ʰh̲alfExp a m̲ax ʰb̲ig @
75  (ʰs̲ign y) × ((k ʰp̲ick lo)_ hi)
76}
77
78⍝│ c_osh x: the hyperbolic cosine, (e^x + e^-x) / 2, never below 1.
79⍝│ Up to |x| = 20 as 1 + m^2 / (2 (m + 1)) with m = e^|x| - 1, past it
80⍝│ e^|x| / 2. Infinite past about 710.47.
81ˡc̲osh ← { x →
82  a ← a̲bs f̲loat x
83  m ← e̲xpm1 a m̲in ʰb̲ig @
84  k ← f̲loat a ≤ ʰb̲ig @
85  lo ← 1.0 + (m × m) ÷ 2.0 × m + 1.0
86  hi ← ʰh̲alfExp a m̲ax ʰb̲ig @
87  ((k ʰp̲ick lo)_ hi)
88}
89
90⍝│ Where the inverses' large-argument form takes over: past it
91⍝│ sqrt(a^2 + 1) is a to the last bit (and a^2 cannot overflow below).
92ʰh̲uge ← { u → 100000000.0 }
93
94⍝│ ln(2a), for a past h:h_uge: the inverses' large-argument form.
95ʰl̲ogTwice ← { a → (l̲og a) + l̲og 2.0 }
96
97⍝│ a_sinh x: the inverse hyperbolic sine, ln(x + sqrt(x^2 + 1)), for
98⍝│ every x. Up to |x| = 1e8 as log1p(a + a^2 / (1 + sqrt(1 + a^2)))
99⍝│ (l_og1p), so no cancellation near 0; past it ln(2a).
100ˡa̲sinh ← { x →
101  y ← f̲loat x
102  a ← a̲bs y
103  b ← a m̲in ʰh̲uge @
104  k ← f̲loat a ≤ ʰh̲uge @
105  lo ← l̲og1p b + (b × b) ÷ 1.0 + s̲qrt 1.0 + b × b
106  hi ← ʰl̲ogTwice a m̲ax ʰh̲uge @
107  (ʰs̲ign y) × ((k ʰp̲ick lo)_ hi)
108}
109
110⍝│ a_cosh x: the inverse hyperbolic cosine, ln(x + sqrt(x^2 - 1)), for
111⍝│ x of 1 or more, from 0 up. Up to 1e8 as log1p(t + sqrt(t) sqrt(t + 2))
112⍝│ with t = x - 1, exact near 1; past it ln(2x). Below 1 the square
113⍝│ root of t (s_qrt) refuses (error[domain]).
114ˡa̲cosh ← { x →
115  y ← f̲loat x
116  t ← (y m̲in ʰh̲uge @) − 1.0
117  k ← f̲loat y ≤ ʰh̲uge @
118  lo ← l̲og1p t + (s̲qrt t) × s̲qrt t + 2.0
119  hi ← ʰl̲ogTwice y m̲ax ʰh̲uge @
120  ((k ʰp̲ick lo)_ hi)
121}
122
123⍝│ a_tanh x: the inverse hyperbolic tangent, ln((1 + x) / (1 - x)) / 2,
124⍝│ for x strictly between -1 and 1, as log1p(2a / (1 - a)) / 2 on the
125⍝│ magnitude a (no cancellation near 0). At -1 or 1 the division
126⍝│ refuses (error[division-by-zero]); beyond them l_og1p does
127⍝│ (error[domain]).
128ˡa̲tanh ← { x →
129  y ← f̲loat x
130  a ← a̲bs y
131  (ʰs̲ign y) × 0.5 × l̲og1p (2.0 × a) ÷ 1.0 − a
132}
133
134⍝━━ Units and lengths
135
136⍝│ d_egrees x: the angle x (radians) in degrees, x * 180 / pi.
137ˡd̲egrees ← { x → (f̲loat x) × 180.0 ÷ p̲i @ }
138
139⍝│ r_adians x: the angle x (degrees) in radians, x * pi / 180.
140ˡr̲adians ← { x → (f̲loat x) × (p̲i @) ÷ 180.0 }
141
142⍝│ y h_ypot x: the length of (x, y), sqrt(x^2 + y^2), without
143⍝│ overflow or underflow: both are divided by the larger magnitude m
144⍝│ before squaring (so the squares are at most 1) and the root is
145⍝│ multiplied by m again; 0 when both are 0.
146ˡh̲ypot ← { y x →
147  a ← a̲bs f̲loat x
148  b ← a̲bs f̲loat y
149  m ← a m̲ax b
150  d ← m + f̲loat m = 0.0
151  m × s̲qrt ((a ÷ d) × a ÷ d) + (b ÷ d) × b ÷ d
152}