sourcelibs/Manifold/src/Manifold.xtl
1⍝│ Manifold: the unit sphere, the Poincare ball, and the pictures of
2⍝│ the circle, the sphere, the torus and the Swiss roll: the distance
3⍝│ measured along the space, putting a point onto it, the tangent part
4⍝│ of a vector at a point, a step from a point along a tangent vector
5⍝│ that stays on the space and the arrow that undoes it, a mean that
6⍝│ stays on the space, how far off it a point is; the points of each
7⍝│ space at given coordinates of its map, for drawing it; and the
8⍝│ distortion of an embedding's distances against those it should keep.
9⍝│ Import it with an alias of your choice: "mf:" u_se< "Manifold".
10⍝│ Put libs/Manifold/src on XETAL_PATH ("just path"); the reference is libs/Manifold/docs.
11⍝│ Names with l: are exported; those under h: are private to this file.
12⍝│
13⍝│ A point is a vector of Floats, in any number of dimensions; a set
14⍝│ of points is a matrix, one point per row. Each space has a name
15⍝│ that starts its functions (s_phere, b_all the Poincare ball,
16⍝│ c_ircle, t_orus, r_oll the Swiss roll), and the word after it is
17⍝│ the same verb for every space: Project (onto it), Distance (along
18⍝│ it), Tangent (the part of a vector along it at a point), Step (from
19⍝│ a point, along a tangent vector, staying on it), Arrow (the tangent
20⍝│ vector a Step would take from one point to another), Mean (of a
21⍝│ set, on it), Off (how far a point is from it); and for drawing,
22⍝│ Points (the points at rows of coordinates of the space's map), Arc
23⍝│ or Geodesic (points along the shortest path between two points).
24⍝│ So a program written for the sphere reads the same on the ball.
25⍝│ Every result is a Float.
26
27"tr:" u̲se< "Trig"
28
29⍝━━ Vectors
30
31⍝│ u d̲ot v: the sum of the products, as a Float.
32ʰd̲ot ← { u v → '+ r̲/ (f̲loat u) × f̲loat v }
33
34⍝│ n̲orm v: the length of v, straight through the space it sits in.
35ʰn̲orm ← { v → s̲qrt v ʰd̲ot v }
36
37⍝│ s v s̲caleRows M: each row of M scaled by the matching item of s.
38ʰs̲caleRows ← { s M → M × s 'l̲eft t̲able o̲ffsets 2 s̲elect s̲hape M }
39
40⍝│ m a̲ngle s: the angle, 0 to pi, between two unit vectors p and
41⍝│ q, from the lengths m = |p - q| and s = |p + q|: 2 atan2(m, s).
42⍝│ Accurate near 0 and near pi alike, where a_cos of the product
43⍝│ p.q is not (near 0 it loses half the digits), because neither
44⍝│ length is a difference of nearly equal numbers.
45ʰa̲ngle ← { m s → 2.0 × (f̲loat m) ᵗʳa̲tan2 f̲loat s }
46
47⍝│ p a̲part q: |p - q| and |p + q|, as a tuple.
48ʰa̲part ← { p q → (ʰn̲orm (f̲loat p) − f̲loat q, ʰn̲orm (f̲loat p) + f̲loat q) }
49
50⍝│ P a̲partRows p: |r - p| and |r + p| for every row r of P.
51ʰa̲partRows ← { P p →
52 Q ← ((f̲irst s̲hape P) r̲eshape 1.0) '× t̲able f̲loat p ⍝ p in every row
53 D ← (f̲loat P) − Q
54 S ← (f̲loat P) + Q
55 (s̲qrt '+ r̲/₂ D × D, s̲qrt '+ r̲/₂ S × S)
56}
57
58⍝│ a̲partPairs P: |r - r'| and |r + r'| for every pair of rows of
59⍝│ P, two matrices, summed one coordinate at a time (so no number is
60⍝│ a difference of nearly equal ones).
61ʰa̲partPairs ← { P →
62 n ← f̲irst s̲hape P
63 zero ← (n c̲at n) r̲eshape 0.0
64 (D, S, _) ← (2 s̲elect s̲hape P) '{ (d, s, k) →
65 c ← f̲loat k s̲elect₂ P
66 (d + (c '− t̲able c) ^ 2.0, s + (c '+ t̲able c) ^ 2.0, k + 1)
67 } p̲ower (zero, zero, 1)
68 (s̲qrt D, s̲qrt S)
69}
70
71⍝━━ The unit sphere
72
73⍝│ s̲phereProject v: v put onto the unit sphere, scaled to length 1;
74⍝│ the direction of v, kept. Given a set (a matrix), every row. The
75⍝│ zero vector has no direction and is refused.
76⍝│ >> "mf:" u_se< "Manifold"
77⍝│ >> mf:s_phereProject 3.0 4.0 0.0
78⍝│ 0.6 0.8 0.0
79⍝│ >> mf:s_phereProject 2 3 r_eshape 3.0 4.0 0.0 0.0 0.0 2.0
80⍝│ 0.6 0.8 0.0
81⍝│ 0.0 0.0 1.0
82ˡs̲phereProject ← { v →
83 2 = t̲ally s̲hape v ? ʰp̲rojectRows v
84 n ← ʰn̲orm v
85 n = 0.0 ? @ p̲anic< "the zero vector has no direction: it cannot be put onto the sphere"
86 (f̲loat v) ÷ n
87}
88
89⍝│ Every row of P put onto the sphere.
90ʰp̲rojectRows ← { P →
91 n ← s̲qrt '+ r̲/₂ (f̲loat P) × f̲loat P
92 '∨ r̲/ n = 0.0 ? @ p̲anic< "a zero row has no direction: it cannot be put onto the sphere"
93 (f̲loat P) ÷ n 'l̲eft t̲able o̲ffsets 2 s̲elect s̲hape P
94}
95
96⍝│ p s̲phereDistance q: the great-circle distance between two points
97⍝│ of the sphere, the angle between them in radians (0 to pi): the
98⍝│ shortest path along the sphere, as against the straight line
99⍝│ through it, which is shorter (its length is 2 sin of half the
100⍝│ angle). p and q must be on the sphere (s_phereProject).
101⍝│ >> "mf:" u_se< "Manifold"
102⍝│ >> 1.0 0.0 0.0 mf:s_phereDistance 0.0 1.0 0.0
103⍝│ 1.5707963267948966
104ˡs̲phereDistance ← { p q → (m, s) ← p ʰa̲part q◆ m ʰa̲ngle s }
105
106⍝│ s̲phereDistances P: the great-circle distance between every pair
107⍝│ of rows of P (points on the sphere): a matrix, 0 on its diagonal.
108⍝│ >> "mf:" u_se< "Manifold"
109⍝│ >> mf:s_phereDistances 2 2 r_eshape 1.0 0.0 0.0 1.0
110⍝│ 0.0 1.5707963267948966
111⍝│ 1.5707963267948966 0.0
112ˡs̲phereDistances ← { P → (M, S) ← ʰa̲partPairs P◆ M ʰa̲ngle S }
113
114⍝│ p s̲phereTangent v: the part of v that lies along the sphere at
115⍝│ p (in the tangent plane there): v less its part along p, the
116⍝│ normal. A tangent vector at p is one with no part along p.
117⍝│ >> "mf:" u_se< "Manifold"
118⍝│ >> 1.0 0.0 0.0 mf:s_phereTangent 2.0 3.0 0.0
119⍝│ 0.0 3.0 0.0
120ˡs̲phereTangent ← { p v → (f̲loat v) − (p ʰd̲ot v) × f̲loat p }
121
122⍝│ p s̲phereStep v: the point reached from p by going along the
123⍝│ sphere in the direction of the tangent vector v, for a distance
124⍝│ equal to v's length (the exponential map): cos t p + sin t v / t
125⍝│ with t = |v|, put back onto the sphere so no rounding takes it
126⍝│ off. A zero v goes nowhere. (A quarter turn from the x axis along
127⍝│ y lands on the y axis; the x left over is cos of a quarter turn,
128⍝│ 0 up to the machine's rounding.)
129⍝│ >> "mf:" u_se< "Manifold"
130⍝│ >> 1.0 0.0 0.0 mf:s_phereStep 0.0 1.5707963267948966 0.0
131⍝│ 0.00000000000000006123233995736766 1.0 0.0
132ˡs̲phereStep ← { p v →
133 t ← ʰn̲orm v
134 t = 0.0 ? f̲loat p
135 ˡs̲phereProject ((c̲os t) × f̲loat p) + (s̲in t) × (f̲loat v) ÷ t
136}
137
138⍝│ p s̲phereArrow q: the tangent vector at p that s_phereStep
139⍝│ carries to q (the logarithmic map): it points along the sphere
140⍝│ toward q and is as long as the distance between them. For q = p
141⍝│ it is zero; for q the opposite point, which every direction
142⍝│ reaches at distance pi, it is zero too (there is no one arrow).
143⍝│ >> "mf:" u_se< "Manifold"
144⍝│ >> 1.0 0.0 0.0 mf:s_phereArrow 0.0 1.0 0.0
145⍝│ 0.0 1.5707963267948966 0.0
146ˡs̲phereArrow ← { p q →
147 u ← p ˡs̲phereTangent q
148 n ← ʰn̲orm u
149 n = 0.0 ? u
150 (p ˡs̲phereDistance q) × u ÷ n
151}
152
153⍝│ The arrows at p to every row of P, one per row (s_phereArrow, over
154⍝│ a set at once): the mean's one step.
155ʰa̲rrows ← { p P →
156 c ← (f̲loat P) '+ '× i̲nner f̲loat p
157 U ← (f̲loat P) − c '× t̲able f̲loat p
158 n ← s̲qrt '+ r̲/₂ U × U
159 (m, s) ← P ʰa̲partRows p
160 ((m ʰa̲ngle s) ÷ n m̲ax 1e-300) ʰs̲caleRows U
161}
162
163⍝│ s̲phereMean P: the mean of the rows of P on the sphere (the
164⍝│ Frechet or Karcher mean): the point whose summed squared distances
165⍝│ along the sphere to the rows is least. The plain mean of the rows
166⍝│ lies inside the sphere; this starts from that mean put onto the
167⍝│ sphere and takes 30 steps, each along the average of the arrows
168⍝│ to the rows, which is a fixed point when the arrows cancel. Points
169⍝│ spread over less than a hemisphere have one mean, and the steps
170⍝│ reach it; rows whose plain mean is the zero vector (say, two
171⍝│ opposite points) are refused.
172⍝│ >> "mf:" u_se< "Manifold"
173⍝│ >> mf:s_phereMean 2 3 r_eshape 1.0 0.0 0.0 0.0 1.0 0.0
174⍝│ 0.7071067811865476 0.7071067811865476 0.0
175ˡs̲phereMean ← { P →
176 n ← f̲loat f̲irst s̲hape P
177 start ← ˡs̲phereProject ('+ r̲/ f̲loat P) ÷ n
178 30 '{ p → p ˡs̲phereStep ('+ r̲/ p ʰa̲rrows P) ÷ n } p̲ower start
179}
180
181⍝│ s̲phereOff p: how far p is from the sphere: the difference
182⍝│ between its length and 1, as a size. 0 for a point on it; what a
183⍝│ test of a constraint reports. Given a set (a matrix), one per row.
184⍝│ >> "mf:" u_se< "Manifold"
185⍝│ >> mf:s_phereOff 0.0 0.5 0.0
186⍝│ 0.5
187⍝│ >> mf:s_phereOff 2 2 r_eshape 0.0 0.5 3.0 4.0
188⍝│ 0.5 4.0
189ˡs̲phereOff ← { p →
190 2 = t̲ally s̲hape p ? a̲bs (s̲qrt '+ r̲/₂ (f̲loat p) × f̲loat p) − 1.0
191 a̲bs (ʰn̲orm p) − 1.0
192}
193
194⍝│ k c̲oordinate p: the k-th coordinate of a point, or the k-th
195⍝│ column of a set of points (a matrix), as Floats.
196ʰc̲oordinate ← { k p → 2 = t̲ally s̲hape p ? f̲loat k s̲elect₂ p◆ f̲loat k s̲elect p }
197
198⍝│ c̲olumns M: the columns of a (rows x n) matrix given as one vector
199⍝│ of n columns laid end to end, back as the matrix: the transpose of
200⍝│ the n x rows reshape. For writing a set of points coordinate by
201⍝│ coordinate.
202ʰc̲olumns ← { n v → o̲\ (n c̲at (t̲ally v) d̲iv n) r̲eshape v }
203
204⍝━━ The sphere, drawn
205
206⍝│ s̲pherePoints LL: the point of the unit sphere at each
207⍝│ (longitude, latitude) row of LL, in radians (longitude round the
208⍝│ z axis from the x axis, latitude up from the equator): one row
209⍝│ (x, y, z) per row of LL. These are the sphere's map: a parallel is
210⍝│ rows of one latitude and many longitudes, a meridian rows of one
211⍝│ longitude. (A quarter turn of latitude is the pole: its x is cos
212⍝│ of a quarter turn, 0 up to the machine's rounding.)
213⍝│ >> "mf:" u_se< "Manifold"
214⍝│ >> mf:s_pherePoints 2 2 r_eshape 0.0 0.0 0.0 1.5707963267948966
215⍝│ 1.0 0.0 0.0
216⍝│ 0.00000000000000006123233995736766 0.0 1.0
217ˡs̲pherePoints ← { LL →
218 lon ← 1 ʰc̲oordinate LL
219 lat ← 2 ʰc̲oordinate LL
220 3 ʰc̲olumns ((c̲os lat) × c̲os lon) c̲at ((c̲os lat) × s̲in lon) c̲at s̲in lat
221}
222
223⍝│ (p, q) s̲phereArc n: n + 1 points along the great circle from p
224⍝│ to q (both on the sphere), p first and q last, evenly spaced along
225⍝│ it: the points s_phereStep reaches from p along 0, 1/n, ..., 1 of
226⍝│ the arrow to q. For q = p, or q the opposite point (no one arc),
227⍝│ every point is p.
228⍝│ >> "mf:" u_se< "Manifold"
229⍝│ >> (1.0 0.0 0.0, 0.0 1.0 0.0) mf:s_phereArc 2
230⍝│ 1.0 0.0 0.0
231⍝│ 0.7071067811865476 0.7071067811865475 0.0
232⍝│ 0.00000000000000006123233995736766 1.0 0.0
233ˡs̲phereArc ← { (p, q) n →
234 v ← p ˡs̲phereArrow q
235 t ← ʰn̲orm v
236 rows ← (n + 1) r̲eshape 1.0
237 t = 0.0 ? rows '× t̲able f̲loat p
238 T ← (t ÷ f̲loat n) × (f̲loat r̲ange n + 1) − 1.0
239 ʰp̲rojectRows ((c̲os T) '× t̲able f̲loat p) + (s̲in T) '× t̲able v ÷ t
240}
241
242⍝━━ The circle
243
244⍝│ c̲irclePoints A: the point of the unit circle at each angle of A
245⍝│ (radians from the x axis, counterclockwise): one row (x, y) per
246⍝│ angle. The circle is the sphere of one dimension, the unit sphere
247⍝│ in two: s_phereProject puts a point onto it, s_phereDistance
248⍝│ measures along it (the angle between two points), s_phereArc draws
249⍝│ between them.
250⍝│ >> "mf:" u_se< "Manifold"
251⍝│ >> mf:c_irclePoints 0.0 3.141592653589793
252⍝│ 1.0 0.0
253⍝│ -1.0 0.00000000000000012246467991473532
254ˡc̲irclePoints ← { A → a ← f̲loat A◆ 2 ʰc̲olumns (c̲os a) c̲at s̲in a }
255
256⍝━━ The torus
257
258⍝│ radii t̲orusPoints UV: the point of the torus of radii R r (R
259⍝│ from the middle of the hole to the middle of the tube, r the
260⍝│ tube's) at each (u, v) row of UV, in radians (u round the hole, v
261⍝│ round the tube): ((R + r cos v) cos u, (R + r cos v) sin u,
262⍝│ r sin v), one row per row of UV. The torus's map is the (u, v)
263⍝│ rectangle, 0 to 2 pi each way, each pair of opposite edges joined.
264⍝│ >> "mf:" u_se< "Manifold"
265⍝│ >> 3.0 1.0 mf:t_orusPoints 2 2 r_eshape 0.0 0.0 0.0 3.141592653589793
266⍝│ 4.0 0.0 0.0
267⍝│ 2.0 0.0 0.00000000000000012246467991473532
268ˡt̲orusPoints ← { radii UV →
269 R ← f̲loat 1 s̲elect radii
270 r ← f̲loat 2 s̲elect radii
271 u ← 1 ʰc̲oordinate UV
272 v ← 2 ʰc̲oordinate UV
273 w ← R + r × c̲os v
274 3 ʰc̲olumns (w × c̲os u) c̲at (w × s̲in u) c̲at r × s̲in v
275}
276
277⍝│ radii t̲orusOff p: how far p is from the torus of radii R r: the
278⍝│ distance from p to the circle through the middle of the tube, less
279⍝│ r, as a size. 0 for a point on it. Given a set (a matrix), one per
280⍝│ row.
281⍝│ >> "mf:" u_se< "Manifold"
282⍝│ >> 3.0 1.0 mf:t_orusOff 0.0 0.0 0.0
283⍝│ 2.0
284⍝│ >> 3.0 1.0 mf:t_orusOff 2 3 r_eshape 4.0 0.0 0.0 3.0 0.0 0.0
285⍝│ 0.0 1.0
286ˡt̲orusOff ← { radii p →
287 R ← f̲loat 1 s̲elect radii
288 r ← f̲loat 2 s̲elect radii
289 x ← 1 ʰc̲oordinate p
290 y ← 2 ʰc̲oordinate p
291 z ← 3 ʰc̲oordinate p
292 rho ← s̲qrt (x × x) + y × y
293 a̲bs (s̲qrt ((rho − R) ^ 2.0) + z × z) − r
294}
295
296⍝━━ The Swiss roll
297
298⍝│ r̲ollPoints TH: the point of the Swiss roll at each (t, h) row of
299⍝│ TH: (t cos t, h, t sin t). The roll is a flat sheet rolled up:
300⍝│ along t it winds round the y axis, farther out with every turn (a
301⍝│ spiral whose radius is t); along h it runs straight up the axis.
302⍝│ Its map is the sheet itself, the rectangle of (r_ollAlong t, h).
303⍝│ >> "mf:" u_se< "Manifold"
304⍝│ >> mf:r_ollPoints 2 2 r_eshape 0.0 0.0 3.141592653589793 1.0
305⍝│ 0.0 0.0 0.0
306⍝│ -3.141592653589793 1.0 0.00000000000000038473413874435795
307ˡr̲ollPoints ← { TH →
308 t ← 1 ʰc̲oordinate TH
309 h ← 2 ʰc̲oordinate TH
310 3 ʰc̲olumns (t × c̲os t) c̲at h c̲at t × s̲in t
311}
312
313⍝│ r̲ollAlong t: the length of the roll's spiral from its middle
314⍝│ (t = 0) out to t, (t sqrt(1 + t^2) + asinh t) / 2: how far along
315⍝│ the unrolled sheet the points at t lie. Its map's first coordinate.
316⍝│ >> "mf:" u_se< "Manifold"
317⍝│ >> mf:r_ollAlong 0.0 1.0
318⍝│ 0.0 1.147793574696319
319ˡr̲ollAlong ← { t → x ← f̲loat t◆ ((x × s̲qrt 1.0 + x × x) + ᵗʳa̲sinh x) ÷ 2.0 }
320
321⍝│ a r̲ollDistance b: the distance along the sheet between the
322⍝│ points of the roll at (t, h) a and b. The sheet is flat, so it is
323⍝│ the plain distance between the two unrolled,
324⍝│ sqrt((along a - along b)^2 + (h_a - h_b)^2): two points a turn
325⍝│ apart are near in three dimensions and far along the sheet.
326⍝│ >> "mf:" u_se< "Manifold"
327⍝│ >> 0.0 0.0 mf:r_ollDistance 1.0 0.0
328⍝│ 1.147793574696319
329ˡr̲ollDistance ← { a b →
330 da ← (ˡr̲ollAlong 1 s̲elect a) − ˡr̲ollAlong 1 s̲elect b
331 dh ← (f̲loat 2 s̲elect a) − f̲loat 2 s̲elect b
332 s̲qrt (da × da) + dh × dh
333}
334
335⍝━━ The Poincare ball
336
337⍝│ u m̲obius v: Mobius addition, the ball's own way of adding (u
338⍝│ then v, each inside the unit ball): ((1 + 2 u.v + |v|^2) u +
339⍝│ (1 - |u|^2) v) / (1 + 2 u.v + |u|^2 |v|^2). The result stays in
340⍝│ the ball; (-u) m_obius v is v seen from u, as a vector at the
341⍝│ middle.
342ʰm̲obius ← { u v →
343 uv ← u ʰd̲ot v
344 uu ← u ʰd̲ot u
345 vv ← v ʰd̲ot v
346 (((1.0 + (2.0 × uv) + vv) × f̲loat u) + (1.0 − uu) × f̲loat v) ÷ (1.0 + (2.0 × uv) + uu × vv)
347}
348
349⍝│ p b̲allDistance q: the distance between two points of the
350⍝│ Poincare ball (the inside of the unit sphere with the hyperbolic
351⍝│ metric; in two dimensions the Poincare disk), 2 atanh of the
352⍝│ length of q seen from p, (-p) m_obius q. Near the middle it is
353⍝│ twice the plain distance; toward the rim it grows without bound,
354⍝│ so points near the rim a little apart in the plane are far apart
355⍝│ along the ball. p and q must be inside the ball (length under 1).
356⍝│ >> "mf:" u_se< "Manifold"
357⍝│ >> 0.0 0.0 mf:b_allDistance 0.5 0.0
358⍝│ 1.0986122886681098
359ˡb̲allDistance ← { p q → 2.0 × ᵗʳa̲tanh ʰn̲orm (n̲eg p) ʰm̲obius q }
360
361⍝│ (p, q) b̲allGeodesic n: n + 1 points along the geodesic (the
362⍝│ shortest path, by the ball's distance) from p to q, p first and q
363⍝│ last, evenly spaced along it: p m_obius the arrow to q scaled by
364⍝│ 0, 1/n, ..., 1 in the ball's way (tanh of that fraction of atanh
365⍝│ its length). Drawn flat, a geodesic is an arc of a circle that
366⍝│ meets the rim at right angles, or a diameter.
367⍝│ >> "mf:" u_se< "Manifold"
368⍝│ >> (0.0 0.0, 0.6 0.0) mf:b_allGeodesic 2
369⍝│ 0.0 0.0
370⍝│ 0.3333333333333333 0.0
371⍝│ 0.6 0.0
372ˡb̲allGeodesic ← { (p, q) n →
373 w ← (n̲eg p) ʰm̲obius q
374 r ← ʰn̲orm w
375 rows ← (n + 1) r̲eshape 1.0
376 r = 0.0 ? rows '× t̲able f̲loat p
377 S ← (t̲anh (ᵗʳa̲tanh r) × ((f̲loat r̲ange n + 1) − 1.0) ÷ f̲loat n) ÷ r
378 d ← t̲ally p
379 (M, _) ← (n + 1) '{ (M, k) → (M c̲at (1 c̲at d) r̲eshape p ʰm̲obius (k s̲elect S) × w, k + 1) } p̲ower ((0 c̲at d) r̲eshape 0.0, 1)
380 M
381}
382
383⍝━━ Distortion
384
385⍝│ D d̲istortion E: how badly the distances E (between every pair of
386⍝│ points of an embedding, a matrix as s_phereDistances gives) keep
387⍝│ the distances D they should have (the same shape: the data's own
388⍝│ distances, say along a tree): for every pair of different points,
389⍝│ |E / D - 1|, the fraction by which the pair's distance is too long
390⍝│ or too short; the mean over the pairs and the worst pair, as a
391⍝│ pair of Floats. 0 0 when every distance is kept; 1 for a pair twice
392⍝│ or half as far as it should be (half: 0.5). Pairs at D = 0 (the
393⍝│ same point twice) are left out; the diagonal always is. The measure
394⍝│ of H1 (docs/research-questions.md), the same whatever space E was
395⍝│ measured along.
396⍝│ >> "mf:" u_se< "Manifold"
397⍝│ >> D := 3 3 r_eshape 0.0 1.0 2.0 1.0 0.0 1.0 2.0 1.0 0.0
398⍝│ >> D mf:d_istortion D
399⍝│ 0.0 0.0
400⍝│ >> D mf:d_istortion 2.0 * D
401⍝│ 1.0 1.0
402⍝│ >> D mf:d_istortion 3 3 r_eshape 0.0 1.0 1.0 1.0 0.0 1.0 1.0 1.0 0.0
403⍝│ 0.16666666666666666 0.5
404ˡd̲istortion ← { D E →
405 n ← f̲irst s̲hape D
406 pair ← ((r̲ange n) '< t̲able r̲ange n) ∧ (f̲loat D) > 0.0 ⍝ each pair once, where there is a distance to keep
407 k ← '+ r̲/ r̲avel pair
408 k = 0 ? 0.0 0.0
409 off ← r̲avel (f̲loat pair) × a̲bs ((f̲loat E) ÷ (f̲loat D) m̲ax 1e-300) − 1.0
410 (('+ r̲/ off) ÷ f̲loat k) c̲at 'm̲ax r̲/ off
411}