(* OCaml *)
let tape_backward last =
  let tp = !tape in
  let vals = tp.vals in
  let grads = tp.grads in
  for cur = last downto 0 do
    let g = Float.Array.get grads cur in
    if g <> 0.0 then begin
      let op = ops.(cur) in
      if op = op_mul then begin
        Float.Array.set grads p1
          (grads p1 +. g *. vals p2);
        Float.Array.set grads p2
          (grads p2 +. g *. vals p1)
      end
    end
  done

let softmax logits =
  let m = Array.fold_left max
    neg_infinity logits in
  let exps = Array.map
    (fun x -> exp (x -. m)) logits in
  let s = Array.fold_left
    (+.) 0.0 exps in
  Array.map (fun x -> x /. s) exps

-- Haskell
randomGaussian :: Double -> Double -> IO Double
randomGaussian μ σ = do
  u1 <- randomIO :: IO Double
  u2 <- randomIO :: IO Double
  let z0 = sqrt (-(2.0 * log u1))
           * cos (2.0 * pi * u2)
  pure $ μ + σ * z0

applyGradients :: Model -> LearningParameters
  -> Int -> UV.Vector Double
  -> Training -> IO ()
applyGradients Model{weightValues}
    params step grads
    Training{mAdam, vAdam} = do
  let n = UMV.length weightValues
  forM_ [0 .. n - 1] $ \i -> do
    let m_i' = beta1 * m_i
              + (1 - beta1) * w_grad
        v_i' = beta2 * v_i
              + (1 - beta2) * w_grad * w_grad
    UMV.write weightValues i
      $ w - lr * m_hat / (sqrt v_hat + eps)

; Scheme
(define (tape:eval-linear last)
  (let ((tp (*tape*)))
    (let loop ((i 0))
      (when (fx<= i last)
        (let ((op-tag (vector-ref
                (tape-op tp) i))
              (vals (tape-val tp)))
          (cond
            ((fx= op-tag op-add)
             (f64vector-set! vals i
               (+ (f64vector-ref vals
                    (vector-ref
                      (tape-p1 tp) i))
                  (f64vector-ref vals
                    (vector-ref
                      (tape-p2 tp) i)))))
            ((fx= op-tag op-mul)
             (f64vector-set! vals i
               (* (f64vector-ref vals
                    (vector-ref
                      (tape-p1 tp) i))
                  (f64vector-ref vals
                    (vector-ref
                      (tape-p2 tp) i)))))))
        (loop (fx+ i 1))))))

\ Forth
: gaussian ( mean stddev -- val )
  1e frnd f- frnd
  flog 2e fnegate f* fsqrt fswap
  2e pi f* f* fcos f* f* f+ ;

: matvec-mul
  ( mat vec rows cols -- res )
  { mat vec rows cols }
  rows cells allocate throw
  { result }
  rows 0 do i { row }
    mat row cols * cells + @
    vec @ *.
    cols 1 do
      mat row cols * i +
        cells + @
      vec i cells + @ *. +.
    loop
    result row cells + !
  loop result ;

: layer-forward
  ( input w b n-in n-out )
  { inp w b n-in n-out }
  n-out cells allocate throw
  { out }
  n-out 0 do
    b i cells + @ { acc }
    n-in 0 do
      inp j cells + @
      w i n-in * j + cells + @
      *. acc +. to acc
    loop
    acc out i cells + !
  loop out ;