{-# 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.NonEmpty (NonEmpty (..))
import qualified Data.List.NonEmpty.Compat as NE
import Diagrams.Core.V
import Diagrams.Parametric
import Linear.V1

-- | Compute the binomial coefficients of degree n.
binomials :: Num n => Int -> NonEmpty n
binomials :: forall n. Num n => Int -> NonEmpty n
binomials Int
n = (Int -> n) -> NonEmpty Int -> NonEmpty n
forall a b. (a -> b) -> NonEmpty a -> NonEmpty b
NE.map Int -> n
forall a b. (Integral a, Num b) => a -> b
fromIntegral (NonEmpty Int -> NonEmpty n) -> NonEmpty Int -> NonEmpty n
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> Int -> [Int] -> NonEmpty Int
forall (f :: * -> *) b a.
Foldable f =>
(b -> a -> b) -> b -> f a -> NonEmpty b
NE.scanl (\Int
x Int
m -> Int
x Int -> Int -> Int
forall a. Num a => a -> a -> a
* (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
m) Int
1 [Int
1 .. Int
n]

data BernsteinPoly n = BernsteinPoly
  { forall n. BernsteinPoly n -> Int
bernsteinDegree :: Int
  , forall n. BernsteinPoly n -> NonEmpty n
bernsteinCoeffs :: NonEmpty n
  }
  deriving (Int -> BernsteinPoly n -> ShowS
[BernsteinPoly n] -> ShowS
BernsteinPoly n -> String
(Int -> BernsteinPoly n -> ShowS)
-> (BernsteinPoly n -> String)
-> ([BernsteinPoly n] -> ShowS)
-> Show (BernsteinPoly n)
forall n. Show n => Int -> BernsteinPoly n -> ShowS
forall n. Show n => [BernsteinPoly n] -> ShowS
forall n. Show n => BernsteinPoly n -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall n. Show n => Int -> BernsteinPoly n -> ShowS
showsPrec :: Int -> BernsteinPoly n -> ShowS
$cshow :: forall n. Show n => BernsteinPoly n -> String
show :: BernsteinPoly n -> String
$cshowList :: forall n. Show n => [BernsteinPoly n] -> ShowS
showList :: [BernsteinPoly n] -> ShowS
Show, (forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b)
-> (forall a b. a -> BernsteinPoly b -> BernsteinPoly a)
-> Functor BernsteinPoly
forall a b. a -> BernsteinPoly b -> BernsteinPoly a
forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
forall (f :: * -> *).
(forall a b. (a -> b) -> f a -> f b)
-> (forall a b. a -> f b -> f a) -> Functor f
$cfmap :: forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
fmap :: forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
$c<$ :: forall a b. a -> BernsteinPoly b -> BernsteinPoly a
<$ :: forall a b. a -> BernsteinPoly b -> BernsteinPoly a
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 :: forall n. Fractional n => [n] -> BernsteinPoly n
listToBernstein [n]
l = case [n] -> Maybe (NonEmpty n)
forall a. [a] -> Maybe (NonEmpty a)
NE.nonEmpty [n]
l of
  Maybe (NonEmpty n)
Nothing -> BernsteinPoly n
0
  Just NonEmpty n
ne -> Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly ([n] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [n]
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) NonEmpty n
ne

-- | Degree elevate a bernstein polynomial a number of times.
degreeElevate :: Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate :: forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
b Int
0 = BernsteinPoly n
b
degreeElevate (BernsteinPoly Int
lp NonEmpty n
p) Int
times =
  BernsteinPoly n -> Int -> BernsteinPoly n
forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate (Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly (Int
lp Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) (NonEmpty n -> n
forall a. NonEmpty a -> a
NE.head NonEmpty n
p n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| NonEmpty n -> n -> [n]
inner NonEmpty n
p n
1)) (Int
times Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
 where
  n :: n
n = Int -> n
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
lp

  inner :: NonEmpty n -> n -> [n]
inner (n
a :| []) n
_ = [n
a]
  inner (n
a :| n
b : [n]
rest) n
i = (n
i n -> n -> n
forall a. Num a => a -> a -> a
* n
a n -> n -> n
forall a. Fractional a => a -> a -> a
/ (n
n n -> n -> n
forall a. Num a => a -> a -> a
+ n
1) n -> n -> n
forall a. Num a => a -> a -> a
+ n
b n -> n -> n
forall a. Num a => a -> a -> a
* (n
1 n -> n -> n
forall a. Num a => a -> a -> a
- n
i n -> n -> n
forall a. Fractional a => a -> a -> a
/ (n
n n -> n -> n
forall a. Num a => a -> a -> a
+ n
1))) n -> [n] -> [n]
forall a. a -> [a] -> [a]
: NonEmpty n -> n -> [n]
inner (n
b n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| [n]
rest) (n
i n -> n -> n
forall a. Num a => a -> a -> a
+ n
1)

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

  go :: n -> n -> n -> n -> [n] -> n
go n
tn n
bc n
tmp n
_ [n
b] = n
tmp n -> n -> n
forall a. Num a => a -> a -> a
+ n
tn n -> n -> n
forall a. Num a => a -> a -> a
* n
bc n -> n -> n
forall a. Num a => a -> a -> a
* n
b
  go n
tn n
bc n
tmp n
i (n
b : [n]
rest) =
    n -> n -> n -> n -> [n] -> n
go
      (n
tn n -> n -> n
forall a. Num a => a -> a -> a
* n
t) -- tn
      (n
bc n -> n -> n
forall a. Num a => a -> a -> a
* (n
n n -> n -> n
forall a. Num a => a -> a -> a
- n
i n -> n -> n
forall a. Num a => a -> a -> a
+ n
1) n -> n -> n
forall a. Fractional a => a -> a -> a
/ n
i) -- bc
      ((n
tmp n -> n -> n
forall a. Num a => a -> a -> a
+ n
tn n -> n -> n
forall a. Num a => a -> a -> a
* n
bc n -> n -> n
forall a. Num a => a -> a -> a
* n
b) n -> n -> n
forall a. Num a => a -> a -> a
* n
u) -- tmp
      (n
i n -> n -> n
forall a. Num a => a -> a -> a
+ n
1) -- i
      [n]
rest
  go n
_ n
_ n
_ n
_ [] = String -> n
forall a. HasCallStack => String -> a
error String
"evaluateBernstein: impossible"

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

-- | Find the derivative of a bernstein polynomial.
bernsteinDeriv :: Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv :: forall n. Fractional n => BernsteinPoly n -> BernsteinPoly n
bernsteinDeriv (BernsteinPoly Int
_ (n
_ :| [])) = BernsteinPoly n
0
bernsteinDeriv (BernsteinPoly Int
lp p :: NonEmpty n
p@(n
_ :| (n
a1 : [n]
as))) =
  -- BernsteinPoly (lp-1) $ map (* fromIntegral lp) $ zipWith (-) (drop 1 p) p
  Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly (Int
lp Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith (\n
a n
b -> (n
a n -> n -> n
forall a. Num a => a -> a -> a
- n
b) n -> n -> n
forall a. Num a => a -> a -> a
* Int -> n
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
lp) (n
a1 n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| [n]
as) NonEmpty n
p

instance Fractional n => Parametric (BernsteinPoly n) where
  atParam :: BernsteinPoly n
-> N (BernsteinPoly n)
-> Codomain (BernsteinPoly n) (N (BernsteinPoly n))
atParam BernsteinPoly n
b = n -> V1 n
forall a. a -> V1 a
V1 (n -> V1 n) -> (n -> n) -> n -> V1 n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. BernsteinPoly n -> n -> n
forall n. Fractional n => BernsteinPoly n -> n -> n
evaluateBernstein BernsteinPoly n
b
instance Num n => DomainBounds (BernsteinPoly n)
instance Fractional n => EndValues (BernsteinPoly n)
instance Fractional n => Sectionable (BernsteinPoly n) where
  splitAtParam :: BernsteinPoly n
-> N (BernsteinPoly n) -> (BernsteinPoly n, BernsteinPoly n)
splitAtParam = BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
BernsteinPoly n
-> N (BernsteinPoly n) -> (BernsteinPoly n, BernsteinPoly n)
forall n.
Num n =>
BernsteinPoly n -> n -> (BernsteinPoly n, BernsteinPoly n)
bernsteinSplit
  reverseDomain :: BernsteinPoly n -> BernsteinPoly n
reverseDomain (BernsteinPoly Int
i NonEmpty n
xs) = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
i (NonEmpty n -> NonEmpty n
forall a. NonEmpty a -> NonEmpty a
NE.reverse NonEmpty n
xs)

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

  -- terp :: NonEmpty n -> [NonEmpty n]
  terp :: NonEmpty n -> [NonEmpty n]
terp l :: NonEmpty n
l@(n
_ :| [n]
as) = case [n] -> Maybe (NonEmpty n)
forall a. [a] -> Maybe (NonEmpty a)
NE.nonEmpty [n]
as of
    Maybe (NonEmpty n)
Nothing -> []
    Just NonEmpty n
as' ->
      let ctrs :: NonEmpty n
ctrs = (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
interp NonEmpty n
l NonEmpty n
as'
       in NonEmpty n
ctrs NonEmpty n -> [NonEmpty n] -> [NonEmpty n]
forall a. a -> [a] -> [a]
: NonEmpty n -> [NonEmpty n]
terp NonEmpty n
ctrs
  controls :: NonEmpty (NonEmpty n)
controls = NonEmpty n
p NonEmpty n -> [NonEmpty n] -> NonEmpty (NonEmpty n)
forall a. a -> [a] -> NonEmpty a
:| NonEmpty n -> [NonEmpty n]
terp NonEmpty n
p

instance Fractional n => Num (BernsteinPoly n) where
  ba :: BernsteinPoly n
ba@(BernsteinPoly Int
la NonEmpty n
a) + :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
+ bb :: BernsteinPoly n
bb@(BernsteinPoly Int
lb NonEmpty n
b)
    | Int
la Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
lb = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
lb (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(+) (BernsteinPoly n -> NonEmpty n
forall n. BernsteinPoly n -> NonEmpty n
bernsteinCoeffs (BernsteinPoly n -> NonEmpty n) -> BernsteinPoly n -> NonEmpty n
forall a b. (a -> b) -> a -> b
$ BernsteinPoly n -> Int -> BernsteinPoly n
forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
ba (Int -> BernsteinPoly n) -> Int -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ Int
lb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
la) NonEmpty n
b
    | Int
la Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
lb = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
la (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(+) NonEmpty n
a (BernsteinPoly n -> NonEmpty n
forall n. BernsteinPoly n -> NonEmpty n
bernsteinCoeffs (BernsteinPoly n -> NonEmpty n) -> BernsteinPoly n -> NonEmpty n
forall a b. (a -> b) -> a -> b
$ BernsteinPoly n -> Int -> BernsteinPoly n
forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
bb (Int -> BernsteinPoly n) -> Int -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ Int
la Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
lb)
    | Bool
otherwise = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
la (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(+) NonEmpty n
a NonEmpty n
b

  ba :: BernsteinPoly n
ba@(BernsteinPoly Int
la NonEmpty n
a) - :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
- bb :: BernsteinPoly n
bb@(BernsteinPoly Int
lb NonEmpty n
b)
    | Int
la Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
lb = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
lb (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith (-) (BernsteinPoly n -> NonEmpty n
forall n. BernsteinPoly n -> NonEmpty n
bernsteinCoeffs (BernsteinPoly n -> NonEmpty n) -> BernsteinPoly n -> NonEmpty n
forall a b. (a -> b) -> a -> b
$ BernsteinPoly n -> Int -> BernsteinPoly n
forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
ba (Int
lb Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
la)) NonEmpty n
b
    | Int
la Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
lb = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
la (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith (-) NonEmpty n
a (BernsteinPoly n -> NonEmpty n
forall n. BernsteinPoly n -> NonEmpty n
bernsteinCoeffs (BernsteinPoly n -> NonEmpty n) -> BernsteinPoly n -> NonEmpty n
forall a b. (a -> b) -> a -> b
$ BernsteinPoly n -> Int -> BernsteinPoly n
forall n. Fractional n => BernsteinPoly n -> Int -> BernsteinPoly n
degreeElevate BernsteinPoly n
bb (Int
la Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
lb))
    | Bool
otherwise = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
la (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$ (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith (-) NonEmpty n
a NonEmpty n
b

  (BernsteinPoly Int
la NonEmpty n
a) * :: BernsteinPoly n -> BernsteinPoly n -> BernsteinPoly n
* (BernsteinPoly Int
lb NonEmpty n
b) =
    Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly (Int
la Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
lb) (NonEmpty n -> BernsteinPoly n) -> NonEmpty n -> BernsteinPoly n
forall a b. (a -> b) -> a -> b
$
      (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith ((n -> n -> n) -> n -> n -> n
forall a b c. (a -> b -> c) -> b -> a -> c
flip n -> n -> n
forall a. Fractional a => a -> a -> a
(/)) (Int -> NonEmpty n
forall n. Num n => Int -> NonEmpty n
binomials (Int
la Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
lb)) (NonEmpty n -> NonEmpty n) -> NonEmpty n -> NonEmpty n
forall a b. (a -> b) -> a -> b
$
        (NonEmpty n -> n) -> NonEmpty (NonEmpty n) -> NonEmpty n
forall a b. (a -> b) -> NonEmpty a -> NonEmpty b
NE.map NonEmpty n -> n
forall a. Num a => NonEmpty a -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (NonEmpty (NonEmpty n) -> NonEmpty n)
-> NonEmpty (NonEmpty n) -> NonEmpty n
forall a b. (a -> b) -> a -> b
$
          NonEmpty (NonEmpty n) -> [NonEmpty n] -> NonEmpty (NonEmpty n)
forall a. NonEmpty a -> [a] -> NonEmpty a
NE.appendList
            ((NonEmpty n -> NonEmpty n)
-> NonEmpty (NonEmpty n) -> NonEmpty (NonEmpty n)
forall a b. (a -> b) -> NonEmpty a -> NonEmpty b
NE.map ((n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(*) NonEmpty n
a') (NonEmpty n -> NonEmpty (NonEmpty n)
forall {a}. NonEmpty a -> NonEmpty (NonEmpty a)
down NonEmpty n
b'))
            ((NonEmpty n -> NonEmpty n) -> [NonEmpty n] -> [NonEmpty n]
forall a b. (a -> b) -> [a] -> [b]
map ((n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(*) (NonEmpty n -> NonEmpty n
forall a. NonEmpty a -> NonEmpty a
NE.reverse NonEmpty n
b')) (NonEmpty (NonEmpty n) -> [NonEmpty n]
forall a. NonEmpty a -> [a]
NE.tail (NonEmpty (NonEmpty n) -> [NonEmpty n])
-> NonEmpty (NonEmpty n) -> [NonEmpty n]
forall a b. (a -> b) -> a -> b
$ NonEmpty n -> NonEmpty (NonEmpty n)
forall {a}. NonEmpty a -> NonEmpty (NonEmpty a)
NE.tails1 NonEmpty n
a'))
   where
    down :: NonEmpty a -> NonEmpty (NonEmpty a)
down (a
x :| [a]
xs) = (NonEmpty a -> a -> NonEmpty a)
-> NonEmpty a -> [a] -> NonEmpty (NonEmpty a)
forall (f :: * -> *) b a.
Foldable f =>
(b -> a -> b) -> b -> f a -> NonEmpty b
NE.scanl ((a -> NonEmpty a -> NonEmpty a) -> NonEmpty a -> a -> NonEmpty a
forall a b c. (a -> b -> c) -> b -> a -> c
flip a -> NonEmpty a -> NonEmpty a
forall a. a -> NonEmpty a -> NonEmpty a
NE.cons) (a
x a -> [a] -> NonEmpty a
forall a. a -> [a] -> NonEmpty a
:| []) [a]
xs -- [[1], [2, 1], [3, 2, 1], ...
    a' :: NonEmpty n
a' = (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(*) NonEmpty n
a (Int -> NonEmpty n
forall n. Num n => Int -> NonEmpty n
binomials Int
la)
    b' :: NonEmpty n
b' = (n -> n -> n) -> NonEmpty n -> NonEmpty n -> NonEmpty n
forall a b c.
(a -> b -> c) -> NonEmpty a -> NonEmpty b -> NonEmpty c
NE.zipWith n -> n -> n
forall a. Num a => a -> a -> a
(*) NonEmpty n
b (Int -> NonEmpty n
forall n. Num n => Int -> NonEmpty n
binomials Int
lb)

  fromInteger :: Integer -> BernsteinPoly n
fromInteger Integer
a = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
0 (Integer -> n
forall a. Num a => Integer -> a
fromInteger Integer
a n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| [])

  signum :: BernsteinPoly n -> BernsteinPoly n
signum (BernsteinPoly Int
_ (n
a :| [n]
_)) = Int -> NonEmpty n -> BernsteinPoly n
forall n. Int -> NonEmpty n -> BernsteinPoly n
BernsteinPoly Int
0 (n -> n
forall a. Num a => a -> a
signum n
a n -> [n] -> NonEmpty n
forall a. a -> [a] -> NonEmpty a
:| [])

  abs :: BernsteinPoly n -> BernsteinPoly n
abs = (n -> n) -> BernsteinPoly n -> BernsteinPoly n
forall a b. (a -> b) -> BernsteinPoly a -> BernsteinPoly b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap n -> n
forall a. Num a => a -> a
abs