    {-# LANGUAGE DeriveFunctor #-}
    {-# LANGUAGE TypeFamilies  #-}
    -----------------------------------------------------------------------------
    -- |
    -- Module      :  Diagrams.TwoD.Segment.Bernstein
    -- Copyright   :  (c) 2014-2015 diagrams-lib team (see LICENSE)
    -- License     :  BSD-style (see LICENSE)
    -- Maintainer  :  diagrams-discuss@googlegroups.com
    --
    -- Bernstein polynomials, used internally by code to find
    -- intersections of paths.  This module is probably not of any
    -- relevance to most users of diagrams.
    -----------------------------------------------------------------------------
    module Diagrams.TwoD.Segment.Bernstein
      ( BernsteinPoly (..)
      , listToBernstein
      , evaluateBernstein

      , degreeElevate
      , bernsteinDeriv
      , evaluateBernsteinDerivs
      ) where

    import           Data.List           (tails)
    import           Diagrams.Core.V
    import           Diagrams.Parametric
    import           Linear.V1

    -- | Compute the binomial coefficients of degree n.
    binomials :: Num n => Int -> [n]
    binomials n = map fromIntegral $ scanl (\x m -> x * (n - m+1) `quot` m) 1 [1..n]

    data BernsteinPoly n = BernsteinPoly
      { bernsteinDegree :: Int
      , bernsteinCoeffs :: [n]
      } deriving (Show, Functor)

    type instance V        (BernsteinPoly n) = V1
    type instance N        (BernsteinPoly n) = n
    type instance Codomain (BernsteinPoly n) = V1

    -- | Create a bernstein polynomial from a list of coëfficients.
    listToBernstein :: Fractional n => [n] -> BernsteinPoly n
    listToBernstein [] = 0
    listToBernstein l  = BernsteinPoly (length l - 1) l

    -- | Degree elevate a bernstein polynomial a number of times.
    degreeElevate :: Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
    degreeElevate b                    0     = b
    degreeElevate (BernsteinPoly lp p) times =
      degreeElevate (BernsteinPoly (lp+1) (head p:inner p 1)) (times-1)
      where
        n = fromIntegral lp

        inner []         _ = [0]
        inner [a]        _ = [a]
        inner (a:b:rest) i = (i*a/(n+1) + b*(1 - i/(n+1))) : inner (b:rest) (i+1)

    -- | Evaluate the bernstein polynomial.
    evaluateBernstein :: Fractional n => BernsteinPoly n -> n -> n
    evaluateBernstein (BernsteinPoly _ [])       _ = 0
    evaluateBernstein (BernsteinPoly _ [b])      _ = b
    evaluateBernstein (BernsteinPoly lp (b':bs)) t = go t n (b'*u) 2 bs
      where
        u = 1-t
        n = fromIntegral lp

        go tn bc tmp _ [b]      = tmp + tn*bc*b
        go tn bc tmp i (b:rest) =
          go (tn*t)              -- tn
             (bc*(n - i+1)/i)    -- bc
             ((tmp + tn*bc*b)*u) -- tmp
             (i+1)               -- i
             rest
        go _ _ _ _ []           = error "evaluateBernstein: impossible"

    -- | Evaluate the bernstein polynomial and its derivatives.
    evaluateBernsteinDerivs :: Fractional n => BernsteinPoly n -> n -> [n]
    evaluateBernsteinDerivs b t
      | bernsteinDegree b == 0 = [evaluateBernstein b t]
      | otherwise              = evaluateBernstein b t : evaluateBernsteinDerivs (bernsteinDeriv b) t

    -- | Find the derivative of a bernstein polynomial.
    bernsteinDeriv :: Fractional n => BernsteinPoly n -> BernsteinPoly n
    bernsteinDeriv (BernsteinPoly 0 _)  = 0
    bernsteinDeriv (BernsteinPoly lp p) =
      -- BernsteinPoly (lp-1) $ map (* fromIntegral lp) $ zipWith (-) (tail p) p
      BernsteinPoly (lp-1) $ zipWith (\a b -> (a - b) * fromIntegral lp) (tail p) p

    instance Fractional n => Parametric (BernsteinPoly n) where
      atParam b = V1 . evaluateBernstein b
    instance Num n        => DomainBounds (BernsteinPoly n)
    instance Fractional n => EndValues    (BernsteinPoly n)
    instance Fractional n => Sectionable  (BernsteinPoly n) where
      splitAtParam  = bernsteinSplit
      reverseDomain (BernsteinPoly i xs) = BernsteinPoly i (reverse xs)

    -- | Split a bernstein polynomial.
    bernsteinSplit :: Num n => BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
    bernsteinSplit (BernsteinPoly lp p) t =
      (BernsteinPoly lp $ map head controls,
       BernsteinPoly lp $ reverse $ map last controls)
      where
        interp a b = (1-t)*a + t*b

        terp [_] = []
        terp l   = let ctrs = zipWith interp l (tail l)
                   in  ctrs : terp ctrs
        controls = p : terp p

    instance Fractional n => Num (BernsteinPoly n) where
      ba@(BernsteinPoly la a) + bb@(BernsteinPoly lb b)
        | la < lb   = BernsteinPoly lb $ zipWith (+) (bernsteinCoeffs $ degreeElevate ba $ lb - la) b
        | la > lb   = BernsteinPoly la $ zipWith (+) a (bernsteinCoeffs $ degreeElevate bb $ la - lb)
        | otherwise = BernsteinPoly la $ zipWith (+) a b

      ba@(BernsteinPoly la a) - bb@(BernsteinPoly lb b)
        | la < lb   = BernsteinPoly lb $ zipWith (-) (bernsteinCoeffs $ degreeElevate ba (lb - la)) b
        | la > lb   = BernsteinPoly la $ zipWith (-) a (bernsteinCoeffs $ degreeElevate bb (la - lb))
        | otherwise = BernsteinPoly la $ zipWith (-) a b

      (BernsteinPoly la a) * (BernsteinPoly lb b) =
        BernsteinPoly (la+lb) $
        zipWith (flip (/)) (binomials (la + lb)) $
                       init $ map sum $
                       map (zipWith (*) a') (down b') ++
                       map (zipWith (*) (reverse b')) (tail $ tails a')
                       -- zipWith (zipWith (*)) (tail $ tails a') (repeat $ reverse b')
        where down l = tail $ scanl (flip (:)) [] l -- [[1], [2, 1], [3, 2, 1], ...
              a' = zipWith (*) a (binomials la)
              b' = zipWith (*) b (binomials lb)

      fromInteger a = BernsteinPoly 0 [fromInteger a]

      signum (BernsteinPoly _ [])    = 0
      signum (BernsteinPoly _ (a:_)) = BernsteinPoly 0 [signum a]

      abs = fmap abs