PowerSeries.hs

{-# LANGUAGE NoImplicitPrelude #-}
module Test.MathObj.PowerSeries where

import qualified MathObj.PowerSeries         as PST
import qualified MathObj.PowerSeries.Core    as PS
import qualified MathObj.PowerSeries.Example as PSE

import qualified Test.QuickCheck.Modifiers as Mod
import Test.NumericPrelude.Utility (equalInfLists, testUnit)
import Test.QuickCheck (quickCheck)
-- import Test.QuickCheck (Property, quickCheck, (==>))
import qualified Test.HUnit as HUnit


import NumericPrelude.Base as P
import NumericPrelude.Numeric as NP


identitiesExplODE, identitiesSeriesFunction, identitiesInverses ::
   [(String, Int, [Rational],[Rational])]

identitiesExplODE =
   ("exp",   500, PSE.expExpl,   PSE.expODE) :
   ("sin",   500, PSE.sinExpl,   PSE.sinODE) :
   ("cos",   500, PSE.cosExpl,   PSE.cosODE) :
   ("tan",    50, PSE.tanExpl,   PSE.tanODE) :
   ("tan",    50, PSE.tanExpl,   PSE.tanExplSieve) :
   ("tan",    50, PSE.tanODE,    PSE.tanODESieve) :
   ("log",   500, PSE.logExpl,   PSE.logODE) :
   ("asin",   50, PSE.asinODE,   snd (PS.inv PSE.sinODE)) :
   ("atan",  500, PSE.atanExpl,  PSE.atanODE) :
   ("sinh",  500, PSE.sinhExpl,  PSE.sinhODE) :
   ("cosh",  500, PSE.coshExpl,  PSE.coshODE) :
   ("atanh", 500, PSE.atanhExpl, PSE.atanhODE) :
   ("sqrt",  100, PSE.sqrtExpl,  PSE.sqrtODE) :
   []

identitiesSeriesFunction =
   ("exp",   500, PSE.expExpl,  PS.exp (\0 -> 1) [0,1]) :
   ("sin",   500, PSE.sinExpl,  PS.sin (\0 -> (0,1)) [0,1]) :
   ("cos",   500, PSE.cosExpl,  PS.cos (\0 -> (0,1)) [0,1]) :
   ("tan",    50, PSE.tanExpl,  PS.tan (\0 -> (0,1)) [0,1]) :
   ("sqrt",   50, PSE.sqrtExpl, PS.sqrt (\1 -> 1) [1,1]) :
   ("power", 500, PSE.powExpl (-1/3), PS.pow (\1 -> 1) (-1/3) [1,1]) :
   ("power",  50, PSE.powExpl (-1/3), PS.exp (\0 -> 1) (PS.scale (-1/3) PSE.log)) :
   ("log",   500, PSE.logExpl, PS.log (\1 -> 0) [1,1]) :
   ("asin",   50, PSE.asin, PS.asin (\1 -> 1) (\0 -> 0) [0,1]) :
 --  ("acos",  50, PSE.acos, PS.acos (\1 -> 1) (\0 -> pi/2) [0,1]) :
   ("atan",  500, PSE.atan, PS.atan (\0 -> 0) [0,1]) :
   []

identitiesInverses =
   ("exp",   100, 1:1:repeat 0, PS.exp  (\0 -> 1) PSE.log) :
   ("log",   100, 0:1:repeat 0, PS.log  (\1 -> 0) PSE.exp) :
   ("tan",    50, 0:1:repeat 0, PS.tan  (\0 -> (0,1)) PSE.atan) :
   ("atan",   50, 0:1:repeat 0, PS.atan (\0 -> 0) PSE.tan) :
   ("sin",    50, 0:1:repeat 0, PS.sin  (\0 -> (0,1)) PSE.asin) :
   ("asin",  100, 0:1:repeat 0, PS.asin (\1 -> 1) (\0 -> 0) PSE.sin) :
   ("sqrt",  500, 1:1:repeat 0, PS.sqrt (\1 -> 1) (PS.mul [1,1] [1,1])) :
   []

identitiesHoles :: [(String, Int, [Rational] -> [Rational], Rational)]
identitiesHoles =
   ("exp",    30, PS.exp  (\0 -> 1), 0) :
   ("log",    30, PS.log  (\1 -> 0), 1) :
   ("tan",    20, PS.tan  (\0 -> (0,1)), 0) :
   ("atan",   20, PS.atan (\0 -> 0), 0) :
   ("sin",    20, PS.sin  (\0 -> (0,1)), 0) :
   ("cos",    20, PS.cos  (\0 -> (0,1)), 0) :
   ("asin",   30, PS.asin (\1 -> 1) (\0 -> 0), 0) :
   ("sqrt",   50, PS.sqrt (\1 -> 1), 1) :
   ("pow13",  30, PS.pow  (\1 -> 1) (1/3), 1) :
   ("pow25",  30, PS.pow  (\1 -> 1) (2/5), 1) :
   []

testSeriesIdentity :: (String, Int, [Rational], [Rational]) -> HUnit.Test
testSeriesIdentity (label, len, x, y) =
   HUnit.test (HUnit.assertBool label (equalInfLists len [x,y]))

testSeriesIdentities ::
   String -> [(String, Int, [Rational], [Rational])] -> HUnit.Test
testSeriesIdentities label ids =
   HUnit.TestLabel label $
     HUnit.TestList $ map testSeriesIdentity ids

_checkSeriesIdentities ::
   [(String, Int, [Rational], [Rational])] -> [(String,Bool)]
_checkSeriesIdentities =
   map (\(label, len, x, y) -> (label, equalInfLists len [x,y]))


holesMultiplicative :: Int -> Int -> Int -> [Rational] -> Bool
holesMultiplicative trunc expon0 expon1 xs =
   let n0 = 1 + mod expon0 10
       n1 = 1 + mod expon1 10
   in  equalInfLists trunc
          [PS.insertHoles n0 $ PS.insertHoles n1 xs,
           PS.insertHoles n1 $ PS.insertHoles n0 xs,
           PS.insertHoles (n0*n1) xs]

testHolesIdentity ::
   (String, Int, [Rational] -> [Rational], Rational) -> HUnit.Test
testHolesIdentity (label, len, f, x0) =
   HUnit.test $ testUnit $ (,) ("holes in " ++ label) $
   quickCheck $ \expon0 xs -> checkHoles len expon0 f x0 xs


checkHoles ::
   Int -> Int -> ([Rational] -> [Rational]) ->
   Rational -> [Rational] -> Bool
checkHoles trunc expon0 f x xs =
   let expon = 1 + mod expon0 10
   in  equalInfLists trunc
          [(f $ PS.insertHoles expon (x:xs)) ++ repeat zero,
           (PS.insertHoles expon $ f $ x:xs) ++ repeat zero]


powerMultSeries :: Int -> Integer -> Mod.Positive Rational -> [Rational] -> Bool
powerMultSeries trunc expon0 xp xs =
   let expon = 1 + mod expon0 10
       x = Mod.getPositive xp
       xt = x:xs
   in  equalInfLists trunc
          [PS.pow
              (const x) (1 % expon)
              (PST.coeffs (PST.fromCoeffs xt ^ expon))
            ++ repeat zero,
           xt ++ repeat zero]

powerMult :: Int -> Rational -> Rational -> Bool
powerMult trunc exp0 exp1 =
   equalInfLists trunc
      [PS.mul (PSE.pow exp0) (PSE.pow exp1), PSE.pow (exp0+exp1)]

powerExplODE :: Int -> Rational -> Bool
powerExplODE trunc expon =
   equalInfLists trunc [PSE.powODE expon, PSE.powExpl expon]

invDiff :: Int -> Rational -> Mod.NonZero Rational -> [Rational] -> Bool
invDiff trunc x0 x1 xs_ =
   let xs = x0 : Mod.getNonZero x1 : xs_
       (y,ys) = PS.inv xs
       (z,zs) = PS.invDiff xs
   in  y==z && equalInfLists trunc [ys, zs]


tests :: HUnit.Test
tests =
   HUnit.TestLabel "power series" $
   HUnit.TestList [
      testSeriesIdentities "explicit vs. ODE solution" identitiesExplODE,
      testSeriesIdentities "transcendent functions of series" identitiesSeriesFunction,
      testSeriesIdentities "inverses of some series" identitiesInverses,
      HUnit.TestLabel "laws" $
      HUnit.TestList $
         map testHolesIdentity identitiesHoles
         ++
         (map testUnit $
            ("multiplicative holes",   quickCheck (holesMultiplicative 100)) :
            ("powers of series",       quickCheck (powerMultSeries 15)) :
            ("products of powers",     quickCheck (powerMult 30)) :
            ("power explicit vs. ODE", quickCheck (powerExplODE 50)) :
            ("inv vs. invDiff",        quickCheck (invDiff 15)) :
            [])
    ]