{-# LANGUAGE BangPatterns, ViewPatterns #-}
-- | Algebra on polynomials in the Bernstein form.  It is based on the
-- paper /Algebraic manipulation in the Bernstein form made simple via
-- convolutions/ by J. Sanchez-Reyes.  It's an efficient
-- implementation using the scaled basis, and using ghc rewrite rules
-- to eliminate intermediate polynomials.

module Math.BernsteinPoly
       (BernsteinPoly(..), bernsteinSubsegment, listToBernstein, zeroPoly,
        (~*), (*~), (~+), (~-), degreeElevate, bernsteinSplit, bernsteinEval,
        bernsteinEvalDeriv, binCoeff, convolve, bernsteinEvalDerivs, bernsteinDeriv)
       where
import Data.Vector.Unboxed as V
import Data.Vector.Unboxed.Mutable as M
import qualified Data.Vector as B

data BernsteinPoly a = BernsteinPoly {
  forall a. BernsteinPoly a -> Vector a
bernsteinCoeffs :: V.Vector a}
                   deriving Int -> BernsteinPoly a -> ShowS
forall a. (Show a, Unbox a) => Int -> BernsteinPoly a -> ShowS
forall a. (Show a, Unbox a) => [BernsteinPoly a] -> ShowS
forall a. (Show a, Unbox a) => BernsteinPoly a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [BernsteinPoly a] -> ShowS
$cshowList :: forall a. (Show a, Unbox a) => [BernsteinPoly a] -> ShowS
show :: BernsteinPoly a -> String
$cshow :: forall a. (Show a, Unbox a) => BernsteinPoly a -> String
showsPrec :: Int -> BernsteinPoly a -> ShowS
$cshowsPrec :: forall a. (Show a, Unbox a) => Int -> BernsteinPoly a -> ShowS
Show

data ScaledPoly a = ScaledPoly {
  forall a. ScaledPoly a -> Vector a
scaledCoeffs :: V.Vector a }
                deriving Int -> ScaledPoly a -> ShowS
forall a. (Show a, Unbox a) => Int -> ScaledPoly a -> ShowS
forall a. (Show a, Unbox a) => [ScaledPoly a] -> ShowS
forall a. (Show a, Unbox a) => ScaledPoly a -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [ScaledPoly a] -> ShowS
$cshowList :: forall a. (Show a, Unbox a) => [ScaledPoly a] -> ShowS
show :: ScaledPoly a -> String
$cshow :: forall a. (Show a, Unbox a) => ScaledPoly a -> String
showsPrec :: Int -> ScaledPoly a -> ShowS
$cshowsPrec :: forall a. (Show a, Unbox a) => Int -> ScaledPoly a -> ShowS
Show
infixl 7 ~*, *~
infixl 6 ~+, ~-

{-# RULES "toScaled/fromScaled" forall a. toScaled (fromScaled a) = a;
  "fromScaled/toScaled" forall a. fromScaled (toScaled a) = a; #-}

toScaled :: (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled :: forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled (BernsteinPoly Vector a
v) =
  forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$
  forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(*) Vector a
v forall a b. (a -> b) -> a -> b
$ forall a. (Num a, Unbox a) => Int -> Vector a
binCoeff forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Num a => a -> a -> a
- Int
1
{-# NOINLINE [2] toScaled #-}

fromScaled :: (Unbox a, Fractional a) => ScaledPoly a -> BernsteinPoly a
fromScaled :: forall a.
(Unbox a, Fractional a) =>
ScaledPoly a -> BernsteinPoly a
fromScaled (ScaledPoly Vector a
v) =
    forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$
    forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Fractional a => a -> a -> a
(/) Vector a
v forall a b. (a -> b) -> a -> b
$ forall a. (Num a, Unbox a) => Int -> Vector a
binCoeff forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Num a => a -> a -> a
- Int
1
{-# NOINLINE [2] fromScaled #-}

-- | Create a bernstein polynomail from a list of coëfficients.
listToBernstein :: (Unbox a, Num a) => [a] -> BernsteinPoly a
listToBernstein :: forall a. (Unbox a, Num a) => [a] -> BernsteinPoly a
listToBernstein [] = forall a. (Num a, Unbox a) => BernsteinPoly a
zeroPoly
listToBernstein [a]
l = forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => [a] -> Vector a
V.fromList [a]
l
{-# INLINE listToBernstein #-}

-- | The constant zero.
zeroPoly :: (Num a, Unbox a) => BernsteinPoly a
zeroPoly :: forall a. (Num a, Unbox a) => BernsteinPoly a
zeroPoly = forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => [a] -> Vector a
V.fromList [a
0]
{-# SPECIALIZE zeroPoly :: BernsteinPoly Double #-}

-- | Return the subsegment between the two parameters.
bernsteinSubsegment :: (Unbox a, Ord a, Fractional a) =>
                       BernsteinPoly a -> a -> a -> BernsteinPoly a
bernsteinSubsegment :: forall a.
(Unbox a, Ord a, Fractional a) =>
BernsteinPoly a -> a -> a -> BernsteinPoly a
bernsteinSubsegment BernsteinPoly a
b a
t1 a
t2 
  | a
t1 forall a. Ord a => a -> a -> Bool
> a
t2   = forall a.
(Unbox a, Ord a, Fractional a) =>
BernsteinPoly a -> a -> a -> BernsteinPoly a
bernsteinSubsegment BernsteinPoly a
b a
t2 a
t1
  | Bool
otherwise = forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a.
(Unbox a, Num a) =>
BernsteinPoly a -> a -> (BernsteinPoly a, BernsteinPoly a)
bernsteinSplit (a
t1forall a. Fractional a => a -> a -> a
/a
t2) forall a b. (a -> b) -> a -> b
$
                forall a b. (a, b) -> a
fst forall a b. (a -> b) -> a -> b
$ forall a.
(Unbox a, Num a) =>
BernsteinPoly a -> a -> (BernsteinPoly a, BernsteinPoly a)
bernsteinSplit BernsteinPoly a
b a
t2
{-# INLINE bernsteinSubsegment #-}                

-- | Calculate the convolution of two vectors.
convolve :: (Unbox a, Num a) => Vector a -> Vector a -> Vector a
convolve :: forall a. (Unbox a, Num a) => Vector a -> Vector a -> Vector a
convolve Vector a
x Vector a
h = forall a. Unbox a => (forall s. ST s (MVector s a)) -> Vector a
V.create forall a b. (a -> b) -> a -> b
$ do
  let xN :: Int
xN = forall a. Unbox a => Vector a -> Int
V.length Vector a
x
  let hN :: Int
hN = forall a. Unbox a => Vector a -> Int
V.length Vector a
h
  let xIndices :: Vector Int
xIndices = forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Int
0 Int
xN
  let hIndices :: Vector Int
hIndices = forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Int
0 Int
hN

  MVector s a
xM <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
V.unsafeThaw Vector a
x
  MVector s a
hM <- forall a (m :: * -> *).
(Unbox a, PrimMonad m) =>
Vector a -> m (MVector (PrimState m) a)
V.unsafeThaw Vector a
h
  MVector s a
yM <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> a -> m (MVector (PrimState m) a)
M.replicate (Int
xN forall a. Num a => a -> a -> a
+ Int
hN forall a. Num a => a -> a -> a
- Int
1) a
0

  forall (m :: * -> *) a b.
(Monad m, Unbox a) =>
Vector a -> (a -> m b) -> m ()
V.forM_ Vector Int
xIndices forall a b. (a -> b) -> a -> b
$ \Int
i -> do
    a
a <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
M.unsafeRead MVector s a
xM Int
i
    forall (m :: * -> *) a b.
(Monad m, Unbox a) =>
Vector a -> (a -> m b) -> m ()
V.forM_ Vector Int
hIndices forall a b. (a -> b) -> a -> b
$ \Int
j -> do
      a
b <- forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> m a
M.unsafeRead MVector s a
hM Int
j
      forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> (a -> a) -> Int -> m ()
M.unsafeModify MVector s a
yM (forall a. Num a => a -> a -> a
+ a
a forall a. Num a => a -> a -> a
* a
b) (Int
i forall a. Num a => a -> a -> a
+ Int
j)

  forall (m :: * -> *) a. Monad m => a -> m a
return MVector s a
yM
{-# SPECIALIZE convolve :: Vector Double -> Vector Double -> Vector Double #-}

-- | Multiply two bernstein polynomials using convolution.  The final
-- degree will be the sum of either degrees.  This operation takes
-- O((n+m)^2) with n and m the degree of the beziers.  Note that
-- convolution can be done in O(n log n) using the FFT, which may be
-- faster for large polynomials.

(~*) :: (Unbox a, Fractional a) =>
        BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a
(forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled -> ScaledPoly a
a) ~* :: forall a.
(Unbox a, Fractional a) =>
BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a
~* (forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled -> ScaledPoly a
b) =
  forall a.
(Unbox a, Fractional a) =>
ScaledPoly a -> BernsteinPoly a
fromScaled forall a b. (a -> b) -> a -> b
$ forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
mulScaled ScaledPoly a
a ScaledPoly a
b
{-# INLINE (~*) #-}  

mulScaled :: (Unbox a, Num a) => ScaledPoly a -> ScaledPoly a -> ScaledPoly a
mulScaled :: forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
mulScaled (ScaledPoly Vector a
a) (ScaledPoly Vector a
b) =
  forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => Vector a -> Vector a -> Vector a
convolve Vector a
a Vector a
b
{-# INLINE mulScaled #-}    

-- | Give the binomial coefficients of degree n.
binCoeff :: (Num a, Unbox a) => Int -> V.Vector a
binCoeff :: forall a. (Num a, Unbox a) => Int -> Vector a
binCoeff Int
n = forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$
             forall a b.
(Unbox a, Unbox b) =>
(a -> b -> a) -> a -> Vector b -> Vector a
V.scanl (\Int
x Int
m -> Int
x forall a. Num a => a -> a -> a
* (Int
nforall a. Num a => a -> a -> a
-Int
mforall a. Num a => a -> a -> a
+Int
1) forall a. Integral a => a -> a -> a
`quot` Int
m)
             Int
1 (forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Int
1 Int
n)
{-# INLINE binCoeff #-}             

-- | Degree elevate a bernstein polynomial a number of times.
degreeElevateScaled :: (Unbox a, Num a)
                       => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled :: forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled b :: ScaledPoly a
b@(ScaledPoly Vector a
p) Int
times
  | Int
times forall a. Ord a => a -> a -> Bool
<= Int
0 = ScaledPoly a
b
  | Bool
otherwise = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => Vector a -> Vector a -> Vector a
convolve (forall a. (Num a, Unbox a) => Int -> Vector a
binCoeff Int
times) Vector a
p
{-# SPECIALIZE degreeElevateScaled :: ScaledPoly Double ->
    Int -> ScaledPoly Double #-}                

degreeElevate :: (Unbox a, Fractional a)
                 => BernsteinPoly a -> Int -> BernsteinPoly a
degreeElevate :: forall a.
(Unbox a, Fractional a) =>
BernsteinPoly a -> Int -> BernsteinPoly a
degreeElevate (forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled -> ScaledPoly a
b) Int
times =
  forall a.
(Unbox a, Fractional a) =>
ScaledPoly a -> BernsteinPoly a
fromScaled (forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled ScaledPoly a
b Int
times)
{-# INLINE degreeElevate #-}  

-- | Evaluate the bernstein polynomial using the horner rule adapted
-- for bernstein polynomials.
bernsteinEval :: (Unbox a, Fractional a)
                 => BernsteinPoly a -> a -> a
bernsteinEval :: forall a. (Unbox a, Fractional a) => BernsteinPoly a -> a -> a
bernsteinEval (BernsteinPoly Vector a
v) a
_
  | forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Eq a => a -> a -> Bool
== Int
0 = a
0
bernsteinEval (BernsteinPoly Vector a
v) a
_
  | forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Eq a => a -> a -> Bool
== Int
1 = forall a. Unbox a => Vector a -> a
V.unsafeHead Vector a
v
bernsteinEval (BernsteinPoly Vector a
v) a
t =
  a -> a -> a -> Int -> a
go a
t (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector a
v Int
0 forall a. Num a => a -> a -> a
* a
u) Int
1
  where u :: a
u = a
1forall a. Num a => a -> a -> a
-a
t
        n :: Int
n = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Num a => a -> a -> a
- Int
1
        go :: a -> a -> a -> Int -> a
go !a
tn !a
bc !a
tmp !Int
i
          | Int
i forall a. Eq a => a -> a -> Bool
== Int
n = a
tmp forall a. Num a => a -> a -> a
+ a
tnforall a. Num a => a -> a -> a
*forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector a
v Int
n
          | Bool
otherwise =
            a -> a -> a -> Int -> a
go (a
tnforall a. Num a => a -> a -> a
*a
t) -- tn
            (a
bcforall a. Num a => a -> a -> a
*forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
nforall a. Num a => a -> a -> a
-Int
i)forall a. Fractional a => a -> a -> a
/(forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i forall a. Num a => a -> a -> a
+ a
1)) -- bc
            ((a
tmp forall a. Num a => a -> a -> a
+ a
tnforall a. Num a => a -> a -> a
*a
bcforall a. Num a => a -> a -> a
*forall a. Unbox a => Vector a -> Int -> a
V.unsafeIndex Vector a
v Int
i)forall a. Num a => a -> a -> a
*a
u) -- tmp
            (Int
iforall a. Num a => a -> a -> a
+Int
1) -- i
{-# SPECIALIZE bernsteinEval :: BernsteinPoly Double -> Double -> Double #-}            
            
-- | Evaluate the bernstein polynomial and first derivative
bernsteinEvalDeriv :: (Unbox t, Fractional t) => BernsteinPoly t -> t -> (t,t)
bernsteinEvalDeriv :: forall t. (Unbox t, Fractional t) => BernsteinPoly t -> t -> (t, t)
bernsteinEvalDeriv b :: BernsteinPoly t
b@(BernsteinPoly Vector t
v) t
t
  | forall a. Unbox a => Vector a -> Int
V.length Vector t
v forall a. Ord a => a -> a -> Bool
<= Int
1 = (forall a. Unbox a => Vector a -> a
V.unsafeHead Vector t
v, t
0)
  | Bool
otherwise = (forall a. (Unbox a, Fractional a) => BernsteinPoly a -> a -> a
bernsteinEval BernsteinPoly t
b t
t, forall a. (Unbox a, Fractional a) => BernsteinPoly a -> a -> a
bernsteinEval (forall a. (Unbox a, Num a) => BernsteinPoly a -> BernsteinPoly a
bernsteinDeriv BernsteinPoly t
b) t
t)
{-# INLINE bernsteinEvalDeriv #-}                

            
-- | Evaluate the bernstein polynomial and its derivatives.
bernsteinEvalDerivs :: (Unbox t, Fractional t) => BernsteinPoly t -> t -> [t]
bernsteinEvalDerivs :: forall t. (Unbox t, Fractional t) => BernsteinPoly t -> t -> [t]
bernsteinEvalDerivs b :: BernsteinPoly t
b@(BernsteinPoly Vector t
v) t
t
  | forall a. Unbox a => Vector a -> Int
V.length Vector t
v forall a. Ord a => a -> a -> Bool
<= Int
1 = [forall a. Unbox a => Vector a -> a
V.unsafeHead Vector t
v, t
0]
  | Bool
otherwise = forall a. (Unbox a, Fractional a) => BernsteinPoly a -> a -> a
bernsteinEval BernsteinPoly t
b t
t forall a. a -> [a] -> [a]
:
                forall t. (Unbox t, Fractional t) => BernsteinPoly t -> t -> [t]
bernsteinEvalDerivs (forall a. (Unbox a, Num a) => BernsteinPoly a -> BernsteinPoly a
bernsteinDeriv BernsteinPoly t
b) t
t
{-# INLINE bernsteinEvalDerivs #-}                

-- | Find the derivative of a bernstein polynomial.
bernsteinDeriv :: (Unbox a, Num a) => BernsteinPoly a -> BernsteinPoly a
bernsteinDeriv :: forall a. (Unbox a, Num a) => BernsteinPoly a -> BernsteinPoly a
bernsteinDeriv (BernsteinPoly Vector a
v)
  | forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Eq a => a -> a -> Bool
== Int
0 = forall a. (Num a, Unbox a) => BernsteinPoly a
zeroPoly
bernsteinDeriv (BernsteinPoly Vector a
v) =
  forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$
  forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a. Num a => a -> a -> a
* forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a. Unbox a => Vector a -> Int
V.length Vector a
v forall a. Num a => a -> a -> a
- Int
1)) forall a b. (a -> b) -> a -> b
$
  forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) (forall a. Unbox a => Vector a -> Vector a
V.tail Vector a
v) Vector a
v
{-# SPECIALIZE bernsteinDeriv :: BernsteinPoly Double ->
    BernsteinPoly Double #-}  

-- | Split a bernstein polynomial
bernsteinSplit :: (Unbox a, Num a) =>
                  BernsteinPoly a -> a -> (BernsteinPoly a, BernsteinPoly a)
bernsteinSplit :: forall a.
(Unbox a, Num a) =>
BernsteinPoly a -> a -> (BernsteinPoly a, BernsteinPoly a)
bernsteinSplit (BernsteinPoly Vector a
v) a
t =
  (forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> Vector a -> Vector b
B.map forall a. Unbox a => Vector a -> a
V.head Vector (Vector a)
interpVecs,
   forall a. Vector a -> BernsteinPoly a
BernsteinPoly forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Vector a
V.reverse forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert forall a b. (a -> b) -> a -> b
$
   forall a b. (a -> b) -> Vector a -> Vector b
B.map forall a. Unbox a => Vector a -> a
V.last forall a b. (a -> b) -> a -> b
$ forall (v :: * -> *) a (w :: * -> *).
(Vector v a, Vector w a) =>
v a -> w a
convert Vector (Vector a)
interpVecs)
  where
    interp :: a -> a -> a
interp a
a a
b = (a
1forall a. Num a => a -> a -> a
-a
t)forall a. Num a => a -> a -> a
*a
a forall a. Num a => a -> a -> a
+ a
tforall a. Num a => a -> a -> a
*a
b
    interpVecs :: Vector (Vector a)
interpVecs = forall a. Int -> (a -> a) -> a -> Vector a
B.iterateN (forall a. Unbox a => Vector a -> Int
V.length Vector a
v) Vector a -> Vector a
interpVec Vector a
v
    interpVec :: Vector a -> Vector a
interpVec Vector a
v2 = forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith a -> a -> a
interp Vector a
v2 (forall a. Unbox a => Vector a -> Vector a
V.tail Vector a
v2)
{-# SPECIALIZE bernsteinSplit :: BernsteinPoly Double -> Double ->
    (BernsteinPoly Double, BernsteinPoly Double) #-}

addScaled :: (Unbox a, Num a) => ScaledPoly a -> ScaledPoly a -> ScaledPoly a
addScaled :: forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
addScaled ba :: ScaledPoly a
ba@(ScaledPoly Vector a
a) bb :: ScaledPoly a
bb@(ScaledPoly Vector a
b)
  | Int
la forall a. Ord a => a -> a -> Bool
< Int
lb = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$
              forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(+) (forall a. ScaledPoly a -> Vector a
scaledCoeffs forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled ScaledPoly a
ba forall a b. (a -> b) -> a -> b
$ Int
lbforall a. Num a => a -> a -> a
-Int
la) Vector a
b
  | Int
la forall a. Ord a => a -> a -> Bool
> Int
lb = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$
              forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(+) Vector a
a (forall a. ScaledPoly a -> Vector a
scaledCoeffs forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled ScaledPoly a
bb forall a b. (a -> b) -> a -> b
$ Int
laforall a. Num a => a -> a -> a
-Int
lb)
  | Bool
otherwise = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith forall a. Num a => a -> a -> a
(+) Vector a
a Vector a
b
  where la :: Int
la = forall a. Unbox a => Vector a -> Int
V.length Vector a
a
        lb :: Int
lb = forall a. Unbox a => Vector a -> Int
V.length Vector a
b
{-# SPECIALIZE addScaled :: ScaledPoly Double -> ScaledPoly Double ->
    ScaledPoly Double #-}        

-- | Sum two bernstein polynomials.  The final degree will be the
-- maximum of either degrees.
(~+) :: (Unbox a, Fractional a) =>
        BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a
BernsteinPoly a
a ~+ :: forall a.
(Unbox a, Fractional a) =>
BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a
~+ BernsteinPoly a
b = forall a.
(Unbox a, Fractional a) =>
ScaledPoly a -> BernsteinPoly a
fromScaled forall a b. (a -> b) -> a -> b
$ forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
addScaled (forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled BernsteinPoly a
a) (forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled BernsteinPoly a
b)
{-# INLINE (~+) #-}

subScaled :: (Unbox a, Num a) => ScaledPoly a -> ScaledPoly a -> ScaledPoly a
subScaled :: forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
subScaled ba :: ScaledPoly a
ba@(ScaledPoly Vector a
a) bb :: ScaledPoly a
bb@(ScaledPoly Vector a
b)
  | Int
la forall a. Ord a => a -> a -> Bool
< Int
lb = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$
              forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) (forall a. ScaledPoly a -> Vector a
scaledCoeffs forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled ScaledPoly a
ba forall a b. (a -> b) -> a -> b
$ Int
lbforall a. Num a => a -> a -> a
-Int
la) Vector a
b
  | Int
la forall a. Ord a => a -> a -> Bool
> Int
lb = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$
              forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) Vector a
a (forall a. ScaledPoly a -> Vector a
scaledCoeffs forall a b. (a -> b) -> a -> b
$ forall a. (Unbox a, Num a) => ScaledPoly a -> Int -> ScaledPoly a
degreeElevateScaled ScaledPoly a
bb forall a b. (a -> b) -> a -> b
$ Int
laforall a. Num a => a -> a -> a
-Int
lb)
  | Bool
otherwise = forall a. Vector a -> ScaledPoly a
ScaledPoly forall a b. (a -> b) -> a -> b
$ forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
V.zipWith (-) Vector a
a Vector a
b
  where la :: Int
la = forall a. Unbox a => Vector a -> Int
V.length Vector a
a
        lb :: Int
lb = forall a. Unbox a => Vector a -> Int
V.length Vector a
b
{-# SPECIALIZE subScaled :: ScaledPoly Double -> ScaledPoly Double ->
    ScaledPoly Double #-}        

-- | Subtract two bernstein polynomials.  The final degree will be the
-- maximum of either degrees.
(~-) :: (Unbox a, Fractional a) =>
        BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a

(forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled -> ScaledPoly a
a) ~- :: forall a.
(Unbox a, Fractional a) =>
BernsteinPoly a -> BernsteinPoly a -> BernsteinPoly a
~- (forall a. (Unbox a, Num a) => BernsteinPoly a -> ScaledPoly a
toScaled -> ScaledPoly a
b) = forall a.
(Unbox a, Fractional a) =>
ScaledPoly a -> BernsteinPoly a
fromScaled forall a b. (a -> b) -> a -> b
$ forall a.
(Unbox a, Num a) =>
ScaledPoly a -> ScaledPoly a -> ScaledPoly a
subScaled ScaledPoly a
a ScaledPoly a
b
{-# INLINE (~-) #-}

-- | Scale a bernstein polynomial by a constant.
(*~) :: (Unbox a, Num a) => a -> BernsteinPoly a -> BernsteinPoly a
a
a *~ :: forall a.
(Unbox a, Num a) =>
a -> BernsteinPoly a -> BernsteinPoly a
*~ (BernsteinPoly Vector a
v) = forall a. Vector a -> BernsteinPoly a
BernsteinPoly (forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (forall a. Num a => a -> a -> a
*a
a) Vector a
v)
{-# INLINE (*~) #-}