{-# 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 :: Accuracy -> f -> f
sineWithAccuracyGuide = Bool -> Accuracy -> f -> f
forall f.
(HasDomain f, CanMinimiseOverDom f (Domain f),
 ConvertibleExactly (MinimumOverDomType f (Domain f)) MPBall,
 ConvertibleExactly (MaximumOverDomType f (Domain f)) MPBall,
 CanMaximiseOverDom f (Domain f), IsBall f,
 CanSub f (CSequence MPBall), HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 CanNeg f,
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer,
 MulType f f ~ f, NegType f ~ f,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 SubType f (CSequence MPBall) ~ f) =>
Bool -> Accuracy -> f -> f
sineCosineWithAccuracyGuide Bool
True

cosineWithAccuracyGuide ::
  _ => Accuracy -> f -> f
cosineWithAccuracyGuide :: Accuracy -> f -> f
cosineWithAccuracyGuide = Bool -> Accuracy -> f -> f
forall f.
(HasDomain f, CanMinimiseOverDom f (Domain f),
 ConvertibleExactly (MinimumOverDomType f (Domain f)) MPBall,
 ConvertibleExactly (MaximumOverDomType f (Domain f)) MPBall,
 CanMaximiseOverDom f (Domain f), IsBall f,
 CanSub f (CSequence MPBall), HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 CanNeg f,
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer,
 MulType f f ~ f, NegType f ~ f,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 SubType f (CSequence MPBall) ~ f) =>
Bool -> Accuracy -> f -> f
sineCosineWithAccuracyGuide Bool
False

sineCosineWithAccuracyGuide ::
  _ => Bool -> Accuracy -> f -> f
sineCosineWithAccuracyGuide :: Bool -> Accuracy -> f -> f
sineCosineWithAccuracyGuide Bool
isSine Accuracy
acGuide f
x =
    maybeTrace
    (
        [Char]
"ChPoly.sineCosine: input:"
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n isSine = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Bool -> [Char]
forall a. Show a => a -> [Char]
show Bool
isSine
        -- ++ "\n xC = " ++ show xC
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n xE = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ ErrorBound -> [Char]
forall a. Show a => a -> [Char]
show ErrorBound
xE
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n xAccuracy = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Accuracy -> [Char]
forall a. Show a => a -> [Char]
show Accuracy
xAccuracy
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n r = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ MPBall -> [Char]
forall a. Show a => a -> [Char]
show MPBall
r
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n k = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
k
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n trM = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ MPBall -> [Char]
forall a. Show a => a -> [Char]
show MPBall
trM
    ) (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$
    maybeTrace
    (
        [Char]
"ChPoly.sineCosine: output:"
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n Taylor series degree = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Integer -> [Char]
forall a. Show a => a -> [Char]
show Integer
n
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n getAccuracy taylorSum = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Accuracy -> [Char]
forall a. Show a => a -> [Char]
show (f -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy f
taylorSum)
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n taylorSumE = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ ErrorBound -> [Char]
forall a. Show a => a -> [Char]
show ErrorBound
taylorSumE
        [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
"\n getAccuracy result = " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Accuracy -> [Char]
forall a. Show a => a -> [Char]
show (f -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy f
res)
    ) (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$
--    xPoly (prec 100) -- dummy
    f
res
    where
    -- showB = show . getApproximate (bits 30)
    -- showAP = show . getApproximate (bits 50) . cheb2Power

    isCosine :: NegType Bool
isCosine = Bool -> NegType Bool
forall t. CanNeg t => t -> NegType t
not Bool
isSine

    -- first separate the centre of the polynomial x from its radius:
    xC :: f
xC = f -> f
forall t. IsBall t => t -> t
centreAsBall f
x
    xE :: ErrorBound
xE = f -> ErrorBound
forall t. IsBall t => t -> ErrorBound
radius f
x
    xAccuracy :: Accuracy
xAccuracy = f -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy f
x

    -- compute (rC+-rE) = range(x):
    dom :: Domain f
dom = f -> Domain f
forall f. HasDomain f => f -> Domain f
getDomain f
x
    -- r = mpBall $ applyApprox xC (getDomain xC)

    r :: MPBall
r = MPBall -> MPBall -> MPBall
forall i.
(IsInterval i, CanMinMaxSameType (IntervalEndpoint i)) =>
i -> i -> i
fromEndpointsAsIntervals (MinimumOverDomType f (Domain f) -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall (MinimumOverDomType f (Domain f) -> MPBall)
-> MinimumOverDomType f (Domain f) -> MPBall
forall a b. (a -> b) -> a -> b
$ f -> Domain f -> MinimumOverDomType f (Domain f)
forall f d.
CanMinimiseOverDom f d =>
f -> d -> MinimumOverDomType f d
minimumOverDom f
x Domain f
dom) (MaximumOverDomType f (Domain f) -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall (MaximumOverDomType f (Domain f) -> MPBall)
-> MaximumOverDomType f (Domain f) -> MPBall
forall a b. (a -> b) -> a -> b
$ f -> Domain f -> MaximumOverDomType f (Domain f)
forall f d.
CanMaximiseOverDom f d =>
f -> d -> MaximumOverDomType f d
maximumOverDom f
x Domain f
dom)
    rC :: MPBall
rC = MPBall -> MPBall
forall t. IsBall t => t -> t
centreAsBall MPBall
r :: MPBall

    -- compute k = round(rC/(pi/2)):
    k :: Integer
k = (Integer, Integer) -> Integer
forall a b. (a, b) -> a
fst ((Integer, Integer) -> Integer) -> (Integer, Integer) -> Integer
forall a b. (a -> b) -> a -> b
$ MPBall -> (Integer, Integer)
forall t. HasIntegerBounds t => t -> (Integer, Integer)
integerBounds (MPBall -> (Integer, Integer)) -> MPBall -> (Integer, Integer)
forall a b. (a -> b) -> a -> b
$ Rational
0.5 Rational -> MPBall -> AddType Rational MPBall
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ (Integer
2Integer -> MPBall -> MulType Integer MPBall
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
*MPBall
rC MPBall -> CSequence MPBall -> DivType MPBall (CSequence MPBall)
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/ CSequence MPBall
pi)

    -- shift xC near 0 using multiples of pi/2:
    txC :: Accuracy -> SubType f (CSequence MPBall)
txC Accuracy
ac = (Accuracy -> f -> f
forall t. (HasPrecision t, CanSetPrecision t) => Accuracy -> t -> t
setPrecisionAtLeastAccuracy (Accuracy
ac) f
xC) f -> CSequence MPBall -> SubType f (CSequence MPBall)
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
k Integer -> CSequence MPBall -> MulType Integer (CSequence MPBall)
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* CSequence MPBall
pi CSequence MPBall -> Integer -> DivType (CSequence MPBall) Integer
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/ Integer
2
    -- work out an absolute range bound for txC:
    (MPBall
_, MPBall
trM) = MPBall -> (MPBall, MPBall)
forall i. IsInterval i => i -> (i, i)
endpointsAsIntervals (MPBall -> (MPBall, MPBall)) -> MPBall -> (MPBall, MPBall)
forall a b. (a -> b) -> a -> b
$ MPBall -> AbsType MPBall
forall t. CanAbs t => t -> AbsType t
abs (MPBall -> AbsType MPBall) -> MPBall -> AbsType MPBall
forall a b. (a -> b) -> a -> b
$ MPBall
r MPBall -> CSequence MPBall -> SubType MPBall (CSequence MPBall)
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
k Integer -> CSequence MPBall -> MulType Integer (CSequence MPBall)
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* CSequence MPBall
pi CSequence MPBall -> Integer -> DivType (CSequence MPBall) Integer
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/ Integer
2

    -- compute sin or cos of txC = xC-k*pi/2 using Taylor series:
    (f
taylorSum, ErrorBound
taylorSumE, Integer
n)
        | Bool
isSine Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
k = (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
forall f.
(HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 MulType f f ~ f,
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer) =>
(Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
sineTaylorSum Accuracy -> f
Accuracy -> SubType f (CSequence MPBall)
txC MPBall
trM Accuracy
acGuide
        | Bool
NegType Bool
isCosine Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer -> Bool
forall a. Integral a => a -> Bool
odd Integer
k = (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
forall f.
(HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 MulType f f ~ f,
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer) =>
(Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
sineTaylorSum Accuracy -> f
Accuracy -> SubType f (CSequence MPBall)
txC MPBall
trM Accuracy
acGuide
        | Bool
otherwise = (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
forall f.
(HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 MulType f f ~ f,
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer) =>
(Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
cosineTaylorSum Accuracy -> f
Accuracy -> SubType f (CSequence MPBall)
txC MPBall
trM Accuracy
acGuide
    -- if k mod 4 = 2 then negate result,
    -- if k mod 4 = 3 then negate result:
    km4 :: ModType Integer Integer
km4 = Integer
k Integer -> Integer -> ModType Integer Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> ModType t1 t2
`mod` Integer
4
    resC :: f
resC
        | Bool
isSine Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer
2 Integer -> Integer -> OrderCompareType Integer Integer
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
<= Integer
ModType Integer Integer
km4 Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer
ModType Integer Integer
km4 Integer -> Integer -> OrderCompareType Integer Integer
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
<= Integer
3 = -f
taylorSum
        | Bool
NegType Bool
isCosine Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer
1 Integer -> Integer -> OrderCompareType Integer Integer
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
<= Integer
ModType Integer Integer
km4 Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer
ModType Integer Integer
km4 Integer -> Integer -> OrderCompareType Integer Integer
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
<= Integer
2 = -f
taylorSum
        | Bool
otherwise = f
taylorSum
    -- add xE to the error bound of the resulting polynomial:
    res :: f
res = (ErrorBound -> ErrorBound) -> f -> f
forall t. IsBall t => (ErrorBound -> ErrorBound) -> t -> t
updateRadius (ErrorBound -> ErrorBound -> AddType ErrorBound ErrorBound
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ (ErrorBound
taylorSumE ErrorBound -> ErrorBound -> AddType ErrorBound ErrorBound
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ ErrorBound
xE)) f
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 :: (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
sineTaylorSum = Bool
-> (Accuracy -> f)
-> MPBall
-> Accuracy
-> (f, ErrorBound, Integer)
forall f.
(HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 MulType f f ~ f,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer) =>
Bool
-> (Accuracy -> f)
-> MPBall
-> Accuracy
-> (f, ErrorBound, Integer)
sineCosineTaylorSum Bool
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 :: (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
cosineTaylorSum = Bool
-> (Accuracy -> f)
-> MPBall
-> Accuracy
-> (f, ErrorBound, Integer)
forall f.
(HasPrecision f, CanSetPrecision f,
 CanDiv (MulType Integer f) Integer, CanMulAsymmetric f f,
 CanMulAsymmetric Integer f, HasAccuracy f,
 HasAccuracy (DivType (MulType Integer f) Integer),
 CanAddAsymmetric (DivType (MulType Integer f) Integer) Integer,
 MulType f f ~ f,
 AddType (DivType (MulType Integer f) Integer) Integer ~ f,
 CanAddAsymmetric
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer),
 AddType
   (DivType (MulType Integer f) Integer)
   (DivType (MulType Integer f) Integer)
 ~ DivType (MulType Integer f) Integer) =>
Bool
-> (Accuracy -> f)
-> MPBall
-> Accuracy
-> (f, ErrorBound, Integer)
sineCosineTaylorSum Bool
False

sineCosineTaylorSum ::
  _ => Bool -> (Accuracy -> f) -> MPBall -> Accuracy -> (f, ErrorBound, Integer)
sineCosineTaylorSum :: Bool
-> (Accuracy -> f)
-> MPBall
-> Accuracy
-> (f, ErrorBound, Integer)
sineCosineTaylorSum Bool
isSine (Accuracy -> f
xAC :: Accuracy -> f) MPBall
xM Accuracy
acGuidePre =
    let
    acGuide :: AddType Accuracy Integer
acGuide = Accuracy
acGuidePre Accuracy -> Integer -> AddType Accuracy Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
4
    _isCosine :: NegType Bool
_isCosine = Bool -> NegType Bool
forall t. CanNeg t => t -> NegType t
not Bool
isSine

    -- Work out the degree of the highest term we need to get the
    -- Lagrange error bound acGuide-accurate:
    n :: SubType Int Integer
n = Map Integer (Integer, MPBall, ErrorBound) -> Int
forall k a. Map k a -> Int
Map.size Map Integer (Integer, MPBall, ErrorBound)
factorialsE Int -> Integer -> SubType Int Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1 -- the last one is used only for the error term
    (Integer
_, (Integer
_,MPBall
_,ErrorBound
termSumEB)) = Map Integer (Integer, MPBall, ErrorBound)
-> (Integer, (Integer, MPBall, ErrorBound))
forall k a. Map k a -> (k, a)
Map.findMax Map Integer (Integer, MPBall, ErrorBound)
factorialsE -- the Lagrange error bound for T_n
    -- At the same time, compute the factorials and keep the Lagrange error bounds:
    factorialsE :: Map Integer (Integer, MPBall, ErrorBound)
factorialsE =
      maybeTrace ("sineCosineTaylorSum: n = " ++ show (Map.size res - 1)) res
      where
      res :: Map Integer (Integer, MPBall, ErrorBound)
res = [(Integer, (Integer, MPBall, ErrorBound))]
-> Map Integer (Integer, MPBall, ErrorBound)
forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList ([(Integer, (Integer, MPBall, ErrorBound))]
 -> Map Integer (Integer, MPBall, ErrorBound))
-> [(Integer, (Integer, MPBall, ErrorBound))]
-> Map Integer (Integer, MPBall, ErrorBound)
forall a b. (a -> b) -> a -> b
$ [(Integer, (Integer, MPBall, ErrorBound))]
-> [(Integer, (Integer, MPBall, ErrorBound))]
takeUntilAccurate ([(Integer, (Integer, MPBall, ErrorBound))]
 -> [(Integer, (Integer, MPBall, ErrorBound))])
-> [(Integer, (Integer, MPBall, ErrorBound))]
-> [(Integer, (Integer, MPBall, ErrorBound))]
forall a b. (a -> b) -> a -> b
$ ((Integer, Integer) -> (Integer, (Integer, MPBall, ErrorBound)))
-> [(Integer, Integer)]
-> [(Integer, (Integer, MPBall, ErrorBound))]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Integer) -> (Integer, (Integer, MPBall, ErrorBound))
addE [(Integer, Integer)]
factorials
      factorials :: [(Integer, Integer)]
factorials = Integer -> Integer -> [(Integer, Integer)]
forall t2 t.
(CanAddAsymmetric t2 Integer, CanMulAsymmetric t t2,
 AddType t2 Integer ~ t2, MulType t t2 ~ t) =>
t2 -> t -> [(t2, t)]
aux Integer
0 Integer
1
        where aux :: t2 -> t -> [(t2, t)]
aux t2
i t
fc_i = (t2
i,t
fc_i) (t2, t) -> [(t2, t)] -> [(t2, t)]
forall a. a -> [a] -> [a]
: t2 -> t -> [(t2, t)]
aux (t2
it2 -> Integer -> AddType t2 Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1) (t
fc_it -> t2 -> MulType t t2
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
*(t2
it2 -> Integer -> AddType t2 Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1))
      addE :: (Integer, Integer) -> (Integer, (Integer, MPBall, ErrorBound))
addE (Integer
i, Integer
fc_i) = (Integer
i, (Integer
fc_i, MPBall
PowType MPBall Integer
xM_i, ErrorBound
e_i))
        where
        e_i :: ErrorBound
e_i = MPBall -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound (MPBall -> ErrorBound) -> MPBall -> ErrorBound
forall a b. (a -> b) -> a -> b
$ MPBall
PowType MPBall Integer
xM_iMPBall -> Integer -> DivType MPBall Integer
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/Integer
fc_i
        xM_i :: PowType MPBall Integer
xM_i = MPBall
xMMPBall -> Integer -> PowType MPBall Integer
forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^Integer
i
      takeUntilAccurate :: [(Integer, (Integer, MPBall, ErrorBound))]
-> [(Integer, (Integer, MPBall, ErrorBound))]
takeUntilAccurate (t_i :: (Integer, (Integer, MPBall, ErrorBound))
t_i@(Integer
i,(Integer
_fc_i, MPBall
_xM_i,ErrorBound
e_i)):[(Integer, (Integer, MPBall, ErrorBound))]
rest)
        | ErrorBound -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy ErrorBound
e_i Accuracy -> Accuracy -> OrderCompareType Accuracy Accuracy
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
> Accuracy
AddType Accuracy Integer
acGuide Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& (Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
i Bool -> Bool -> EqCompareType Bool Bool
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Bool
isSine) = [(Integer, (Integer, MPBall, ErrorBound))
t_i]
        | Bool
otherwise = (Integer, (Integer, MPBall, ErrorBound))
t_i (Integer, (Integer, MPBall, ErrorBound))
-> [(Integer, (Integer, MPBall, ErrorBound))]
-> [(Integer, (Integer, MPBall, ErrorBound))]
forall a. a -> [a] -> [a]
: [(Integer, (Integer, MPBall, ErrorBound))]
-> [(Integer, (Integer, MPBall, ErrorBound))]
takeUntilAccurate [(Integer, (Integer, MPBall, ErrorBound))]
rest
      takeUntilAccurate [] = [Char] -> [(Integer, (Integer, MPBall, ErrorBound))]
forall a. HasCallStack => [Char] -> a
error [Char]
"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 :: Map Integer Accuracy
powerAccuracies0 =
      -- maybeTrace ("sineCosineTaylorSum: powerAccuracies0 = " ++ show res)
      Map Integer Accuracy
res
      where
      res :: Map Integer Accuracy
res = ((Integer, MPBall, ErrorBound) -> Accuracy)
-> Map Integer (Integer, MPBall, ErrorBound)
-> Map Integer Accuracy
forall a b k. (a -> b) -> Map k a -> Map k b
Map.map (Integer, MPBall, ErrorBound) -> Accuracy
(Integer, MPBall, ErrorBound) -> AddType Accuracy Accuracy
aux Map Integer (Integer, MPBall, ErrorBound)
factorialsE
      aux :: (Integer, MPBall, ErrorBound) -> AddType Accuracy Accuracy
aux (Integer
fc_i,MPBall
_xM_i,ErrorBound
_e_i) =
        -- the accuracy needed of the power to give a sufficiently accurate term:
        Accuracy
AddType Accuracy Integer
acGuide Accuracy -> Integer -> AddType Accuracy Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
1 Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ (NormLog -> Accuracy
forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits (NormLog -> Accuracy) -> NormLog -> Accuracy
forall a b. (a -> b) -> a -> b
$ Integer -> NormLog
forall a. HasNorm a => a -> NormLog
getNormLog Integer
fc_i)
    -- Ensure the accuracies in powers are sufficient
    -- to compute accurate higher powers by their multiplications:
    powerAccuracies :: Map Integer Accuracy
powerAccuracies =
      -- maybeTrace ("sineCosineTaylorSum: powerAccuracies = " ++ show res)
      Map Integer Accuracy
res
      where
      res :: Map Integer Accuracy
res =
        (Map Integer Accuracy
 -> (Integer, Accuracy) -> Map Integer Accuracy)
-> Map Integer Accuracy
-> [(Integer, Accuracy)]
-> Map Integer Accuracy
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Map Integer Accuracy -> (Integer, Accuracy) -> Map Integer Accuracy
updateAccuracies Map Integer Accuracy
powerAccuracies0 ([(Integer, Accuracy)] -> Map Integer Accuracy)
-> [(Integer, Accuracy)] -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
          Integer -> [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall n a. CanBeInteger n => n -> [a] -> [a]
drop Integer
1 ([(Integer, Accuracy)] -> [(Integer, Accuracy)])
-> [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall a b. (a -> b) -> a -> b
$ [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall a. [a] -> [a]
reverse ([(Integer, Accuracy)] -> [(Integer, Accuracy)])
-> [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall a b. (a -> b) -> a -> b
$ -- from second-highest down to the lowest
            Integer -> [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall n a. CanBeInteger n => n -> [a] -> [a]
drop Integer
2 ([(Integer, Accuracy)] -> [(Integer, Accuracy)])
-> [(Integer, Accuracy)] -> [(Integer, Accuracy)]
forall a b. (a -> b) -> a -> b
$ Map Integer Accuracy -> [(Integer, Accuracy)]
forall k a. Map k a -> [(k, a)]
Map.toAscList Map Integer Accuracy
powerAccuracies0 -- the 2 lowest are computed directly
      updateAccuracies :: Map Integer Accuracy -> (Integer, Accuracy) -> Map Integer Accuracy
updateAccuracies Map Integer Accuracy
powerACs (Integer
i, Accuracy
ac_i)
        | Integer -> Bool
forall a. Integral a => a -> Bool
odd Integer
i Bool -> Bool -> AndOrType Bool Bool
forall a b. CanAndOrAsymmetric a b => a -> b -> AndOrType a b
&& Integer -> Bool
forall a. Integral a => a -> Bool
odd Integer
DivIType Integer Integer
j =  -- pw_(2j+1) = x * pw_j * pw_j
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC Integer
DivIType Integer Integer
j (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_j Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_x) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC Integer
1 (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_j Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_j) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Map Integer Accuracy
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
        | Integer -> Bool
forall a. Integral a => a -> Bool
odd Integer
i  = -- pw_(2j+1) = x * pw_(j-1) * pw_(j+1)
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC (Integer
DivIType Integer Integer
jInteger -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1) (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jU Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_x) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC (Integer
DivIType Integer Integer
jInteger -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1) (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jD Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_x) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC Integer
1 (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jU Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jD) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Map Integer Accuracy
powerACs
        | Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
DivIType Integer Integer
j = -- pw_(2j) = (power j) * (power j)
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC Integer
DivIType Integer Integer
j (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_j) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Map Integer Accuracy
powerACs
        | Bool
otherwise = -- pw_(2j) = (power (j-1)) * (power (j+1))
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC (Integer
DivIType Integer Integer
jInteger -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1) (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jU) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Integer -> Accuracy -> Map Integer Accuracy -> Map Integer Accuracy
forall k t1 t2.
(Ord k, CanMinMaxAsymmetric t1 t2, MinMaxType t1 t2 ~ t2) =>
k -> t1 -> Map k t2 -> Map k t2
updateAC (Integer
DivIType Integer Integer
jInteger -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1) (Accuracy
ac_i Accuracy -> Accuracy -> AddType Accuracy Accuracy
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Accuracy
log_pw_jD) (Map Integer Accuracy -> Map Integer Accuracy)
-> Map Integer Accuracy -> Map Integer Accuracy
forall a b. (a -> b) -> a -> b
$
            Map Integer Accuracy
powerACs
        where
        updateAC :: k -> t1 -> Map k t2 -> Map k t2
updateAC k
k t1
ac_k = (t2 -> t2) -> k -> Map k t2 -> Map k t2
forall k a. Ord k => (a -> a) -> k -> Map k a -> Map k a
Map.adjust (t1 -> t2 -> MinMaxType t1 t2
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max t1
ac_k) k
k
        j :: DivIType Integer Integer
j = Integer
i Integer -> Integer -> DivIType Integer Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> DivIType t1 t2
`divI` Integer
2
        log_x :: Accuracy
log_x = Integer -> Accuracy
getLogXM Integer
1
        log_pw_j :: Accuracy
log_pw_j = Integer -> Accuracy
getLogXM Integer
DivIType Integer Integer
j
        log_pw_jU :: Accuracy
log_pw_jU = Integer -> Accuracy
getLogXM (Integer
DivIType Integer Integer
jInteger -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1)
        log_pw_jD :: Accuracy
log_pw_jD = Integer -> Accuracy
getLogXM (Integer
DivIType Integer Integer
jInteger -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1)
        getLogXM :: Integer -> Accuracy
getLogXM Integer
k =
          case (Integer
-> Map Integer (Integer, MPBall, ErrorBound)
-> Maybe (Integer, MPBall, ErrorBound)
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
k Map Integer (Integer, MPBall, ErrorBound)
factorialsE) of
            Just (Integer
_fc_k,MPBall
xM_k,ErrorBound
_e_k) -> NormLog -> Accuracy
forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits (NormLog -> Accuracy) -> NormLog -> Accuracy
forall a b. (a -> b) -> a -> b
$ MPBall -> NormLog
forall a. HasNorm a => a -> NormLog
getNormLog MPBall
xM_k
            Maybe (Integer, MPBall, ErrorBound)
_ -> [Char] -> Accuracy
forall a. HasCallStack => [Char] -> a
error [Char]
"sineCosineTaylorSum: internal error"

    x :: f
x = case Integer -> Map Integer Accuracy -> Maybe Accuracy
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
1 Map Integer Accuracy
powerAccuracies of
      Just Accuracy
ac1 -> Accuracy -> f
xAC Accuracy
ac1
      Maybe Accuracy
_ -> [Char] -> f
forall a. HasCallStack => [Char] -> a
error [Char]
"sineCosineTaylorSum: internal error"

    -- Compute the powers needed for the terms, reducing their size while
    -- respecting the required accuracy:
    powers :: Map Integer f
powers
      | Bool
isSine = Map Integer f
powersSine
      | Bool
otherwise = Map Integer f
powersCosine
      where
      powersSine :: Map Integer f
powersSine =
        -- maybeTrace ("sineCosineTaylorSum: powerSine accuracies:\n"
        --   ++ (showPowerAccuracies res))
        Map Integer f
res
        where
        res :: Map Integer f
res = (Map Integer f -> (Integer, Integer) -> Map Integer f)
-> Map Integer f -> [(Integer, Integer)] -> Map Integer f
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Map Integer f -> (Integer, Integer) -> Map Integer f
addPower Map Integer f
initPowers ([(Integer, Integer)] -> Map Integer f)
-> [(Integer, Integer)] -> Map Integer f
forall a b. (a -> b) -> a -> b
$ [Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
1..] [Integer
3,Integer
5..Integer
SubType Int Integer
n]
        initPowers :: Map Integer f
initPowers = [(Integer, f)] -> Map Integer f
forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList [(Integer
1, f
x)]
        addPower :: Map Integer f -> (Integer, Integer) -> Map Integer f
addPower Map Integer f
prevPowers (Integer
j,Integer
i) =
          maybeTrace (showPowerDebug i rpw_i) $
          Integer -> f -> Map Integer f -> Map Integer f
forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert Integer
i f
rpw_i Map Integer f
prevPowers
          where
          rpw_i :: f
rpw_i = Integer -> f -> f
reduce Integer
i f
pw_i
          pw_i :: f
pw_i
            | Integer -> Bool
forall a. Integral a => a -> Bool
odd Integer
j = f
x f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr Integer
j f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr Integer
j
            | Bool
otherwise = f
x f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr (Integer
jInteger -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1) f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr (Integer
jInteger -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1)
          pwr :: Integer -> f
pwr Integer
k = case Integer -> Map Integer f -> Maybe f
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
k Map Integer f
prevPowers of
            Just f
r -> f
r
            Maybe f
_ -> [Char] -> f
forall a. HasCallStack => [Char] -> a
error [Char]
"sineCosineTaylorSum: internal error (powersSine: pwr k)"
      powersCosine :: Map Integer f
powersCosine =
        -- maybeTrace ("sineCosineTaylorSum: powerCosine accuracies:\n"
        --   ++ (showPowerAccuracies res))
        Map Integer f
res
        where
        res :: Map Integer f
res = (Map Integer f -> (Integer, Integer) -> Map Integer f)
-> Map Integer f -> [(Integer, Integer)] -> Map Integer f
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl Map Integer f -> (Integer, Integer) -> Map Integer f
addPower Map Integer f
initPowers ([(Integer, Integer)] -> Map Integer f)
-> [(Integer, Integer)] -> Map Integer f
forall a b. (a -> b) -> a -> b
$ [Integer] -> [Integer] -> [(Integer, Integer)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
2..] [Integer
4,Integer
6..Integer
SubType Int Integer
n]
        initPowers :: Map Integer f
initPowers = [(Integer, f)] -> Map Integer f
forall k a. Eq k => [(k, a)] -> Map k a
Map.fromAscList [(Integer
2, f
xxR)]
        xxR :: f
xxR = Integer -> f -> f
reduce Integer
2 (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$ f
xf -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
*f
x
        addPower :: Map Integer f -> (Integer, Integer) -> Map Integer f
addPower Map Integer f
prevPowers (Integer
j,Integer
i) =
          maybeTrace (showPowerDebug i rpw_i) $
          Integer -> f -> Map Integer f -> Map Integer f
forall k a. Ord k => k -> a -> Map k a -> Map k a
Map.insert Integer
i f
rpw_i Map Integer f
prevPowers
          where
          rpw_i :: f
rpw_i = Integer -> f -> f
reduce Integer
i f
pw_i
          pw_i :: f
pw_i
            | Integer -> Bool
forall a. Integral a => a -> Bool
even Integer
j = Integer -> f
pwr Integer
j f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr Integer
j
            | Bool
otherwise = Integer -> f
pwr (Integer
jInteger -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1) f -> f -> MulType f f
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* Integer -> f
pwr (Integer
jInteger -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
1)
          pwr :: Integer -> f
pwr Integer
k = case Integer -> Map Integer f -> Maybe f
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
k Map Integer f
prevPowers of
            Just f
r -> f
r
            Maybe f
_ -> [Char] -> f
forall a. HasCallStack => [Char] -> a
error [Char]
"sineCosineTaylorSum: internal error (powersCosine: pwr k)"
      showPowerDebug :: Integer -> f -> String
      showPowerDebug :: Integer -> f -> [Char]
showPowerDebug Integer
i f
rpw_i =
        [Char] -> Integer -> [Char] -> [Char] -> [Char]
forall r. PrintfType r => [Char] -> r
printf [Char]
"power %d: accuracy req: %s, actual accuracy: %s" -- , degree: %d"
          Integer
i (Accuracy -> [Char]
forall a. Show a => a -> [Char]
show Accuracy
pa) (Accuracy -> [Char]
forall a. Show a => a -> [Char]
show (Accuracy -> [Char]) -> Accuracy -> [Char]
forall a b. (a -> b) -> a -> b
$ f -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getAccuracy f
rpw_i) -- (terms_degree $  poly_coeffs $ chPoly_poly p)
          where
          Just Accuracy
pa = Integer -> Map Integer Accuracy -> Maybe Accuracy
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
i Map Integer Accuracy
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 :: Integer -> f -> f
reduce Integer
i = Accuracy -> f -> f
forall t. (HasPrecision t, CanSetPrecision t) => Accuracy -> t -> t
setPrecisionAtLeastAccuracy (Accuracy
ac_i Accuracy -> Integer -> AddType Accuracy Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
10)
        where
        ac_i :: Accuracy
ac_i = case Integer -> Map Integer Accuracy -> Maybe Accuracy
forall k a. Ord k => k -> Map k a -> Maybe a
Map.lookup Integer
i Map Integer Accuracy
powerAccuracies of
          Just Accuracy
ac -> Accuracy
ac
          Maybe Accuracy
_ -> [Char] -> Accuracy
forall a. HasCallStack => [Char] -> a
error [Char]
"sineCosineTaylorSum: internal error"
    termSum :: f
termSum =
      maybeTrace ("sineCosineTaylorSum: term accuracies = "
        ++ (show (map (\(i,t) -> (i,getAccuracy t)) terms))) (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$
      maybeTrace ("sineCosineTaylorSum: term partial sum accuracies = "
        ++ (show (map (getAccuracy . sumP) (tail $ List.inits (map snd terms))))) (f -> f) -> f -> f
forall a b. (a -> b) -> a -> b
$
      -- maybeTrace ("sineCosineTaylorSum: terms = " ++ (show terms)) $
      -- maybeTrace ("sineCosineTaylorSum: term partial sums = " 
      --   ++ (show (map sumP (tail $ List.inits (map snd terms))))) $
      [DivType (MulType Integer f) Integer]
-> DivType (MulType Integer f) Integer
sumP (((Integer, DivType (MulType Integer f) Integer)
 -> DivType (MulType Integer f) Integer)
-> [(Integer, DivType (MulType Integer f) Integer)]
-> [DivType (MulType Integer f) Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, DivType (MulType Integer f) Integer)
-> DivType (MulType Integer f) Integer
forall a b. (a, b) -> b
snd [(Integer, DivType (MulType Integer f) Integer)]
terms) DivType (MulType Integer f) Integer
-> Integer -> AddType (DivType (MulType Integer f) Integer) Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
initNum
      where
      sumP :: [DivType (MulType Integer f) Integer]
-> DivType (MulType Integer f) Integer
sumP = (DivType (MulType Integer f) Integer
 -> DivType (MulType Integer f) Integer
 -> DivType (MulType Integer f) Integer)
-> [DivType (MulType Integer f) Integer]
-> DivType (MulType Integer f) Integer
forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 DivType (MulType Integer f) Integer
-> DivType (MulType Integer f) Integer
-> DivType (MulType Integer f) Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
(+)
      terms :: [(Integer, DivType (MulType Integer f) Integer)]
terms =
        Map Integer (DivType (MulType Integer f) Integer)
-> [(Integer, DivType (MulType Integer f) Integer)]
forall k a. Map k a -> [(k, a)]
Map.toAscList (Map Integer (DivType (MulType Integer f) Integer)
 -> [(Integer, DivType (MulType Integer f) Integer)])
-> Map Integer (DivType (MulType Integer f) Integer)
-> [(Integer, DivType (MulType Integer f) Integer)]
forall a b. (a -> b) -> a -> b
$ (Integer
 -> f
 -> (Integer, MPBall, ErrorBound)
 -> DivType (MulType Integer f) Integer)
-> Map Integer f
-> Map Integer (Integer, MPBall, ErrorBound)
-> Map Integer (DivType (MulType Integer f) Integer)
forall k a b c.
Ord k =>
(k -> a -> b -> c) -> Map k a -> Map k b -> Map k c
Map.intersectionWithKey Integer
-> f
-> (Integer, MPBall, ErrorBound)
-> DivType (MulType Integer f) Integer
forall t1 t2 t2 b c.
(Integral (DivIType t1 Integer), CanDivIMod t1 Integer,
 CanDiv (MulType Integer t2) t2, CanMulAsymmetric Integer t2) =>
t1 -> t2 -> (t2, b, c) -> DivType (MulType Integer t2) t2
makeTerm Map Integer f
powers Map Integer (Integer, MPBall, ErrorBound)
factorialsE
      initNum :: Integer
initNum | Bool
isSine = Integer
0
              | Bool
otherwise = Integer
1
    makeTerm :: t1 -> t2 -> (t2, b, c) -> DivType (MulType Integer t2) t2
makeTerm t1
i t2
pwr (t2
fact,b
_,c
_e) =
      Integer
IfThenElseType Bool Integer
sign Integer -> t2 -> MulType Integer t2
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* t2
pwrMulType Integer t2 -> t2 -> DivType (MulType Integer t2) t2
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/t2
fact -- alternate signs
      where
      sign :: IfThenElseType Bool Integer
sign = if (DivIType t1 Integer -> Bool
forall a. Integral a => a -> Bool
even (DivIType t1 Integer -> Bool) -> DivIType t1 Integer -> Bool
forall a b. (a -> b) -> a -> b
$ t1
i t1 -> Integer -> DivIType t1 Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> DivIType t1 t2
`divI` Integer
2) then Integer
1 else -Integer
1
    in
    (f
termSum, ErrorBound
termSumEB, Integer
SubType Int Integer
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"