packages feed

aern2-fun (empty) → 0.2.9.0

raw patch · 7 files changed

+1034/−0 lines, 7 filesdep +QuickCheckdep +aern2-mpdep +aern2-real

Dependencies added: QuickCheck, aern2-mp, aern2-real, base, collect-errors, containers, hspec, mixed-types-num, psqueues

Files

+ LICENSE view
@@ -0,0 +1,12 @@+Copyright (c) 2015-2020 Michal Konecny+All rights reserved.++Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:++1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.++2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++3. Neither the name of the software nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ aern2-fun.cabal view
@@ -0,0 +1,60 @@+cabal-version: 1.12++-- This file has been generated from package.yaml by hpack version 0.34.4.+--+-- see: https://github.com/sol/hpack++name:           aern2-fun+version:        0.2.9.0+synopsis:       Generic operations for real functions+description:    Please see the README on GitHub at <https://github.com/michalkonecny/aern2/#readme>+category:       Math+homepage:       https://github.com/michalkonecny/aern2#readme+bug-reports:    https://github.com/michalkonecny/aern2/issues+author:         Michal Konecny, Eike Neumann+maintainer:     mikkonecny@gmail.com+copyright:      2015-2021 Michal Konecny, Eike Neumann+license:        BSD3+license-file:   LICENSE+build-type:     Simple++source-repository head+  type: git+  location: https://github.com/michalkonecny/aern2++library+  exposed-modules:+      AERN2.Interval+      AERN2.PQueue+      AERN2.RealFun.Operations+      AERN2.RealFun.SineCosine+      AERN2.RealFun.Tests+  other-modules:+      Paths_aern2_fun+  hs-source-dirs:+      src+  default-extensions:+      RebindableSyntax,+      ScopedTypeVariables,+      TypeFamilies,+      TypeOperators,+      ConstraintKinds,+      DefaultSignatures,+      MultiParamTypeClasses,+      FlexibleContexts,+      FlexibleInstances,+      UndecidableInstances+  other-extensions:+      TemplateHaskell+  ghc-options: -Wall+  build-depends:+      QuickCheck+    , aern2-mp >=0.2.9+    , aern2-real >=0.2.9+    , base ==4.*+    , collect-errors >=0.1.5+    , containers+    , hspec+    , mixed-types-num >=0.5.9+    , psqueues+  default-language: Haskell2010
+ src/AERN2/Interval.hs view
@@ -0,0 +1,260 @@+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE DeriveGeneric #-}+{-|+    Module      :  AERN2.Interval+    Description :  Intervals for use as function domains+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Intervals for use as function domains+-}+module AERN2.Interval+(+  Interval(..), singleton+  , width, split+  , arbitraryNonEmptyInterval+  , arbitraryNonEmptySmallInterval+  , intersect, intersects+  , DyadicInterval, CanBeDyadicInterval, dyadicInterval+  , RealInterval, CanBeRealInterval, realInterval+)+where++import MixedTypesNumPrelude+-- import Numeric.CollectErrors (NumErrors, CanTakeErrors(..))+import qualified Numeric.CollectErrors as CN++import qualified Prelude as P+import Text.Printf+-- import Text.Regex.TDFA++-- import Data.Maybe+++import GHC.Generics+import Data.Typeable++-- import qualified Data.List as List++-- import Test.Hspec+import Test.QuickCheck++import AERN2.MP.Enclosure+import AERN2.MP.Dyadic+import AERN2.MP.Ball hiding (intersect)+-- import qualified AERN2.MP.Ball as MPBall++import AERN2.Real++{- type -}++data Interval l r = Interval l r+  deriving (P.Eq, Generic)++instance (Show l, Show r) => Show (Interval l r) where+    show (Interval l r) =+      printf "Interval (%s) (%s)" (show l) (show r)+      -- printf "[%s,%s]" (show l) (show r)++instance (Read l, Read r) => Read (Interval l r) where+  readsPrec _pr intervalS+    | prefix1 == "Interval (" =+      case reads afterP1 of+        [(l,afterL)] ->+          if prefix2 == ") ("+            then+              case reads afterP2 of+                [(r,')':rest)] -> [(Interval l r, rest)]+                _ -> []+            else []+            where+            (prefix2, afterP2) = splitAt (length ") (") afterL+        _ -> []+    | otherwise = []+    where+    (prefix1, afterP1) = splitAt (length "Interval (") intervalS++singleton :: a -> Interval a a+singleton a = Interval a a++instance IsInterval (Interval e e) where+  type IntervalEndpoint (Interval e e) = e+  fromEndpoints l r = Interval l r+  endpoints (Interval l r) = (l,r)++width :: (CanSub r l) => Interval l r -> SubType r l+width (Interval l r) = r - l++split ::+  (CanAddSameType t, CanMulBy t Dyadic)+  =>+  (Interval t t) -> (Interval t t, Interval t t)+split (Interval l r) = (Interval l m, Interval m r)+  where+  m = (l + r)*(dyadic 0.5)++instance+  (Arbitrary l, Arbitrary r, HasOrderCertainlyAsymmetric l r)+  =>+  Arbitrary (Interval l r)+  where+  arbitrary =+    do+    l <- arbitrary+    r <- arbitrary+    if l !<=! r then return (Interval l r) else arbitrary++arbitraryNonEmptyInterval ::+  (Arbitrary l, Arbitrary r, HasOrderCertainlyAsymmetric l r)+  =>+  Gen (Interval l r)+arbitraryNonEmptyInterval =+  do+  l <- arbitrary+  r <- arbitrary+  if l !<! r then return (Interval l r) else arbitraryNonEmptyInterval++arbitraryNonEmptySmallInterval ::+  (Arbitrary e, CanAddThis e Integer)+  =>+  Gen (Interval e e)+arbitraryNonEmptySmallInterval =+  do+  l <- arbitrary+  w <- growingElements [1..10]+  return (Interval l (l+w))++{- containment -}++instance+  (HasOrderAsymmetric l l',  OrderCompareType l l' ~ Bool,+  HasOrderAsymmetric r' r,  OrderCompareType r' r ~ Bool)+  =>+  CanTestContains (Interval l r) (Interval l' r')+  where+  contains (Interval l r) (Interval l' r') =+    l <= l' && r' <= r++$(declForTypes+  [[t| Integer |], [t| Int |], [t| Rational |], [t| Dyadic |]]+  (\ t -> [d|+    instance+      (HasOrderAsymmetric l $t,  OrderCompareType l $t ~ Bool,+      HasOrderAsymmetric $t r,  OrderCompareType $t r ~ Bool)+      =>+      CanTestContains (Interval l r) $t+      where+      contains (Interval l r) e = l <= e && e <= r+  |]))++instance+  (CanSubSameType e, CanAddSubMulBy t e+  , HasIntegerBounds t, CanSubThis t Integer, CanDivBy t Integer)+  =>+  CanMapInside (Interval e e) t+  where+  mapInside (Interval l r) x =+    l + xUnit * (r - l)+    where+    xUnit = (x - xL) / (max 1 $ xU - xL)+    (xL,xU) = integerBounds x++{- intersection -}++instance+  (CanMinMaxSameType l, CanMinMaxSameType r, HasOrderCertainly l r)+  =>+  CanIntersectAsymmetric (Interval l r) (Interval l r)+  where+  type IntersectionType (Interval l r) (Interval l r) = CN (Interval l r)+  intersect (Interval l1 r1) (Interval l2 r2)+    | l !<=! r = pure (Interval l r)+    | l !>! r = CN.noValueNumErrorCertain err+    | otherwise = CN.prependErrorPotential err $ cn (Interval l r)+    where+    l = l1 `max` l2+    r = r1 `min` r2+    err = CN.NumError "empty intersection"++intersects ::+  (CanMinMaxSameType l, CanMinMaxSameType r, HasOrderCertainly l r)+  =>+  Interval l r -> Interval l r -> Bool+intersects i1 i2 = not . CN.hasError $ intersect i1 i2++{- comparison -}++instance+  (HasEqAsymmetric l1 l2, HasEqAsymmetric r1 r2+  , EqCompareType l1 l2 ~ EqCompareType r1 r2+  , CanAndOrSameType (EqCompareType l1 l2))+  =>+  HasEqAsymmetric (Interval l1 r1) (Interval l2 r2)+  where+  type EqCompareType (Interval l1 r1) (Interval l2 r2) = EqCompareType l1 l2+  equalTo (Interval l1 r1) (Interval l2 r2) =+    (l1 == l2) && (r1 == r2)++{- Dyadic intervals -}++type DyadicInterval = Interval Dyadic Dyadic+type CanBeDyadicInterval t = ConvertibleExactly t DyadicInterval++dyadicInterval :: (CanBeDyadicInterval t) => t -> DyadicInterval+dyadicInterval = convertExactly++instance+  (CanBeDyadic l, CanBeDyadic r, HasOrderCertainly l r, Show l, Show r,+   Typeable l, Typeable r)+  =>+  ConvertibleExactly (l, r) DyadicInterval where+  safeConvertExactly (l,r)+    | l !<=! r = Right $ Interval (dyadic l) (dyadic r)+    | otherwise = convError "endpoints are not in the correct order" (l,r)++instance ConvertibleExactly Dyadic DyadicInterval where+  safeConvertExactly d =+    Right $ Interval d d++instance ConvertibleExactly Integer DyadicInterval where+  safeConvertExactly n =+    do+    nD <- safeConvertExactly n+    Right $ Interval nD nD++instance ConvertibleExactly Rational DyadicInterval where+  safeConvertExactly r =+    do+    rD <- safeConvertExactly r+    Right $ Interval rD rD++instance ConvertibleExactly MPBall DyadicInterval where+  safeConvertExactly ball =+    Right $ Interval (centre l) (centre r)+    where+    (l,r) = endpointsAsIntervals ball++instance ConvertibleExactly DyadicInterval MPBall where+  safeConvertExactly (Interval lD rD) =+    Right $ fromEndpointsAsIntervals (mpBall lD) (mpBall rD)++{- CauchyReal intervals -}++type RealInterval = Interval CReal CReal+type CanBeRealInterval t = ConvertibleExactly t RealInterval++realInterval :: (CanBeRealInterval t) => t -> RealInterval+realInterval = convertExactly++instance+  (CanBeCReal l, CanBeCReal r, HasOrderCertainly l r, Show l, Show r,+   Typeable l, Typeable r)+  =>+  ConvertibleExactly (l, r) RealInterval where+  safeConvertExactly (l,r)+    | l !<=! r = Right $ Interval (creal l) (creal r)+    | otherwise = convError "endpoints are not in the correct order" (l,r)
+ src/AERN2/PQueue.hs view
@@ -0,0 +1,63 @@+{-|+    Module      :  AERN2.PQueue+    Description :  IntPSQ hiding keys and values+    Copyright   :  2016 (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    IntPSQ with simplified API, hiding the keys+-}++module AERN2.PQueue+(+  -- * Type+  PQueue+  -- * Query+  , null+  , size+  -- * Construction+  , singleton+  , empty+  -- * Insertion+  , insert+  -- * Deletion+  , minView+)+where++import Prelude hiding (null)+-- import qualified Prelude as P++import qualified Data.IntPSQ as Q++data PQueue p = PQueue { pq_queue :: Q.IntPSQ p (), pq_nextKey :: Int }++null :: PQueue p -> Bool+null = Q.null . pq_queue++size :: PQueue p -> Int+size = Q.size . pq_queue++singleton :: (Ord p) => p -> PQueue p+singleton p = PQueue (Q.singleton 0 p ()) 1++empty :: (Ord p) => PQueue p+empty = PQueue Q.empty 0++insert :: (Ord p) => p -> PQueue p -> PQueue p+insert p q+  | nextKey >= 0 =+    PQueue (Q.insert nextKey p () (pq_queue q)) (nextKey + 1)+  | otherwise =+    error "PQueue insert: key overflow"+  where+  nextKey = pq_nextKey q++minView :: (Ord p) => PQueue p -> Maybe (p, PQueue p)+minView q =+  case Q.minView (pq_queue q) of+    Just (_,p,_,qq) -> Just (p, q { pq_queue = qq })+    Nothing -> Nothing
+ src/AERN2/RealFun/Operations.hs view
@@ -0,0 +1,204 @@+{-# LANGUAGE CPP #-}+-- #define DEBUG+{-# LANGUAGE PartialTypeSignatures #-}+{-# OPTIONS_GHC -Wno-partial-type-signatures #-}+{-|+    Module      :  AERN2.RealFun.Operations+    Description :  Classes for real number function operations+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Classes for real number function operations+-}+module AERN2.RealFun.Operations+(+  HasDomain(..)+  , SameDomFnPair(..), ArbitraryWithDom(..)+  , CanApply(..)+  , CanApplyApprox(..), sampledRange+  , HasFnConstructorInfo(..)+  , HasConstFunctions, constFn, specEvalConstFn+  , HasVars(..), specEvalUnaryVarFn+  , CanMaximiseOverDom(..), CanMinimiseOverDom(..)+  , specCanMaximiseOverDom+  , CanIntegrateOverDom(..)+)+where++#ifdef DEBUG+import Debug.Trace (trace)+#define maybeTrace trace+#define maybeTraceIO putStrLn+#else+#define maybeTrace (\ (_ :: String) t -> t)+#define maybeTraceIO (\ (_ :: String) -> return ())+#endif++import MixedTypesNumPrelude+-- import qualified Prelude as P+import Text.Printf++-- import Data.Typeable++-- import qualified Data.List as List++import Test.Hspec+import Test.QuickCheck++import AERN2.Interval++import AERN2.MP.Dyadic+import AERN2.MP.Enclosure++{- domain -}++class HasDomain f where+  type Domain f+  getDomain :: f -> Domain f++data SameDomFnPair f = SameDomFnPair (f,f) deriving Show++instance (ArbitraryWithDom f, Arbitrary f) => (Arbitrary (SameDomFnPair f)) where+  arbitrary =+    do+    f1 <- arbitrary+    f2 <- arbitraryWithDom (getDomain f1)+    return $ SameDomFnPair (f1,f2)++class (HasDomain f) => ArbitraryWithDom f where+  arbitraryWithDom :: (Domain f) -> Gen f++{- evaluation -}++class CanApply f x where+  type ApplyType f x+  {-| compute @f(x)@  -}+  apply :: f {-^ @f@ -} -> x {-^ @x@ -} -> ApplyType f x++{-|+  Give an unsafe etimate of the function's range which is fast to compute.+  Intended to be used in optimisation heuristics.+-}+class CanApplyApprox f x where+  type ApplyApproxType f x+  {-| compute a cheap and unsafe approximation of @f(x)@  -}+  applyApprox :: f {-^ @f@ -} -> x {-^ @x@ -} -> ApplyApproxType f x++{-|+  Evaluate a function on a regular grid of the given size and return+  the largerst and smallest values found.  Useful for making instances+  of class 'CanApplyApprox'.+-}+sampledRange ::+  (CanApply f t, ApplyType f t ~ t,+   CanMinMaxSameType t, ConvertibleExactly Dyadic t, Show t)+  =>+  DyadicInterval -> Integer -> f -> Interval t t+sampledRange (Interval l r) depth f =+    maybeTrace+    ( "sampledRange:"+    ++ "\n samplePointsT = " ++ (show samplePointsT)+    ++ "\n samples = " ++ show samples+    ) $+    Interval minValue maxValue+    where+    minValue = foldl1 min samples+    maxValue = foldl1 max samples+    samples = map (apply f) samplePointsT+    samplePointsT = map convertExactly samplePoints+    _ = minValue : samplePointsT+    samplePoints :: [Dyadic]+    samplePoints = [(l*i + r*(size - i))*(dyadic (1/size)) | i <- [0..size]]+    size = 2^depth+++{- constructing basic functions -}++class HasFnConstructorInfo f where+  type FnConstructorInfo f+  getFnConstructorInfo :: f -> FnConstructorInfo f++type HasConstFunctions t f =+  (ConvertibleExactly (FnConstructorInfo f, t) f)++constFn :: (HasConstFunctions t f) => (FnConstructorInfo f) -> t -> f+constFn = curry convertExactly++specEvalConstFn ::+  _ => T c-> T f -> T x -> Spec+specEvalConstFn (T cName :: T c) (T fName :: T f) (T xName :: T x) =+  it (printf "Evaluating %s-constant functions %s on %s" cName fName xName) $+    property $+      \ (c :: c) (constrInfo :: FnConstructorInfo f) (xPres :: [x]) ->+        let f = constFn constrInfo c :: f in+        let dom = getDomain f in+        and $ flip map xPres $ \xPre ->+          apply f (mapInside dom xPre) ?==? c++class HasVars f where+  type Var f+  {-| the function @x@, ie the function that project the domain to the given variable @x@  -}+  varFn ::+    FnConstructorInfo f {-^ eg domain and/or accuracy guide -}->+    Var f {-^ @x@ -} ->+    f++specEvalUnaryVarFn ::+  _ => T f -> T x -> Spec+specEvalUnaryVarFn (T fName :: T f) (T xName :: T x) =+  it (printf "Evaluating variable functions %s on %s" fName xName) $ property $+    \ (constrInfo :: FnConstructorInfo f) (xPres :: [x]) ->+      and $ flip map xPres $ \xPre ->+        let f = varFn constrInfo () :: f in+        let x = mapInside (getDomain f) xPre in+        apply f x ?==? x+++{- range computation -}++class CanMaximiseOverDom f d where+  type MaximumOverDomType f d+  maximumOverDom :: f -> d -> MaximumOverDomType f d++class CanMinimiseOverDom f d where+  type MinimumOverDomType f d+  minimumOverDom :: f -> d -> MinimumOverDomType f d++-- specCanMaximiseOverDom ::+--   _ => (T f) -> Spec+-- specCanMaximiseOverDom (T fName :: T f) =+--   describe ("CanMaximiseOverDom " ++ fName) $ do+--     it "is consistent over a split domain" $ property $+--         \ (f :: f) ->+--           let dom = getDomain f in+--           let (dom1, dom2) = split dom in+--           let maxOnDom = maximumOverDom f dom in+--           let maxOnDom1 = maximumOverDom f dom1 in+--           let maxOnDom2 = maximumOverDom f dom2 in+--           maxOnDom ?>=? maxOnDom1+--           &&+--           maxOnDom ?>=? maxOnDom2++specCanMaximiseOverDom ::+  _ => (T f) -> (T x) -> Spec+specCanMaximiseOverDom (T fName :: T f) (T _xName :: T x) =+  describe ("CanMaximiseOverDom " ++ fName) $ do+    it "is consistent with evaluation" $ property $+      \ (f :: f) (xPres :: [x]) ->+        let dom = getDomain f in+        let maxOnDom = maximumOverDom f dom in+        and $ flip map xPres $ \xPre ->+          let x = mapInside dom xPre in+          let v1 = apply f x in+          maxOnDom ?>=? v1+++{- integration -}++class CanIntegrateOverDom f bounds where+  type IntegralOverDomType f bounds+  integrateOverDom :: f -> bounds -> IntegralOverDomType f bounds
+ src/AERN2/RealFun/SineCosine.hs view
@@ -0,0 +1,331 @@+{-# LANGUAGE CPP #-}+-- #define DEBUG+{-# LANGUAGE PartialTypeSignatures #-}+{-# OPTIONS_GHC -Wno-partial-type-signatures #-}+{-|+    Module      :  AERN2.RealFun.SineCosine+    Description :  Pointwise sine and cosine for functions+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Pointwise sine and cosine for functions+-}+module AERN2.RealFun.SineCosine+-- (+-- )+where++#ifdef DEBUG+import Debug.Trace (trace)+#define maybeTrace trace+#else+#define maybeTrace (flip const)+#endif++import MixedTypesNumPrelude+-- import qualified Prelude as P+import Text.Printf++import qualified Data.Map as Map+import qualified Data.List as List++-- import Test.Hspec+-- import Test.QuickCheck++import AERN2.MP+-- import qualified AERN2.MP.Ball as MPBall+-- import AERN2.MP.Dyadic++import AERN2.Real++-- import AERN2.Interval+import AERN2.RealFun.Operations++{-+    To compute sin(xC+-xE):++    * compute (rC+-rE) = range(xC)+    * compute k = round(rC/(pi/2))+    * compute sin or cos of txC = xC-k*pi/2 using Taylor series+      * use sin for even k and cos for odd k+      * which degree to use?+        * keep trying higher and higher degrees until+            * the accuracy of the result worsens+            * OR the accuracy of the result is 8x higher than xE+    * if k mod 4 = 2 then negate result+    * if k mod 4 = 3 then negate result+    * add xE to the error bound of the resulting polynomial+-}++sineWithAccuracyGuide ::+  _ => Accuracy -> f -> f+sineWithAccuracyGuide = sineCosineWithAccuracyGuide True++cosineWithAccuracyGuide ::+  _ => Accuracy -> f -> f+cosineWithAccuracyGuide = sineCosineWithAccuracyGuide False++sineCosineWithAccuracyGuide ::+  _ => Bool -> Accuracy -> f -> f+sineCosineWithAccuracyGuide isSine acGuide x =+    maybeTrace+    (+        "ChPoly.sineCosine: input:"+        ++ "\n isSine = " ++ show isSine+        -- ++ "\n xC = " ++ show xC+        ++ "\n xE = " ++ show xE+        ++ "\n xAccuracy = " ++ show xAccuracy+        ++ "\n r = " ++ show r+        ++ "\n k = " ++ show k+        ++ "\n trM = " ++ show trM+    ) $+    maybeTrace+    (+        "ChPoly.sineCosine: output:"+        ++ "\n Taylor series degree = " ++ show n+        ++ "\n getAccuracy taylorSum = " ++ show (getAccuracy taylorSum)+        ++ "\n taylorSumE = " ++ show taylorSumE+        ++ "\n getAccuracy result = " ++ show (getAccuracy res)+    ) $+--    xPoly (prec 100) -- dummy+    res+    where+    -- showB = show . getApproximate (bits 30)+    -- showAP = show . getApproximate (bits 50) . cheb2Power++    isCosine = not isSine++    -- first separate the centre of the polynomial x from its radius:+    xC = centreAsBall x+    xE = radius x+    xAccuracy = getAccuracy x++    -- compute (rC+-rE) = range(x):+    dom = getDomain x+    -- r = mpBall $ applyApprox xC (getDomain xC)++    r = fromEndpointsAsIntervals (mpBall $ minimumOverDom x dom) (mpBall $ maximumOverDom x dom)+    rC = centreAsBall r :: MPBall++    -- compute k = round(rC/(pi/2)):+    k = fst $ integerBounds $ 0.5 + (2*rC / pi)++    -- shift xC near 0 using multiples of pi/2:+    txC ac = (setPrecisionAtLeastAccuracy (ac) xC) - k * pi / 2+    -- work out an absolute range bound for txC:+    (_, trM) = endpointsAsIntervals $ abs $ r - k * pi / 2++    -- compute sin or cos of txC = xC-k*pi/2 using Taylor series:+    (taylorSum, taylorSumE, n)+        | isSine && even k = sineTaylorSum txC trM acGuide+        | isCosine && odd k = sineTaylorSum txC trM acGuide+        | otherwise = cosineTaylorSum txC trM acGuide+    -- if k mod 4 = 2 then negate result,+    -- if k mod 4 = 3 then negate result:+    km4 = k `mod` 4+    resC+        | isSine && 2 <= km4 && km4 <= 3 = -taylorSum+        | isCosine && 1 <= km4 && km4 <= 2 = -taylorSum+        | otherwise = taylorSum+    -- add xE to the error bound of the resulting polynomial:+    res = updateRadius (+ (taylorSumE + xE)) resC+++{-|+    For a given polynomial @p@, compute a partial Taylor sum of @cos(p)@ and return+    it together with its error bound @e@ and the degree of the polynomial @n@.+-}+sineTaylorSum ::+  _ => (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)+sineTaylorSum = sineCosineTaylorSum True++{-|+    For a given polynomial @p@, compute a partial Taylor sum of @cos(p)@ and return+    it together with its error bound @e@ and the degree of the polynomial @n@.+-}+cosineTaylorSum ::+  _ => (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)+cosineTaylorSum = sineCosineTaylorSum False++sineCosineTaylorSum ::+  _ => Bool -> (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)+sineCosineTaylorSum isSine (xAC :: Accuracy -> f) xM acGuidePre =+    let+    acGuide = acGuidePre + 4+    _isCosine = not isSine++    -- Work out the degree of the highest term we need to get the+    -- Lagrange error bound acGuide-accurate:+    n = Map.size factorialsE - 1 -- the last one is used only for the error term+    (_, (_,_,termSumEB)) = Map.findMax factorialsE -- the Lagrange error bound for T_n+    -- At the same time, compute the factorials and keep the Lagrange error bounds:+    factorialsE =+      maybeTrace ("sineCosineTaylorSum: n = " ++ show (Map.size res - 1)) res+      where+      res = Map.fromAscList $ takeUntilAccurate $ map addE factorials+      factorials = aux 0 1+        where aux i fc_i = (i,fc_i) : aux (i+1) (fc_i*(i+1))+      addE (i, fc_i) = (i, (fc_i, xM_i, e_i))+        where+        e_i = errorBound $ xM_i/fc_i+        xM_i = xM^i+      takeUntilAccurate (t_i@(i,(_fc_i, _xM_i,e_i)):rest)+        | getAccuracy e_i > acGuide && (even i == isSine) = [t_i]+        | otherwise = t_i : takeUntilAccurate rest+      takeUntilAccurate [] = error "sineCosineTaylorSum: internal error"++    -- Work out accuracy needed for each power x^n, given that x^n/n! should have+    -- accuracy around acGuide + 1:+    --    -log_2 (\eps/n!) ~ acGuide + 1+    --    -log_2 (\eps) ~ acGuide + 1 - (-log_2(1/n!))+    powerAccuracies0 =+      -- maybeTrace ("sineCosineTaylorSum: powerAccuracies0 = " ++ show res)+      res+      where+      res = Map.map aux factorialsE+      aux (fc_i,_xM_i,_e_i) =+        -- the accuracy needed of the power to give a sufficiently accurate term:+        acGuide + 1 + (bits $ getNormLog fc_i)+    -- Ensure the accuracies in powers are sufficient+    -- to compute accurate higher powers by their multiplications:+    powerAccuracies =+      -- maybeTrace ("sineCosineTaylorSum: powerAccuracies = " ++ show res)+      res+      where+      res =+        foldl updateAccuracies powerAccuracies0 $+          drop 1 $ reverse $ -- from second-highest down to the lowest+            drop 2 $ Map.toAscList powerAccuracies0 -- the 2 lowest are computed directly+      updateAccuracies powerACs (i, ac_i)+        | odd i && odd j =  -- pw_(2j+1) = x * pw_j * pw_j+            updateAC j (ac_i + log_pw_j + log_x) $+            updateAC 1 (ac_i + log_pw_j + log_pw_j) $+            powerACs+            -- pw_(2j+1) + e_pw2j1 =  (x+e_x) * (pw_j + e_pwj) * (pw_j + e_pwj)+            -- = e_x * e_pwj * e_pwj+            -- ...+            -- + x * e_pwj * pw_j -- assume this term puts most constraint on the size of e_pwj+            -- + e_x * pw_j * pw_j -- assume this term puts most constraint on the size of e_x+            -- ...+            -- + x*pw_j*pw_j+        | odd i  = -- pw_(2j+1) = x * pw_(j-1) * pw_(j+1)+            updateAC (j-1) (ac_i + log_pw_jU + log_x) $+            updateAC (j+1) (ac_i + log_pw_jD + log_x) $+            updateAC 1 (ac_i + log_pw_jU + log_pw_jD) $+            powerACs+        | even j = -- pw_(2j) = (power j) * (power j)+            updateAC j (ac_i + log_pw_j) $+            powerACs+        | otherwise = -- pw_(2j) = (power (j-1)) * (power (j+1))+            updateAC (j-1) (ac_i + log_pw_jU) $+            updateAC (j+1) (ac_i + log_pw_jD) $+            powerACs+        where+        updateAC k ac_k = Map.adjust (max ac_k) k+        j = i `divI` 2+        log_x = getLogXM 1+        log_pw_j = getLogXM j+        log_pw_jU = getLogXM (j+1)+        log_pw_jD = getLogXM (j-1)+        getLogXM k =+          case (Map.lookup k factorialsE) of+            Just (_fc_k,xM_k,_e_k) -> bits $ getNormLog xM_k+            _ -> error "sineCosineTaylorSum: internal error"++    x = case Map.lookup 1 powerAccuracies of+      Just ac1 -> xAC ac1+      _ -> error "sineCosineTaylorSum: internal error"++    -- Compute the powers needed for the terms, reducing their size while+    -- respecting the required accuracy:+    powers+      | isSine = powersSine+      | otherwise = powersCosine+      where+      powersSine =+        -- maybeTrace ("sineCosineTaylorSum: powerSine accuracies:\n"+        --   ++ (showPowerAccuracies res))+        res+        where+        res = foldl addPower initPowers $ zip [1..] [3,5..n]+        initPowers = Map.fromAscList [(1, x)]+        addPower prevPowers (j,i) =+          maybeTrace (showPowerDebug i rpw_i) $+          Map.insert i rpw_i prevPowers+          where+          rpw_i = reduce i pw_i+          pw_i+            | odd j = x * pwr j * pwr j+            | otherwise = x * pwr (j-1) * pwr (j+1)+          pwr k = case Map.lookup k prevPowers of+            Just r -> r+            _ -> error "sineCosineTaylorSum: internal error (powersSine: pwr k)"+      powersCosine =+        -- maybeTrace ("sineCosineTaylorSum: powerCosine accuracies:\n"+        --   ++ (showPowerAccuracies res))+        res+        where+        res = foldl addPower initPowers $ zip [2..] [4,6..n]+        initPowers = Map.fromAscList [(2, xxR)]+        xxR = reduce 2 $ x*x+        addPower prevPowers (j,i) =+          maybeTrace (showPowerDebug i rpw_i) $+          Map.insert i rpw_i prevPowers+          where+          rpw_i = reduce i pw_i+          pw_i+            | even j = pwr j * pwr j+            | otherwise = pwr (j-1) * pwr (j+1)+          pwr k = case Map.lookup k prevPowers of+            Just r -> r+            _ -> error "sineCosineTaylorSum: internal error (powersCosine: pwr k)"+      showPowerDebug :: Integer -> f -> String+      showPowerDebug i rpw_i =+        printf "power %d: accuracy req: %s, actual accuracy: %s" -- , degree: %d"+          i (show pa) (show $ getAccuracy rpw_i) -- (terms_degree $  poly_coeffs $ chPoly_poly p)+          where+          Just pa = Map.lookup i powerAccuracies+      -- showPowerAccuracies pwrs =+      --   unlines $ map showAAA $ Map.toAscList $+      --     Map.intersectionWith (\p (pa0, pa) -> (pa0,pa, p)) pwrs $+      --       Map.intersectionWith (,) powerAccuracies0 powerAccuracies+      --   where+      --   showAAA (i,(pa0,pa,p)) =+      --     printf "power %d: accuracy req 0: %s, accuracy req: %s, actual accuracy: %s" -- , degree: %d"+      --       i (show pa0) (show pa) (show $ getAccuracy p) -- (terms_degree $  poly_coeffs $ chPoly_poly p)+      reduce i = setPrecisionAtLeastAccuracy (ac_i + 10)+        where+        ac_i = case Map.lookup i powerAccuracies of+          Just ac -> ac+          _ -> error "sineCosineTaylorSum: internal error"+    termSum =+      maybeTrace ("sineCosineTaylorSum: term accuracies = "+        ++ (show (map (\(i,t) -> (i,getAccuracy t)) terms))) $+      maybeTrace ("sineCosineTaylorSum: term partial sum accuracies = "+        ++ (show (map (getAccuracy . sumP) (tail $ List.inits (map snd terms))))) $+      -- maybeTrace ("sineCosineTaylorSum: terms = " ++ (show terms)) $+      -- maybeTrace ("sineCosineTaylorSum: term partial sums = " +      --   ++ (show (map sumP (tail $ List.inits (map snd terms))))) $+      sumP (map snd terms) + initNum+      where+      sumP = foldl1 (+)+      terms =+        Map.toAscList $ Map.intersectionWithKey makeTerm powers factorialsE+      initNum | isSine = 0+              | otherwise = 1+    makeTerm i pwr (fact,_,_e) =+      sign * pwr/fact -- alternate signs+      where+      sign = if (even $ i `divI` 2) then 1 else -1+    in+    (termSum, termSumEB, n)+--+-- lookupForce :: P.Ord k => k -> Map.Map k a -> a+-- lookupForce j amap =+--     case Map.lookup j amap of+--         Just t -> t+--         Nothing -> error "internal error in SineCosine.lookupForce"
+ src/AERN2/RealFun/Tests.hs view
@@ -0,0 +1,104 @@+{-# LANGUAGE CPP #-}+-- #define DEBUG+{-# LANGUAGE PartialTypeSignatures #-}+{-# OPTIONS_GHC -Wno-partial-type-signatures #-}+{-|+    Module      :  AERN2.RealFun.Tests+    Description :  Test support for real number function operations+    Copyright   :  (c) Michal Konecny+    License     :  BSD3++    Maintainer  :  mikkonecny@gmail.com+    Stability   :  experimental+    Portability :  portable++    Test support for real number function operations+-}+module AERN2.RealFun.Tests+(+  FnAndDescr(..)+  , specFnPointwiseOp1, specFnPointwiseOp2+)+where++#ifdef DEBUG+import Debug.Trace (trace)+#define maybeTrace trace+#else+#define maybeTrace (flip const)+#endif++import MixedTypesNumPrelude+-- import qualified Prelude as P+-- import Text.Printf++-- import Data.Typeable++-- import qualified Data.List as List++import Test.Hspec+import Test.QuickCheck++-- import AERN2.Interval++-- import AERN2.MP.Dyadic+import AERN2.MP.Enclosure++import AERN2.RealFun.Operations++data FnAndDescr f = FnAndDescr f String++instance Show f => Show (FnAndDescr f) where+  show (FnAndDescr f descr) =+    show f ++ "[" ++ descr ++ "]"++instance (HasDomain f) => HasDomain (FnAndDescr f) where+  type Domain (FnAndDescr f) = Domain f+  getDomain (FnAndDescr f _) = getDomain f++specFnPointwiseOp2 ::+  _ =>+  (T f) -> (T x) ->+  String ->+  (f -> f -> f) ->+  (v -> v -> v) ->+  (FnAndDescr f -> FnAndDescr f) ->+  (FnAndDescr f -> FnAndDescr f) ->+  Spec+specFnPointwiseOp2+    (T fName :: T f) (T _xName :: T x)+    opName opFn (opVal :: v -> v -> v) reshapeFn1 reshapeFn2+  =+  it ("pointwise " ++ opName ++ " on " ++ fName +++      " corresponds to " ++ opName ++ " on values") $ property $+      \ (SameDomFnPair (f1Pre,f2Pre) :: SameDomFnPair (FnAndDescr f)) (xPres :: [x]) ->+          let FnAndDescr f1 _d1 = reshapeFn1 f1Pre in+          let FnAndDescr f2 _d2 = reshapeFn2 f2Pre in+          and $ flip map xPres $ \xPre ->+            let x = mapInside (getDomain f1) xPre in+            let v1 = apply f1 x in+            let v2 = apply f2 x in+            let vr = opVal v1 v2 in+            apply (opFn f1 f2) x ?==? (vr :: v)++specFnPointwiseOp1 ::+  _ =>+  (T f) -> (T x) ->+  String ->+  (f -> f) ->+  (v -> v) ->+  (FnAndDescr f -> FnAndDescr f) ->+  Spec+specFnPointwiseOp1+    (T fName :: T f) (T _xName :: T x)+    opName opFn (opVal :: v -> v) reshapeFn1+  =+  it ("pointwise " ++ opName ++ " on " ++ fName +++      " corresponds to " ++ opName ++ " on values") $ property $+      \ (f1Pre :: FnAndDescr f) (xPres :: [x]) ->+          let FnAndDescr f1 _d1 = reshapeFn1 f1Pre in+          and $ flip map xPres $ \xPre ->+            let x = mapInside (getDomain f1) xPre in+            let v1 = apply f1 x in+            let vr = opVal v1 in+            apply (opFn f1) x ?==? (vr :: v)