{-# LANGUAGE CPP                    #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs                  #-}
module Math.Regression.Simple (
    -- * Regressions
    linear,
    linearWithErrors,
    quadratic,
    quadraticWithErrors,
    quadraticAndLinear,
    quadraticAndLinearWithErrors,
    -- * Operations
    Add (..),
    Eye (..),
    Mult (..),
    Det (..),
    Inv (..),
    -- * Zeros
    zerosLin,
    zerosQuad,
    optimaQuad,
    -- * Two dimensions
    V2 (..),
    M22 (..),
    -- * Three dimensions
    V3 (..),
    M33 (..),
    -- * Auxiliary classes
    Foldable' (..),
    IsDoublePair (..),
    ) where

import Data.Complex (Complex (..))

import qualified Data.List           as L
import qualified Data.Vector         as V
import qualified Data.Vector.Unboxed as U

-- $setup
-- >>> :set -XTypeApplications
--
-- >>> import Numeric (showFFloat)
--
-- Don't show too much decimal digits
--
-- >>> showDouble x = showFFloat @Double (Just (min 5 (5 - ceiling (logBase 10 x)))) x
-- >>> showDouble 123.456 ""
-- "123.46"
--
-- >>> showDouble 1234567 ""
-- "1234567"
--
-- >>> showDouble 123.4567890123456789 ""
-- "123.46"
--
-- >>> showDouble 0.0000000000000012345 ""
-- "0.00000"
-- 
-- >>> newtype PP a = PP a
-- >>> class Show' a where showsPrec' :: Int -> a -> ShowS
-- >>> instance Show' a => Show (PP a) where showsPrec d (PP x) = showsPrec' d x
-- >>> instance Show' Double where showsPrec' d x = if x < 0 then showParen (d > 6) (showChar '-' . showDouble (negate x)) else showDouble x
-- >>> instance (Show' a, Show' b) => Show' (a, b) where showsPrec' d (x, y) = showParen True $ showsPrec' 0 x . showString ", " . showsPrec' 0 y
-- >>> instance Show' V2 where showsPrec' d (V2 x y)   = showParen (d > 10) $ showString "V2 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y
-- >>> instance Show' V3 where showsPrec' d (V3 x y z) = showParen (d > 10) $ showString "V3 " . showsPrec' 11 x . showChar ' ' . showsPrec' 11 y . showChar ' ' . showsPrec' 11 z
--
-- Inputs:
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
--

-------------------------------------------------------------------------------
-- Classes
-------------------------------------------------------------------------------

-- | Addition
class Add a where
    zero :: a
    add  :: a -> a -> a

-- | Identity
class Eye a where
    eye :: a

-- | Multiplication of different things.
class Eye a => Mult a b c | a b -> c where
    mult :: a -> b -> c

-- | Determinant
class Eye a => Det a where
    det :: a -> Double

-- | Inverse
class Det a => Inv a where
    inv :: a -> a

infixl 6 `add`
infixl 7 `mult`

instance Eye Double where
    eye :: Double
eye = Double
1

instance Add Double where
    zero :: Double
zero = Double
0
    add :: Double -> Double -> Double
add = Double -> Double -> Double
forall a. Num a => a -> a -> a
(+)

instance Det Double where
    det :: Double -> Double
det = Double -> Double
forall a. a -> a
id

instance Inv Double where
    inv :: Double -> Double
inv = Double -> Double
forall a. Fractional a => a -> a
recip

-------------------------------------------------------------------------------
-- Zeros
-------------------------------------------------------------------------------

-- | Solve linear equation.
--
-- >>> zerosLin (V2 1 2)
-- -2.0
--
zerosLin :: V2 -> Double
zerosLin :: V2 -> Double
zerosLin (V2 Double
a Double
b) = Double -> Double
forall a. Num a => a -> a
negate (Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a)

-- | Solve quadratic equation.
--
-- >>> zerosQuad (V3 2 0 (-1))
-- Right (-0.7071067811865476,0.7071067811865476)
--
-- >>> zerosQuad (V3 2 0 1)
-- Left ((-0.0) :+ (-0.7071067811865476),(-0.0) :+ 0.7071067811865476)
--
-- Double root is not treated separately:
--
-- >>> zerosQuad (V3 1 0 0)
-- Right (-0.0,0.0)
--
-- >>> zerosQuad (V3 1 (-2) 1)
-- Right (1.0,1.0)
--
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad :: V3 -> Either (Complex Double, Complex Double) (Double, Double)
zerosQuad (V3 Double
a Double
b Double
c)
    | Double
delta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = (Complex Double, Complex Double)
-> Either (Complex Double, Complex Double) (Double, Double)
forall a b. a -> Either a b
Left ((-Double
bDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (-Double
sqrtNDeltaDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da), (-Double
bDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ (Double
sqrtNDeltaDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
da))
    | Bool
otherwise = (Double, Double)
-> Either (Complex Double, Complex Double) (Double, Double)
forall a b. b -> Either a b
Right ((- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sqrtDelta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
da, (-Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
sqrtDelta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
da)
  where
    delta :: Double
delta = Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
    sqrtDelta :: Double
sqrtDelta = Double -> Double
forall a. Floating a => a -> a
sqrt Double
delta
    sqrtNDelta :: Double
sqrtNDelta = Double -> Double
forall a. Floating a => a -> a
sqrt (- Double
delta)
    da :: Double
da = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a

-- | Find an optima point.
--
-- >>> optimaQuad (V3 1 (-2) 0)
-- 1.0
--
-- compare to
--
-- >>> zerosQuad (V3 1 (-2) 0)
-- Right (0.0,2.0)
--
optimaQuad :: V3 -> Double
optimaQuad :: V3 -> Double
optimaQuad (V3 Double
a Double
b Double
_) = V2 -> Double
zerosLin (Double -> Double -> V2
V2 (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) Double
b)

-------------------------------------------------------------------------------
-- 2 dimensions
-------------------------------------------------------------------------------

-- | 2d vector. Strict pair of 'Double's.
--
-- Also used to represent linear polynomial: @V2 a b@  \(= a x + b\).
--
data V2 = V2 !Double !Double
  deriving (V2 -> V2 -> Bool
(V2 -> V2 -> Bool) -> (V2 -> V2 -> Bool) -> Eq V2
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: V2 -> V2 -> Bool
$c/= :: V2 -> V2 -> Bool
== :: V2 -> V2 -> Bool
$c== :: V2 -> V2 -> Bool
Eq, Int -> V2 -> ShowS
[V2] -> ShowS
V2 -> String
(Int -> V2 -> ShowS)
-> (V2 -> String) -> ([V2] -> ShowS) -> Show V2
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [V2] -> ShowS
$cshowList :: [V2] -> ShowS
show :: V2 -> String
$cshow :: V2 -> String
showsPrec :: Int -> V2 -> ShowS
$cshowsPrec :: Int -> V2 -> ShowS
Show)

instance Add V2 where
    zero :: V2
zero = Double -> Double -> V2
V2 Double
0 Double
0
    add :: V2 -> V2 -> V2
add (V2 Double
x Double
y) (V2 Double
x' Double
y') = Double -> Double -> V2
V2 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x') (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Mult Double V2 V2 where
    mult :: Double -> V2 -> V2
mult Double
k (V2 Double
x Double
y) = Double -> Double -> V2
V2 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
    {-# INLINE mult #-}

-- | 2×2 matrix.
data M22 = M22 !Double !Double !Double !Double
  deriving (M22 -> M22 -> Bool
(M22 -> M22 -> Bool) -> (M22 -> M22 -> Bool) -> Eq M22
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: M22 -> M22 -> Bool
$c/= :: M22 -> M22 -> Bool
== :: M22 -> M22 -> Bool
$c== :: M22 -> M22 -> Bool
Eq, Int -> M22 -> ShowS
[M22] -> ShowS
M22 -> String
(Int -> M22 -> ShowS)
-> (M22 -> String) -> ([M22] -> ShowS) -> Show M22
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [M22] -> ShowS
$cshowList :: [M22] -> ShowS
show :: M22 -> String
$cshow :: M22 -> String
showsPrec :: Int -> M22 -> ShowS
$cshowsPrec :: Int -> M22 -> ShowS
Show)

instance Add M22 where
    zero :: M22
zero = Double -> Double -> Double -> Double -> M22
M22 Double
0 Double
0 Double
0 Double
0
    add :: M22 -> M22 -> M22
add (M22 Double
a Double
b Double
c Double
d) (M22 Double
a' Double
b' Double
c' Double
d') = Double -> Double -> Double -> Double -> M22
M22 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a') (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b') (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c') (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye M22 where
    eye :: M22
eye = Double -> Double -> Double -> Double -> M22
M22 Double
1 Double
0 Double
0 Double
1
    {-# INLINE eye #-}

instance Det M22 where det :: M22 -> Double
det = M22 -> Double
det2
instance Inv M22 where inv :: M22 -> M22
inv = M22 -> M22
inv2

instance Mult Double M22 M22 where
    mult :: Double -> M22 -> M22
mult Double
k (M22 Double
a Double
b Double
c Double
d) = Double -> Double -> Double -> Double -> M22
M22 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d)
    {-# INLINE mult #-}

instance Mult M22 V2 V2 where
    mult :: M22 -> V2 -> V2
mult (M22 Double
a Double
b Double
c Double
d) (V2 Double
u Double
v) = Double -> Double -> V2
V2 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v) (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v)
    {-# INLINE mult #-}

-- | >>> M22 1 2 3 4 `mult` eye @M22
-- M22 1.0 2.0 3.0 4.0
instance Mult M22 M22 M22 where
    mult :: M22 -> M22 -> M22
mult (M22 Double
a Double
b Double
c Double
d) (M22 Double
x Double
y Double
z Double
w) = Double -> Double -> Double -> Double -> M22
M22
        (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
        (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z) (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
    {-# INLINE mult #-}

det2 :: M22 -> Double
det2 :: M22 -> Double
det2 (M22 Double
a Double
b Double
c Double
d) = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c
{-# INLINE det2 #-}

inv2 :: M22 -> M22
inv2 :: M22 -> M22
inv2 m :: M22
m@(M22 Double
a Double
b Double
c Double
d) = Double -> Double -> Double -> Double -> M22
M22
    (  Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det) (- Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det)
    (- Double
c Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det) (  Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det)
  where
    det :: Double
det = M22 -> Double
det2 M22
m
{-# INLINE inv2 #-}

-------------------------------------------------------------------------------
-- 3 dimensions
-------------------------------------------------------------------------------

-- | 3d vector. Strict triple of 'Double's.
--
-- Also used to represent quadratic polynomial: @V3 a b c@  \(= a x^2 + b x + c\).
data V3 = V3 !Double !Double !Double
  deriving (V3 -> V3 -> Bool
(V3 -> V3 -> Bool) -> (V3 -> V3 -> Bool) -> Eq V3
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: V3 -> V3 -> Bool
$c/= :: V3 -> V3 -> Bool
== :: V3 -> V3 -> Bool
$c== :: V3 -> V3 -> Bool
Eq, Int -> V3 -> ShowS
[V3] -> ShowS
V3 -> String
(Int -> V3 -> ShowS)
-> (V3 -> String) -> ([V3] -> ShowS) -> Show V3
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [V3] -> ShowS
$cshowList :: [V3] -> ShowS
show :: V3 -> String
$cshow :: V3 -> String
showsPrec :: Int -> V3 -> ShowS
$cshowsPrec :: Int -> V3 -> ShowS
Show)

instance Add V3 where
    zero :: V3
zero = Double -> Double -> Double -> V3
V3 Double
0 Double
0 Double
0
    add :: V3 -> V3 -> V3
add (V3 Double
x Double
y Double
z) (V3 Double
x' Double
y' Double
z') = Double -> Double -> Double -> V3
V3 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x') (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y') (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
z')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Mult Double V3 V3 where
    mult :: Double -> V3 -> V3
mult Double
k (V3 Double
x Double
y Double
z) = Double -> Double -> Double -> V3
V3 (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
z)
    {-# INLINE mult #-}

-- | 3×3 matrix.
data M33 = M33
    !Double !Double !Double
    !Double !Double !Double
    !Double !Double !Double
  deriving (M33 -> M33 -> Bool
(M33 -> M33 -> Bool) -> (M33 -> M33 -> Bool) -> Eq M33
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: M33 -> M33 -> Bool
$c/= :: M33 -> M33 -> Bool
== :: M33 -> M33 -> Bool
$c== :: M33 -> M33 -> Bool
Eq, Int -> M33 -> ShowS
[M33] -> ShowS
M33 -> String
(Int -> M33 -> ShowS)
-> (M33 -> String) -> ([M33] -> ShowS) -> Show M33
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [M33] -> ShowS
$cshowList :: [M33] -> ShowS
show :: M33 -> String
$cshow :: M33 -> String
showsPrec :: Int -> M33 -> ShowS
$cshowsPrec :: Int -> M33 -> ShowS
Show)

instance Add M33 where
    zero :: M33
zero = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0 Double
0

    add :: M33 -> M33 -> M33
add (M33 Double
a Double
b Double
c Double
d Double
e Double
f Double
g Double
h Double
i) (M33 Double
a' Double
b' Double
c' Double
d' Double
e' Double
f' Double
g' Double
h' Double
i') = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33
        (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a') (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b') (Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c')
        (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d') (Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e') (Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f')
        (Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g') (Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h') (Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
i')
    {-# INLINE zero #-}
    {-# INLINE add #-}

instance Eye M33 where
    eye :: M33
eye = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
1 Double
0 Double
0 Double
0 Double
1 Double
0 Double
0 Double
0 Double
1
    {-# INLINE eye #-}

instance Det M33 where det :: M33 -> Double
det = M33 -> Double
det3
instance Inv M33 where inv :: M33 -> M33
inv = M33 -> M33
inv3

instance Mult Double M33 M33 where
    mult :: Double -> M33 -> M33
mult Double
k (M33 Double
a Double
b Double
c Double
d Double
e Double
f Double
g Double
h Double
i) = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33
        (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
c)
        (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
e) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
f)
        (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
g) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
i)
    {-# INLINE mult #-}

instance Mult M33 V3 V3 where
    mult :: M33 -> V3 -> V3
mult (M33 Double
a Double
b Double
c
           Double
d Double
e Double
f
           Double
g Double
h Double
i) (V3 Double
u Double
v Double
w) = Double -> Double -> Double -> V3
V3
        (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
        (Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
f Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
        (Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
    {-# INLINE mult #-}

-- TODO: instance Mult M33 M33 M33 where

det3 :: M33 -> Double
det3 :: M33 -> Double
det3 (M33 Double
a Double
b Double
c
          Double
d Double
e Double
f
          Double
g Double
h Double
i)
    = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
eDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
iDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
fDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
iDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
h) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
bDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
fDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
cDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
e)
{-# INLINE det3 #-}

inv3 :: M33 -> M33
inv3 :: M33 -> M33
inv3 m :: M33
m@(M33 Double
a Double
b Double
c
            Double
d Double
e Double
f
            Double
g Double
h Double
i)
    = Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
a' Double
b' Double
c'
          Double
d' Double
e' Double
f'
          Double
g' Double
h' Double
i'
  where
    a' :: Double
a' = Double -> Double -> Double -> Double -> Double
cofactor Double
e Double
f Double
h Double
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    b' :: Double
b' = Double -> Double -> Double -> Double -> Double
cofactor Double
c Double
b Double
i Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    c' :: Double
c' = Double -> Double -> Double -> Double -> Double
cofactor Double
b Double
c Double
e Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    d' :: Double
d' = Double -> Double -> Double -> Double -> Double
cofactor Double
f Double
d Double
i Double
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    e' :: Double
e' = Double -> Double -> Double -> Double -> Double
cofactor Double
a Double
c Double
g Double
i Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    f' :: Double
f' = Double -> Double -> Double -> Double -> Double
cofactor Double
c Double
a Double
f Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    g' :: Double
g' = Double -> Double -> Double -> Double -> Double
cofactor Double
d Double
e Double
g Double
h Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    h' :: Double
h' = Double -> Double -> Double -> Double -> Double
cofactor Double
b Double
a Double
h Double
g Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    i' :: Double
i' = Double -> Double -> Double -> Double -> Double
cofactor Double
a Double
b Double
d Double
e Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
det
    cofactor :: Double -> Double -> Double -> Double -> Double
cofactor Double
q Double
r Double
s Double
t = M22 -> Double
det2 (Double -> Double -> Double -> Double -> M22
M22 Double
q Double
r Double
s Double
t)
    det :: Double
det = M33 -> Double
det3 M33
m
{-# INLINE inv3 #-}

-------------------------------------------------------------------------------
-- Regressions
-------------------------------------------------------------------------------

-- | Linear regression.
--
-- The type is
--
-- @
-- 'linear' :: [('Double', 'Double')] -> 'V2'
-- @
--
-- but overloaded to work with boxed and unboxed 'Vector's.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> PP $ linear input1
-- V2 2.0000 1.00000
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ linear input2
-- V2 2.0063 0.88685
--
linear :: (Foldable' xs x, IsDoublePair x) => xs -> V2
linear :: xs -> V2
linear = (V2, V2) -> V2
forall a b. (a, b) -> a
fst ((V2, V2) -> V2) -> (xs -> (V2, V2)) -> xs -> V2
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> (V2, V2)
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> (V2, V2)
linearWithErrors

-- | Like 'linear' but also return parameters' standard errors.
--
-- To get confidence intervals you should multiply the errors
-- by @quantile (studentT (n - 2)) ci'@ from @statistics@ package
-- or similar.
-- For big @n@ using value 1 gives 68% interval and using value 2 gives 95% confidence interval.
-- See https://en.wikipedia.org/wiki/Student%27s_t-distribution#Table_of_selected_values
-- (@quantile@ calculates one-sided values, you need two-sided, thus adjust @ci@ value).
--
-- The first input is perfect fit:
--
-- >>> PP $ linearWithErrors input1
-- (V2 2.0000 1.00000, V2 0.00000 0.00000)
--
-- The second input is quite good:
--
-- >>> PP $ linearWithErrors input2
-- (V2 2.0063 0.88685, V2 0.09550 0.23826)
--
-- But the third input isn't so much,
-- standard error of a slope argument is 20%.
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ linearWithErrors input3
-- (V2 3.0000 1.00000, V2 0.63246 1.1832)
--
-- @since 0.1.1
--
linearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V2, V2)
linearWithErrors :: xs -> (V2, V2)
linearWithErrors = Kahan2 -> (V2, V2)
linearImpl (Kahan2 -> (V2, V2)) -> (xs -> Kahan2) -> xs -> (V2, V2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> Kahan2
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan2
kahan2

linearImpl :: Kahan2 -> (V2, V2)
linearImpl :: Kahan2 -> (V2, V2)
linearImpl (K2 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
y Double
_) (V2 Double
y2 Double
_) (V2 Double
xy Double
_)) =
    (V2
params, V2
errors)
  where
    n :: Double
    n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'

    matrix :: M22
matrix@(M22 Double
a11 Double
_ Double
_ Double
a22) = M22 -> M22
inv2 (Double -> Double -> Double -> Double -> M22
M22 Double
x2 Double
x Double
x Double
n)
    params :: V2
params@(V2 Double
a Double
b)          = M22 -> V2 -> V2
forall a b c. Mult a b c => a -> b -> c
mult M22
matrix (Double -> Double -> V2
V2 Double
xy Double
y)

    errors :: V2
errors = Double -> Double -> V2
V2 Double
sa Double
sb

    -- ensure that error is always non-negative.
    -- Due rounding errors, in perfect fit situations it can be slightly negative.
    err :: Double
err = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
    sa :: Double
sa  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a11 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2))
    sb :: Double
sb  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a22 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2))

-- | Quadratic regression.
--
-- The type is
--
-- @
-- 'quadratic' :: [('Double', 'Double')] -> 'V3'
-- @
--
-- but overloaded to work with boxed and unboxed 'Vector's.
--
-- >>> let input1 = [(0, 1), (1, 3), (2, 5)]
-- >>> quadratic input1
-- V3 0.0 2.0 1.0
--
-- >>> let input2 = [(0.1, 1.2), (1.3, 3.1), (1.9, 4.9), (3.0, 7.1), (4.1, 9.0)]
-- >>> PP $ quadratic input2
-- V3 (-0.00589) 2.0313 0.87155
--
-- >>> let input3 = [(0, 2), (1, 3), (2, 6), (3, 11)]
-- >>> PP $ quadratic input3
-- V3 1.00000 0.00000 2.0000
--
quadratic :: (Foldable' xs x, IsDoublePair x) => xs -> V3
quadratic :: xs -> V3
quadratic xs
data_ = M33 -> V3 -> V3
forall a b c. Mult a b c => a -> b -> c
mult (M33 -> M33
inv3 (Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
x2 Double
x Double
n Double
x3 Double
x2 Double
x Double
x4 Double
x3 Double
x2)) (Double -> Double -> Double -> V3
V3 Double
y Double
xy Double
x2y)
  where
    K3 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
x3 Double
_) (V2 Double
x4 Double
_) (V2 Double
y Double
_) (V2 Double
xy Double
_) (V2 Double
x2y Double
_) (V2 Double
_y2 Double
_) = xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 xs
data_

    n :: Double
    n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'

-- | Like 'quadratic' but also return parameters' standard errors.
--
-- >>> PP $ quadraticWithErrors input2
-- (V3 (-0.00589) 2.0313 0.87155, V3 0.09281 0.41070 0.37841)
--
-- >>> PP $ quadraticWithErrors input3
-- (V3 1.00000 0.00000 2.0000, V3 0.00000 0.00000 0.00000)
--
-- @since 0.1.1
--
quadraticWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V3)
quadraticWithErrors :: xs -> (V3, V3)
quadraticWithErrors = Kahan3 -> (V3, V3)
quadraticImpl (Kahan3 -> (V3, V3)) -> (xs -> Kahan3) -> xs -> (V3, V3)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3

quadraticImpl :: Kahan3 -> (V3, V3)
quadraticImpl :: Kahan3 -> (V3, V3)
quadraticImpl (K3 Int
n' (V2 Double
x Double
_) (V2 Double
x2 Double
_) (V2 Double
x3 Double
_) (V2 Double
x4 Double
_) (V2 Double
y Double
_) (V2 Double
xy Double
_) (V2 Double
x2y Double
_) (V2 Double
y2 Double
_)) =
    (V3
params, V3
errors)
  where
    n :: Double
    n :: Double
n = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n'

    matrix :: M33
matrix@(M33 Double
a11 Double
_   Double
_
                Double
_   Double
a22 Double
_
                Double
_   Double
_   Double
a33) = M33 -> M33
inv3 (Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> Double
-> M33
M33 Double
x4 Double
x3 Double
x2 Double
x3 Double
x2 Double
x Double
x2 Double
x Double
n)

    params :: V3
params@(V3 Double
a Double
b Double
c) = M33 -> V3 -> V3
forall a b c. Mult a b c => a -> b -> c
mult M33
matrix (Double -> Double -> Double -> V3
V3 Double
x2y Double
xy Double
y)

    errors :: V3
errors = Double -> Double -> Double -> V3
V3 Double
sa Double
sb Double
sc

    err :: Double
err = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double
y2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x2y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xy Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)
    sa :: Double
sa  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a11 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))
    sb :: Double
sb  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a22 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))
    sc :: Double
sc  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
a33 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3))

-- | Do both linear and quadratic regression in one data scan.
--
-- >>> PP $ quadraticAndLinear input2
-- (V3 (-0.00589) 2.0313 0.87155, V2 2.0063 0.88685)
--
quadraticAndLinear :: (Foldable' xs x, IsDoublePair x) => xs -> (V3, V2)
quadraticAndLinear :: xs -> (V3, V2)
quadraticAndLinear = ((V3, V2), (V3, V2)) -> (V3, V2)
forall a b. (a, b) -> a
fst (((V3, V2), (V3, V2)) -> (V3, V2))
-> (xs -> ((V3, V2), (V3, V2))) -> xs -> (V3, V2)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. xs -> ((V3, V2), (V3, V2))
forall xs x.
(Foldable' xs x, IsDoublePair x) =>
xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors

-- | Like 'quadraticAndLinear' but also return parameters' standard errors
--
-- >>> PP $ quadraticAndLinearWithErrors input2
-- ((V3 (-0.00589) 2.0313 0.87155, V2 2.0063 0.88685), (V3 0.09281 0.41070 0.37841, V2 0.09550 0.23826))
--
-- @since 0.1.1
--
quadraticAndLinearWithErrors :: (Foldable' xs x, IsDoublePair x) => xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors :: xs -> ((V3, V2), (V3, V2))
quadraticAndLinearWithErrors xs
data_ =
    ((V3
paramsQ, V2
paramsL), (V3
errorsQ, V2
errorsL))
  where
    k3 :: Kahan3
k3@(K3 Int
n V2
x V2
x2 V2
x3 V2
x4 V2
y V2
xy V2
x2y V2
y2) = xs -> Kahan3
forall xs x. (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 xs
data_
    k2 :: Kahan2
k2 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2 Int
n V2
x V2
x2 V2
y V2
y2 V2
xy

    (V2
paramsL, V2
errorsL) = Kahan2 -> (V2, V2)
linearImpl Kahan2
k2
    (V3
paramsQ, V3
errorsQ) = Kahan3 -> (V3, V3)
quadraticImpl Kahan3
k3

-------------------------------------------------------------------------------
-- Input
-------------------------------------------------------------------------------

-- | Like 'Foldable' but with element in the class definition.
class Foldable' xs x | xs -> x where
    foldl' :: (b -> x -> b) -> b -> xs -> b

instance              Foldable' [a]          a where foldl' :: (b -> a -> b) -> b -> [a] -> b
foldl' = (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
L.foldl'
instance              Foldable' (V.Vector a) a where foldl' :: (b -> a -> b) -> b -> Vector a -> b
foldl' = (b -> a -> b) -> b -> Vector a -> b
forall a b. (a -> b -> a) -> a -> Vector b -> a
V.foldl'
instance U.Unbox a => Foldable' (U.Vector a) a where foldl' :: (b -> a -> b) -> b -> Vector a -> b
foldl' = (b -> a -> b) -> b -> Vector a -> b
forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
U.foldl'

-- | Class witnessing that @dp@ has a pair of 'Double's.
class IsDoublePair dp where
    withDP :: dp -> (Double -> Double -> r) -> r
    makeDP :: Double -> Double -> dp

instance IsDoublePair V2 where
    withDP :: V2 -> (Double -> Double -> r) -> r
withDP (V2 Double
x Double
y) Double -> Double -> r
k = Double -> Double -> r
k Double
x Double
y
    makeDP :: Double -> Double -> V2
makeDP = Double -> Double -> V2
V2

instance (a ~ Double, b ~ Double) => IsDoublePair (a, b) where
    withDP :: (a, b) -> (Double -> Double -> r) -> r
withDP ~(a
x, b
y) Double -> Double -> r
k = Double -> Double -> r
k a
Double
x b
Double
y
    makeDP :: Double -> Double -> (a, b)
makeDP = (,)

-------------------------------------------------------------------------------
-- Kahan2
-------------------------------------------------------------------------------

data Kahan2 = K2
    { Kahan2 -> Int
k2n  :: {-# UNPACK #-} !Int
    , Kahan2 -> V2
k2x  :: {-# UNPACK #-} !V2
    , Kahan2 -> V2
k2x2 :: {-# UNPACK #-} !V2
    , Kahan2 -> V2
k2y  :: {-# UNPACK #-} !V2
    , Kahan2 -> V2
k2y2 :: {-# UNPACK #-} !V2
    , Kahan2 -> V2
k2xy :: {-# UNPACK #-} !V2
    }

zeroKahan2 :: Kahan2
zeroKahan2 :: Kahan2
zeroKahan2 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2 Int
0 V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero

-- | https://en.wikipedia.org/wiki/Kahan_summation_algorithm
addKahan :: V2 -> Double -> V2
addKahan :: V2 -> Double -> V2
addKahan (V2 Double
acc Double
c) Double
i =
    let y :: Double
y = Double
i Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
c
        t :: Double
t = Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y
    in Double -> Double -> V2
V2 Double
t ((Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
acc) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)

kahan2 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan2
kahan2 :: xs -> Kahan2
kahan2 = (Kahan2 -> x -> Kahan2) -> Kahan2 -> xs -> Kahan2
forall xs x b. Foldable' xs x => (b -> x -> b) -> b -> xs -> b
foldl' Kahan2 -> x -> Kahan2
forall dp. IsDoublePair dp => Kahan2 -> dp -> Kahan2
f Kahan2
zeroKahan2 where
    f :: Kahan2 -> dp -> Kahan2
f (K2 Int
n V2
x V2
x2 V2
y V2
y2 V2
xy) dp
uv = dp -> (Double -> Double -> Kahan2) -> Kahan2
forall dp r. IsDoublePair dp => dp -> (Double -> Double -> r) -> r
withDP dp
uv ((Double -> Double -> Kahan2) -> Kahan2)
-> (Double -> Double -> Kahan2) -> Kahan2
forall a b. (a -> b) -> a -> b
$ \Double
u Double
v -> Int -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan2
K2
        (Int -> Int
forall a. Enum a => a -> a
succ Int
n)
        (V2 -> Double -> V2
addKahan V2
x Double
u)
        (V2 -> Double -> V2
addKahan V2
x2 (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u))
        (V2 -> Double -> V2
addKahan V2
y Double
v)
        (V2 -> Double -> V2
addKahan V2
y2 (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
        (V2 -> Double -> V2
addKahan V2
xy (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))

-------------------------------------------------------------------------------
-- Kahan3
-------------------------------------------------------------------------------

data Kahan3 = K3
    { Kahan3 -> Int
k3n   :: {-# UNPACK #-} !Int
    , Kahan3 -> V2
k3x   :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3x2  :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3x3  :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3x4  :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3y   :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3xy  :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
k3x2y :: {-# UNPACK #-} !V2
    , Kahan3 -> V2
ky2   :: {-# UNPACK #-} !V2
    }

zeroKahan3 :: Kahan3
zeroKahan3 :: Kahan3
zeroKahan3 = Int -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan3
K3 Int
0 V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero V2
forall a. Add a => a
zero

kahan3 :: (Foldable' xs x, IsDoublePair x) => xs -> Kahan3
kahan3 :: xs -> Kahan3
kahan3 = (Kahan3 -> x -> Kahan3) -> Kahan3 -> xs -> Kahan3
forall xs x b. Foldable' xs x => (b -> x -> b) -> b -> xs -> b
foldl' Kahan3 -> x -> Kahan3
forall dp. IsDoublePair dp => Kahan3 -> dp -> Kahan3
f Kahan3
zeroKahan3 where
    f :: Kahan3 -> dp -> Kahan3
f (K3 Int
n V2
x V2
x2 V2
x3 V2
x4 V2
y V2
xy V2
x2y V2
y2) dp
uv = dp -> (Double -> Double -> Kahan3) -> Kahan3
forall dp r. IsDoublePair dp => dp -> (Double -> Double -> r) -> r
withDP dp
uv ((Double -> Double -> Kahan3) -> Kahan3)
-> (Double -> Double -> Kahan3) -> Kahan3
forall a b. (a -> b) -> a -> b
$ \Double
u Double
v ->
        let u2 :: Double
u2 = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u
        in Int -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> V2 -> Kahan3
K3
            (Int -> Int
forall a. Enum a => a -> a
succ Int
n)
            (V2 -> Double -> V2
addKahan V2
x Double
u)
            (V2 -> Double -> V2
addKahan V2
x2 Double
u2)
            (V2 -> Double -> V2
addKahan V2
x3 (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u2))
            (V2 -> Double -> V2
addKahan V2
x4 (Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u2))
            (V2 -> Double -> V2
addKahan V2
y Double
v)
            (V2 -> Double -> V2
addKahan V2
xy (Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
            (V2 -> Double -> V2
addKahan V2
x2y (Double
u2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))
            (V2 -> Double -> V2
addKahan V2
y2 (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v))