blog

Elliptic Fourier series are a good example to kick the tires on array programming systems; J and Python, however, are both dynamically typed.

A modern statically approach makes things significantly easier on the programmer; here I will use my Apple language, which aspires to the interactivity and conciseness of the APL family.

One has a principled yet safe version of Haskell's zip or J's " (re-rank on a dyadic verb):

 > :ty \f.\x.\y. f`x y
((a → (b → c)) → (Arr (i `Cons` Nil) a → (Arr (i `Cons` Nil) b → Arr (i `Cons` Nil) c)))

This avoids the pitfalls of zip-ing different-length lists (one cannot write eqList = zipWith (==) in Haskell) and avoids deferring failures to runtime as in J.

One can inspect the type signature of cons:

 > :ty (<|)
(a → (Arr (i `Cons` Nil) a → Arr ((i + 1) `Cons` Nil) a))

There is a safe last:

 > :ty }.
(Arr ((i + 1) `Cons` Nil) a → a)

These work together:

 > :ty \x.\xs. }.(x<|xs)
(b → (Arr (i `Cons` Nil) b → b))

I believe this combination of interactivity and type inference in array programming will be quite productive.

To address the elliptic Fourier series, we focus on the offsets $A_0$ and $C_0$. This is straightforward to express with mathematical notation:

$\displaystyle A_0 = \frac{1}{T} \sum_{p=1}^{K} \left(\frac{\Delta x_p}{2\Delta t_p}(t_p^2-t_{p-1}^2) + \xi_p(t_p-t_{p-1}) \right) \\ \displaystyle C_0 = \frac{1}{T} \sum_{p=1}^{K} \left(\frac{\Delta y_p}{2\Delta t_p}(t_p^2-t_{p-1}^2) + \delta_p(t_p-t_{p-1}) \right)$

where

$\xi_p = \displaystyle\sum_{j=1}^{p-1}\Delta x_j - \frac{\Delta x_p}{\Delta t_p}\sum_{j=1}^{p-1}\Delta t_j \\ \delta_p = \displaystyle\sum_{j=1}^{p-1}\Delta y_j - \frac{\Delta y_p}{\Delta t_p}\sum_{j=1}^{p-1}\Delta t_j$

In Apple:

λxs.λys.
  { sum ⇐ [(+)/1 0 x]
  ; dxs ⟜ (-)\~xs; dys ⟜ (-)\~ys
  ; dts ⟜ [√(x^2+y^2)]`dxs dys
  ; pxs ← (+)Λ 0 dxs; pys ← (+)Λ 0 dys; pts ⟜ (+)Λ 0 dts
  ; dtss ⟜ 0<|(-)\~((^2)'1 pts)
  ; T ⟜}.pts
  ; 𝜉 ← (-)`pxs ((*)`((%)`dxs dts) pts)
  ; 𝛿 ← (-)`pys ((*)`((%)`dys dts) pts)
  ; A ← ((sum ((*)`((%)`dxs dts) dtss))%2 + (sum ((*)`𝜉 dts)))%T
  ; C ← ((sum ((*)`((%)`dys dts) dtss))%2 + (sum ((*)`𝛿 dts)))%T
  ; (A,C)
  }

This all typechecks! So we can be sure that there will not be off-by-one errors.