--
-- Leibniz formula (Fourier series coefficients)
--

def f (x: MathExpr) : MathExpr := x

def multSd (x: MathExpr) (f: MathExpr) (G: MathExpr) : MathExpr :=
  let F := Sd x f
   in F * G - Sd x (f * d/d G x)

-- Test multSd with various trigonometric functions
assertEqual "multSd cos x"
  (multSd x (cos x) (f x))
  (x * sin x + cos x - sin x)

assertEqual "multSd cos 2x"
  (multSd x (cos (2 * x)) (f x))
  (x * sin (2 * x) / 2 + cos (2 * x) / 4 - sin (2 * x) / 2)

assertEqual "multSd sin x"
  (multSd x (sin x) (f x))
  (- x * cos x + sin x + cos x)

-- Fourier coefficients for f(x) = x
def coeffAs : [MathExpr] :=
  map
    (\n ->
      let F := multSd x (cos (n * x)) (f x)
       in (substitute [(x, π)] F - substitute [(x, - π)] F) / π)
    nats

-- First 10 coefficients should all be 0 (f(x) = x is an odd function)
assertEqual "first 10 Fourier cosine coefficients"
  (take 10 coeffAs)
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

def bs : [MathExpr] :=
  map
    (\n ->
      let F := multSd x (sin (n * x)) (f x)
       in (substitute [(x, π)] F - substitute [(x, - π)] F) / π)
    (take 10 nats)

assertEqual "first 10 Fourier sine coefficients"
  (take 10 bs)
  [2, -1, 2/3, -1/2, 2/5, -1/3, 2/7, -1/4, 2/9, -1/5]

def f' : [MathExpr] := map (\(k, b) -> b * sin (k * x)) (zip nats bs)

-- Fourier series terms
assertEqual "first 10 Fourier series terms"
  (take 10 f')
  [2 * sin x, - sin (2 * x), 2 * sin (3 * x) / 3, - sin (4 * x) / 2,
   2 * sin (5 * x) / 5, - sin (6 * x) / 3, 2 * sin (7 * x) / 7, - sin (8 * x) / 4,
   2 * sin (9 * x) / 9, - sin (10 * x) / 5]