diff --git a/AERN-Real-Interval.cabal b/AERN-Real-Interval.cabal new file mode 100644 --- /dev/null +++ b/AERN-Real-Interval.cabal @@ -0,0 +1,66 @@ +Name: AERN-Real-Interval +Version: 2011.1 +Cabal-Version: >= 1.2 +Build-Type: Simple +License: BSD3 +License-File: LICENCE +Author: Michal Konecny (Aston University) +Copyright: (c) 2011 Michal Konecny, Jan Duracz +Maintainer: mikkonecny@gmail.com +Homepage: http://code.google.com/p/aern/ +Stability: experimental +Category: Data, Math +Synopsis: arbitrary precision real interval arithmetic +Tested-with: GHC ==6.12.3 +Description: + A concrete implementation for the refinement-order type classes from AERN-Real in the form + of interval arithmetic. The arithmetic supports also anti-consistent + intervals (ie @[l,r]@ with @l >= r@) and maintains monotonicity in the + refinement order so that both inner and outer approximations of exact interval + expressions can be safely computed. + . + The package AERN-Real-Double makes it possible to use ordinary + machine 'Double' values as endpoints, although limited to its fixed granularity (ie precision). + . + A package AERN-Real-MPFR will be provided in future to facilitate (via package hmpfr) + the use of the MPFR arbitrary granularity (ie precision) floating point numbers + as endpoints. + +Extra-Source-Files: + CHANGES + +Library + hs-source-dirs: src + ghc-options: -O2 + Build-Depends: + base >= 4 && < 5, + QuickCheck >= 2.1 && < 3, + test-framework >= 0.2 && < 0.4, test-framework-quickcheck2 >= 0.2 && < 0.4, + deepseq >= 1.1 && < 2.0, + AERN-Basics == 2011.1, AERN-Real == 2011.1 + Exposed-modules: + Numeric.AERN.Basics.Interval + Numeric.AERN.RealArithmetic.Interval + Numeric.AERN.RealArithmetic.Interval.Mutable + Numeric.AERN.RealArithmetic.Interval.ElementaryFromBasis + Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps + Numeric.AERN.RealArithmetic.Interval.Mutable.ElementaryFromFieldOps + + Other-modules: + Numeric.AERN.Basics.Interval.Basics + Numeric.AERN.Basics.Interval.Consistency + Numeric.AERN.Basics.Interval.NumericOrder + Numeric.AERN.Basics.Interval.RefinementOrder + Numeric.AERN.Basics.Interval.Mutable + Numeric.AERN.RealArithmetic.Interval.ExactOps + Numeric.AERN.RealArithmetic.Interval.Measures + Numeric.AERN.RealArithmetic.Interval.Conversion + Numeric.AERN.RealArithmetic.Interval.FieldOps + Numeric.AERN.RealArithmetic.Interval.MixedFieldOps + Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps + Numeric.AERN.RealArithmetic.Interval.Mutable.FieldOps + Numeric.AERN.RealArithmetic.Interval.Mutable.MixedFieldOps + Numeric.AERN.RealArithmetic.Interval.SpecialConst + Numeric.AERN.RealArithmetic.Interval.Floating + Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps.Sqrt + diff --git a/CHANGES b/CHANGES new file mode 100644 --- /dev/null +++ b/CHANGES @@ -0,0 +1,3 @@ +2011.1: 6th May 2011 + * initial release of AERN-Real-Interval + diff --git a/LICENCE b/LICENCE new file mode 100644 --- /dev/null +++ b/LICENCE @@ -0,0 +1,30 @@ +Copyright (c) 2010 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 author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE 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 AUTHORS 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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/src/Numeric/AERN/Basics/Interval.hs b/src/Numeric/AERN/Basics/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval.hs @@ -0,0 +1,29 @@ +{-# LANGUAGE TypeFamilies #-} +{-| + Module : Numeric.AERN.Basics.Interval + Description : a minimal interval datatype + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + A minimal interval datatype and its instances. +-} +module Numeric.AERN.Basics.Interval +( + module Numeric.AERN.Basics.Interval.Basics, + module Numeric.AERN.Basics.Interval.Consistency, + module Numeric.AERN.Basics.Interval.NumericOrder, + module Numeric.AERN.Basics.Interval.RefinementOrder, + module Numeric.AERN.Basics.Interval.Mutable +) +where + +import Numeric.AERN.Basics.Interval.Basics +import Numeric.AERN.Basics.Interval.Consistency +import Numeric.AERN.Basics.Interval.NumericOrder +import Numeric.AERN.Basics.Interval.RefinementOrder +import Numeric.AERN.Basics.Interval.Mutable + diff --git a/src/Numeric/AERN/Basics/Interval/Basics.hs b/src/Numeric/AERN/Basics/Interval/Basics.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval/Basics.hs @@ -0,0 +1,95 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-| + Module : Numeric.AERN.Basics.Interval.Basics + Description : interval datatype and its basic instances + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Interval datatype and its basic instances. + + This is a hidden module reexported via its parent. +-} +module Numeric.AERN.Basics.Interval.Basics +( + Interval(..), + getEndpoints, fromEndpoints, mapBothEndpoints, mapEachEndpoint, mapEndpointPair +) +where + +import Prelude hiding (EQ, LT, GT) + +import Numeric.AERN.Basics.ShowInternals +import Numeric.AERN.Basics.PartialOrdering + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import Numeric.AERN.Basics.NumericOrder.OpsImplicitEffort + +import Control.DeepSeq + +{-| + Pairs of endpoints. An end user should not use this type directly + but use the classes of which this is an instance. +-} +data Interval e = + Interval + { + leftEndpoint :: ! e, + rightEndpoint :: ! e + } + +instance (ShowInternals e, NumOrd.PartialComparison e) => (ShowInternals (Interval e)) + where + type ShowInternalsIndicator (Interval e) = ShowInternalsIndicator e + defaultShowIndicator (Interval l r) = defaultShowIndicator l + showInternals indicator (Interval l r) = + case (NumOrd.pCompareEff (NumOrd.pCompareDefaultEffort l) l r) of + Just EQ -> "<" ++ showL ++ ">" + Just LT -> showConsistent + Just LEE -> showConsistent + Just GT -> showAnticonsistent + Just GEE -> showAnticonsistent + _ -> showUnknown + where + showL = showInternals indicator l + showR = showInternals indicator r + showConsistent = "[." ++ showL ++ "," ++ showR ++ "^]" + showAnticonsistent = "[^" ++ showL ++ "," ++ showR ++ ".]" + showUnknown = "[?" ++ showL ++ "," ++ showR ++ "?]" + +instance (ShowInternals e, NumOrd.PartialComparison e) => (Show (Interval e)) + where + show i = showInternals (defaultShowIndicator i) i + +instance (NFData e) => NFData (Interval e) where + rnf (Interval l r) = + rnf l `seq` rnf r `seq` () +-- l `seq` r `seq` () + +-- | Given an argument interval 'i' 'getEndpoints' returns the endpoint pair +-- ('leftEndpoint' 'i','rightEndpoint' 'i'). +getEndpoints :: Interval t -> (t,t) +getEndpoints (Interval l r) = (l, r) + +-- | Constructs an interval from an endpoint pair. +fromEndpoints :: (t,t) -> Interval t +fromEndpoints (l,r) = Interval l r + +mapBothEndpoints f (Interval l r) = Interval (f l) (f r) + +mapEachEndpoint fl fh (Interval l r) = Interval (fl l) (fh r) + +mapEndpointPair f (Interval l r) = + Interval fl fr + where + (fl, fr) = f (l, r) + + + + + diff --git a/src/Numeric/AERN/Basics/Interval/Consistency.hs b/src/Numeric/AERN/Basics/Interval/Consistency.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval/Consistency.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.Basics.Interval.Basics + Description : consistency instances for intervals + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Consistency instances for intervals. + + This is a hidden module reexported via its parent. +-} + +module Numeric.AERN.Basics.Interval.Consistency +( + ConsistentInterval(..), + AntiConsistentInterval(..), + ConsistentOrACInterval(..) +) +where + +import Prelude hiding (LT,EQ) +import Numeric.AERN.Basics.PartialOrdering + +import Numeric.AERN.Basics.Interval.Basics + +import Numeric.AERN.Basics.Consistency + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import Numeric.AERN.Basics.NumericOrder.OpsImplicitEffort + + +import Test.QuickCheck +import Numeric.AERN.Misc.QuickCheck + +import Test.Framework (testGroup, Test) +import Test.Framework.Providers.QuickCheck2 (testProperty) + + +instance (NumOrd.PartialComparison e) => HasConsistency (Interval e) + where + type ConsistencyEffortIndicator (Interval e) = + NumOrd.PartialCompareEffortIndicator e + consistencyDefaultEffort (Interval l r) = + NumOrd.pCompareDefaultEffort l + isConsistentEff effort (Interval l r) = + NumOrd.pLeqEff effort l r + +instance (NumOrd.PartialComparison e) => HasAntiConsistency (Interval e) + where + isAntiConsistentEff effort (Interval l r) = + NumOrd.pLeqEff effort r l + flipConsistency (Interval l r) = Interval r l + +instance HasThinRepresentative (Interval e) + where + getThinRepresentative (Interval l r) = Interval r r + +-- random generation of intervals with no guarantee of consistency: +instance (NumOrd.ArbitraryOrderedTuple e) => Arbitrary (Interval e) + where + arbitrary = + do + (NumOrd.UniformlyOrderedPair (l,r)) <- arbitrary + return $ Interval l r + +{-| type for random generation of consistent intervals -} +data ConsistentInterval e = ConsistentInterval (Interval e) deriving (Show) + +instance (NumOrd.ArbitraryOrderedTuple e) => (Arbitrary (ConsistentInterval e)) + where + arbitrary = + case NumOrd.arbitraryPairRelatedBy LT of + Just gen -> + do + (l,r) <- gen + shouldBeSingleton <- arbitraryBoolRatio 1 10 + case shouldBeSingleton of + True -> return $ ConsistentInterval (Interval l l) + False -> return $ ConsistentInterval (Interval l r) + + +{-| type for random generation of anti-consistent intervals -} +data AntiConsistentInterval e = AntiConsistentInterval (Interval e) deriving (Show) + +instance (NumOrd.ArbitraryOrderedTuple e) => (Arbitrary (AntiConsistentInterval e)) + where + arbitrary = + case NumOrd.arbitraryPairRelatedBy LT of + Just gen -> + do + (l,r) <- gen + shouldBeSingleton <- arbitraryBoolRatio 1 10 + case shouldBeSingleton of + True -> return $ AntiConsistentInterval (Interval l l) + False -> return $ AntiConsistentInterval (Interval r l) + +{-| type for random generation of consistent and anti-consistent intervals + with the same probability -} +data ConsistentOrACInterval e = ConsistentOrACInterval (Interval e) deriving (Show) + +instance (NumOrd.ArbitraryOrderedTuple e) => (Arbitrary (ConsistentOrACInterval e)) + where + arbitrary = + do + consistent <- arbitrary + case consistent of + True -> + do + (ConsistentInterval i) <- arbitrary + return $ ConsistentOrACInterval i + False -> + do + (AntiConsistentInterval i) <- arbitrary + return $ ConsistentOrACInterval i + diff --git a/src/Numeric/AERN/Basics/Interval/Mutable.hs b/src/Numeric/AERN/Basics/Interval/Mutable.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval/Mutable.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE TypeFamilies #-} +{-| + Module : Numeric.AERN.Basics.Interval.Mutable + Description : mutable version of Interval + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Mutable version of Interval. + + This is a private module reexported publicly via its parent. +-} + +module Numeric.AERN.Basics.Interval.Mutable where + +import Numeric.AERN.Basics.Interval.Basics +import Numeric.AERN.Basics.Mutable + +import Numeric.AERN.RealArithmetic.ExactOps + +instance (CanBeMutable e) => CanBeMutable (Interval e) where + data Mutable (Interval e) s = + MInterval { mIntervalLeft :: Mutable e s, mIntervalRight :: Mutable e s } + makeMutable (Interval l r) = + do + lM <- makeMutable l + rM <- makeMutable r + return $ MInterval lM rM + unsafeMakeMutable (Interval l r) = + do + lM <- unsafeMakeMutable l + rM <- unsafeMakeMutable r + return $ MInterval lM rM + writeMutable (MInterval lM rM) (Interval l r) = + do + writeMutable lM l + writeMutable rM r + unsafeWriteMutable (MInterval lM rM) (Interval l r) = + do + unsafeWriteMutable lM l + unsafeWriteMutable rM r + readMutable (MInterval lM rM) = + do + l <- readMutable lM + r <- readMutable rM + return $ Interval l r + unsafeReadMutable (MInterval lM rM) = + do + l <- unsafeReadMutable lM + r <- unsafeReadMutable rM + return $ Interval l r diff --git a/src/Numeric/AERN/Basics/Interval/NumericOrder.hs b/src/Numeric/AERN/Basics/Interval/NumericOrder.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval/NumericOrder.hs @@ -0,0 +1,292 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-| + Module : Numeric.AERN.Basics.Interval.NumericOrder + Description : interval instances of numeric-ordered structures + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Interval instances of numeric-ordered structures. + + This is a hidden module reexported via its parent. +-} +module Numeric.AERN.Basics.Interval.NumericOrder where + +import Prelude hiding (EQ, LT, GT) + +import Numeric.AERN.Basics.Consistency +import Numeric.AERN.Basics.PartialOrdering + +import Numeric.AERN.Basics.Interval.Basics +import Numeric.AERN.Basics.Interval.Consistency +import Numeric.AERN.Basics.Interval.Mutable + +import Numeric.AERN.Misc.List + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd + +import Test.QuickCheck + +import Data.Maybe + +instance + (NumOrd.PartialComparison e) => + (NumOrd.PartialComparison (Interval e)) + where + type NumOrd.PartialCompareEffortIndicator (Interval e) = + NumOrd.PartialCompareEffortIndicator e + pCompareDefaultEffort (Interval l r) = NumOrd.pCompareDefaultEffort l + pCompareEff effort i1@(Interval l1 r1) i2@(Interval l2 r2) = + case (isConsistentEff effort i1, isConsistentEff effort i2) of + (Just True, Just True) -> + case (r1 `leq` l2, r1 `less` l2, r2 `leq` l1, r2 `less` l1) of + (Just True, _, Just True, _) -> Just EQ + (_, Just True, _, _) -> Just LT + (Just True, _, _, _) -> Just LEE + (_, _, _, Just True) -> Just GT + (_, _, Just True, _) -> Just GEE + _ -> Nothing + _ -> Nothing + where + leq = NumOrd.pLeqEff effort + less = NumOrd.pLessEff effort + + +instance + (NumOrd.RoundedLatticeEffort e) => + (NumOrd.OuterRoundedLatticeEffort (Interval e)) + where + type NumOrd.MinmaxOuterEffortIndicator (Interval e) = NumOrd.MinmaxEffortIndicator e + minmaxOuterDefaultEffort (Interval l r) = NumOrd.minmaxDefaultEffort l + +instance + (NumOrd.RoundedLattice e) => + (NumOrd.OuterRoundedLattice (Interval e)) + where + minOutEff effort (Interval l1 r1) (Interval l2 r2) = + Interval (NumOrd.minDnEff effort l1 l2) (NumOrd.minUpEff effort r1 r2) + maxOutEff effort (Interval l1 r1) (Interval l2 r2) = + Interval (NumOrd.maxDnEff effort l1 l2) (NumOrd.maxUpEff effort r1 r2) + +instance + (NumOrd.RoundedLatticeEffort e) => + (NumOrd.InnerRoundedLatticeEffort (Interval e)) + where + type NumOrd.MinmaxInnerEffortIndicator (Interval e) = NumOrd.MinmaxEffortIndicator e + minmaxInnerDefaultEffort (Interval l r) = NumOrd.minmaxDefaultEffort l +instance + (NumOrd.RoundedLattice e) => + (NumOrd.InnerRoundedLattice (Interval e)) + where + minInEff effort (Interval l1 r1) (Interval l2 r2) = + Interval (NumOrd.minUpEff effort l1 l2) (NumOrd.minDnEff effort r1 r2) + maxInEff effort (Interval l1 r1) (Interval l2 r2) = + Interval (NumOrd.maxUpEff effort l1 l2) (NumOrd.maxDnEff effort r1 r2) + +instance + (NumOrd.RoundedLattice e) => + (NumOrd.RefinementRoundedLattice (Interval e)) + +instance + (NumOrd.RoundedLatticeInPlace e) => + (NumOrd.OuterRoundedLatticeInPlace (Interval e)) + where + minOutInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.minDnInPlaceEff effort resLM l1M l2M + NumOrd.minUpInPlaceEff effort resRM r1M r2M + maxOutInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxDnInPlaceEff effort resLM l1M l2M + NumOrd.maxUpInPlaceEff effort resRM r1M r2M + + +instance + (NumOrd.RoundedLatticeInPlace e) => + (NumOrd.InnerRoundedLatticeInPlace (Interval e)) + where + minInInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.minUpInPlaceEff effort resLM l1M l2M + NumOrd.minDnInPlaceEff effort resRM r1M r2M + maxInInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxUpInPlaceEff effort resLM l1M l2M + NumOrd.maxDnInPlaceEff effort resRM r1M r2M + +instance + (NumOrd.RoundedLatticeInPlace e) => + (NumOrd.RefinementRoundedLatticeInPlace (Interval e)) + +instance (NumOrd.HasLeast e) => (NumOrd.HasLeast (Interval e)) + where + least = Interval NumOrd.least NumOrd.least + +instance (NumOrd.HasGreatest e) => (NumOrd.HasGreatest (Interval e)) + where + greatest = Interval NumOrd.greatest NumOrd.greatest + +instance (NumOrd.HasExtrema e) => (NumOrd.HasExtrema (Interval e)) + +instance (NumOrd.ArbitraryOrderedTuple e) => NumOrd.ArbitraryOrderedTuple (Interval e) where + type NumOrd.Area (Interval e) = NumOrd.Area e + areaWhole (Interval l r) = NumOrd.areaWhole l + arbitraryTupleInAreaRelatedBy area = + arbitraryIntervalTupleInAreaNumericallyRelatedBy (Just area) + arbitraryTupleRelatedBy = + arbitraryIntervalTupleInAreaNumericallyRelatedBy Nothing + +arbitraryIntervalTupleInAreaNumericallyRelatedBy maybeArea indices constraints = + case endpointGens of + [] -> Nothing + _ -> Just $ + do + gen <- elements endpointGens + endpointTuple <- gen + return $ endpointsToIntervals endpointTuple + where + endpointGens = + case maybeArea of + (Just area) -> + catMaybes $ + map (NumOrd.arbitraryTupleInAreaRelatedBy area endpointIndices) + endpointConstraintsVersions + Nothing -> + catMaybes $ + map (NumOrd.arbitraryTupleRelatedBy endpointIndices) + endpointConstraintsVersions + endpointIndices = + concat $ map (\ix -> [(ix,-1), (ix,1)]) indices + endpointsToIntervals [] = [] + endpointsToIntervals (l : r : rest) = + (Interval l r) : (endpointsToIntervals rest) + endpointConstraintsVersions = +-- unsafePrintReturn +-- ("arbitraryIntervalTupleRelatedBy:" +-- ++ "\n indices = " ++ show indices +-- ++ "\n constraints = " ++ show constraints +-- ++ "\n endpointIndices = " ++ show endpointIndices +-- ++ "\n endpointConstraintsVersions = " +-- ) $ + map concat $ combinations $ map intervalConstraintsToEndpointConstraints constraints + intervalConstraintsToEndpointConstraints :: + ((ix, ix), [PartialOrdering]) -> [[(((ix,Int), (ix,Int)), [PartialOrdering])]] + intervalConstraintsToEndpointConstraints ((ix1, ix2),rels) = + concat $ map forEachRel rels + where + endpoints1Comparable = [(((ix1,-1),(ix1, 1)), [EQ,LT,LEE,GT,GEE])] + endpoints2Comparable = [(((ix2,-1),(ix2, 1)), [EQ,LT,LEE,GT,GEE])] + endpointsComparable = endpoints1Comparable ++ endpoints2Comparable + forEachRel EQ = -- both must be thin and equal + [[(((ix1,-1),(ix1,1)), [EQ]), (((ix1,1),(ix2,1)), [EQ]), (((ix2,-1),(ix2,1)), [EQ])]] + ++ -- some cases where the order is not decided: + -- or the interval ix1 is indide ix2 + ix1 does not coincide with ix2's endpoint + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,-1)), [EQ, GT, GEE])] ++ + [(((ix1,1),(ix2,-1)), [GT, GEE])] ++ + [(((ix1,-1),(ix2,1)), [LT, LEE])] ++ + [(((ix1,1),(ix2,1)), [EQ, LT, LEE])] + ] + ++ + -- or the interval ix2 is indide ix1 + ix2 does not coincide with ix1's endpoint + [ + endpointsComparable ++ + [(((ix2,-1),(ix1,-1)), [EQ, GT, GEE])] ++ + [(((ix2,1),(ix1,-1)), [GT, GEE])] ++ + [(((ix2,-1),(ix1,1)), [LT, LEE])] ++ + [(((ix2,1),(ix1,1)), [EQ, LT, LEE])] + ] + forEachRel LT = -- both endpoints of ix1 must be less than both endpoints of ix2 + [ + endpointsComparable ++ + [(((ix1,side1),(ix2,side2)), [LT]) | side1 <- [-1,1], side2 <- [-1,1]] + ] + ++ -- some undecidable cases: + -- or the interval ix1 overlaps ix2 and ix1 is slightly to the left of ix2 + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,-1)), [EQ, LT, LEE]), + (((ix1,1),(ix2,1)), [EQ,LT,LEE]), + (((ix1,-1),(ix2,1)), [LT, LEE]), + (((ix1,1),(ix2,-1)), [GT, GEE])] + ] + forEachRel GT = -- both endpoints of ix1 must be greater than both endpoints of ix2 + [ + endpointsComparable ++ + [(((ix1,side1),(ix2,side2)), [GT]) | side1 <- [-1,1], side2 <- [-1,1]] + ] + ++ -- some undecidable cases: + -- or the interval ix1 overlaps ix2 and ix1 is slightly to the right of ix2 + [ + endpointsComparable ++ + [(((ix2,-1),(ix1,-1)), [EQ, LT, LEE]), + (((ix2,1),(ix1,1)), [EQ,LT,LEE]), + (((ix2,-1),(ix1,1)), [LT, LEE]), + (((ix2,1),(ix1,-1)), [GT, GEE])] + ] + forEachRel LEE = + [ + endpointsComparable ++ + [(((ix1,1),(ix2,-1)), [EQ]), + (((ix1,1),(ix2,1)), [LT,LEE,EQ]), (((ix1,-1),(ix2,-1)), [LT,LEE,EQ]), + (((ix1,-1),(ix2,1)), [LT,LEE])] + ] + forEachRel GEE = + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,1)), [EQ]), + (((ix1,1),(ix2,1)), [GT,GEE,EQ]), (((ix1,-1),(ix2,-1)), [GT,GEE,EQ]), + (((ix1,1),(ix2,-1)), [GT,GEE])] + ] + forEachRel NC = + -- either some pair of endpoints is NC: + [ endpointsComparable ++ [(((ix1,side1), (ix2, side2)),[NC])] + | side1 <- [-1,1], side2 <- [-1,1] + ] +-- ++ +-- -- or the interval ix1 is indide ix2 + ix1 does not coincide with ix2's endpoint +-- [ +-- endpointsComparable ++ +-- [(((ix1,-1),(ix2,-1)), [EQ, GT, GEE])] ++ +-- [(((ix1,1),(ix2,-1)), [GT, GEE])] ++ +-- [(((ix1,-1),(ix2,1)), [LT, LEE])] ++ +-- [(((ix1,1),(ix2,1)), [EQ, LT, LEE])] +-- ] +-- ++ +-- -- or the interval ix2 is indide ix1 + ix2 does not coincide with ix1's endpoint +-- [ +-- endpointsComparable ++ +-- [(((ix2,-1),(ix1,-1)), [EQ, GT, GEE])] ++ +-- [(((ix2,1),(ix1,-1)), [GT, GEE])] ++ +-- [(((ix2,-1),(ix1,1)), [LT, LEE])] ++ +-- [(((ix2,1),(ix1,1)), [EQ, LT, LEE])] +-- ] +-- ++ +-- -- or the interval ix1 overlaps ix2 and ix1 is slightly to the left of ix2 +-- [ +-- endpointsComparable ++ +-- [(((ix1,-1),(ix2,-1)), [EQ, LT, LEE]), +-- (((ix1,1),(ix2,1)), [EQ,LT,LEE]), +-- (((ix1,-1),(ix2,1)), [LT, LEE]), +-- (((ix1,1),(ix2,-1)), [GT, GEE])] +-- ] +-- ++ +-- -- or the interval ix1 overlaps ix2 and ix1 is slightly to the right of ix2 +-- [ +-- endpointsComparable ++ +-- [(((ix2,-1),(ix1,-1)), [EQ, LT, LEE]), +-- (((ix2,1),(ix1,1)), [EQ,LT,LEE]), +-- (((ix2,-1),(ix1,1)), [LT, LEE]), +-- (((ix2,1),(ix1,-1)), [GT, GEE])] +-- ] + diff --git a/src/Numeric/AERN/Basics/Interval/RefinementOrder.hs b/src/Numeric/AERN/Basics/Interval/RefinementOrder.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/Basics/Interval/RefinementOrder.hs @@ -0,0 +1,341 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-| + Module : Numeric.AERN.Basics.Interval.RefinementOrder + Description : interval instances of refinement-ordered structures + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Interval instances of refinement-ordered structures. + + This is a hidden module reexported via its parent. +-} +module Numeric.AERN.Basics.Interval.RefinementOrder where + +import Prelude hiding (EQ, LT, GT) + +import Numeric.AERN.Basics.Effort +import Numeric.AERN.Basics.PartialOrdering + +import Numeric.AERN.Basics.Interval.Basics +import Numeric.AERN.Basics.Interval.Mutable +import Numeric.AERN.Basics.Interval.NumericOrder + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd + +import Numeric.AERN.Basics.Mutable + +import Numeric.AERN.Misc.List + +import Test.QuickCheck + +import Data.Maybe + +instance + (NumOrd.PartialComparison e) => + (RefOrd.PartialComparison (Interval e)) + where + type RefOrd.PartialCompareEffortIndicator (Interval e) = + NumOrd.PartialCompareEffortIndicator e + pCompareDefaultEffort (Interval l r) = + NumOrd.pCompareDefaultEffort l + pCompareEff effort (Interval l1 r1) (Interval l2 r2) = + case (c l1 l2, c r1 r2) of + (Just EQ, Just EQ) -> Just EQ + (Just LT, Just GT) -> Just LT + (Just LT, Just EQ) -> Just LT + (Just EQ, Just GT) -> Just LT + (Just GT, Just LT) -> Just GT + (Just GT, Just EQ) -> Just GT + (Just EQ, Just LT) -> Just GT + (Just _, Just _) -> Just NC + _ -> Nothing + where + c = NumOrd.pCompareEff effort + + +instance (NumOrd.HasExtrema e) => (RefOrd.HasTop (Interval e)) + where + top = Interval NumOrd.greatest NumOrd.least + +instance (NumOrd.HasExtrema e) => (RefOrd.HasBottom (Interval e)) + where + bottom = Interval NumOrd.least NumOrd.greatest + +instance (NumOrd.HasExtrema e) => (RefOrd.HasExtrema (Interval e)) + +instance (NumOrd.RoundedLatticeEffort e, NumOrd.PartialComparison e) + => RefOrd.OuterRoundedBasisEffort (Interval e) + where + type RefOrd.PartialJoinOutEffortIndicator (Interval e) = + (NumOrd.MinmaxEffortIndicator e, NumOrd.PartialCompareEffortIndicator e) + partialJoinOutDefaultEffort (Interval l r) = + (NumOrd.minmaxDefaultEffort l, NumOrd.pCompareDefaultEffort l) + +instance + (NumOrd.RoundedLattice e, NumOrd.PartialComparison e) => + RefOrd.OuterRoundedBasis (Interval e) + where + partialJoinOutEff (effortMinmax, effortComp) (Interval l1 r1) (Interval l2 r2) = + case l <=? r of + Just True -> Just $ Interval l r + _ -> Nothing + where + (<=?) = NumOrd.pLeqEff effortComp + l = NumOrd.maxDnEff effortMinmax l1 l2 + r = NumOrd.minUpEff effortMinmax r1 r2 + +instance (NumOrd.RoundedLatticeEffort e, NumOrd.PartialComparison e) + => RefOrd.InnerRoundedBasisEffort (Interval e) + where + type RefOrd.PartialJoinInEffortIndicator (Interval e) = + (NumOrd.MinmaxEffortIndicator e, NumOrd.PartialCompareEffortIndicator e) + partialJoinInDefaultEffort (Interval l r) = + (NumOrd.minmaxDefaultEffort l, NumOrd.pCompareDefaultEffort l) + +instance + (NumOrd.RoundedLattice e, NumOrd.PartialComparison e) => + RefOrd.InnerRoundedBasis (Interval e) + where + partialJoinInEff (effortMinmax, effortComp) (Interval l1 r1) (Interval l2 r2) = + case l <=? r of + Just True -> Just $ Interval l r + _ -> Nothing + where + (<=?) = NumOrd.pLeqEff effortComp + l = NumOrd.maxUpEff effortMinmax l1 l2 + r = NumOrd.minDnEff effortMinmax r1 r2 + +instance + (NumOrd.RoundedLattice e, NumOrd.PartialComparison e) + => + (RefOrd.RoundedBasis (Interval e)) + +instance + (NumOrd.RoundedLatticeInPlace e, NumOrd.PartialComparison e) => + (RefOrd.OuterRoundedBasisInPlace (Interval e)) + where + partialJoinOutInPlaceEff (effortMinmax, effortComp) + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxDnInPlaceEff effortMinmax resLM l1M l2M + NumOrd.minUpInPlaceEff effortMinmax resRM r1M r2M + l <- unsafeReadMutable resLM + r <- unsafeReadMutable resRM + let (<=?) = NumOrd.pLeqEff effortComp + case l <=? r of + Just True -> return True + _ -> return False + +instance + (NumOrd.RoundedLatticeInPlace e, NumOrd.PartialComparison e) => + (RefOrd.InnerRoundedBasisInPlace (Interval e)) + where + partialJoinInInPlaceEff (effortMinmax, effortComp) + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxUpInPlaceEff effortMinmax resLM l1M l2M + NumOrd.minDnInPlaceEff effortMinmax resRM r1M r2M + l <- unsafeReadMutable resLM + r <- unsafeReadMutable resRM + let (<=?) = NumOrd.pLeqEff effortComp + case l <=? r of + Just True -> return True + _ -> return False + +instance + (NumOrd.RoundedLatticeInPlace e, NumOrd.PartialComparison e) + => + (RefOrd.RoundedBasisInPlace (Interval e)) + + +instance + (NumOrd.RoundedLatticeEffort e) => + (RefOrd.OuterRoundedLatticeEffort (Interval e)) + where + type RefOrd.JoinMeetOutEffortIndicator (Interval e) = + NumOrd.MinmaxEffortIndicator e + joinmeetOutDefaultEffort (Interval l r) = + NumOrd.minmaxDefaultEffort l + +instance + (NumOrd.RoundedLattice e) => + (RefOrd.OuterRoundedLattice (Interval e)) + where + joinOutEff effort (Interval l1 r1) (Interval l2 r2) = + Interval l r + where + l = NumOrd.maxDnEff effort l1 l2 + r = NumOrd.minUpEff effort r1 r2 + meetOutEff effort (Interval l1 r1) (Interval l2 r2) = + Interval l r + where + l = NumOrd.minDnEff effort l1 l2 + r = NumOrd.maxUpEff effort r1 r2 + +instance + (NumOrd.RoundedLatticeEffort e) => + (RefOrd.InnerRoundedLatticeEffort (Interval e)) + where + type RefOrd.JoinMeetInEffortIndicator (Interval e) = + NumOrd.MinmaxEffortIndicator e + joinmeetInDefaultEffort (Interval l r) = + NumOrd.minmaxDefaultEffort l + +instance + (NumOrd.RoundedLattice e) => + (RefOrd.InnerRoundedLattice (Interval e)) + where + joinInEff effort (Interval l1 r1) (Interval l2 r2) = + Interval l r + where + l = NumOrd.maxUpEff effort l1 l2 + r = NumOrd.minDnEff effort r1 r2 + meetInEff effort (Interval l1 r1) (Interval l2 r2) = + Interval l r + where + l = NumOrd.minUpEff effort l1 l2 + r = NumOrd.maxDnEff effort r1 r2 + +instance + (NumOrd.RoundedLattice e, NumOrd.PartialComparison e) => + (RefOrd.RoundedLattice (Interval e)) + +instance + (NumOrd.RoundedLatticeInPlace e) => + (RefOrd.OuterRoundedLatticeInPlace (Interval e)) + where + joinOutInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxDnInPlaceEff effort resLM l1M l2M + NumOrd.minUpInPlaceEff effort resRM r1M r2M + meetOutInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.minDnInPlaceEff effort resLM l1M l2M + NumOrd.maxUpInPlaceEff effort resRM r1M r2M + +instance + (NumOrd.RoundedLatticeInPlace e) => + (RefOrd.InnerRoundedLatticeInPlace (Interval e)) + where + joinInInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.maxUpInPlaceEff effort resLM l1M l2M + NumOrd.minDnInPlaceEff effort resRM r1M r2M + meetInInPlaceEff effort + (MInterval resLM resRM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + NumOrd.minUpInPlaceEff effort resLM l1M l2M + NumOrd.maxDnInPlaceEff effort resRM r1M r2M + +instance (NumOrd.ArbitraryOrderedTuple e) => RefOrd.ArbitraryOrderedTuple (Interval e) where + type RefOrd.Area (Interval e) = NumOrd.Area e + areaWhole (Interval l r) = NumOrd.areaWhole l + arbitraryTupleInAreaRelatedBy area = + arbitraryIntervalTupleInAreaRefinementRelatedBy (Just area) + arbitraryTupleRelatedBy = + arbitraryIntervalTupleInAreaRefinementRelatedBy Nothing + +arbitraryIntervalTupleInAreaRefinementRelatedBy maybeArea indices thinIndices constraints = + case endpointGens of + [] -> Nothing + _ -> Just $ + do + gen <- elements endpointGens + endpointTuple <- gen + return $ endpointsToIntervals endpointTuple + where + endpointGens = + case maybeArea of + (Just area) -> + catMaybes $ + map (NumOrd.arbitraryTupleInAreaRelatedBy area endpointIndices) + endpointConstraintsVersions + Nothing -> + catMaybes $ + map (NumOrd.arbitraryTupleRelatedBy endpointIndices) + endpointConstraintsVersions + endpointIndices = + concat $ map (\ix -> [(ix,-1), (ix,1)]) indices + endpointsToIntervals [] = [] + endpointsToIntervals (l : r : rest) = + (Interval l r) : (endpointsToIntervals rest) + endpointConstraintsVersions = +-- unsafePrintReturn +-- ("arbitraryIntervalTupleRelatedBy:" +-- ++ "\n indices = " ++ show indices +-- ++ "\n constraints = " ++ show constraints +-- ++ "\n endpointIndices = " ++ show endpointIndices +-- ++ "\n endpointConstraintsVersions = " +-- ) $ + map ((++ thinnessConstraints) . concat) $ + combinations $ map intervalConstraintsToEndpointConstraints constraints + thinnessConstraints = map (\ix -> (((ix,-1),(ix,1)),[EQ])) thinIndices + intervalConstraintsToEndpointConstraints :: + ((ix, ix), [PartialOrdering]) -> [[(((ix,Int), (ix,Int)), [PartialOrdering])]] + intervalConstraintsToEndpointConstraints ((ix1, ix2),rels) = + concat $ map forEachRel rels + where + endpoints1Comparable = [(((ix1,-1),(ix1, 1)), [EQ,LT,LEE,GT,GEE])] + endpoints2Comparable = [(((ix2,-1),(ix2, 1)), [EQ,LT,LEE,GT,GEE])] + endpointsComparable = endpoints1Comparable ++ endpoints2Comparable + forEachRel EQ = -- both endpoints agree + [[(((ix1,-1),(ix2,-1)), [EQ]), (((ix1,1),(ix2,1)), [EQ])]] + forEachRel GT = + -- the interval ix1 is indide ix2, but not equal + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,-1)), [GT, GEE])] ++ + [(((ix1,1),(ix2,1)), [EQ, LT, LEE])] + ] + ++ + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,-1)), [EQ, GT, GEE])] ++ + [(((ix1,1),(ix2,1)), [LT, LEE])] + ] + forEachRel LT = + -- the interval ix2 is indide ix1, but not equal + [ + endpointsComparable ++ + [(((ix2,-1),(ix1,-1)), [GT, GEE])] ++ + [(((ix2,1),(ix1,1)), [EQ, LT, LEE])] + ] + ++ + [ + endpointsComparable ++ + [(((ix2,-1),(ix1,-1)), [EQ, GT, GEE])] ++ + [(((ix2,1),(ix1,1)), [LT, LEE])] + ] + forEachRel NC = + -- either some pair of endpoints is NC: + [ endpointsComparable ++ [(((ix1,side1), (ix2, side2)),[NC])] + | side1 <- [-1,1], side2 <- [-1,1] + ] + ++ + -- or the interval ix1 is to the left of ix2 + [ + endpointsComparable ++ + [(((ix1,-1),(ix2,-1)), [LT, LEE]), + (((ix1,1),(ix2,1)), [LT,LEE])] + ] + ++ + -- or the interval ix2 is to the left of ix1 + [ + endpointsComparable ++ + [(((ix2,-1),(ix1,-1)), [LT, LEE]), + (((ix2,1),(ix1,1)), [LT,LEE])] + ] + forEachRel _ = [] + + + diff --git a/src/Numeric/AERN/RealArithmetic/Interval.hs b/src/Numeric/AERN/RealArithmetic/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval.hs @@ -0,0 +1,34 @@ +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval + Description : instances of arithmetic classes for Intervals + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Instances of arithmetic classes for Intervals. +-} + +module Numeric.AERN.RealArithmetic.Interval +() +where + +import Numeric.AERN.RealArithmetic.Interval.ExactOps +import Numeric.AERN.RealArithmetic.Interval.Measures +import Numeric.AERN.RealArithmetic.Interval.Conversion +import Numeric.AERN.RealArithmetic.Interval.FieldOps +import Numeric.AERN.RealArithmetic.Interval.MixedFieldOps +import Numeric.AERN.RealArithmetic.Interval.SpecialConst +import Numeric.AERN.RealArithmetic.Interval.Floating + +import Numeric.AERN.Basics.Exception + +import Numeric.AERN.Basics.Interval + +instance (HasLegalValues e) => HasLegalValues (Interval e) where + isLegal (Interval l r) = isLegal l && isLegal r + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Conversion.hs b/src/Numeric/AERN/RealArithmetic/Interval/Conversion.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Conversion.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE ImplicitParams #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Conversion + Description : conversions between intervals and standard numeric types + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Conversion between intervals and standard numeric types. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.Conversion +() +where + +import Numeric.AERN.RealArithmetic.RefinementOrderRounding + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn + +import Numeric.AERN.Basics.Exception +import Numeric.AERN.Basics.Interval + +import Control.Exception + +instance (ArithUpDn.Convertible t e, Show t) => Convertible t (Interval e) where + type ConvertEffortIndicator t (Interval e) = + ArithUpDn.ConvertEffortIndicator t e + convertDefaultEffort i (Interval l r) = ArithUpDn.convertDefaultEffort i l + convertInEff effort x = + Interval xUp xDn + where + xUp = convertUpEffException effort x + xDn = convertDnEffException effort x + convertOutEff effort x = + Interval xDn xUp + where + xUp = convertUpEffException effort x + xDn = convertDnEffException effort x + +convertUpEffException effort x = + case ArithUpDn.convertUpEff effort x of + Just xUp -> xUp + _ -> throw $ AERNException $ + "failed to convert to interval: x = " ++ show x + +convertDnEffException effort x = + case ArithUpDn.convertDnEff effort x of + Just xDn -> xDn + _ -> throw $ AERNException $ + "failed to convert to interval: x = " ++ show x + +instance (ArithUpDn.Convertible e t) => + ArithUpDn.Convertible (Interval e) t where + type ArithUpDn.ConvertEffortIndicator (Interval e) t = + ArithUpDn.ConvertEffortIndicator e t + convertDefaultEffort (Interval l r) i = ArithUpDn.convertDefaultEffort l i + convertUpEff effort (Interval l r) = ArithUpDn.convertUpEff effort r + convertDnEff effort (Interval l r) = ArithUpDn.convertDnEff effort l + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromBasis.hs b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromBasis.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromBasis.hs @@ -0,0 +1,50 @@ +{-# LANGUAGE TypeFamilies #-} +--{-# LANGUAGE FlexibleContexts #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.ElementaryFromBasis + Description : elementary operations using basis-level operations + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Elementary operations using basis-level operations. +-} + +module Numeric.AERN.RealArithmetic.Interval.ElementaryFromBasis where + +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn + +import Numeric.AERN.Basics.Interval + +instance (ArithUpDn.RoundedExponentiationEffort e) => + (ArithInOut.RoundedExponentiationEffort (Interval e)) + where + type ArithInOut.ExpEffortIndicator (Interval e) = ArithUpDn.ExpEffortIndicator e + expDefaultEffort (Interval l r) = ArithUpDn.expDefaultEffort l + +instance (ArithUpDn.RoundedExponentiation e) => + (ArithInOut.RoundedExponentiation (Interval e)) + where + expInEff effort (Interval l r) = + Interval (ArithUpDn.expUpEff effort l) (ArithUpDn.expDnEff effort r) + expOutEff effort (Interval l r) = + Interval (ArithUpDn.expDnEff effort l) (ArithUpDn.expUpEff effort r) + +instance (ArithUpDn.RoundedSquareRootEffort e) => + (ArithInOut.RoundedSquareRootEffort (Interval e)) + where + type ArithInOut.SqrtEffortIndicator (Interval e) = ArithUpDn.SqrtEffortIndicator e + sqrtDefaultEffort (Interval l r) = ArithUpDn.sqrtDefaultEffort l + +instance (ArithUpDn.RoundedSquareRoot e) => + (ArithInOut.RoundedSquareRoot (Interval e)) + where + sqrtInEff effort (Interval l r) = + Interval (ArithUpDn.sqrtUpEff effort l) (ArithUpDn.sqrtDnEff effort r) + sqrtOutEff effort (Interval l r) = + Interval (ArithUpDn.sqrtDnEff effort l) (ArithUpDn.sqrtUpEff effort r) + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps.hs @@ -0,0 +1,287 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps + Description : elementary operations using generic direct implementation + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Elementary operations using generic direct implementation. +-} + +module Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps +(expDefaultEffortWithIters, expOutIters, expInIters, + sqrtDefaultEffortWithIters, sqrtOutIters, sqrtInIters) +where + +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.ElementaryFromFieldOps.Exponentiation + +import Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps.Sqrt + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.OpsImplicitEffort +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd +import Numeric.AERN.Basics.RefinementOrder.OpsImplicitEffort + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval + +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Consistency +import Numeric.AERN.Basics.Effort + +instance + (ArithInOut.RoundedMixedField (Interval e) Int, + ArithInOut.RoundedField (Interval e), + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e), + NumOrd.PartialComparison e, + RefOrd.OuterRoundedLatticeEffort (Interval e)) + => + (ArithInOut.RoundedExponentiationEffort (Interval e)) + where + type ArithInOut.ExpEffortIndicator (Interval e) = + ((ArithInOut.FieldOpsEffortIndicator (Interval e), + ArithInOut.MixedFieldOpsEffortIndicator (Interval e) Int) + , + Int1To10 + , + ((RefOrd.JoinMeetOutEffortIndicator (Interval e), + NumOrd.PartialCompareEffortIndicator e), + (ArithUpDn.ConvertEffortIndicator (Interval e) Int, + ArithInOut.ConvertEffortIndicator Double (Interval e))) + ) + expDefaultEffort i@(Interval l r) = + ((ArithInOut.fieldOpsDefaultEffort i, + ArithInOut.mixedFieldOpsDefaultEffort i sampleI) + , + Int1To10 10 + , + ((RefOrd.joinmeetOutDefaultEffort i, + NumOrd.pCompareDefaultEffort l), + (ArithUpDn.convertDefaultEffort i sampleI, + ArithInOut.convertDefaultEffort sampleD i)) + ) + where + sampleI = 1 :: Int + sampleD = 1 :: Double + +expDefaultEffortWithIters :: + (NumOrd.PartialComparison e, + RefOrd.OuterRoundedLatticeEffort (Interval e), + ArithInOut.RoundedFieldEffort (Interval e), + ArithInOut.RoundedMixedFieldEffort (Interval e) Int, + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e)) + => + (Interval e) -> + Int -> + ArithInOut.ExpEffortIndicator (Interval e) +expDefaultEffortWithIters i@(Interval l r) n = + ((ArithInOut.fieldOpsDefaultEffort i, + ArithInOut.mixedFieldOpsDefaultEffort i sampleI) + , + Int1To10 n + , + ((RefOrd.joinmeetOutDefaultEffort i, + NumOrd.pCompareDefaultEffort l), + (ArithUpDn.convertDefaultEffort i sampleI, + ArithInOut.convertDefaultEffort sampleD i)) + ) + where + sampleI = 1 :: Int + sampleD = 1 :: Double + + +instance + (ArithInOut.RoundedMixedField (Interval e) Int, + ArithInOut.RoundedField (Interval e), + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e), + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + RefOrd.OuterRoundedLattice (Interval e)) + => + (ArithInOut.RoundedExponentiation (Interval e)) + where + expOutEff + ((effortField, effortMixedField), + (Int1To10 effortTaylor), + ((effortMeet, effortComp), effortConv)) + (Interval l r) = Interval resL resR + where + Interval resL _ = + expOutThinArg + effortField effortMixedField + effortMeet effortComp effortComp effortConv + effortTaylor + (Interval l l) + Interval _ resR = + expOutThinArg + effortField effortMixedField + effortMeet effortComp effortComp effortConv + effortTaylor + (Interval r r) + expInEff + ((effortField, effortMixedField), + (Int1To10 effortTaylor), + ((effortMeet, effortComp), effortConv)) + (Interval l r) = Interval resL resR + where + Interval _ resL = + expOutThinArg + effortField effortMixedField + effortMeet effortComp effortComp effortConv + effortTaylor + (Interval l l) + Interval resR _ = + expOutThinArg + effortField effortMixedField + effortMeet effortComp effortComp effortConv + effortTaylor + (Interval r r) + +expOutIters, expInIters :: + (ArithInOut.RoundedMixedField (Interval e) Int, + ArithInOut.RoundedField (Interval e), + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e), + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + RefOrd.OuterRoundedLattice (Interval e)) + => + Int -> (Interval e) -> (Interval e) +expOutIters n i = ArithInOut.expOutEff (expDefaultEffortWithIters i n) i +expInIters n i = ArithInOut.expInEff (expDefaultEffortWithIters i n) i + +instance + (ArithUpDn.RoundedMixedFieldEffort e Int, + ArithUpDn.RoundedFieldEffort e, + ArithUpDn.Convertible e Double, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) + => + (ArithInOut.RoundedSquareRootEffort (Interval e)) + where + type ArithInOut.SqrtEffortIndicator (Interval e) = + ((ArithUpDn.FieldOpsEffortIndicator e, + ArithUpDn.MixedFieldOpsEffortIndicator e Int) + , + Int1To10 + , + ((NumOrd.MinmaxEffortIndicator e, NumOrd.PartialCompareEffortIndicator e), + ArithUpDn.ConvertEffortIndicator e Double) + ) + + sqrtDefaultEffort i@(Interval l r) = + ((ArithUpDn.fieldOpsDefaultEffort l, + ArithUpDn.mixedFieldOpsDefaultEffort l sampleI) + , + Int1To10 10 + , + ((NumOrd.minmaxDefaultEffort l, NumOrd.pCompareDefaultEffort l), + ArithUpDn.convertDefaultEffort l sampleD) + ) + where + sampleI = 1 :: Int + sampleD = 1 :: Double + +sqrtDefaultEffortWithIters :: + (ArithUpDn.RoundedMixedFieldEffort e Int, + ArithUpDn.RoundedFieldEffort e, + ArithUpDn.Convertible e Double, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) + => + (Interval e) -> + Int -> + ArithInOut.SqrtEffortIndicator (Interval e) +sqrtDefaultEffortWithIters i@(Interval l r) n = + ((ArithUpDn.fieldOpsDefaultEffort l, + ArithUpDn.mixedFieldOpsDefaultEffort l sampleI) + , + Int1To10 n + , + ((NumOrd.minmaxDefaultEffort l, NumOrd.pCompareDefaultEffort l), + ArithUpDn.convertDefaultEffort l sampleD) + ) + where + sampleI = 1 :: Int + sampleD = 1 :: Double + + +instance + (ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedField e, + ArithUpDn.Convertible e Double, + HasZero e, HasOne e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + Show e) + => + (ArithInOut.RoundedSquareRoot (Interval e)) + where + sqrtOutEff + ((effortField, effortMixedField), + (Int1To10 effortNewton), + ((effortMinmax, effortComp), effortConv)) + (Interval l r) = + case NumOrd.pEqualEff effortComp l r of + Just True -> sqrtL + _ -> Interval sqrtLL sqrtRR + + where + sqrtL@(Interval sqrtLL _) = sqrt l + sqrtR@(Interval _ sqrtRR) = sqrt r + sqrt = + sqrtOutThinArg + effortField + effortMixedField + effortMinmax + effortComp + effortConv + effortNewton + sqrtInEff + ((effortField, effortMixedField), + (Int1To10 effortNewton), + ((effortMinmax, effortComp), effortConv)) + (Interval l r) = + case NumOrd.pEqualEff effortComp l r of + Just True -> Interval sqrtLR sqrtLL -- invert + _ -> Interval sqrtLR sqrtRL + + where + (Interval sqrtLL sqrtLR) = sqrt l + (Interval sqrtRL sqrtRR) = sqrt r + sqrt = + sqrtOutThinArg + effortField + effortMixedField + effortMinmax + effortComp + effortConv + effortNewton + +sqrtOutIters, sqrtInIters :: + (ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedField e, + ArithUpDn.Convertible e Double, + HasZero e, HasOne e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + Show e) + => + Int -> (Interval e) -> (Interval e) +sqrtOutIters n i = ArithInOut.sqrtOutEff (sqrtDefaultEffortWithIters i n) i +sqrtInIters n i = ArithInOut.sqrtInEff (sqrtDefaultEffortWithIters i n) i + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps/Sqrt.hs b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps/Sqrt.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/ElementaryFromFieldOps/Sqrt.hs @@ -0,0 +1,387 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps.Sqrt + Description : an interval-specific implementation of sqrt + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + An interval-specific implementation of sqrt. +-} + +module Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps.Sqrt where + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import Numeric.AERN.RealArithmetic.NumericOrderRounding.OpsImplicitEffort +import Numeric.AERN.RealArithmetic.NumericOrderRounding.InPlace.OpsImplicitEffort +--import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.OpsImplicitEffort +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd +import Numeric.AERN.Basics.RefinementOrder.OpsImplicitEffort + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval + +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Consistency +import Numeric.AERN.Basics.Effort +import Numeric.AERN.Basics.Mutable +import Numeric.AERN.Basics.Exception + +import Control.Exception (throw) +import Control.Monad.ST (ST) + +sqrtOutThinArg :: + (HasZero e, HasOne e, Show e, + NumOrd.RoundedLattice e, + NumOrd.PartialComparison e, + ArithUpDn.Convertible e Double, + ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedField e) => + ArithUpDn.FieldOpsEffortIndicator e -> + ArithUpDn.MixedFieldOpsEffortIndicator e Int -> + NumOrd.MinmaxEffortIndicator e -> + NumOrd.PartialCompareEffortIndicator e -> + ArithUpDn.ConvertEffortIndicator e Double -> + Int {-^ the highest number of iterations of Newton method to make -} -> + e {-^ @x@ as a singleton interval -} -> + (Interval e) {-^ @sqrt(x)@ -} +sqrtOutThinArg + effortField + effortMixedField + effortMinmax + effortCompare + effortToDouble + maxIters + x + | sureIsZero x = zero + | not (sureAbove0 x) = + case (sureAbove0 (neg x)) of + True -> + throw $ AERNDomViolationException $ + "sqrtOutThinArg: applied to a negative argument " ++ show x + _ -> + throw $ AERNMaybeDomViolationException $ + "sqrtOutThinArg: cannot check that sqrt is applied to a positive argument " ++ show x + | xRecipSqrtDownInFastRegion = +-- unsafePrint ("AERN: sqrtOutThinArg: lower bound in fast region") $ + Interval + (x *. xRecipSqrtDown) + (x *^ xRecipSqrtUp) -- best upper bound estimate based on an error estimate of the lower bound + | sureAbove0 xRecipSqrtDown = +-- unsafePrint ("AERN: sqrtOutThinArg: lower bound NOT in fast region, using division") $ + Interval + (x *. xRecipSqrtDown) + (recipUp xRecipSqrtDown) + -- an upper bound using division - introduces a fairly large error; used when iteration has not reached the fast region + | otherwise = +-- unsafePrint ("AERN: sqrtOutThinArg: lower bound too close to zero, using dummy upper bound") $ + Interval + (x *. xRecipSqrtDown) + (NumOrd.maxUpEff effortMinmax x one) + -- a dummy fallback upper bound where lower bound is too close to 0 + where + (xRecipSqrtDownPrev, xRecipSqrtDown) = recipSqrtDown + xRecipSqrtDownInFastRegion = + case ArithUpDn.convertDnEff effortToDouble t of + Just lowerBound -> lowerBound > (0.381966012 :: Double) -- (3 - sqrt 5)/2 + Nothing -> False + where + t = (xRecipSqrtDownPrev *. xRecipSqrtDownPrev) *. x + xRecipSqrtUp = + -- only valid in "fast" region, ie where the error is smaller + -- than the gap between the results of the last two iterations + (xRecipSqrtLastUp +^ xRecipSqrtLastUp) +^ (neg xRecipSqrtDownPrev) + xRecipSqrtLastUp = + (xRecipSqrtDownPrev /^| 2 ) + *^ + (3 |+^ (neg $ x *. (xRecipSqrtDownPrev *. xRecipSqrtDownPrev))) + + sureAbove0 t = + case ArithUpDn.convertDnEff effortToDouble t of + Just lowerBound -> lowerBound > (0 :: Double) + Nothing -> False + + sureIsZero t = + case NumOrd.pEqualEff effortCompare t zero of + Just True -> True + _ -> False + + x1 +^ x2 = ArithUpDn.addUpEff effortAdd x1 x2 + x1 *^ x2 = ArithUpDn.multUpEff effortMult x1 x2 + x1 *. x2 = ArithUpDn.multDnEff effortMult x1 x2 + recipUp x = ArithUpDn.recipUpEff effortDiv x + recipDn x = ArithUpDn.recipDnEff effortDiv x + n |+^ x = ArithUpDn.mixedAddUpEff effortAddInt x (n :: Int) + n |+. x = ArithUpDn.mixedAddDnEff effortAddInt x (n :: Int) + x /^| n = ArithUpDn.mixedDivUpEff effortDivInt x (n :: Int) + x /.| n = ArithUpDn.mixedDivDnEff effortDivInt x (n :: Int) + + effortAdd = ArithUpDn.fldEffortAdd x effortField + effortMult = ArithUpDn.fldEffortMult x effortField + effortDiv = ArithUpDn.fldEffortDiv x effortField + effortAddInt = ArithUpDn.mxfldEffortAdd x (0::Int) effortMixedField + effortDivInt = ArithUpDn.mxfldEffortDiv x (0::Int) effortMixedField + + recipSqrtDown + | q0OK = -- computed an approximation in the stable region: + iterRecipSqrt maxIters zero q0 -- then iterate! + | otherwise = (zero, zero) -- zero is an always correct lower approximation + where + (q0OK, q0) = + (sureAbove0 xPlusOneUp && sureAbove0 babylon2, + recipDn babylon2) + where + -- babylon2 = (x+5)/4 - 1/(x+1) rounded upwards + -- ie two Babylonian iterations + -- \ t -> (t + x/t)/2 + -- starting with t_0 = x: + -- t_1 = (x + 1)/2 + -- t_2 = ((x + 1)/2 + x/((x + 1)/2))/2 = + -- ((x + 1)/2 + 2x/(x + 1))/2 = + -- ((x + 1)^2 + 4x)/4(x+1) = + -- (x^2 + 6x + 1)/4(x+1) = + -- (x^2 + 6x + 5 - 4)/4(x+1) = + -- ((x + 5)(x + 1) - 4)/4(x+1) = + -- (x + 5) - 1/(x+1) + babylon2 = xPlus5Div4Up +^ (neg xPlusOneRecipDn) + xPlus5Div4Up = ((5::Int) |+^ x) /^| (4::Int) + xPlusOneRecipDn = recipDn xPlusOneUp + xPlusOneUp = (1::Int) |+^ x + -- iteratively improve q, a lower bound on sqrt(1/x) + -- using the formula q_{n+1} = (q_n / 2) * (3 - x * q_n * q_n) + -- quoted eg in http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Iterative_methods_for_reciprocal_square_roots + iterRecipSqrt maxIters qNm2 qNm1 + | maxIters > 0 && sureAbove0 qNm1 = +-- unsafePrint ("AERN: sqrtOutThinArg: recipSqrtDown: iterRecipSqrt: maxIters = " ++ show maxIters) $ + iterRecipSqrt (maxIters - 1) qNm1 qN + | otherwise = (qNm2, qNm1) + where + qN = + (qNm1 /.| (2::Int)) + *. + ((3::Int) |+. (neg $ x *^ (qNm1 *^ qNm1))) + + +sqrtOutThinArgInPlace :: + (CanBeMutable e, + HasZero e, HasOne e, Show e, + NumOrd.RoundedLattice e, + NumOrd.RoundedLatticeInPlace e, + NumOrd.PartialComparison e, + ArithUpDn.Convertible e Double, + ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedMixedFieldInPlace e Int, + ArithUpDn.RoundedField e, + ArithUpDn.RoundedFieldInPlace e) => + ArithUpDn.FieldOpsEffortIndicator e -> + ArithUpDn.MixedFieldOpsEffortIndicator e Int -> + NumOrd.MinmaxEffortIndicator e -> + NumOrd.PartialCompareEffortIndicator e -> + ArithUpDn.ConvertEffortIndicator e Double -> + Mutable (Interval e) s {-^ where to write the result @sqrt(x)@ -} -> + Int {-^ the highest number of iterations of Newton method to make -} -> + Mutable e s {-^ @x@ viewed as a singleton interval -} -> + ST s () +sqrtOutThinArgInPlace + effortField + effortMixedField + effortMinmax + effortCompare + effortToDouble + resM@(MInterval resLM resRM) + maxIters + xM + = + do + -- we need x - a pure version of xM for branching conditions: + x <- unsafeReadMutable xM + -- unsafe is OK because we do not write into xM and we write to resM only as the last thing + -- and the value of x does not escape beyond this function + let ?addUpDnEffort = ArithUpDn.fldEffortAdd x effortField + let ?multUpDnEffort = ArithUpDn.fldEffortMult x effortField + let ?mixedAddUpDnEffort = ArithUpDn.mxfldEffortAdd x (0::Int) effortMixedField + let ?mixedMultUpDnEffort = ArithUpDn.mxfldEffortMult x (0::Int) effortMixedField + let ?mixedDivUpDnEffort = ArithUpDn.mxfldEffortDiv x (0::Int) effortMixedField + computeSqrt xM x + where + computeSqrt xM x = + do + continue + where + continue + | sureIsZero x = + writeMutable resM zero + | not (sureAbove0 x) = + case (sureAbove0 (neg x)) of + True -> + throw $ AERNDomViolationException $ + "sqrtOutThinArgInPlace: applied to a negative argument " ++ show x + _ -> + throw $ AERNMaybeDomViolationException $ + "sqrtOutThinArgInPlace: cannot check that sqrt is applied to a positive argument " ++ show x + | otherwise = + do + -- declare some variables: + temp1M <- makeMutable zero + temp2M <- makeMutable zero + -- iterate using Newton's method, assign results of last two iterations to the above vars: + prevFirst <- recipSqrtDown temp1M temp2M + case prevFirst of + True -> assignBounds temp1M temp2M + False -> assignBounds temp2M temp1M + where + assignBounds xRecipSqrtDownPrevM xRecipSqrtDownM = + do + xRecipSqrtDown <- unsafeReadMutable xRecipSqrtDownM + xRecipSqrtDownPrev <- unsafeReadMutable xRecipSqrtDownPrevM + -- assign lower bound resL := x *. xRecipSqrtDown: + ArithUpDn.multDnInPlaceEff effortMult resLM xM xRecipSqrtDownM + -- assign upper bound the best applicable method out of three methods: + constructUpperBound xRecipSqrtDownM xRecipSqrtDown xRecipSqrtDownPrevM xRecipSqrtDownPrev + constructUpperBound xRecipSqrtDownM xRecipSqrtDown xRecipSqrtDownPrevM xRecipSqrtDownPrev + | xRecipSqrtDownInFastRegion = + do + -- in fast region, use the difference between the last two approx: + + -- first, we need an upwards-rounded version of xRecipSqrtDown: + -- xRecipSqrtLastUp := newton... xRecipSqrtDownPrev: + let xRecipSqrtLastUpM = xRecipSqrtDownM -- safely reuse variable + newtonIterateUp xRecipSqrtLastUpM xRecipSqrtDownPrevM + + -- now compute and use the difference: + -- xRecipSqrtUp := 2*^xRecipSqrtLastUp -^ xRecipSqrtDownPrev: + -- only valid in "fast" region, ie where the error is smaller + -- than the gap between the results of the last two iterations + let xRecipSqrtUpM = xRecipSqrtLastUpM -- safely reuse variable + xRecipSqrtLastUpM *^|= (2 :: Int) + xRecipSqrtLastUpM -^= xRecipSqrtDownPrevM + + -- assign upper bound resR := x *^ xRecipSqrtUp: + ArithUpDn.multUpInPlaceEff effortMult resRM xM xRecipSqrtUpM + | sureAbove0 xRecipSqrtDown = + do + -- compute upper bound resR := 1 /^ xRecipSqrtDown: + -- introduces a fairly large error; + -- used when iteration has not reached the fast region + ArithUpDn.recipUpInPlaceEff effortDiv resRM xRecipSqrtDownM + | otherwise = + do + -- a dummy fallback upper bound where lower bound is too close to 0: + unsafeWriteMutable resRM $ + NumOrd.maxUpEff effortMinmax x one + where + xRecipSqrtDownInFastRegion = + case ArithUpDn.convertDnEff effortToDouble t of + Just lowerBound -> lowerBound > (0.381966012 :: Double) -- (3 - sqrt 5)/2 + Nothing -> False + where + t = + (xRecipSqrtDownPrev *. xRecipSqrtDownPrev) *. x + newtonIterateUp resM tM = + -- assumes no aliasing between resM and tM, does not change tM + do + -- res := + -- (t /^| 2 ) + -- *^ + -- (3 |+^ (neg $ x *. (t *. t))) + ArithUpDn.multUpInPlaceEff effortMult resM tM tM + resM *^= xM + negInPlace resM resM + resM +.|= (3 :: Int) + resM *.= tM + resM /.|= (2 :: Int) + + newtonIterateDn resM tM = + -- assumes no aliasing between resM and tM, does not change tM + do + -- res := + -- (t /.| 2 ) + -- *. + -- (3 |+. (neg $ x *^ (t *^ t))) + ArithUpDn.multDnInPlaceEff effortMult resM tM tM + resM *.= xM + negInPlace resM resM + resM +^|= (3 :: Int) + resM *^= tM + resM /^|= (2 :: Int) + + sureAbove0 t = + case ArithUpDn.convertDnEff effortToDouble t of + Just lowerBound -> lowerBound > (0 :: Double) + Nothing -> False + + sureIsZero t = + case NumOrd.pEqualEff effortCompare t zero of + Just True -> True + _ -> False + + effortAdd = ArithUpDn.fldEffortAdd x effortField + effortMult = ArithUpDn.fldEffortMult x effortField + effortDiv = ArithUpDn.fldEffortDiv x effortField + effortMultInt = ArithUpDn.mxfldEffortMult x (0::Int) effortMixedField + effortAddInt = ArithUpDn.mxfldEffortAdd x (0::Int) effortMixedField + effortDivInt = ArithUpDn.mxfldEffortDiv x (0::Int) effortMixedField + + recipSqrtDown aM bM + | q0OK = -- computed an approximation in the stable region: + do + writeMutable aM zero + writeMutable bM q0 + iterRecipSqrt maxIters True aM bM -- then iterate! + | otherwise = + do + writeMutable aM zero -- zero is an always correct lower approximation + writeMutable bM zero + return True + where + (q0OK, q0) = + (sureAbove0 xPlusOneUp && sureAbove0 babylon2, + recipDn babylon2) + where + -- babylon2 = (x+5)/4 - 1/(x+1) rounded upwards + -- ie two Babylonian iterations + -- \ t -> (t + x/t)/2 + -- starting with t_0 = x: + -- t_1 = (x + 1)/2 + -- t_2 = ((x + 1)/2 + x/((x + 1)/2))/2 = + -- ((x + 1)/2 + 2x/(x + 1))/2 = + -- ((x + 1)^2 + 4x)/4(x+1) = + -- (x^2 + 6x + 1)/4(x+1) = + -- (x^2 + 6x + 5 - 4)/4(x+1) = + -- ((x + 5)(x + 1) - 4)/4(x+1) = + -- (x + 5) - 1/(x+1) + babylon2 + = xPlus5Div4Up +^ (neg xPlusOneRecipDn) + xPlus5Div4Up + = ((5::Int) |+^ x) /^| (4::Int) + xPlusOneRecipDn = recipDn xPlusOneUp + xPlusOneUp + = (1::Int) |+^ x + recipDn = ArithUpDn.recipDnEff effortDiv + -- iteratively improve q, a lower bound on sqrt(1/x) + -- using the formula q_{n+1} = (q_n / 2) * (3 - x * q_n * q_n) + -- quoted eg in http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Iterative_methods_for_reciprocal_square_roots + iterRecipSqrt maxIters prevFirst aM bM = + -- assuming aM and bM do not alias + do + qNm1 <- unsafeReadMutable qNm1M + case maxIters > 0 && sureAbove0 qNm1 of + False -> -- should not or cannot continue iterating + return prevFirst -- indicate which of the two variables has the older result + True -> + do + newtonIterateDn qN qNm1M + iterRecipSqrt (maxIters - 1) (not prevFirst) bM aM -- swap the variables + where + (qNm2M, qNm1M) + | prevFirst = (aM,bM) + | otherwise = (bM,aM) + qN = qNm2M diff --git a/src/Numeric/AERN/RealArithmetic/Interval/ExactOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/ExactOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/ExactOps.hs @@ -0,0 +1,38 @@ +{-| + Module : Numeric.AERN.RealArithmetic.Interval.ExactOps + Description : exact zero, one and neg for intervals + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Exact zero, one and neg for intervals. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.ExactOps +() +where + +import Numeric.AERN.Basics.Interval +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import Numeric.AERN.RealArithmetic.ExactOps + +instance (HasZero e, NumOrd.PartialComparison e) => HasZero (Interval e) where + zero = Interval zero zero + +instance (HasOne e) => HasOne (Interval e) where + one = Interval one one + +instance (HasInfinities e) => HasInfinities (Interval e) where + plusInfinity = Interval plusInfinity plusInfinity + minusInfinity = Interval minusInfinity minusInfinity + excludesPlusInfinity (Interval l r) = excludesPlusInfinity r + excludesMinusInfinity (Interval l r) = excludesMinusInfinity l + +instance (Neg e) => Neg (Interval e) where + neg (Interval l r) = Interval (neg r) (neg l) + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/FieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/FieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/FieldOps.hs @@ -0,0 +1,438 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-# LANGUAGE UndecidableInstances #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.FieldOps + Description : refinement rounded basic operations for intervals + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Refinement rounded basic operations for intervals. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.FieldOps +() +where + +import Numeric.AERN.Basics.Interval + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval.ExactOps + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import Numeric.AERN.RealArithmetic.RefinementOrderRounding + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd + + +instance (ArithUpDn.RoundedAddEffort e) => RoundedAddEffort (Interval e) where + type AddEffortIndicator (Interval e) = ArithUpDn.AddEffortIndicator e + addDefaultEffort (Interval l r) = ArithUpDn.addDefaultEffort l + +instance (ArithUpDn.RoundedAdd e) => RoundedAdd (Interval e) where + addInEff effort (Interval l1 r1) (Interval l2 r2) = + Interval + (ArithUpDn.addUpEff effort l1 l2) + (ArithUpDn.addDnEff effort r1 r2) + addOutEff effort (Interval l1 r1) (Interval l2 r2) = + Interval + (ArithUpDn.addDnEff effort l1 l2) + (ArithUpDn.addUpEff effort r1 r2) + +instance (ArithUpDn.RoundedAdd e, Neg e) => RoundedSubtr (Interval e) + +instance + (NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedAbsEffort (Interval e) + where + type AbsEffortIndicator (Interval e) = + (NumOrd.PartialCompareEffortIndicator e, NumOrd.MinmaxEffortIndicator e) + absDefaultEffort (Interval l r) = + (NumOrd.pCompareDefaultEffort l, NumOrd.minmaxDefaultEffort l) + +instance + (ArithUpDn.RoundedAbs e, + HasZero e, Neg e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedAbs (Interval e) + where + absOutEff = absOutUsingCompMax + absInEff = absInUsingCompMax + + +instance + (ArithUpDn.RoundedMultiplyEffort e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedMultiplyEffort (Interval e) + where + type MultEffortIndicator (Interval e) = + (NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e, + ArithUpDn.MultEffortIndicator e) + multDefaultEffort (Interval l r) = + (NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l, + ArithUpDn.multDefaultEffort l) + +instance + (ArithUpDn.RoundedMultiply e, + HasZero e, Neg e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedMultiply (Interval e) + where + multOutEff (effortComp, effortMinmax, effortMult) i1 i2 = + fromEndpoints $ + multiplyIntervals + (pNonnegNonposEff effortComp) + (ArithUpDn.multDnEff effortMult) (ArithUpDn.multUpEff effortMult) + (NumOrd.minDnEff effortMinmax) -- minL + (NumOrd.minUpEff effortMinmax) -- minR + (NumOrd.maxDnEff effortMinmax) -- maxL + (NumOrd.maxUpEff effortMinmax) -- maxR + (NumOrd.minDnEff effortMinmax) + (NumOrd.maxUpEff effortMinmax) + i1 i2 + multInEff (effortComp, effortMinmax, effortMult) i1 i2 = + fromEndpoints $ + multiplyIntervals + (pNonnegNonposEff effortComp) + (ArithUpDn.multUpEff effortMult) (ArithUpDn.multDnEff effortMult) + (NumOrd.minUpEff effortMinmax) -- minL + (NumOrd.minDnEff effortMinmax) -- minR + (NumOrd.maxUpEff effortMinmax) -- maxL + (NumOrd.maxDnEff effortMinmax) -- maxR + (NumOrd.maxUpEff effortMinmax) + (NumOrd.minDnEff effortMinmax) + i1 i2 + +multiplyIntervals + pNonnegNonpos timesL timesR minL minR maxL maxR + combineL combineR + (Interval l1 r1) (Interval l2 r2) = + let _ = [minL, maxR, combineL, combineR] in + case (pNonnegNonpos l1, -- sign of l1 + pNonnegNonpos r1, -- sign of r1 + pNonnegNonpos l2, -- sign of l2 + pNonnegNonpos r2 -- sign of r2 + ) of + + ----------------------------------------------------------- + -- cases where i1 or i2 is known to be positive or negative + ----------------------------------------------------------- + -- i1 negative, i2 positive + ((_, Just True), (_, Just True), (Just True, _), (Just True, _)) -> + (l1 `timesL` r2, r1 `timesR` l2) + -- i1 negative, i2 negative + ((_, Just True), (_, Just True), (_, Just True), (_, Just True)) -> + (r1 `timesL` r2, l1 `timesR` l2) + -- i1 negative, i2 consistent and containing zero + ((_, Just True), (_, Just True), (_, Just True), (Just True, _)) -> + (l1 `timesL` r2, l1 `timesR` l2) + -- i1 negative, i2 anti-consistent and anti-containing zero + ((_, Just True), (_, Just True), (Just True, _), (_, Just True)) -> + (r1 `timesL` r2, r1 `timesR` l2) + -- i1 negative, nothing known about i2: + ((_, Just True), (_, Just True), _, _) -> + ((r1 `timesL` r2) `combineL` (l1 `timesL` r2), + (r1 `timesR` l2) `combineR` (l1 `timesR` l2)) + + -- i1 positive, i2 positive + ((Just True, _), (Just True, _), (Just True, _), (Just True, _)) -> + (l1 `timesL` l2, r1 `timesR` r2) + -- i1 positive, i2 negative + ((Just True, _), (Just True, _), (_, Just True), (_, Just True)) -> + (r1 `timesL` l2, l1 `timesR` r2) + -- i1 positive, i2 consistent and containing zero + ((Just True, _), (Just True, _), (_, Just True), (Just True, _)) -> + (r1 `timesL` l2, r1 `timesR` r2) + -- i1 positive, i2 anti-consistent and anti-containing zero + ((Just True, _), (Just True, _), (Just True, _), (_, Just True)) -> + (l1 `timesL` l2, l1 `timesR` r2) + + -- i1 positive, nothing known about i2: + ((Just True, _), (Just True, _), _, _) -> + ((r1 `timesL` l2) `combineL` (l1 `timesL` l2), + (r1 `timesR` r2) `combineR` (l1 `timesR` r2)) + + + -- i1 consistent and containing zero, i2 positive + ((_, Just True), (Just True, _), (Just True, _), (Just True, _)) -> + (l1 `timesL` r2, r1 `timesR` r2) + -- i1 anti-consistent and anti-containing zero, i2 positive + ((Just True, _), (_, Just True), (Just True, _), (Just True, _)) -> + (l1 `timesL` l2, r1 `timesR` l2) + -- nothing known about i1, i2 positive + (_, _, (Just True, _), (Just True, _)) -> + ((l1 `timesL` r2) `combineL` (l1 `timesL` l2), + (r1 `timesR` r2) `combineR` (r1 `timesR` l2)) + + -- i1 consistent and containing zero, i2 negative + ((_, Just True), (Just True, _), (_, Just True), (_, Just True)) -> + (r1 `timesL` l2, l1 `timesR` l2) + -- i1 anti-consistent and anti-containing zero, i2 negative + ((Just True, _), (_, Just True), (_, Just True), (_, Just True)) -> + (r1 `timesL` r2, l1 `timesR` r2) + -- nothing known about i1, i2 negative + (_, _, (_, Just True), (_, Just True)) -> + ((r1 `timesL` r2) `combineL` (r1 `timesL` l2), + (l1 `timesR` r2) `combineR` (l1 `timesR` l2)) + + ----------------------------------------------------------- + -- cases where both i1 or i2 are around zero + ----------------------------------------------------------- + + -- i1 consistent and containing zero, i2 consistent and containing zero + ((_, Just True), (Just True, _), (_, Just True), (Just True, _)) -> + ((l1 `timesL` r2) `minL` (r1 `timesL` l2), + (l1 `timesR` l2) `maxR` (r1 `timesR` r2)) + -- i1 consistent and containing zero, i2 anti-consistent and anti-containing zero + ((_, Just True), (Just True, _), (Just True, _), (_, Just True)) -> + (zero, zero) + -- i1 consistent and containing zero, i2 unknown + ((_, Just True), (Just True, _), _, _) -> + (((l1 `timesL` r2) `combineL` (r1 `timesL` l2)) `combineL` zero, + ((l1 `timesR` l2) `combineR` (r1 `timesR` r2)) `combineR` zero) + + -- i1 anti-consistent and anti-containing zero, i2 consistent and containing zero + ((Just True, _), (_, Just True), (_, Just True), (Just True, _)) -> + (zero, zero) + -- i1 anti-consistent and anti-containing zero, i2 anti-consistent and anti-containing zero + ((Just True, _), (_, Just True), (Just True, _), (_, Just True)) -> + ((l1 `timesL` l2) `maxL` (r1 `timesL` r2), + (l1 `timesR` r2) `minR` (r1 `timesR` l2)) + -- i1 anti-consistent and anti-containing zero, i2 unknown + ((Just True, _), (_, Just True), _, _) -> + ((l1 `timesL` l2) `combineL` (r1 `timesL` r2) `combineL` zero, + (l1 `timesR` r2) `combineR` (r1 `timesR` l2) `combineR` zero) + + -- i1 unknown, i2 anti-consistent and anti-containing zero + (_, _, (Just True, _), (_, Just True)) -> + ((l1 `timesL` l2) `combineL` (r1 `timesL` r2) `combineL` zero, + (l1 `timesR` r2) `combineR` (r1 `timesR` l2) `combineR` zero) + + -- i1 unknown, i2 consistent and containing zero + (_, _, (_, Just True), (Just True, _)) -> + ((l1 `timesL` r2) `combineL` (r1 `timesL` l2) `combineL` zero, + (l1 `timesR` l2) `combineR` (r1 `timesR` r2) `combineR` zero) + + -- both i1 and i2 unknown sign + _ -> + (foldl1 combineL [l1 `timesL` r2, r1 `timesL` l2, l1 `timesL` l2, r1 `timesL` r2], + foldl1 combineR [l1 `timesR` r2, r1 `timesR` l2, l1 `timesR` l2, r1 `timesR` r2]) + +instance + (ArithUpDn.RoundedPowerNonnegToNonnegIntEffort e, + ArithUpDn.RoundedMultiplyEffort e, + NumOrd.PartialComparison e, NumOrd.RoundedLatticeEffort e + ) => + RoundedPowerToNonnegIntEffort (Interval e) + where + type PowerToNonnegIntEffortIndicator (Interval e) = + (ArithUpDn.PowerNonnegToNonnegIntEffortIndicator e, + NumOrd.PartialCompareEffortIndicator e, + PowerToNonnegIntEffortIndicatorFromMult (Interval e)) + powerToNonnegIntDefaultEffort i@(Interval l r) = + (ArithUpDn.powerNonnegToNonnegIntDefaultEffort l, + NumOrd.pCompareDefaultEffort l, + powerToNonnegIntDefaultEffortFromMult i) + +instance + (ArithUpDn.RoundedPowerNonnegToNonnegInt e, + ArithUpDn.RoundedMultiply e, + HasZero e, HasOne e, Neg e, + NumOrd.PartialComparison e, NumOrd.RoundedLattice e + ) => + RoundedPowerToNonnegInt (Interval e) + where + powerToNonnegIntInEff + (effPowerEndpt, effComp, effPowerFromMult@(_,effMinMax,_)) + i@(Interval l r) n = + case (pNonnegNonposEff effComp l, pNonnegNonposEff effComp r) of + ((Just True, _), (Just True, _)) -> -- both non-negative + Interval lPowerUp hPowerDn + ((_, Just True), (_, Just True)) -> -- both non-positive + case even n of + True -> Interval hNegPowerUp lNegPowerDn -- switching sign! + False -> Interval lNegNegPowerUp hNegNegPowerDn + _ -> -- may involve crossing zero, revert to the default: + case even n of + True -> + NumOrd.maxInEff effMinMax zero iPowerFromMult + -- take advantage of the fact that the result is non-negative + False -> iPowerFromMult + where + lPowerUp = ArithUpDn.powerNonnegToNonnegIntUpEff effPowerEndpt l n + hPowerDn = ArithUpDn.powerNonnegToNonnegIntDnEff effPowerEndpt r n + lNegPowerDn = ArithUpDn.powerNonnegToNonnegIntDnEff effPowerEndpt (neg l) n + hNegPowerUp = ArithUpDn.powerNonnegToNonnegIntUpEff effPowerEndpt (neg r) n + lNegNegPowerUp = neg lNegPowerDn + hNegNegPowerDn = neg hNegPowerUp + iPowerFromMult = powerToNonnegIntInEffFromMult effPowerFromMult i n + powerToNonnegIntOutEff + (effPowerEndpt, effComp, effPowerFromMult@(_,effMinMax,_)) + i@(Interval l r) n = + case (pNonnegNonposEff effComp l, pNonnegNonposEff effComp r) of + ((Just True, _), (Just True, _)) -> -- both non-negative + Interval lPowerDn hPowerUp + ((_, Just True), (_, Just True)) -> -- both non-positive + case even n of + True -> Interval hNegPowerDn lNegPowerUp -- switching sign! + False -> Interval lNegNegPowerDn hNegNegPowerUp + _ -> -- may involve crossing zero, revert to the default: + case even n of + True -> + NumOrd.maxOutEff effMinMax zero iPowerFromMult + -- take advantage of the fact that the result is non-negative + False -> iPowerFromMult + where + lPowerDn = ArithUpDn.powerNonnegToNonnegIntDnEff effPowerEndpt l n + hPowerUp = ArithUpDn.powerNonnegToNonnegIntUpEff effPowerEndpt r n + lNegPowerUp = ArithUpDn.powerNonnegToNonnegIntUpEff effPowerEndpt (neg l) n + hNegPowerDn = ArithUpDn.powerNonnegToNonnegIntDnEff effPowerEndpt (neg r) n + lNegNegPowerDn = neg lNegPowerUp + hNegNegPowerUp = neg hNegPowerDn + iPowerFromMult = powerToNonnegIntOutEffFromMult effPowerFromMult i n + +instance + (ArithUpDn.RoundedMultiplyEffort e, ArithUpDn.RoundedDivideEffort e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedDivideEffort (Interval e) + where + type DivEffortIndicator (Interval e) = + (NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e, + (ArithUpDn.MultEffortIndicator e, + ArithUpDn.DivEffortIndicator e)) + divDefaultEffort (Interval l r) = + (NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l, + (ArithUpDn.multDefaultEffort l, + ArithUpDn.divDefaultEffort l)) + +instance + (ArithUpDn.RoundedMultiply e, ArithUpDn.RoundedDivide e, + HasZero e, Neg e, HasOne e, NumOrd.HasExtrema e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedDivide (Interval e) + where + divOutEff (effortComp, effortMinmax, (effortMult, effortDiv)) i1 i2 = + multOutEff (effortComp, effortMinmax, effortMult) i1 $ + recipInterval + (pPosNonnegNegNonposEff effortComp) + (ArithUpDn.divDnEff effortDiv) + (ArithUpDn.divUpEff effortDiv) + RefOrd.bottom + i2 + divInEff (effortComp, effortMinmax, (effortMult, effortDiv)) i1 i2 = + multInEff (effortComp, effortMinmax, effortMult) i1 $ + recipInterval + (pPosNonnegNegNonposEff effortComp) + (ArithUpDn.divUpEff effortDiv) + (ArithUpDn.divDnEff effortDiv) + RefOrd.top + i2 + + +recipInterval pPosNonnegNegNonpos divL divR fallback (Interval l r) = + case (pPosNonnegNegNonpos l, pPosNonnegNegNonpos r) of + -- positive: + ((Just True, _, _, _), (Just True, _, _, _)) -> + Interval (divL one r) (divR one l) + -- negative: + ((_, _, Just True, _), (_, _, Just True, _)) -> + Interval (divL one r) (divR one l) + -- consistent around zero: + ((_, _, _, Just True), (_, Just True, _, _)) -> + RefOrd.bottom + -- anti-consistent around zero: + ((_, Just True, _, _), (_,_,_, Just True)) -> + RefOrd.top + -- unknown: + _ -> + fallback + +instance + (ArithUpDn.RoundedRingEffort e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedRingEffort (Interval e) + where + type RingOpsEffortIndicator (Interval e) = + (ArithUpDn.RingOpsEffortIndicator e, + NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e) + ringOpsDefaultEffort (Interval l r) = + (ArithUpDn.ringOpsDefaultEffort l, + NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l) + ringEffortAdd (Interval l r) (effortRing, effortComp, effortMinmax) = + ArithUpDn.ringEffortAdd l effortRing + ringEffortMult (Interval l r) (effortRing, effortComp, effortMinmax) = + (effortComp, effortMinmax, + ArithUpDn.ringEffortMult l effortRing) + ringEffortPow i@(Interval l r) e@(effortRing, effortComp, effortMinmax) = + (ArithUpDn.ringEffortPow l effortRing, + effortComp, + ringEffortMult i e) + +instance + (ArithUpDn.RoundedRing e, + ArithUpDn.RoundedPowerNonnegToNonnegInt e, + HasOne e, HasZero e, Neg e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedRing (Interval e) + + +instance + (ArithUpDn.RoundedFieldEffort e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedFieldEffort (Interval e) + where + type FieldOpsEffortIndicator (Interval e) = + (ArithUpDn.FieldOpsEffortIndicator e, + NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e) + fieldOpsDefaultEffort (Interval l r) = + (ArithUpDn.fieldOpsDefaultEffort l, + NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l) + fldEffortAdd (Interval l r) (effortField, effortComp, effortMinmax) = + ArithUpDn.fldEffortAdd l effortField + fldEffortMult (Interval l r) (effortField, effortComp, effortMinmax) = + (effortComp, effortMinmax, + ArithUpDn.fldEffortMult l effortField) + fldEffortPow i@(Interval l r) e@(effortField, effortComp, effortMinmax) = + (ArithUpDn.fldEffortPow l effortField, + effortComp, + fldEffortMult i e) + fldEffortDiv (Interval l r) (effortField, effortComp, effortMinmax) = + (effortComp, effortMinmax, + (ArithUpDn.fldEffortMult l effortField, + ArithUpDn.fldEffortDiv l effortField)) + + +instance + (ArithUpDn.RoundedField e, + ArithUpDn.RoundedPowerNonnegToNonnegInt e, + HasZero e, Neg e, HasOne e, + NumOrd.HasExtrema e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedField (Interval e) + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Floating.hs b/src/Numeric/AERN/RealArithmetic/Interval/Floating.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Floating.hs @@ -0,0 +1,144 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Floating + Description : Floating instance using outer rounding + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Floating instance using outer rounding. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.Floating where + +import Prelude hiding (EQ, LT, GT) +import qualified Prelude +import Numeric.AERN.Basics.PartialOrdering + +import Numeric.AERN.RealArithmetic.Interval.ExactOps +import Numeric.AERN.RealArithmetic.Interval.Measures +import Numeric.AERN.RealArithmetic.Interval.Conversion +import Numeric.AERN.RealArithmetic.Interval.FieldOps +import Numeric.AERN.RealArithmetic.Interval.MixedFieldOps +import Numeric.AERN.RealArithmetic.Interval.SpecialConst + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.OpsDefaultEffort +import Numeric.AERN.RealArithmetic.ExactOps + +import Numeric.AERN.Basics.RefinementOrder.OpsDefaultEffort + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import Numeric.AERN.Basics.NumericOrder.OpsDefaultEffort + +import Numeric.AERN.Basics.ShowInternals +import Numeric.AERN.Basics.Interval + +import Numeric.AERN.Basics.Exception +import Control.Exception + + +instance + (NumOrd.PartialComparison e, ShowInternals e) => + Eq (Interval e) + where + i1 == i2 = + case i1 |==? i2 of + Just r -> r + _ -> throw $ + AERNException $ + "equality cannot be decided for: " + ++ show i1 ++ " == " + ++ show i2 + ++ "\n consider replacing == with NumericOrder.OpsDefaultEffort.|==?" + ++ "\n or with NumericOrder.OpsImplicitEffort.|==?" + ++ "\n or with NumericOrder.pCompareEff" + +instance + (NumOrd.PartialComparison e, ShowInternals e, NumOrd.RoundedLattice e) => + Ord (Interval e) + where + compare i1 i2 = + case NumOrd.pCompareEff (NumOrd.pCompareDefaultEffort i1) i1 i2 of + Just EQ -> Prelude.EQ + Just LT -> Prelude.LT + Just GT -> Prelude.GT + _ -> throw $ + AERNException $ + "comparison cannot be decided for: compare " + ++ show i1 ++ " " + ++ show i2 + ++ "\n consider replacing with ops defined at NumericOrder.OpsDefaultEffort" + ++ "\n or at NumericOrder.OpsImplicitEffort" + ++ "\n or with NumericOrder.pCompareEff" + max i1 i2 = maxOut i1 i2 + min i1 i2 = minOut i1 i2 + +instance + (ArithUpDn.Convertible Integer e, + ShowInternals e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + HasZero e, + ArithUpDn.RoundedRing e, + ArithUpDn.RoundedAbs e) => + Num (Interval e) + where + negate = neg + (+) = (<+>) + (*) = (<*>) + abs = absOut + fromInteger n = + result + where + result = + ArithInOut.convertOutEff (ArithInOut.convertDefaultEffort n result) n + signum a = + error $ "signum not implemented for Interval" + +instance + (ArithUpDn.Convertible Integer e, + ArithUpDn.Convertible Rational e, + Eq e, ShowInternals e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + HasZero e, HasOne e, NumOrd.HasExtrema e, + ArithUpDn.RoundedField e, + ArithUpDn.RoundedAbs e) => + Fractional (Interval e) + where + (/) = () + fromRational r = + result + where + result = + ArithInOut.convertOutEff (ArithInOut.convertDefaultEffort r result) r + +instance + (ArithUpDn.Convertible Integer e, + ArithUpDn.Convertible Rational e, + Eq e, + ShowInternals e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + HasZero e, + ArithUpDn.RoundedField e, + NumOrd.HasExtrema e, + ArithUpDn.RoundedAbs e, + ArithUpDn.RoundedSpecialConst e, + ArithInOut.RoundedExponentiation (Interval e), + ArithInOut.RoundedSquareRoot (Interval e)) => + Floating (Interval e) + where + pi = piOut + exp = expOut + sqrt = sqrtOut diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Measures.hs b/src/Numeric/AERN/RealArithmetic/Interval/Measures.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Measures.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Measures + Description : distance and imprecision for intervals + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Distance and imprecision for intervals. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.Measures +() +where + +import Numeric.AERN.RealArithmetic.Measures +import Numeric.AERN.RealArithmetic.ExactOps + +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Consistency + +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.OpsImplicitEffort + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd +import Numeric.AERN.Basics.RefinementOrder.OpsImplicitEffort + + +instance (HasDistance e, ArithInOut.RoundedAdd (Distance e)) => + HasDistance (Interval e) where + type Distance (Interval e) = Distance e + type DistanceEffortIndicator (Interval e) = + (DistanceEffortIndicator e, ArithInOut.AddEffortIndicator (Distance e)) + distanceDefaultEffort (Interval l r) = + (effortDist, effortAdd) + where + effortDist = distanceDefaultEffort l + effortAdd = ArithInOut.addDefaultEffort d + d = distanceBetweenEff effortDist l r + distanceBetweenEff (effortDist, effortAdd) (Interval l1 r1) (Interval l2 r2) = + let ?addInOutEffort = effortAdd in + distL <+> distR + where + distL = distanceBetweenEff effortDist l1 l2 + distR = distanceBetweenEff effortDist r1 r2 + +instance + (HasDistance e, RefOrd.OuterRoundedLattice (Distance e), Neg (Distance e), + NumOrd.PartialComparison e) => + HasImprecision (Interval e) + where + type Imprecision (Interval e) = Distance e + type ImprecisionEffortIndicator (Interval e) = + (DistanceEffortIndicator e, + RefOrd.JoinMeetOutEffortIndicator (Distance e), + ConsistencyEffortIndicator (Interval e)) + imprecisionDefaultEffort i@(Interval l r) = + (effortDist, effortMeet, consistencyDefaultEffort i) + where + effortDist = distanceDefaultEffort l + effortMeet = RefOrd.joinmeetOutDefaultEffort d + d = distanceBetweenEff effortDist l r + imprecisionOfEff (effortDist, effortMeet, effortConsistency) i@(Interval l r) = + let + ?joinmeetOutEffort = effortMeet + in + case (isConsistentEff effortConsistency i) of + Just True -> dist + Just False -> neg dist + Nothing -> dist <⊓> (neg dist) + where + dist = distanceBetweenEff effortDist l r + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/MixedFieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/MixedFieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/MixedFieldOps.hs @@ -0,0 +1,214 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.MixedFieldOps + Description : rounded basic arithmetic operations mixing 2 types + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Rounded basic arithmetical operations mixing an interval and another type. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.MixedFieldOps +() +where + +import Numeric.AERN.Basics.Interval + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval.ExactOps + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import Numeric.AERN.RealArithmetic.RefinementOrderRounding + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd + +instance (ArithUpDn.RoundedMixedAddEffort e tn) => + RoundedMixedAddEffort (Interval e) tn + where + type MixedAddEffortIndicator (Interval e) tn = ArithUpDn.MixedAddEffortIndicator e tn + mixedAddDefaultEffort (Interval l r) n = ArithUpDn.mixedAddDefaultEffort l n + +instance (ArithUpDn.RoundedMixedAdd e tn) => + RoundedMixedAdd (Interval e) tn + where + mixedAddInEff effort (Interval l2 r2) n = + Interval + (ArithUpDn.mixedAddUpEff effort l2 n) + (ArithUpDn.mixedAddDnEff effort r2 n) + mixedAddOutEff effort (Interval l2 r2) n = + Interval + (ArithUpDn.mixedAddDnEff effort l2 n) + (ArithUpDn.mixedAddUpEff effort r2 n) + +instance (ArithUpDn.RoundedMixedMultiplyEffort e tn, + NumOrd.PartialComparison tn, NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e) => + RoundedMixedMultiplyEffort (Interval e) tn + where + type MixedMultEffortIndicator (Interval e) tn = + ((NumOrd.PartialCompareEffortIndicator tn, + NumOrd.PartialCompareEffortIndicator e), + NumOrd.MinmaxEffortIndicator e, + ArithUpDn.MixedMultEffortIndicator e tn) + mixedMultDefaultEffort (Interval l r) n = + ((NumOrd.pCompareDefaultEffort n, + NumOrd.pCompareDefaultEffort l), + NumOrd.minmaxDefaultEffort l, + ArithUpDn.mixedMultDefaultEffort l n) + +instance (ArithUpDn.RoundedMixedMultiply e tn, + HasZero tn, HasZero e, + NumOrd.PartialComparison tn, NumOrd.PartialComparison e, + NumOrd.RoundedLattice e) => + RoundedMixedMultiply (Interval e) tn + where + mixedMultInEff ((effortCompS, effortCompE), effortMinmax, effortMult) i n = + fromEndpoints $ + multiplySingletonWithInterval + (pNonnegNonposEff effortCompS) + (pNonnegNonposEff effortCompE) + (flip $ ArithUpDn.mixedMultUpEff effortMult) + (flip $ ArithUpDn.mixedMultDnEff effortMult) + (NumOrd.maxUpEff effortMinmax) + (NumOrd.minDnEff effortMinmax) + n i + mixedMultOutEff ((effortCompS, effortCompE), effortMinmax, effortMult) i n = + fromEndpoints $ + multiplySingletonWithInterval + (pNonnegNonposEff effortCompS) + (pNonnegNonposEff effortCompE) + (flip $ ArithUpDn.mixedMultDnEff effortMult) + (flip $ ArithUpDn.mixedMultUpEff effortMult) + (NumOrd.minDnEff effortMinmax) + (NumOrd.maxUpEff effortMinmax) + n i + +multiplySingletonWithInterval + sNonnegNonpos iNonnegNonpos timesL timesR + combineL combineR + s1 (Interval l2 r2) = + let _ = [combineL, combineR] in + case (sNonnegNonpos s1, -- sign of s1 + iNonnegNonpos l2, -- sign of l2 + iNonnegNonpos r2 -- sign of r2 + ) of + + -- s1 is zero + ((Just True, Just True), _, _) -> + (zero, zero) + + -- s1 non negative + ((Just True, _), _, _) -> + (s1 `timesL` l2, s1 `timesR` r2) + + -- s1 non positive + ((_, Just True), _, _) -> + (s1 `timesL` r2, s1 `timesR` l2) + + -- nothing known about s1, i2 positive + (_, (Just True, _), (Just True, _)) -> + ((s1 `timesL` r2) `combineL` (s1 `timesL` l2), + (s1 `timesR` r2) `combineR` (s1 `timesR` l2)) + + -- nothing known about s1, i2 negative + (_, (_, Just True), (_, Just True)) -> + ((s1 `timesL` r2) `combineL` (s1 `timesL` l2), + (s1 `timesR` r2) `combineR` (s1 `timesR` l2)) + + -- both s1 and i2 are around zero + _ -> + ((s1 `timesL` l2) `combineL` (s1 `timesL` r2) `combineL` zero, + (s1 `timesR` l2) `combineR` (s1 `timesR` r2) `combineR` zero) + -- need to include zero to account for + -- consistent vs anti-consistent cases giving constant 0 + +instance (RoundedDivideEffort (Interval e), + Convertible tn (Interval e)) => + RoundedMixedDivideEffort (Interval e) tn + where + type MixedDivEffortIndicator (Interval e) tn = + (DivEffortIndicator (Interval e), + ConvertEffortIndicator tn (Interval e)) + mixedDivDefaultEffort = mixedDivDefaultEffortByConversion + +instance (RoundedDivide (Interval e), + Convertible tn (Interval e)) => + RoundedMixedDivide (Interval e) tn + where + mixedDivInEff = mixedDivInEffByConversion + mixedDivOutEff = mixedDivOutEffByConversion + +instance (RoundedMixedAddEffort (Interval e) tn, + RoundedMixedMultiplyEffort (Interval e) tn, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e, + NumOrd.PartialComparison tn, + ArithUpDn.RoundedMixedRingEffort e tn) => + RoundedMixedRingEffort (Interval e) tn + where + type MixedRingOpsEffortIndicator (Interval e) tn = + (ArithUpDn.MixedRingOpsEffortIndicator e tn, + (NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e, + NumOrd.PartialCompareEffortIndicator tn)) + mixedRingOpsDefaultEffort i@(Interval l r) n = + (ArithUpDn.mixedRingOpsDefaultEffort l n, + (NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l, + NumOrd.pCompareDefaultEffort n)) + mxringEffortAdd (Interval l r) n (effortRing, _) = + ArithUpDn.mxringEffortAdd l n effortRing + mxringEffortMult (Interval l r) n (effortRing, (effortCompEpt, effortMinmax, effortCompS)) = + ((effortCompS, effortCompEpt), + effortMinmax, + ArithUpDn.mxringEffortMult l n effortRing) + +instance (RoundedMixedAdd (Interval e) tn, + RoundedMixedMultiply (Interval e) tn, + RoundedMixedRingEffort (Interval e) tn) => + RoundedMixedRing (Interval e) tn + +instance (RoundedMixedRingEffort (Interval e) tn, + RoundedMixedDivideEffort (Interval e) tn, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeEffort e, + NumOrd.PartialComparison tn, + ArithUpDn.RoundedMixedFieldEffort e tn) => + RoundedMixedFieldEffort (Interval e) tn + where + type MixedFieldOpsEffortIndicator (Interval e) tn = + (ArithUpDn.MixedFieldOpsEffortIndicator e tn, + (NumOrd.PartialCompareEffortIndicator e, + NumOrd.MinmaxEffortIndicator e, + NumOrd.PartialCompareEffortIndicator tn), + MixedDivEffortIndicator (Interval e) tn) + mixedFieldOpsDefaultEffort i@(Interval l r) n = + (ArithUpDn.mixedFieldOpsDefaultEffort l n, + (NumOrd.pCompareDefaultEffort l, + NumOrd.minmaxDefaultEffort l, + NumOrd.pCompareDefaultEffort n), + mixedDivDefaultEffort i n) + mxfldEffortAdd (Interval l r) n (effortFld, _, _) = + ArithUpDn.mxfldEffortAdd l n effortFld + mxfldEffortMult (Interval l r) n (effortFld, (effortCompEpt, effortMinmax, effortCompS), _) = + ((effortCompS, effortCompEpt), + effortMinmax, + ArithUpDn.mxfldEffortMult l n effortFld) + mxfldEffortDiv _ _ (_, _, effortDiv) = effortDiv + +instance (RoundedMixedRing (Interval e) tn, + RoundedMixedDivide (Interval e) tn, + RoundedMixedFieldEffort (Interval e) tn) => + RoundedMixedField (Interval e) tn + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Mutable.hs b/src/Numeric/AERN/RealArithmetic/Interval/Mutable.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Mutable.hs @@ -0,0 +1,25 @@ +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Mutable + Description : instances of mutable arithmetic classes for MIntervals + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Instances of mutable arithmetic classes for MIntervals. +-} + +module Numeric.AERN.RealArithmetic.Interval.Mutable +( + module Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps, + module Numeric.AERN.RealArithmetic.Interval.Mutable.FieldOps, + module Numeric.AERN.RealArithmetic.Interval.Mutable.MixedFieldOps +) +where + +import Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps +import Numeric.AERN.RealArithmetic.Interval.Mutable.FieldOps +import Numeric.AERN.RealArithmetic.Interval.Mutable.MixedFieldOps + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ElementaryFromFieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ElementaryFromFieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ElementaryFromFieldOps.hs @@ -0,0 +1,209 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Mutable.ElementaryFromFieldOps + Description : elementary in-place operations using generic direct implementation + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Elementary in-place operations using generic implementation directly from + field operations. +-} + +module Numeric.AERN.RealArithmetic.Interval.Mutable.ElementaryFromFieldOps + (expOutInPlaceIters, expInInPlaceIters, sqrtOutInPlaceIters, sqrtInInPlaceIters) +where + +import Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps + +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.ElementaryFromFieldOps.Exponentiation + +import Numeric.AERN.RealArithmetic.Interval.ElementaryFromFieldOps.Sqrt + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut +import Numeric.AERN.RealArithmetic.RefinementOrderRounding.OpsImplicitEffort +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd +import Numeric.AERN.Basics.RefinementOrder.OpsImplicitEffort + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval + +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Consistency +import Numeric.AERN.Basics.Effort + +import Numeric.AERN.Basics.Mutable + +instance + (CanBeMutable e, + ArithInOut.RoundedFieldInPlace (Interval e), + ArithInOut.RoundedMixedFieldInPlace (Interval e) Int, + ArithInOut.RoundedPowerToNonnegIntInPlace (Interval e), + ArithInOut.RoundedMixedField (Interval e) Int, + ArithInOut.RoundedField (Interval e), + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e), + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + RefOrd.OuterRoundedLattice (Interval e)) + => + (ArithInOut.RoundedExponentiationInPlace (Interval e)) + where + expOutInPlaceEff + ((effortField, effortMixedField), + (Int1To10 effortTaylor), + ((effortMeet, effortComp), effortConv)) + (MInterval resL resR) + (MInterval lM rM) = + do + (MInterval forgetMeL forgetMeR) <- makeMutable zero + expOutThinArgInPlace + effortField effortMixedField + effortMeet effortComp effortComp effortConv + (MInterval resL forgetMeR) + effortTaylor + (MInterval lM lM) + expOutThinArgInPlace + effortField effortMixedField + effortMeet effortComp effortComp effortConv + (MInterval forgetMeL resR) + effortTaylor + (MInterval rM rM) + expInInPlaceEff + ((effortField, effortMixedField), + (Int1To10 effortTaylor), + ((effortMeet, effortComp), effortConv)) + (MInterval resL resR) + (MInterval lM rM) = + do + (MInterval forgetMeL forgetMeR) <- makeMutable zero + expOutThinArgInPlace + effortField effortMixedField + effortMeet effortComp effortComp effortConv + (MInterval forgetMeL resL) + effortTaylor + (MInterval lM lM) + expOutThinArgInPlace + effortField effortMixedField + effortMeet effortComp effortComp effortConv + (MInterval resR forgetMeR) + effortTaylor + (MInterval rM rM) + +expOutInPlaceIters, expInInPlaceIters :: + (CanBeMutable e, + ArithInOut.RoundedFieldInPlace (Interval e), + ArithInOut.RoundedMixedFieldInPlace (Interval e) Int, + ArithInOut.RoundedPowerToNonnegIntInPlace (Interval e), + ArithInOut.RoundedMixedField (Interval e) Int, + ArithInOut.RoundedField (Interval e), + ArithUpDn.Convertible (Interval e) Int, + ArithInOut.Convertible Double (Interval e), + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + RefOrd.OuterRoundedLattice (Interval e)) + => + Int -> OpMutable1 (Interval e) s +expOutInPlaceIters n resM iM = + do + i <- unsafeReadMutable iM + ArithInOut.expOutInPlaceEff (expDefaultEffortWithIters i n) resM iM +expInInPlaceIters n resM iM = + do + i <- unsafeReadMutable iM + ArithInOut.expInInPlaceEff (expDefaultEffortWithIters i n) resM iM + +instance + (CanBeMutable e, Show e, + ArithUpDn.RoundedFieldInPlace e, + ArithUpDn.RoundedMixedFieldInPlace e Int, + ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedField e, + ArithUpDn.Convertible e Double, + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + NumOrd.RoundedLatticeInPlace e) + => + (ArithInOut.RoundedSquareRootInPlace (Interval e)) + where + sqrtOutInPlaceEff + = pureToMutable1Eff ArithInOut.sqrtOutEff + +-- the following is a proper in-place version - but it currently fails the test + +-- ((effortField, effortMixedField), +-- (Int1To10 effortNewton), +-- ((effortMeet, effortComp), effortConv)) +-- (MInterval resL resR) +-- (MInterval lM rM) = +-- do +-- (MInterval forgetMeL forgetMeR) <- makeMutable zero +-- sqrtOutThinArgInPlace +-- effortField effortMixedField +-- effortMeet effortComp effortConv +-- (MInterval resL forgetMeR) +-- effortNewton +-- lM +-- sqrtOutThinArgInPlace +-- effortField effortMixedField +-- effortMeet effortComp effortConv +-- (MInterval forgetMeL resR) +-- effortNewton +-- rM + sqrtInInPlaceEff + = pureToMutable1Eff ArithInOut.sqrtInEff +-- ((effortField, effortMixedField), +-- (Int1To10 effortNewton), +-- ((effortMeet, effortComp), effortConv)) +-- (MInterval resL resR) +-- (MInterval lM rM) = +-- do +-- (MInterval forgetMeL forgetMeR) <- makeMutable zero +-- sqrtOutThinArgInPlace +-- effortField effortMixedField +-- effortMeet effortComp effortConv +-- (MInterval forgetMeL resL) +-- effortNewton +-- lM +-- sqrtOutThinArgInPlace +-- effortField effortMixedField +-- effortMeet effortComp effortConv +-- (MInterval resR forgetMeR) +-- effortNewton +-- rM + +sqrtOutInPlaceIters, sqrtInInPlaceIters :: + (CanBeMutable e, Show e, + ArithUpDn.RoundedFieldInPlace e, + ArithUpDn.RoundedMixedFieldInPlace e Int, + ArithUpDn.RoundedMixedField e Int, + ArithUpDn.RoundedField e, + ArithUpDn.Convertible e Double, + HasZero e, HasOne e, + HasInfinities e, + NumOrd.PartialComparison e, + NumOrd.RoundedLattice e, + NumOrd.RoundedLatticeInPlace e) + => + Int -> OpMutable1 (Interval e) s +sqrtOutInPlaceIters n resM iM = + do + i <- unsafeReadMutable iM + ArithInOut.sqrtOutInPlaceEff (sqrtDefaultEffortWithIters i n) resM iM +sqrtInInPlaceIters n resM iM = + do + i <- unsafeReadMutable iM + ArithInOut.sqrtInInPlaceEff (sqrtDefaultEffortWithIters i n) resM iM + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ExactOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ExactOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/ExactOps.hs @@ -0,0 +1,31 @@ +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps + Description : neg for mutable intervals + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Exact neg for mutable intervals. + + This module is hidden and reexported via its parent Interval.Mutable. +-} + +module Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps where + +import Numeric.AERN.Basics.Mutable +import Numeric.AERN.Basics.Interval + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval.ExactOps + +instance (NegInPlace e, Neg e) => NegInPlace (Interval e) + where + negInPlace (MInterval lRes hRes) (MInterval lM hM) = + do + temp <- cloneMutable hRes + negInPlace temp lM -- mind potential aliasing hRes - hM + negInPlace lRes hM + assignMutable hRes temp diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Mutable/FieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/FieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/FieldOps.hs @@ -0,0 +1,591 @@ +{-# LANGUAGE FlexibleContexts, UndecidableInstances #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Mutable.FieldOps + Description : field operations for mutable intervals + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Field operations for mutable intervals. + + This module is hidden and reexported via its parent Interval.Mutable. +-} + +module Numeric.AERN.RealArithmetic.Interval.Mutable.FieldOps() where + +import Numeric.AERN.Basics.Mutable +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Interval.Mutable + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import Numeric.AERN.RealArithmetic.RefinementOrderRounding +import Numeric.AERN.RealArithmetic.Interval.FieldOps + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd + +import Control.Monad.ST (ST) + + +instance (ArithUpDn.RoundedAddInPlace e, CanBeMutable e) => + RoundedAddInPlace (Interval e) + where + addInInPlaceEff eff (MInterval resL resR) (MInterval aL aR) (MInterval bL bR) = + do + ArithUpDn.addUpInPlaceEff eff resL aL bL + ArithUpDn.addDnInPlaceEff eff resR aR bR + addOutInPlaceEff eff (MInterval resL resR) (MInterval aL aR) (MInterval bL bR) = + do + ArithUpDn.addDnInPlaceEff eff resL aL bL + ArithUpDn.addUpInPlaceEff eff resR aR bR + +instance + (ArithUpDn.RoundedAddInPlace e, + CanBeMutable e, + NegInPlace e) => + RoundedSubtrInPlace (Interval e) + +instance (RoundedAbs (Interval e), CanBeMutable (Interval e)) => + RoundedAbsInPlace (Interval e) + +instance + (ArithUpDn.RoundedMultiplyInPlace e, + NumOrd.RoundedLatticeInPlace e, + HasZero e, NumOrd.PartialComparison e, + CanBeMutable e) => + RoundedMultiplyInPlace (Interval e) + where + multOutInPlaceEff (effortComp, effortMinmax, effortMult) r i1 i2 = + multiplyIntervalsInPlace + (pNonnegNonposEff effortComp) + (ArithUpDn.multDnInPlaceEff effortMult) + (ArithUpDn.multUpInPlaceEff effortMult) + (NumOrd.minDnInPlaceEff effortMinmax) -- minL + (NumOrd.minUpInPlaceEff effortMinmax) -- minR + (NumOrd.maxDnInPlaceEff effortMinmax) -- maxL + (NumOrd.maxUpInPlaceEff effortMinmax) -- maxR + (NumOrd.minDnInPlaceEff effortMinmax) + (NumOrd.maxUpInPlaceEff effortMinmax) + r i1 i2 + multInInPlaceEff (effortComp, effortMinmax, effortMult) r i1 i2 = + multiplyIntervalsInPlace + (pNonnegNonposEff effortComp) + (ArithUpDn.multUpInPlaceEff effortMult) + (ArithUpDn.multDnInPlaceEff effortMult) + (NumOrd.minUpInPlaceEff effortMinmax) -- minL + (NumOrd.minDnInPlaceEff effortMinmax) -- minR + (NumOrd.maxUpInPlaceEff effortMinmax) -- maxL + (NumOrd.maxDnInPlaceEff effortMinmax) -- maxR + (NumOrd.maxUpInPlaceEff effortMinmax) + (NumOrd.minDnInPlaceEff effortMinmax) + r i1 i2 + +multiplyIntervalsInPlace :: + (CanBeMutable e, HasZero e) => + (e -> (Maybe Bool, Maybe Bool)) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (Mutable (Interval e) s) -> + (Mutable (Interval e) s) -> + (Mutable (Interval e) s) -> + ST s () +multiplyIntervalsInPlace + pNonnegNonpos timesLInPlace timesRInPlace + minLInPlace minRInPlace maxLInPlace maxRInPlace + combineLInPlace combineRInPlace + (MInterval lResM rResM) (MInterval l1M r1M) (MInterval l2M r2M) = + do + let _ = [minLInPlace, maxRInPlace, combineLInPlace, combineRInPlace] + l1 <- readMutable l1M + r1 <- readMutable r1M + l2 <- readMutable l2M + r2 <- readMutable r2M + case (pNonnegNonpos l1, -- sign of l1 + pNonnegNonpos r1, -- sign of r1 + pNonnegNonpos l2, -- sign of l2 + pNonnegNonpos r2 -- sign of r2 + ) of + + ----------------------------------------------------------- + -- cases where i1 or i2 is known to be positive or negative + ----------------------------------------------------------- + -- i1 negative, i2 positive + ((_, Just True), (_, Just True), (Just True, _), (Just True, _)) -> +-- (l1 `timesL` r2, r1 `timesR` l2) + assignResEndpointsUsingTimesLR l1M r2M r1M l2M + -- i1 negative, i2 negative + ((_, Just True), (_, Just True), (_, Just True), (_, Just True)) -> +-- (r1 `timesL` r2, l1 `timesR` l2) + assignResEndpointsUsingTimesLR r1M r2M l1M l2M + -- i1 negative, i2 consistent and containing zero + ((_, Just True), (_, Just True), (_, Just True), (Just True, _)) -> +-- (l1 `timesL` r2, l1 `timesR` l2) + assignResEndpointsUsingTimesLR l1M r2M l1M l2M + -- i1 negative, i2 anti-consistent and anti-containing zero + ((_, Just True), (_, Just True), (Just True, _), (_, Just True)) -> +-- (r1 `timesL` r2, r1 `timesR` l2) + assignResEndpointsUsingTimesLR r1M r2M r1M l2M + -- i1 negative, nothing known about i2: + ((_, Just True), (_, Just True), _, _) -> +-- ((r1 `timesL` r2) `combineL` (l1 `timesL` r2), +-- (r1 `timesR` l2) `combineR` (l1 `timesR` l2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 r1M r2M + timesLInPlace temp2 l1M r2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 r1M l2M + timesRInPlace temp2 l1M l2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + + -- i1 positive, i2 positive + ((Just True, _), (Just True, _), (Just True, _), (Just True, _)) -> +-- (l1 `timesL` l2, r1 `timesR` r2) + do + timesLInPlace lResM l1M l2M + timesRInPlace rResM r1M r2M + -- i1 positive, i2 negative + ((Just True, _), (Just True, _), (_, Just True), (_, Just True)) -> +-- (r1 `timesL` l2, l1 `timesR` r2) + assignResEndpointsUsingTimesLR r1M l2M l1M r2M + -- i1 positive, i2 consistent and containing zero + ((Just True, _), (Just True, _), (_, Just True), (Just True, _)) -> +-- (r1 `timesL` l2, r1 `timesR` r2) + assignResEndpointsUsingTimesLR r1M l2M r1M r2M + -- i1 positive, i2 anti-consistent and anti-containing zero + ((Just True, _), (Just True, _), (Just True, _), (_, Just True)) -> +-- (l1 `timesL` l2, l1 `timesR` r2) + assignResEndpointsUsingTimesLR l1M l2M l1M r2M + -- i1 positive, nothing known about i2: + ((Just True, _), (Just True, _), _, _) -> +-- ((r1 `timesL` l2) `combineL` (l1 `timesL` l2), +-- (r1 `timesR` r2) `combineR` (l1 `timesR` r2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 r1M l2M + timesLInPlace temp2 l1M l2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 r1M r2M + timesRInPlace temp2 l1M r2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + + + -- i1 consistent and containing zero, i2 positive + ((_, Just True), (Just True, _), (Just True, _), (Just True, _)) -> +-- (l1 `timesL` r2, r1 `timesR` r2) + assignResEndpointsUsingTimesLR l1M r2M r1M r2M + -- i1 anti-consistent and anti-containing zero, i2 positive + ((Just True, _), (_, Just True), (Just True, _), (Just True, _)) -> +-- (l1 `timesL` l2, r1 `timesR` l2) + assignResEndpointsUsingTimesLR l1M l2M r1M l2M + -- nothing known about i1, i2 positive + (_, _, (Just True, _), (Just True, _)) -> +-- ((l1 `timesL` r2) `combineL` (l1 `timesL` l2), +-- (r1 `timesR` r2) `combineR` (r1 `timesR` l2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M r2M + timesLInPlace temp2 l1M l2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 r1M r2M + timesRInPlace temp2 r1M l2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + + -- i1 consistent and containing zero, i2 negative + ((_, Just True), (Just True, _), (_, Just True), (_, Just True)) -> +-- (r1 `timesL` l2, l1 `timesR` l2) + assignResEndpointsUsingTimesLR r1M l2M l1M l2M + -- i1 anti-consistent and anti-containing zero, i2 negative + ((Just True, _), (_, Just True), (_, Just True), (_, Just True)) -> +-- (r1 `timesL` r2, l1 `timesR` r2) + assignResEndpointsUsingTimesLR r1M r2M l1M r2M + -- nothing known about i1, i2 negative + (_, _, (_, Just True), (_, Just True)) -> +-- ((r1 `timesL` r2) `combineL` (r1 `timesL` l2), +-- (l1 `timesR` r2) `combineR` (l1 `timesR` l2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 r1M r2M + timesLInPlace temp2 r1M l2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M r2M + timesRInPlace temp2 l1M l2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + + ----------------------------------------------------------- + -- cases where both i1 or i2 are around zero + ----------------------------------------------------------- + + -- i1 consistent and containing zero, i2 consistent and containing zero + ((_, Just True), (Just True, _), (_, Just True), (Just True, _)) -> +-- ((l1 `timesL` r2) `minL` (r1 `timesL` l2), +-- (l1 `timesR` l2) `maxR` (r1 `timesR` r2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M r2M + timesLInPlace temp2 r1M l2M + minLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M l2M + timesRInPlace temp2 r1M r2M + maxRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + -- i1 consistent and containing zero, i2 anti-consistent and anti-containing zero + ((_, Just True), (Just True, _), (Just True, _), (_, Just True)) -> +-- (zero, zero) + do + let z = zero + let _ = [z,l1] + writeMutable lResM z + writeMutable rResM z + -- i1 consistent and containing zero, i2 unknown + ((_, Just True), (Just True, _), _, _) -> +-- (((l1 `timesL` r2) `combineL` (r1 `timesL` l2)) `combineL` zero, +-- ((l1 `timesR` l2) `combineR` (r1 `timesR` r2)) `combineR` zero) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M r2M + timesLInPlace temp2 r1M l2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M l2M + timesRInPlace temp2 r1M r2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + let z = zero + let _ = [z,l1] + writeMutable temp1 z + combineLInPlace lResM lResM temp1 + combineRInPlace rResM rResM temp1 + + -- i1 anti-consistent and anti-containing zero, i2 consistent and containing zero + ((Just True, _), (_, Just True), (_, Just True), (Just True, _)) -> +-- (zero, zero) + do + let z = zero + let _ = [z,l1] + writeMutable lResM z + writeMutable rResM z + -- i1 anti-consistent and anti-containing zero, i2 anti-consistent and anti-containing zero + ((Just True, _), (_, Just True), (Just True, _), (_, Just True)) -> +-- ((l1 `timesL` l2) `maxL` (r1 `timesL` r2), +-- (l1 `timesR` r2) `minR` (r1 `timesR` l2)) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M l2M + timesLInPlace temp2 r1M r2M + maxLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M r2M + timesRInPlace temp2 r1M l2M + minRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + -- i1 anti-consistent and anti-containing zero, i2 unknown + ((Just True, _), (_, Just True), _, _) -> +-- ((l1 `timesL` l2) `combineL` (r1 `timesL` r2) `combineL` zero, +-- (l1 `timesR` r2) `combineR` (r1 `timesR` l2) `combineR` zero) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M l2M + timesLInPlace temp2 r1M r2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M r2M + timesRInPlace temp2 r1M l2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + let z = zero + let _ = [z,l1] + writeMutable temp1 z + combineLInPlace lResM lResM temp1 + combineRInPlace rResM rResM temp1 + + -- i1 unknown, i2 anti-consistent and anti-containing zero + (_, _, (Just True, _), (_, Just True)) -> +-- ((l1 `timesL` l2) `combineL` (r1 `timesL` r2) `combineL` zero, +-- (l1 `timesR` r2) `combineR` (r1 `timesR` l2) `combineR` zero) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 l1M l2M + timesLInPlace temp2 r1M r2M + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 l1M r2M + timesRInPlace temp2 r1M l2M + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + let z = zero + let _ = [z,l1] + writeMutable temp1 z + combineLInPlace lResM lResM temp1 + combineRInPlace rResM rResM temp1 + + -- i1 unknown, i2 consistent and containing zero + (_, _, (_, Just True), (Just True, _)) -> +-- ((l1 `timesL` r2) `combineL` (r1 `timesL` l2) `combineL` zero, +-- (l1 `timesR` l2) `combineR` (r1 `timesR` r2) `combineR` zero) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + + timesLInPlace temp1 l1M r2M + timesLInPlace temp2 r1M l2M + combineLInPlace temp3 temp1 temp2 + + timesRInPlace temp1 l1M l2M + timesRInPlace temp2 r1M r2M + combineRInPlace rResM temp1 temp2 + + assignMutable lResM temp3 + + let z = zero + let _ = [z,l1] + writeMutable temp1 z + combineLInPlace lResM lResM temp1 + combineRInPlace rResM rResM temp1 + + -- both i1 and i2 unknown sign + _ -> +-- (foldl1 combineL [l1 `timesL` r2, r1 `timesL` l2, l1 `timesL` l2, r1 `timesL` r2], +-- foldl1 combineR [l1 `timesR` r2, r1 `timesR` l2, l1 `timesR` l2, r1 `timesR` r2]) + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + + timesLInPlace temp1 l1M r2M + timesLInPlace temp2 r1M l2M + combineLInPlace temp1 temp1 temp2 + timesLInPlace temp2 l1M l2M + combineLInPlace temp1 temp1 temp2 + timesLInPlace temp2 r1M r2M + combineLInPlace temp3 temp1 temp2 + + timesRInPlace temp1 l1M r2M + timesRInPlace temp2 r1M l2M + combineRInPlace temp1 temp1 temp2 + timesRInPlace temp2 l1M l2M + combineRInPlace temp1 temp1 temp2 + timesRInPlace temp2 r1M r2M + combineRInPlace rResM temp1 temp2 + + assignMutable lResM temp3 + where + assignResEndpointsUsingTimesLR l1M l2M r1M r2M = + do + temp1 <- makeMutable zero + timesLInPlace temp1 l1M l2M -- beware of aliasing between res and param + timesRInPlace rResM r1M r2M + assignMutable lResM temp1 + + +instance + (RoundedSubtrInPlace (Interval e), + RoundedMultiplyInPlace (Interval e), + RoundedRingEffort (Interval e)) => + RoundedRingInPlace (Interval e) + +instance + (ArithUpDn.RoundedPowerNonnegToNonnegIntInPlace e, + RoundedPowerToNonnegInt (Interval e), + RoundedMultiplyInPlace (Interval e), + HasOne e, HasZero e, NegInPlace e, + NumOrd.PartialComparison e, NumOrd.RoundedLatticeInPlace e, + CanBeMutable e + ) => + RoundedPowerToNonnegIntInPlace (Interval e) + where + powerToNonnegIntInInPlaceEff + (effPowerEndpt, effComp, effPowerFromMult@(_,effMinMax,_)) + res@(MInterval resL resR) a@(MInterval aL aR) n = + do + l <- readMutable aL + r <- readMutable aR + case (pNonnegNonposEff effComp l, pNonnegNonposEff effComp r) of + ((Just True, _), (Just True, _)) -> -- both non-negative + do + ArithUpDn.powerNonnegToNonnegIntUpInPlaceEff + effPowerEndpt resL aL n + ArithUpDn.powerNonnegToNonnegIntDnInPlaceEff + effPowerEndpt resR aR n + ((_, Just True), (_, Just True)) -> -- both non-positive + do + -- negate the parameters, use the result as a temp space (may alias!): + negInPlace res a + -- compute the power of the positive interval: + ArithUpDn.powerNonnegToNonnegIntUpInPlaceEff + effPowerEndpt resL resL n + ArithUpDn.powerNonnegToNonnegIntDnInPlaceEff + effPowerEndpt resR resR n + case even n of + True -> + return () -- keep result positive + False -> + negInPlace res res -- back to the original sign + _ -> + do + powerToNonnegIntInInPlaceEffFromMult effPowerFromMult res a n + case even n of + True -> + do + let zeroI = zero + let _ = [zeroI, Interval l r] + zeroM <- unsafeMakeMutable zeroI + NumOrd.maxInInPlaceEff effMinMax res res zeroM + False -> return () + powerToNonnegIntOutInPlaceEff + (effPowerEndpt, effComp, effPowerFromMult@(_,effMinMax,_)) + res@(MInterval resL resR) a@(MInterval aL aR) n = + do + l <- readMutable aL + r <- readMutable aR + case (pNonnegNonposEff effComp l, pNonnegNonposEff effComp r) of + ((Just True, _), (Just True, _)) -> -- both non-negative + do + ArithUpDn.powerNonnegToNonnegIntDnInPlaceEff + effPowerEndpt resL aL n + ArithUpDn.powerNonnegToNonnegIntUpInPlaceEff + effPowerEndpt resR aR n + ((_, Just True), (_, Just True)) -> -- both non-positive + do + -- negate the parameters, use the result as a temp space (may alias!): + negInPlace res a + -- compute the power of the positive interval: + ArithUpDn.powerNonnegToNonnegIntDnInPlaceEff + effPowerEndpt resL resL n + ArithUpDn.powerNonnegToNonnegIntUpInPlaceEff + effPowerEndpt resR resR n + case even n of + True -> + return () -- keep result positive + False -> + negInPlace res res -- back to the original sign + _ -> + do + powerToNonnegIntOutInPlaceEffFromMult effPowerFromMult res a n + case even n of + True -> + do + let zeroI = zero + let _ = [zeroI, Interval l r] + zeroM <- unsafeMakeMutable zeroI + NumOrd.maxOutInPlaceEff effMinMax res res zeroM + False -> return () + +instance + (ArithUpDn.RoundedDivideInPlace e, + ArithUpDn.RoundedMultiplyInPlace e, + NumOrd.RoundedLatticeInPlace e, + HasZero e, HasOne e, Neg e, + NumOrd.PartialComparison e, NumOrd.HasExtrema e, + CanBeMutable e) => + RoundedDivideInPlace (Interval e) + where + divOutInPlaceEff + (effortComp, effortMinmax, (effortMult, effortDiv)) + res@(MInterval resL resR) a@(MInterval aL aR) b@(MInterval bL bR) = + do + temp <- makeMutable zero + recipIntervalInPlace + (pPosNonnegNegNonposEff effortComp) + divDn + divUp + bottom + temp b + multOutInPlaceEff (effortComp, effortMinmax, effortMult) res a temp + where + bottom = RefOrd.bottom + sampleI = getDummySample res + _ = [bottom, sampleI] + divUp = ArithUpDn.divUpInPlaceEff effortDiv + divDn = ArithUpDn.divDnInPlaceEff effortDiv + divInInPlaceEff + (effortComp, effortMinmax, (effortMult, effortDiv)) + res@(MInterval resL resR) a@(MInterval aL aR) b@(MInterval bL bR) = + do + temp <- makeMutable zero + recipIntervalInPlace + (pPosNonnegNegNonposEff effortComp) + divUp + divDn + top + temp b + multInInPlaceEff (effortComp, effortMinmax, effortMult) res a temp + where + top = RefOrd.top + sampleI = getDummySample res + _ = [top, sampleI] + divUp = ArithUpDn.divUpInPlaceEff effortDiv + divDn = ArithUpDn.divDnInPlaceEff effortDiv + +recipIntervalInPlace pPosNonnegNegNonpos divL divR fallback + res@(MInterval resL resR) a@(MInterval aL aR) = + do + let oneP = one + let top = RefOrd.top + let bottom = RefOrd.bottom + let _ = [top, bottom, fallback] + oneM <- unsafeMakeMutable oneP + l <- readMutable aL + r <- readMutable aR + let _ = [l, r, oneP] + case (pPosNonnegNegNonpos l, pPosNonnegNegNonpos r) of + -- positive: + ((Just True, _, _, _), (Just True, _, _, _)) -> + do + divL resL oneM aR + divR resR oneM aL + -- negative: + ((_, _, Just True, _), (_, _, Just True, _)) -> + do + divL resL oneM aR + divR resR oneM aL + -- consistent around zero: + ((_, _, _, Just True), (_, Just True, _, _)) -> + writeMutable res bottom + -- anti-consistent around zero: + ((_, Just True, _, _), (_,_,_, Just True)) -> + writeMutable res top + -- unknown: + _ -> + writeMutable res fallback + +instance + (ArithUpDn.RoundedFieldInPlace e, + ArithUpDn.RoundedPowerNonnegToNonnegIntInPlace e, + HasZero e, NegInPlace e, HasOne e, + NumOrd.HasExtrema e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeInPlace e) => + RoundedFieldInPlace (Interval e) + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/Mutable/MixedFieldOps.hs b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/MixedFieldOps.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/Mutable/MixedFieldOps.hs @@ -0,0 +1,205 @@ +{-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances, MultiParamTypeClasses #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.Mutable.MixedFieldOps + Description : mixed field operations for mutable intervals + Copyright : (c) Michal Konecny, Jan Duracz + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Mixed field operations for mutable intervals. + + This module is hidden and reexported via its parent Interval.Mutable. +-} + +module Numeric.AERN.RealArithmetic.Interval.Mutable.MixedFieldOps() where + +import Numeric.AERN.Basics.Mutable +import Numeric.AERN.Basics.Interval +import Numeric.AERN.Basics.Interval.Mutable + +import Numeric.AERN.RealArithmetic.ExactOps +import Numeric.AERN.RealArithmetic.Interval.Mutable.ExactOps + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import Numeric.AERN.RealArithmetic.RefinementOrderRounding +-- import Numeric.AERN.RealArithmetic.Interval.FieldOps +import Numeric.AERN.RealArithmetic.Interval.MixedFieldOps + +import qualified Numeric.AERN.Basics.NumericOrder as NumOrd +import qualified Numeric.AERN.Basics.RefinementOrder as RefOrd + +import Control.Monad.ST (ST) + + +instance (ArithUpDn.RoundedMixedAddInPlace t tn, CanBeMutable t) => + RoundedMixedAddInPlace (Interval t) tn + where + mixedAddInInPlaceEff eff (MInterval resL resR) (MInterval aL aR) n = + do + ArithUpDn.mixedAddUpInPlaceEff eff resL aL n + ArithUpDn.mixedAddDnInPlaceEff eff resR aR n + mixedAddOutInPlaceEff eff (MInterval resL resR) (MInterval aL aR) n = + do + ArithUpDn.mixedAddDnInPlaceEff eff resL aL n + ArithUpDn.mixedAddUpInPlaceEff eff resR aR n + +instance + (ArithUpDn.RoundedMixedMultiplyInPlace e tn, + NumOrd.RoundedLatticeInPlace e, + HasZero e, NumOrd.PartialComparison e, + HasZero tn, NumOrd.PartialComparison tn, + CanBeMutable e) => + RoundedMixedMultiplyInPlace (Interval e) tn + where + mixedMultInInPlaceEff + ((effortCompS,effortCompE), effortMinmax, effortMult) + r i1 s = + multiplySingletonAndIntervalInPlace + (pNonnegNonposEff effortCompS) + (pNonnegNonposEff effortCompE) + (ArithUpDn.mixedMultUpInPlaceEff effortMult) + (ArithUpDn.mixedMultDnInPlaceEff effortMult) + (NumOrd.maxUpInPlaceEff effortMinmax) + (NumOrd.minDnInPlaceEff effortMinmax) + r s i1 + mixedMultOutInPlaceEff + ((effortCompS,effortCompE), effortMinmax, effortMult) + r i1 s = + multiplySingletonAndIntervalInPlace + (pNonnegNonposEff effortCompS) + (pNonnegNonposEff effortCompE) + (ArithUpDn.mixedMultDnInPlaceEff effortMult) + (ArithUpDn.mixedMultUpInPlaceEff effortMult) + (NumOrd.minDnInPlaceEff effortMinmax) + (NumOrd.maxUpInPlaceEff effortMinmax) + r s i1 + +multiplySingletonAndIntervalInPlace :: + (CanBeMutable e, HasZero e) => + (tn -> (Maybe Bool, Maybe Bool)) -> + (e -> (Maybe Bool, Maybe Bool)) -> + (OpMutableNonmut e tn s) -> + (OpMutableNonmut e tn s) -> + (OpMutable2 e s) -> + (OpMutable2 e s) -> + (Mutable (Interval e) s) -> + tn -> + (Mutable (Interval e) s) -> + ST s () +multiplySingletonAndIntervalInPlace + sNonnegNonpos iNonnegNonpos + timesLInPlace timesRInPlace + combineLInPlace combineRInPlace + (MInterval lResM rResM) s1 (MInterval l2M r2M) = + do + let _ = [combineLInPlace, combineRInPlace] + l2 <- readMutable l2M + r2 <- readMutable r2M + let _ = [l2,r2] + case (sNonnegNonpos s1, -- sign of s1 + iNonnegNonpos l2, -- sign of l2 + iNonnegNonpos r2 -- sign of r2 + ) of + + -- s1 is zero + ((Just True, Just True), _, _) -> +-- (zero, zero) + do + let z = zero + let _ = [z,l2] + writeMutable lResM z + writeMutable rResM z + + -- s1 non negative + ((Just True, _), _, _) -> +-- (s1 `timesL` l2, s1 `timesR` r2) + do + assignResEndpointsUsingTimesLR l2M r2M + + -- s1 non positive + ((_, Just True), _, _) -> +-- (s1 `timesL` r2, s1 `timesR` l2) + do + assignResEndpointsUsingTimesLR r2M l2M + + -- nothing known about s1, i2 positive + (_, (Just True, _), (Just True, _)) -> +-- ((s1 `timesL` r2) `combineL` (s1 `timesL` l2), +-- (s1 `timesR` r2) `combineR` (s1 `timesR` l2)) + do + assignResEndpointsUsingBothOptions + return () + + -- nothing known about s1, i2 negative + (_, (_, Just True), (_, Just True)) -> +-- ((s1 `timesL` r2) `combineL` (s1 `timesL` l2), +-- (s1 `timesR` r2) `combineR` (s1 `timesR` l2)) + do + assignResEndpointsUsingBothOptions + return () + + + -- both s1 and i2 are around zero + _ -> +-- ((s1 `timesL` l2) `combineL` (s1 `timesL` r2) `combineL` zero, +-- (s1 `timesR` l2) `combineR` (s1 `timesR` r2) `combineR` zero) +-- -- need to include zero to account for +-- -- consistent vs anti-consistent cases giving constant 0 + do + temp1 <- assignResEndpointsUsingBothOptions + let z = zero + let _ = [z,l2] + writeMutable temp1 z + combineLInPlace lResM lResM temp1 + combineRInPlace rResM rResM temp1 + where + assignResEndpointsUsingTimesLR l2M r2M = + do + temp1 <- makeMutable zero + timesLInPlace temp1 l2M s1 -- beware of aliasing between res and param + timesRInPlace rResM r2M s1 + assignMutable lResM temp1 + assignResEndpointsUsingBothOptions = + do + temp1 <- makeMutable zero + temp2 <- makeMutable zero + temp3 <- makeMutable zero + timesLInPlace temp1 r2M s1 + timesLInPlace temp2 l2M s1 + combineLInPlace temp3 temp1 temp2 + timesRInPlace temp1 r2M s1 + timesRInPlace temp2 l2M s1 + combineRInPlace rResM temp1 temp2 + assignMutable lResM temp3 + return temp1 + +instance (RoundedDivideInPlace (Interval e), + Convertible tn (Interval e)) => + RoundedMixedDivideInPlace (Interval e) tn + where + mixedDivInInPlaceEff = mixedDivInInPlaceEffByConversion + mixedDivOutInPlaceEff = mixedDivOutInPlaceEffByConversion + +instance + (ArithUpDn.RoundedMixedRingInPlace e tn, + NumOrd.PartialComparison tn, + HasZero tn, + HasZero e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeInPlace e) => + RoundedMixedRingInPlace (Interval e) tn + +instance + (ArithUpDn.RoundedMixedFieldInPlace e tn, + RoundedDivideInPlace (Interval e), + Convertible tn (Interval e), + NumOrd.PartialComparison tn, + HasZero tn, + HasZero e, + NumOrd.PartialComparison e, + NumOrd.RoundedLatticeInPlace e) => + RoundedMixedFieldInPlace (Interval e) tn + diff --git a/src/Numeric/AERN/RealArithmetic/Interval/SpecialConst.hs b/src/Numeric/AERN/RealArithmetic/Interval/SpecialConst.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AERN/RealArithmetic/Interval/SpecialConst.hs @@ -0,0 +1,48 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ImplicitParams #-} +{-| + Module : Numeric.AERN.RealArithmetic.Interval.SpecialConst + Description : common constants such as pi + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mikkonecny@gmail.com + Stability : experimental + Portability : portable + + Common constants such as pi. + + This module is hidden and reexported via its parent Interval. +-} + +module Numeric.AERN.RealArithmetic.Interval.SpecialConst +() +where + +import qualified Numeric.AERN.RealArithmetic.NumericOrderRounding as ArithUpDn +import qualified Numeric.AERN.RealArithmetic.RefinementOrderRounding as ArithInOut + +import Numeric.AERN.Basics.Interval + +instance (ArithUpDn.RoundedSpecialConstEffort e) => + (ArithInOut.RoundedSpecialConstEffort (Interval e)) + where + type ArithInOut.SpecialConstEffortIndicator (Interval e) = + ArithUpDn.SpecialConstEffortIndicator e + specialConstDefaultEffort (Interval l r) = + ArithUpDn.specialConstDefaultEffort l + +instance (ArithUpDn.RoundedSpecialConst e) => + (ArithInOut.RoundedSpecialConst (Interval e)) + where + piInEff effort = + Interval (ArithUpDn.piUpEff effort) (ArithUpDn.piDnEff effort) + piOutEff effort = + Interval (ArithUpDn.piDnEff effort) (ArithUpDn.piUpEff effort) + eInEff effort = + Interval (ArithUpDn.eUpEff effort) (ArithUpDn.eDnEff effort) + eOutEff effort = + Interval (ArithUpDn.eDnEff effort) (ArithUpDn.eUpEff effort) +