let epsilon = 1e-9

let sin_taylor x =
  let rec accumulate s y k =
    if abs_float y < epsilon then s else
    accumulate (s +. y) (-. y *. x *. x /. (k +. 1.) /. (k +. 2.)) (k +. 2.)
  in accumulate 0. x 1.

let sin_cube x =
  let rec reduce k y =
    if abs_float y < epsilon then restore k y else
    reduce (k + 1) (y /. 3.)
  and restore k y =
    if k = 0 then y else
    restore (k - 1) (y *. (3. -. 4. *. y *. y))
  in reduce 0 x

let half_pi = 1.570796326794896619231321691638

let reduce_quadrant a =
  let q = floor (a /. half_pi) in
  (truncate q mod 4 + 4) mod 4, a -. half_pi *. q

let sine loop a =
  let (i, x) = reduce_quadrant a in
  let y = loop (if i mod 2 = 0 then x else half_pi -. x) in
  if i < 2 then y else -. y