{-# LANGUAGE TemplateHaskell #-}
{-# OPTIONS_GHC -Wno-orphans #-}
{-|
    Module      :  AERN2.MP.Ball.Elementary
    Description :  Elementary operations on arbitrary precision dyadic balls
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

    Maintainer  :  mikkonecny@gmail.com
    Stability   :  experimental
    Portability :  portable

    Elementary operations on arbitrary precision dyadic balls
-}
module AERN2.MP.Ball.Elementary
(
  -- * Ball operations (see also instances)
  piBallP
  -- * Helpers for constructing ball functions
  , fromApproxWithLipschitz
)
where

import MixedTypesNumPrelude
import qualified Prelude as P

import AERN2.Normalize

import AERN2.MP.Dyadic (Dyadic)
import qualified AERN2.MP.Float as MPFloat
import AERN2.MP.Float (MPFloat, mpFloat, ceduCentreErr)
-- import AERN2.MP.Float.Operators
import AERN2.MP.Precision
-- import qualified AERN2.MP.ErrorBound as EB
import AERN2.MP.ErrorBound (errorBound)

import AERN2.MP.Ball.Type
import AERN2.MP.Ball.Conversions ()
import AERN2.MP.Ball.Comparisons ()
import AERN2.MP.Ball.Field (mulByEndpoints)


{- trigonometrics -}

piBallP :: Precision -> MPBall
piBallP :: Precision -> MPBall
piBallP Precision
p = MPFloat -> ErrorBound -> MPBall
MPBall MPFloat
piC (MPFloat -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound MPFloat
piErr)
  where
  (MPFloat
piC, MPFloat
piErr) = BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a. BoundsCEDU a -> (a, a)
MPFloat.ceduCentreErr (BoundsCEDU MPFloat -> (MPFloat, MPFloat))
-> BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a b. (a -> b) -> a -> b
$ Precision -> BoundsCEDU MPFloat
MPFloat.piCEDU Precision
p

instance CanSinCos MPBall where
  sin :: MPBall -> SinCosType MPBall
sin = Integer -> MPBall -> MPBall
sinB Integer
1
  cos :: MPBall -> SinCosType MPBall
cos = Integer -> MPBall -> MPBall
cosB Integer
1

sinB :: Integer -> MPBall -> MPBall
sinB :: Integer -> MPBall -> MPBall
sinB Integer
i MPBall
x =
    -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.sinDown MPFloat.sinUp lip) x
    Integer -> MPBall -> MinMaxType Integer MPBall
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max (-Integer
1) (MPBall -> MinMaxType Integer MPBall)
-> MPBall -> MinMaxType Integer MPBall
forall a b. (a -> b) -> a -> b
$ Integer -> MPBall -> MinMaxType Integer MPBall
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
min Integer
1 (MPBall -> MinMaxType Integer MPBall)
-> MPBall -> MinMaxType Integer MPBall
forall a b. (a -> b) -> a -> b
$
    (MPFloat -> BoundsCEDU MPFloat) -> MPFloat -> MPBall -> MPBall
fromApproxWithLipschitz MPFloat -> BoundsCEDU MPFloat
MPFloat.sinCEDU MPFloat
lip MPBall
x
    where
    lip :: MPFloat
lip
        | Integer
i Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
0 = Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
1
        | Bool
otherwise = MPBall -> IntervalEndpoint MPBall
forall i. IsInterval i => i -> IntervalEndpoint i
endpointR (MPBall -> IntervalEndpoint MPBall)
-> MPBall -> IntervalEndpoint 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
$ Integer -> MPBall -> MPBall
cosB (Integer
i Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1) MPBall
x

cosB :: Integer -> MPBall -> MPBall
cosB :: Integer -> MPBall -> MPBall
cosB Integer
i MPBall
x =
    -- increasingPrecisionUntilNotImproving (fromApproxWithLipschitz MPFloat.cosDown MPFloat.cosUp lip) x
    Integer -> MPBall -> MinMaxType Integer MPBall
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max (-Integer
1) (MPBall -> MinMaxType Integer MPBall)
-> MPBall -> MinMaxType Integer MPBall
forall a b. (a -> b) -> a -> b
$ Integer -> MPBall -> MinMaxType Integer MPBall
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
min Integer
1 (MPBall -> MinMaxType Integer MPBall)
-> MPBall -> MinMaxType Integer MPBall
forall a b. (a -> b) -> a -> b
$
    (MPFloat -> BoundsCEDU MPFloat) -> MPFloat -> MPBall -> MPBall
fromApproxWithLipschitz MPFloat -> BoundsCEDU MPFloat
MPFloat.cosCEDU MPFloat
lip MPBall
x
    where
    lip :: MPFloat
lip
        | Integer
i Integer -> Integer -> EqCompareType Integer Integer
forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
0 = Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
1
        | Bool
otherwise = MPBall -> IntervalEndpoint MPBall
forall i. IsInterval i => i -> IntervalEndpoint i
endpointR (MPBall -> IntervalEndpoint MPBall)
-> MPBall -> IntervalEndpoint 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
$ Integer -> MPBall -> MPBall
sinB (Integer
i Integer -> Integer -> SubType Integer Integer
forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
- Integer
1) MPBall
x

-- increasingPrecisionUntilNotImproving :: (MPBall -> MPBall) -> (MPBall -> MPBall)
-- increasingPrecisionUntilNotImproving f x =
--   waitUntilNotImproving $ map aux (precisions xPrec (xPrec*2))
--   where
--   xPrec = getPrecision x
--   precisions p1 p2 = p1 : (precisions p2 (p1 + p2))
--   aux p = f $ setPrecision p x
--   waitUntilNotImproving xx@(x1:_) = aux2 (getAccuracy x1) xx
--   waitUntilNotImproving _ = error "AERN2.MP.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving"
--   aux2 x1AC (x1:x2:rest)
--     | x1AC < x2AC = aux2 x2AC (x2:rest)
--     | otherwise = x1
--     where
--     x2AC = getAccuracy x2
--   aux2 _ _ = error "AERN2.MP.Ball.Elementary: internal error in increasingPrecisionUntilNotImproving"

{- exp, log, power -}

instance CanExp MPBall where
  exp :: MPBall -> ExpType MPBall
exp = (IntervalEndpoint MPBall -> IntervalEndpoint MPBall)
-> (IntervalEndpoint MPBall -> IntervalEndpoint MPBall)
-> MPBall
-> MPBall
forall t.
IsInterval t =>
(IntervalEndpoint t -> IntervalEndpoint t)
-> (IntervalEndpoint t -> IntervalEndpoint t) -> t -> t
intervalFunctionByEndpointsUpDown MPFloat -> MPFloat
IntervalEndpoint MPBall -> IntervalEndpoint MPBall
MPFloat.expDown MPFloat -> MPFloat
IntervalEndpoint MPBall -> IntervalEndpoint MPBall
MPFloat.expUp

instance CanLog MPBall where
  type LogType MPBall = MPBall
  log :: MPBall -> LogType MPBall
log MPBall
x
    | MPBall
x_MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
!>! Integer
1 =
        Precision -> MPBall -> MPBall
forall t. CanSetPrecision t => Precision -> t -> t
setPrecision Precision
p (MPBall -> MPBall) -> MPBall -> MPBall
forall a b. (a -> b) -> a -> b
$ (MPBall -> MPBall) -> (MPBall -> ErrorBound) -> MPBall -> MPBall
forall t.
(IsBall t, HasEqCertainly t t) =>
(t -> t) -> (t -> ErrorBound) -> t -> t
ballFunctionUsingLipschitz MPBall -> MPBall
log_ MPBall -> ErrorBound
forall t2.
(Convertible (DivType Integer t2) ErrorBound, CanDiv Integer t2) =>
t2 -> ErrorBound
logLip MPBall
x_
    | MPBall
x_MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
!>! Integer
0 =
        Precision -> MPBall -> MPBall
forall t. CanSetPrecision t => Precision -> t -> t
setPrecision Precision
p (MPBall -> MPBall) -> MPBall -> MPBall
forall a b. (a -> b) -> a -> b
$ (MPBall -> MPBall) -> MPBall -> MPBall
forall t.
(IsInterval t, CanMinMaxSameType (IntervalEndpoint t),
 HasEqCertainly t t) =>
(t -> t) -> t -> t
intervalFunctionByEndpoints MPBall -> MPBall
log_ MPBall
x_
    | MPBall
x MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
!>! Integer
0 =
        (MPBall -> MPBall) -> MPBall -> MPBall
forall t.
(IsInterval t, CanMinMaxSameType (IntervalEndpoint t),
 HasEqCertainly t t) =>
(t -> t) -> t -> t
intervalFunctionByEndpoints MPBall -> MPBall
log_ MPBall
x
    | Bool
otherwise = LogType MPBall
MPBall
err
    where
    p :: Precision
p = MPBall -> Precision
forall t. HasPrecision t => t -> Precision
getPrecision MPBall
x
    x_ :: MPBall
x_ = MPBall -> MPBall
reducePrecionIfInaccurate MPBall
x
    err :: MPBall
err = [Char] -> MPBall
forall a. HasCallStack => [Char] -> a
error ([Char] -> MPBall) -> [Char] -> MPBall
forall a b. (a -> b) -> a -> b
$ [Char]
"log: argument must be > 0: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ MPBall -> [Char]
forall a. Show a => a -> [Char]
show MPBall
x
    log_ :: MPBall -> MPBall
log_ (MPBall MPFloat
c ErrorBound
e) = MPFloat -> ErrorBound -> MPBall
MPBall MPFloat
lc (ErrorBound
e ErrorBound -> ErrorBound -> AddType ErrorBound ErrorBound
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ (MPFloat -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound MPFloat
le))
      where
      (MPFloat
lc, MPFloat
le) = BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a. BoundsCEDU a -> (a, a)
ceduCentreErr (BoundsCEDU MPFloat -> (MPFloat, MPFloat))
-> BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a b. (a -> b) -> a -> b
$ MPFloat -> BoundsCEDU MPFloat
MPFloat.logCEDU MPFloat
c
    logLip :: t2 -> ErrorBound
logLip t2
y = DivType Integer t2 -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound (DivType Integer t2 -> ErrorBound)
-> DivType Integer t2 -> ErrorBound
forall a b. (a -> b) -> a -> b
$ (Integer
1Integer -> t2 -> DivType Integer t2
forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/t2
y)

instance CanPow MPBall MPBall where
  pow :: MPBall -> MPBall -> PowType MPBall MPBall
pow = MPBall
-> (MPBall -> MPBall -> MPBall)
-> (MPBall -> MPBall)
-> MPBall
-> MPBall
-> MPBall
forall t.
(CanLogSameType t, CanExpSameType t, CanTestInteger t,
 CanTestZero t) =>
t -> (t -> t -> t) -> (t -> t) -> t -> t -> t
powUsingExpLog (Integer -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall Integer
1) MPBall -> MPBall -> MPBall
mulByEndpoints MPBall -> MPBall
forall t. CanRecip t => t -> DivType Integer t
recip

instance CanPow MPBall Dyadic where
  pow :: MPBall -> Dyadic -> PowType MPBall Dyadic
pow MPBall
b Dyadic
e = MPBall -> MPBall -> PowType MPBall MPBall
forall b e. CanPow b e => b -> e -> PowType b e
pow MPBall
b (Dyadic -> MPBall
forall t. CanBeMPBall t => t -> MPBall
mpBall Dyadic
e)

instance CanPow MPBall Rational where
  pow :: MPBall -> Rational -> PowType MPBall Rational
pow MPBall
b Rational
e = MPBall -> MPBall -> PowType MPBall MPBall
forall b e. CanPow b e => b -> e -> PowType b e
pow MPBall
b (Precision -> Rational -> MPBall
forall t. CanBeMPBallP t => Precision -> t -> MPBall
mpBallP (MPBall -> Precision
forall t. HasPrecision t => t -> Precision
getPrecision MPBall
b) Rational
e)

instance CanSqrt MPBall where
  type SqrtType MPBall = MPBall
  sqrt :: MPBall -> SqrtType MPBall
sqrt MPBall
x
    | MPBall
x MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
!>=! Integer
0 = MPBall -> MPBall
aux MPBall
x
    | MPBall
x MPBall -> Integer -> Bool
forall a b. HasOrderCertainlyAsymmetric a b => a -> b -> Bool
?>=? Integer
0 = MPBall -> MPBall
aux (MPBall -> MPBall) -> MPBall -> MPBall
forall a b. (a -> b) -> a -> b
$ Integer -> MPBall -> MinMaxType Integer MPBall
forall t1 t2.
CanMinMaxAsymmetric t1 t2 =>
t1 -> t2 -> MinMaxType t1 t2
max Integer
0 MPBall
x
    | Bool
otherwise = SqrtType MPBall
MPBall
err
    where
    aux :: MPBall -> MPBall
aux =
      (IntervalEndpoint MPBall -> IntervalEndpoint MPBall)
-> (IntervalEndpoint MPBall -> IntervalEndpoint MPBall)
-> MPBall
-> MPBall
forall t.
IsInterval t =>
(IntervalEndpoint t -> IntervalEndpoint t)
-> (IntervalEndpoint t -> IntervalEndpoint t) -> t -> t
intervalFunctionByEndpointsUpDown
        (\ IntervalEndpoint MPBall
e -> MPFloat -> MPFloat
MPFloat.sqrtDown (MPFloat -> MPFloat -> MPFloat
forall a. Ord a => a -> a -> a
P.max (Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
0) MPFloat
IntervalEndpoint MPBall
e))
        (\ IntervalEndpoint MPBall
e -> MPFloat -> MPFloat
MPFloat.sqrtUp (MPFloat -> MPFloat -> MPFloat
forall a. Ord a => a -> a -> a
P.max (Integer -> MPFloat
forall t. CanBeMPFloat t => t -> MPFloat
mpFloat Integer
0) MPFloat
IntervalEndpoint MPBall
e))
    err :: MPBall
err = [Char] -> MPBall
forall a. HasCallStack => [Char] -> a
error ([Char] -> MPBall) -> [Char] -> MPBall
forall a b. (a -> b) -> a -> b
$ [Char]
"sqrt: argument must be >= 0: " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ MPBall -> [Char]
forall a. Show a => a -> [Char]
show MPBall
x

{- generic methods for computing real functions from MPFR-approximations -}

{-|
    Computes a real function @f@ from correctly rounded MPFR-approximations and a number @lip@ which is a
    Lipschitz constant for @f@, i.e. @|f(x) - f(y)| <= lip * |x - y|@ for all @x@,@y@.
-}
fromApproxWithLipschitz ::
    (MPFloat -> MPFloat.BoundsCEDU MPFloat) {-^ @fCEDU@: a version of @f@ on MPFloat returning rigorous bounds -} ->
    MPFloat {-^ @lip@ a Lipschitz constant for @f@, @lip > 0@ -} ->
    (MPBall -> MPBall) {-^ @f@ on MPBall rounding *outwards* -}
fromApproxWithLipschitz :: (MPFloat -> BoundsCEDU MPFloat) -> MPFloat -> MPBall -> MPBall
fromApproxWithLipschitz MPFloat -> BoundsCEDU MPFloat
fCEDU MPFloat
lip _x :: MPBall
_x@(MPBall MPFloat
xc ErrorBound
xe) =
    MPBall -> MPBall
forall t. CanNormalize t => t -> t
normalize (MPBall -> MPBall) -> MPBall -> MPBall
forall a b. (a -> b) -> a -> b
$ MPFloat -> ErrorBound -> MPBall
MPBall MPFloat
fxCP AddType ErrorBound ErrorBound
ErrorBound
err
    where
    (MPFloat
fxC, MPFloat
fxErr) = BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a. BoundsCEDU a -> (a, a)
MPFloat.ceduCentreErr (BoundsCEDU MPFloat -> (MPFloat, MPFloat))
-> BoundsCEDU MPFloat -> (MPFloat, MPFloat)
forall a b. (a -> b) -> a -> b
$ MPFloat -> BoundsCEDU MPFloat
fCEDU MPFloat
xc
    (MPBall MPFloat
fxCP ErrorBound
fxe) =
      Precision -> MPBall -> MPBall
forall t. CanSetPrecision t => Precision -> t -> t
setPrecision (MPFloat -> Precision
forall t. HasPrecision t => t -> Precision
getPrecision MPFloat
xc) (MPBall -> MPBall) -> MPBall -> MPBall
forall a b. (a -> b) -> a -> b
$ -- beware, some MPFloat functions may increase precision, eg sine and cosine
        (MPFloat -> ErrorBound -> MPBall
MPBall MPFloat
fxC (MPFloat -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound MPFloat
fxErr))
    err :: AddType ErrorBound ErrorBound
err = (MPFloat -> ErrorBound
forall t. CanBeErrorBound t => t -> ErrorBound
errorBound MPFloat
lip) ErrorBound -> ErrorBound -> MulType ErrorBound ErrorBound
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* ErrorBound
xe  ErrorBound -> ErrorBound -> AddType ErrorBound ErrorBound
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+  ErrorBound
fxe

$(declForTypes
  [[t| Integer |], [t| Int |], [t| Rational |]]
  (\ b -> [d|

  instance 
    CanPow $b MPBall 
    where
    type PowType $b MPBall = MPBall
    pow x e = pow (mpBallP (getPrecision e) x) e
  |]))