diff --git a/AERN-Real.cabal b/AERN-Real.cabal --- a/AERN-Real.cabal +++ b/AERN-Real.cabal @@ -1,5 +1,5 @@ Name: AERN-Real -Version: 0.9.9 +Version: 0.10.0 Cabal-Version: >= 1.2 Build-Type: Simple License: BSD3 @@ -7,7 +7,7 @@ Author: Michal Konecny (Aston University) Copyright: (c) 2007-2008 Michal Konecny, Amin Farjudian, Jan Duracz Maintainer: mik@konecny.aow.cz -Stability: experimental +Stability: beta Category: Data, Math Synopsis: arbitrary precision interval arithmetic for approximating exact real numbers Tested-with: GHC ==6.10.1 @@ -25,10 +25,14 @@ . For an architectural overview, see module "Data.Number.ER.Real". . - Simple examples of usage can be found in module @Demo.hs@ in folder @tests@. -Extra-source-files: - tests/Demo.hs tests/Pi.hs tests/Matrix.hs -Data-files: + Simple examples of usage can be found in folder @demos@. + . + There is a built-in test suite and it can be evoked using + the module in the folder @tests@. + +Extra-Source-Files: + demos/Demo.hs demos/Pi.hs demos/Matrix.hs + tests/RunERIntervalTests.hs ChangeLog Flag use-hmpfr @@ -38,36 +42,43 @@ hs-source-dirs: src if flag(use-hmpfr) Build-Depends: - base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time, hmpfr == 0.1.3 + base >= 3, base < 4, containers, binary, html >= 1.0, regex-compat >= 0.71, stm, time, QuickCheck == 1.2.0.0, filepath, directory, hmpfr == 0.2 cpp-options: -DUSE_MPFR else Build-Depends: - base >= 3, base < 4, containers, binary >= 0.4, html >= 1.0, regex-compat >= 0.71, stm, time + base >= 3, base < 4, containers, binary, html >= 1.0, regex-compat >= 0.71, stm, time, QuickCheck == 1.2.0.0, filepath, directory Exposed-modules: - Data.Number.ER, - Data.Number.ER.Real, - Data.Number.ER.Real.DefaultRepr, - Data.Number.ER.Real.Base.MachineDouble, - Data.Number.ER.Real.Base.CombinedMachineAP, - Data.Number.ER.Real.Base.Rational, - Data.Number.ER.Real.Base.Float, - Data.Number.ER.Real.Base.MPFR, - Data.Number.ER.Real.Base, - Data.Number.ER.Real.Arithmetic.Elementary, - Data.Number.ER.Real.Arithmetic.Integration, - Data.Number.ER.Real.Arithmetic.LinearSolver, - Data.Number.ER.Real.Arithmetic.Taylor, - Data.Number.ER.Real.Arithmetic.Newton, - Data.Number.ER.Real.Approx.Sequence, - Data.Number.ER.Real.Approx.Elementary, - Data.Number.ER.Real.Approx.Interval, - Data.Number.ER.Real.Approx, - Data.Number.ER.Real.DomainBox, - Data.Number.ER.Real.DomainBox.IntMap, - Data.Number.ER.ShowHTML, - Data.Number.ER.PlusMinus, - Data.Number.ER.BasicTypes, - Data.Number.ER.Misc, - Data.Number.ER.MiscSTM, - Data.Number.ER.ExtendedInteger - + Data.Number.ER, + Data.Number.ER.BasicTypes, + Data.Number.ER.BasicTypes.DomainBox, + Data.Number.ER.BasicTypes.DomainBox.IntMap, + Data.Number.ER.BasicTypes.ExtendedInteger, + Data.Number.ER.BasicTypes.PlusMinus, + Data.Number.ER.BasicTypes.Tests.Generate, + Data.Number.ER.Misc, + Data.Number.ER.Misc.STM, + Data.Number.ER.Misc.Tests, + Data.Number.ER.Real, + Data.Number.ER.Real.Approx, + Data.Number.ER.Real.Approx.Elementary, + Data.Number.ER.Real.Approx.Interval, + Data.Number.ER.Real.Approx.OI, + Data.Number.ER.Real.Approx.Sequence, + Data.Number.ER.Real.Approx.Tests.Generate, + Data.Number.ER.Real.Approx.Tests.Properties, + Data.Number.ER.Real.Approx.Tests.Reporting, + Data.Number.ER.Real.Approx.Tests.Run, + Data.Number.ER.Real.Arithmetic.Elementary, + Data.Number.ER.Real.Arithmetic.Integration, + Data.Number.ER.Real.Arithmetic.LinearSolver, + Data.Number.ER.Real.Arithmetic.Newton, + Data.Number.ER.Real.Arithmetic.Taylor, + Data.Number.ER.Real.Base, + Data.Number.ER.Real.Base.CombinedMachineAP, + Data.Number.ER.Real.Base.Float, + Data.Number.ER.Real.Base.MachineDouble, + Data.Number.ER.Real.Base.MPFR, + Data.Number.ER.Real.Base.Rational, + Data.Number.ER.Real.Base.Tests.Generate, + Data.Number.ER.Real.DefaultRepr, + Data.Number.ER.ShowHTML diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ +0.9.10: 28 July 2009 + * switching to beta status + * new QuickCheck test suite covering most functionality + * new support for anti-consistent intervals (eg [2,0]) + (also called directed or improper; using Kaucher arithmetic) + * new support for inner-rounded interval arithemtic + * fixed errors in some elementary functions for extreme values + * fixed performance bug in arctan + * improved hierarchy of auxiliary modules + 0.9.9: 23 February 2009 * Small changes needed in other AERN packages: * New operation for domain boxes: get its dimension. diff --git a/demos/Demo.hs b/demos/Demo.hs new file mode 100644 --- /dev/null +++ b/demos/Demo.hs @@ -0,0 +1,149 @@ +{-# LANGUAGE CPP #-} +{-| + Module : Main + Description : simple examples of using AERN-Real + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Simple examples of using AERN-Real +-} +module Main where + +import qualified Data.Number.ER.Real as AERN +import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq) + +#ifdef USE_MPFR +--type B = AERN.BAP -- use pure Haskell floats +type B = AERN.BMAP -- use combination of double and pure Haskell floats +--type B = AERN.BMPFR -- use MPFR floats +#else +--type B = AERN.BAP -- use pure Haskell floats +type B = AERN.BMAP -- use combination of double and pure Haskell floats +#endif +type RA = AERN.RA B +type IRA = AERN.IRA B +type R = ConvergRealSeq IRA + +one :: R +one = 1 + +two :: R +two = 2 + +piSeq :: R +piSeq = ConvergRealSeq $ AERN.pi + +seqExp = convertFuncRA2Seq $ AERN.exp +seqSine = convertFuncRA2Seq $ AERN.sin +seqCosine = convertFuncRA2Seq $ AERN.cos + +main = + do + AERN.initialiseBaseArithmetic (0 :: RA) + putStrLn "****************************" + putStrLn "Testing interval arithmetic:" + putStrLn "****************************" + putStrLn "**** Fractions:" + putStrLn $ + "(default granularity, show internals) 1/3 =\n " ++ + AERN.showApprox 30 True True (1/3 :: RA) + putStrLn $ + "(granularity 50, show internals) 1/3 =\n " ++ + AERN.showApprox 30 True True ((AERN.setGranularityOuter 50 1/3) :: RA) + putStrLn $ + "(granularity 100, show internals) 1/3 =\n " ++ + AERN.showApprox 40 True True ((AERN.setGranularityOuter 100 1/3) :: RA) + putStrLn $ + "(granularity 100, do not show internals) 1/3 =\n " ++ + AERN.showApprox 40 True False ((AERN.setGranularityOuter 100 1/3) :: RA) + putStrLn $ + "(granularity 100, default show) 1/3 =\n " ++ + show ((AERN.setGranularityOuter 100 1/3) :: RA) + putStrLn "**** Exp:" + putStrLn $ + "(effort 5, granularity 50) exp 1 =\n " ++ + (show $ AERN.exp 5 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn $ + "(effort 10, granularity 50) exp 1 =\n " ++ + (show $ AERN.exp 10 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn $ + "(effort 10, granularity 100) exp 1 =\n " ++ + (show $ AERN.exp 10 (AERN.setGranularityOuter 100 (1::RA))) + putStrLn $ + "(effort 20, granularity 50) exp 1 =\n " ++ + (show $ AERN.exp 20 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn $ + "(effort 20, granularity 100) exp 1 =\n " ++ + (show $ AERN.exp 20 (AERN.setGranularityOuter 100 (1::RA))) + putStrLn "**** Pi:" + putStrLn $ + "(effort 10) pi =\n " ++ + (show $ (AERN.pi 10 :: RA)) + putStrLn $ + "(effort 50) pi =\n " ++ + (AERN.showApprox 20 True False $ (AERN.pi 50 :: RA)) + putStrLn $ + "(effort 100) pi =\n " ++ + (AERN.showApprox 35 True False $ (AERN.pi 100 :: RA)) + putStrLn $ + "(effort 200) pi =\n " ++ + (AERN.showApprox 65 True False $ (AERN.pi 200 :: RA)) + putStrLn $ + "(effort 400) pi =\n " ++ + (AERN.showApprox 125 True False $ (AERN.pi 400 :: RA)) + putStrLn "**** Sine:" + putStrLn $ + "(effort 10, granularity 50) sin 1 =\n " ++ + (show $ AERN.sin 10 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn $ + "(effort 10, granularity 100) sin 1 =\n " ++ + (show $ AERN.sin 10 (AERN.setGranularityOuter 100 (1::RA))) + putStrLn "**** Integration:" + putStrLn $ + "(effort 10, granularity 50) integrate exp 0 1 =\n " ++ + (show $ AERN.integrateContAdapt_R AERN.exp 10 0 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn $ + "(effort 20, granularity 50) integrate exp 0 1 =\n " ++ + (show $ AERN.integrateContAdapt_R AERN.exp 20 0 (AERN.setGranularityOuter 50 (1::RA))) +-- putStrLn $ +-- "(effort 30, granularity 50) integrate exp 0 1 =\n " ++ +-- (show $ AERN.integrateContAdapt_R AERN.exp 30 0 (AERN.setGranularityOuter 50 (1::RA))) + putStrLn "*****************************" + putStrLn "Testing convergent sequences:" + putStrLn "*****************************" +-- putStrLn $ "1 =\n " ++ show one +-- putStrLn $ "1 + 2 =\n " ++ (show $ one + two) + putStrLn "**** Fractions:" + putStrLn $ + "(precision 20) 1/3 =\n " ++ + (AERN.showConvergRealSeqAuto 20 $ one / 3) + putStrLn $ + "(precision 20) 100000000001/300000000000 =\n " ++ + (AERN.showConvergRealSeqAuto 20 $ (one + 100000000000)/300000000000 ) + putStrLn $ + "100000000001/300000000000 =? 1/3:\n " ++ + (show $ one/3 == 100000000001/300000000000) +-- putStrLn $ "abs -1 = " ++ (show $ abs (- one)) +-- putStrLn $ "neg 2 = " ++ (show $ negate two) +-- putStrLn $ "1 + 2 = " ++ (show $ one + 2) + putStrLn "**** Elementary:" + putStrLn $ + "(precision 30) exp 1 =\n " ++ + (AERN.showConvergRealSeqAuto 30 $ seqExp one) + putStrLn $ + "(precision 500) pi =\n " ++ + (AERN.showConvergRealSeqAuto 500 $ piSeq) + putStrLn $ + "(precision 30) cosine(1) =\n " ++ + (AERN.showConvergRealSeqAuto 30 $ seqCosine one) + putStrLn $ + "(precision 30) sine(1) =\n " ++ + (AERN.showConvergRealSeqAuto 30 $ seqSine one) + putStrLn "**** Integration:" + putStrLn $ -- very slow for precision > 4 + "(precision 3) integrate exp 0 1 =\n " ++ + (AERN.showConvergRealSeqAuto 3 $ AERN.integrateCont AERN.exp 0 one) diff --git a/demos/Matrix.hs b/demos/Matrix.hs new file mode 100644 --- /dev/null +++ b/demos/Matrix.hs @@ -0,0 +1,385 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Main + +where + +import qualified Data.Number.ER.Real as AERN +import Data.Number.ER.BasicTypes +import Data.Number.ER.Misc + +import Data.Maybe +import qualified Data.List as List +import qualified Data.Map as Map + +import qualified Data.Array.IArray as IAr +import qualified Data.Array.MArray as MAr +import qualified Data.Array.ST as STAr +import qualified Data.Ix as Ix +import qualified Data.Array.Base as BAr + +import Control.Monad.ST +import GHC.Arr + +#ifdef USE_MPFR +type B = AERN.BAP -- use pure Haskell floats +--type B = AERN.BMPFR -- use MPFR floats +#else +type B = AERN.BAP -- use pure Haskell floats +#endif +type RA = AERN.RA B +type IRA = AERN.IRA B + +testMatrixN = 100 +incrementGran = (+) 50 + +-- Hilbert 100x100 matrix: +addOneDiag = False +targetPrec = 167 -- approx 50 decimal digits after the point +initialGran = 2050 -- 100x100 +--initialGran = 2388 -- 100x100 Norbert's +--initialGran = 750 -- 50x50 +--initialGran = 300 -- 10x10 + +--targetPrec = 34 -- approx 10 decimal digits after the point +--initialGran = 1350 +--initialGran = 50 -- 50x50 + +-- Hilbert matrix + 1: +--addOneDiag = True +--targetPrec = 167 -- approx 50 decimal digits after the point +--initialGran = 200 + +--targetPrec = 34 -- approx 10 decimal digits after the point +--initialGran = 50 + +main = + do + AERN.initialiseBaseArithmetic (0 :: RA) + putStrLn $ + "Inverting the " ++ show n ++ "x" ++ show n ++ " Hilbert matrix " + ++ "with target binary precision " ++ show targetPrec ++ "..." +-- putStrLn $ +-- "sorted matrix elements = \n" ++ (unlines $ map show elemsSortedByPrec) + putStrLn $ + "sum of all elements in inverted matrix = " ++ show (sum elems) +-- putStrLn $ show (Matrix n n rarr) + where + n = testMatrixN + elems = IAr.elems rarr + elemsSortedByPrec = + List.sortBy comparePrec elems + where + comparePrec a b = + compare aPrecLO bPrecLO + where + aPrecLO = fst $ AERN.bounds $ aHI - aLO + (aLO, aHI) = AERN.bounds a + bPrecLO = fst $ AERN.bounds $ bHI - bLO + (bLO, bHI) = AERN.bounds b + rarr = + STAr.runSTArray $ + do + mInv@(Matrix _ _ rowsInv) <- + invert testMatrix +-- m <- testMatrix initialGran +-- mUnit@(Matrix _ _ rowsUnit) <- multM m mInv + return rowsInv + + +testMatrix :: + Granularity -> + ST s (STMatrix s IRA) +testMatrix gran = + do + marr <- MAr.newArray ((1,1),(n,n)) 0 + mapM (updateCell marr) assocsGran + return $ Matrix n n marr + where + assocsGran = map (mapSnd $ AERN.setMinGranularityOuter gran) assocs + assocs = +-- assocsMini + assocsHilbert gran n + assocsMini = + [((1,1),1), + ((1,2),3), + ((2,1),2), + ((2,2),0) + ] + n = testMatrixN + updateCell marr (ix, el) = + do + unsafeMatrixWrite marr n ix el + +assocsHilbert gran n = + [((i,j), coeff i j)| i <- [1..n], j <- [1..n]] + where + coeff i j + | addOneDiag && i == j = + 1 + oneOverIplusJ + | otherwise = + oneOverIplusJ + where + oneOverIplusJ = + recip $ (AERN.setMinGranularityOuter gran $ iRA + jRA + 1) + iRA = fromInteger $ toInteger i + jRA = fromInteger $ toInteger j + + +--invert :: +-- Precision -> +-- () -> +invert getMatrix = + do + gaussElim getMatrixI + where + n = testMatrixN + getMatrixI gran = + do + m <- getMatrix gran + mI <- addIdentity m + return mI + +gaussElim getMatrix = + elimWithMinGran initialGran + where + elimWithMinGran workingGran = + do + mI@(Matrix colN rowN _) <- getMatrix workingGran + idPerm <- MAr.newListArray (1,rowN) [1..rowN] + elimAtRow mI 1 idPerm + where + elimAtRow mI@(Matrix colN rowN mIarr) i perm = + do + success <- ensureNonZeroDiag -- make sure (i,i) is non-zero by permuting + case success of + False -> -- failed - all elements contain 0 -> try larger granularity + unsafePrint ("failed to divide at granularity " ++ show workingGran) $ + elimWithMinGran (incrementGran workingGran) + True -> + do + normaliseRow + eliminateColumn + case i == rowN of + True -> + do + mInv <- permuteRowsDropCols perm testMatrixN mI + mPrec <- getMatrixPrecision mInv + case mPrec >= targetPrec of + False -> -- resulting precision insufficient + unsafePrint + ("insufficient precision " ++ show mPrec ++ + " at granularity " ++ show workingGran) $ + elimWithMinGran (incrementGran workingGran) + True -> + unsafePrint + ("precision " ++ show mPrec ++ + " succeeded at granularity " ++ show workingGran) + return mInv + False -> elimAtRow mI (i+1) perm + where + ensureNonZeroDiag = + do + maybeNonZeroIx <- findNonZeroRow + case maybeNonZeroIx of + Nothing -> + return False + Just ii -> + do + case ii > 0 of + True -> swap i (i + ii) perm + False -> return () + return True + findNonZeroRow = + do + elems <- mapM getElemPerm [(i,rowIx) | rowIx <- [i..rowN]] + return $ List.findIndex (\e -> not $ 0 `AERN.refines` e) elems + getElemPerm (colIx,rowIx) = + do + rowIxPerm <- unsafePermRead perm rowIx + unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) + + normaliseRow = + do + rowIxPerm <- unsafePermRead perm i + e <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) + unsafeMatrixWrite mIarr rowN (i, rowIxPerm) 1 + mapM (divideCellBy e rowIxPerm) [(i+1)..colN] + divideCellBy e rowIxPerm colIx = + do + e2 <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) + unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (e2/e) + + eliminateColumn = + do + iRowPerm <- unsafePermRead perm i + mapM (eliminateColumnRow iRowPerm) $ [1..(i-1)] ++ [(i+1)..rowN] + eliminateColumnRow iRowPerm rowIx = + do + rowIxPerm <- unsafePermRead perm rowIx + c <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) -- remember old element for scaling i'th row + unsafeMatrixWrite mIarr rowN (i,rowIxPerm) 0 -- at column i we set 0 + mapM (eliminateColumnRowColumn iRowPerm rowIxPerm c) [(i+1)..colN] + eliminateColumnRowColumn iRowPerm rowIxPerm c colIx = + do + ei <- unsafeMatrixRead mIarr rowN (colIx, iRowPerm) -- at i'th row + er <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) -- at current row + unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (er - c * ei) -- eliminate by i'th row + + +swap :: + Int -> + Int -> + (STAr.STUArray s Int Int) -> + ST s () +swap i1 i2 perm = + do + a1 <- unsafePermRead perm i1 + a2 <- unsafePermRead perm i2 + unsafePermWrite perm i1 a2 + unsafePermWrite perm i2 a1 + + +unsafePermWrite permArr i e = + do + BAr.unsafeWrite permArr (i - 1) e + +unsafePermRead permArr i = + do + BAr.unsafeRead permArr (i - 1) + + +addIdentity :: + (STMatrix s IRA) -> + ST s (STMatrix s IRA) +addIdentity (Matrix colN rowN marr) = + do +-- (_, (colN,rowN)) <- MAr.getBounds marr + mElems <- MAr.getElems marr + mIarr <- MAr.newListArray ((1,1),(colN+rowN,rowN)) $ mElems ++ (idElems rowN) + return $ Matrix (colN + rowN) rowN mIarr + where + idElems m = + 1 : (concat $ replicate (m-1) $ (replicate m 0) ++ [1]) + + +data Matrix marr el = + Matrix + { + mxRowN :: Int, + mxColN :: Int, + mxRows :: marr (ColIx,RowIx) el + } + +type ColIx = Int +type RowIx = Int + +type IMatrix el = + Matrix Array el + +type STMatrix s el = + Matrix (STArray s) el + +instance + (IAr.IArray marr el,-- IAr.IArray marr (marr Int el), + Show el) => + Show (Matrix marr el) + where + show (Matrix colN rowN rows) = + "\nMatrix:\n" ++ + (concat $ map showCol [1..colN]) + where +-- (_,(colN,rowN)) = IAr.bounds rows + showCol colIx = + unlines $ + map showCell [(colIx, rowIx) | rowIx <- [1..rowN]] + showCell ix@(colIx, rowIx) = + (show ix) ++ + (replicate colIx '.') ++ + (show $ (IAr.!) rows ix) + +getMatrixPrecision (Matrix _ _ marr) = + do + elems <- MAr.getElems marr + return $ foldl1 min $ map AERN.getPrecision elems + +unsafeMatrixWrite marr rowN (i,j) e = + do + BAr.unsafeWrite marr (rowN*(i-1) + j-1) e +-- MAr.writeArray marr (i,j) e + +unsafeMatrixRead marr rowN (i,j) = + do + BAr.unsafeRead marr (rowN*(i-1) + j-1) +-- MAr.readArray marr (i,j) + +permuteRowsDropCols :: + (STAr.STUArray s Int Int) -> + Int {-^ drop this many first columns -} -> + (STMatrix s IRA) -> + ST s (STMatrix s IRA) +permuteRowsDropCols perm dropN (Matrix colN rowN marr) = + do +-- (_, (colN,rowN)) <- MAr.getBounds marr + (_, permN) <- MAr.getBounds perm + rarr <- MAr.newArray ((1,1),(colN - dropN, permN)) 0 + mapM (copyElem marr rarr rowN) [(colIx, rowIx) | colIx <- [1..colN - dropN], rowIx <- [1..permN]] + return (Matrix (colN - dropN) permN rarr) + where + copyElem marr rarr rowN (colIx, rowIx) = + do + permRowIx <- unsafePermRead perm rowIx + e <- unsafeMatrixRead marr rowN (colIx + dropN, permRowIx) + unsafeMatrixWrite rarr rowN (colIx, rowIx) e + + +addM m1 m2 + | mxColN m1 == mxColN m2 && mxRowN m1 == mxRowN m2 = + do + marr <- MAr.newArray ((1,1),(colN, rowN)) 0 + mapM (addCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]] + return (Matrix colN rowN marr) + | otherwise = + error "Matrix: addM mismatch" + where + colN = mxColN m1 + rowN = mxRowN m1 + marr1 = mxRows m1 + marr2 = mxRows m2 + addCell marr (colIx, rowIx) = + do + elem1 <- unsafeMatrixRead marr1 rowN (colIx, rowIx) + elem2 <- unsafeMatrixRead marr2 rowN (colIx, rowIx) + unsafeMatrixWrite marr rowN (colIx, rowIx) (elem1 + elem2) + +multM m1 m2 + | colN1 == rowN2 = + do + marr <- MAr.newArray ((1,1),(colN, rowN)) 0 + mapM (multCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]] + return (Matrix colN rowN marr) + | otherwise = + error "Matrix: multM mismatch" + where + colN1 = mxColN m1 + rowN1 = mxRowN m1 + colN2 = mxColN m2 + rowN2 = mxRowN m2 + colN = colN2 + rowN = rowN1 + marr1 = mxRows m1 + marr2 = mxRows m2 + multCell marr (colIx, rowIx) = + do + elems1 <- mapM (getCell1 rowIx) [1..colN1] + elems2 <- mapM (getCell2 colIx) [1..rowN2] + unsafeMatrixWrite marr rowN (colIx, rowIx) (sum $ zipWith (*) elems1 elems2) + getCell1 rowIx colIx = + do + unsafeMatrixRead marr1 rowN1 (colIx, rowIx) + getCell2 rowIx colIx = + do + unsafeMatrixRead marr2 rowN2 (colIx, rowIx) + diff --git a/demos/Pi.hs b/demos/Pi.hs new file mode 100644 --- /dev/null +++ b/demos/Pi.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE UndecidableInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE DeriveDataTypeable #-} +module Main + +where + +import qualified Data.Number.ER.Real as AERN +import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq) +import Data.Number.ER.BasicTypes +import Data.Number.ER.Misc + +import Data.Maybe + +#ifdef USE_MPFR +--type B = AERN.BMPFR -- use MPFR floats +type B = AERN.BAP -- use pure Haskell floats +#else +type B = AERN.BAP -- use pure Haskell floats +--type B = AERN.BMAP -- use combination of double and pure Haskell floats +#endif +type RA = AERN.RA B +type IRA = AERN.IRA B + + +decimalPrec = 1000 +binaryPrec = + fromInteger $ toInteger $ + snd $ AERN.integerBounds $ + (fromInteger decimalPrec :: RA) * (AERN.log 100 10)/(AERN.log 100 2) + +main = + do + AERN.initialiseBaseArithmetic (0 :: RA) + putStrLn $ + show decimalPrec + ++ " decimal digits of pi = \n" + ++ (AERN.showConvergRealSeqAuto binaryPrec pi) + where + pi :: ConvergRealSeq RA + pi = ConvergRealSeq AERN.pi + diff --git a/src/Data/Number/ER/BasicTypes.hs b/src/Data/Number/ER/BasicTypes.hs --- a/src/Data/Number/ER/BasicTypes.hs +++ b/src/Data/Number/ER/BasicTypes.hs @@ -1,6 +1,6 @@ {-| Module : Data.Number.ER.BasicTypes - Description : generic types for exact real number processing + Description : auxiliary types for exact real number processing Copyright : (c) Michal Konecny License : BSD3 @@ -8,12 +8,12 @@ Stability : experimental Portability : portable - generic types for exact real number processing + auxiliary types for exact real number processing -} module Data.Number.ER.BasicTypes where -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI {-| Precision represents an upper bound on the measure of diff --git a/src/Data/Number/ER/BasicTypes/DomainBox.hs b/src/Data/Number/ER/BasicTypes/DomainBox.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes/DomainBox.hs @@ -0,0 +1,192 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-| + Module : Data.Number.ER.BasicTypes.DomainBox + Description : portions of many-dimensional domains + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Abstractions of the 'Box' datatype, often used to represent + sections of multi-dimensional function domains. + + To be imported qualified, usually with prefix DBox. + + VariableID(..) and DomainBox + are usually imported separately and not qualified. +-} +module Data.Number.ER.BasicTypes.DomainBox +( + VariableID(..), + getNVars, + DomainBox(..), + DomainBoxMappable(..), + DomainIntBox(..) +) +where + +import Data.Number.ER.BasicTypes + +import qualified Data.Set as Set +import qualified Data.Map as Map + +import Prelude hiding (lookup) + + +{-| + A class abstracting a type of variable identifiers + for axes in function domains, polynomials etc. +-} +class (Ord varid) => VariableID varid + where + newVarID :: Set.Set varid -> varid + defaultVar :: varid + defaultVar = newVarID Set.empty + showVar :: varid -> String + +getNVars :: (VariableID varid) => Int -> [varid] +getNVars n = + aux (Set.empty) n + where + aux prevVars n + | n > 0 = + aux (Set.insert (newVarID prevVars) prevVars) (n - 1) + | n == 0 = + Set.toAscList $ prevVars + +{-| + A class abstracting a type of many-dimensional points, intervals + or anything indexed by a subset of dimensions. + + More generally, this class abstracts most of 'Data.Map.Map'. +-} +class (VariableID varid) => DomainBox box varid val + | box -> varid val, varid val -> box + where + noinfo :: box + isNoinfo :: box -> Bool + size :: box -> Int + {-| constructor using 'defaultVar' -} + unary :: val -> box + singleton :: varid -> val -> box + toList :: box -> [(varid, val)] + fromList :: [(varid, val)] -> box + toAscList :: box -> [(varid, val)] + fromAscList :: [(varid, val)] -> box +-- toMap :: box -> Map.Map varid val +-- fromMap :: Map.Map varid val -> box + compare :: (val -> val -> Ordering) -> box -> box -> Ordering + adjust :: (val -> val) -> varid -> box -> box + insert :: varid -> val -> box -> box + insertWith :: (val -> val -> val) -> varid -> val -> box -> box + delete :: varid -> box -> box + member :: varid -> box -> Bool + notMember :: varid -> box -> Bool + union :: box -> box -> box + unionWith :: (val -> val -> val) -> box -> box -> box + keys :: box -> [varid] + elems :: box -> [val] + filter :: (val -> Bool) -> box -> box + fold :: (val -> a -> a) -> a -> box -> a + foldWithKey :: (varid -> val -> a -> a) -> a -> box -> a + {-| + for all variables that appear in both boxes, + apply the function and add the result to the list + -} + zipWith :: (val -> val -> a) -> box -> box -> [(varid, a)] + {-| + For all variables that appear in either of the two boxes, + apply the function and add the result to the list. + + Supply the default value when the variable is missing from either box. + -} + zipWithDefault :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] + {-| + For all variables that appear in the first box, + apply the function and add the result to the list. + + Supply the default value when the variable is missing from the second box. + -} + zipWithDefaultSecond :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] + findWithDefault :: val -> varid -> box -> val + {-| + Pick the extents of a single variable in a domain box. + If there is no information for this variable, assume the + variable ranges over the whole real line. + -} + lookup :: + String {-^ identification of caller location to use in error messages -} -> + varid -> + box -> + val + +{-| + A class linking two domain box types that share the + index type so that boxes of the two types can be + converted etc. +-} +class (DomainBox box1 varid val1, DomainBox box2 varid val2) => + DomainBoxMappable box1 box2 varid val1 val2 + where + map :: (val1 -> val2) -> box1 -> box2 + mapWithKey :: (varid -> val1 -> val2) -> box1 -> box2 + intersectionWith :: (val1 -> val2 -> val1) -> box1 -> box2 -> box1 + difference :: box1 -> box2 -> box1 + +{-| + A class abstracting a type of many-dimensional intervals. +-} +class (DomainBox box varid ira) => DomainIntBox box varid ira + | box -> varid ira, varid ira -> box + where + {-| + Check whether the two domains specify the same + interval for each variable that they share. + -} + compatible :: + box -> + box -> + Bool + {-| + Assuming that two domains are compatible, take the + most information from both of the domains about the + ranges of variables. + -} + unify :: + String {-^ identification of caller location to use in error messages -} -> + box -> + box -> + box + {-| + Find the variable with the largest interval + and return it together with the default splitting point + in its domain. + -} + bestSplit :: + box {-^ box considered for splitting -} -> + (varid, (ira, ira)) + {-^ variable with widest domain, its domain and default split point -} + split :: + box {-^ box to split -} -> + varid {-^ direction to split in -} -> + Maybe ira {-^ point to split the domain of variable @varid@ at, if absent use default -} -> + (box, box) + classifyPosition :: + box {-^ domain @d1@ -} -> + box {-^ domain @d2@ -} -> + (Bool, Bool, Bool, Bool) + {-^ + Answers to these (mutually exclusive) questions: + + * is @d1@ outside and /not/ touching @d2@? + + * is @d1@ outside and touching @d2@? + + * is @d1@ intersecting and not inside @d2@? + + * is @d1@ inside @d2@? + -} + diff --git a/src/Data/Number/ER/BasicTypes/DomainBox/IntMap.hs b/src/Data/Number/ER/BasicTypes/DomainBox/IntMap.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes/DomainBox/IntMap.hs @@ -0,0 +1,207 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeSynonymInstances #-} +{-| + Module : Data.Number.ER.BasicTypes.DomainBox.IntMap + Description : implementation of DomainBox based on Data.IntMap + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + A simple implementation of the 'VariableID' and 'DomainBox' classes. +-} +module Data.Number.ER.BasicTypes.DomainBox.IntMap +( + VarID, Box +) +where + +import qualified Data.Number.ER.Real.Approx as RA +import qualified Data.Number.ER.BasicTypes.DomainBox as DBox +import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) + +import Data.Number.ER.Misc + +import qualified Data.IntMap as IMap +import qualified Data.Set as Set + +type VarID = Int +type Box ira = IMap.IntMap ira + +instance VariableID VarID + where + newVarID prevVars + | Set.null prevVars = 0 + | otherwise = + 1 + (Set.findMax prevVars) + showVar v + | v == 0 = "x" + | otherwise = "x" ++ show v + +instance (Show val) => (DomainBox (Box val) VarID val) + where + noinfo = IMap.empty + isNoinfo = IMap.null + size = IMap.size + unary r = IMap.singleton defaultVar r + singleton = IMap.singleton + toList = IMap.toList + fromList = IMap.fromList + toAscList = IMap.toAscList + fromAscList = IMap.fromAscList +-- toMap = id +-- fromMap = id + compare compareVals b1 b2 = + compareListsWith comparePairs (IMap.toList b1) (IMap.toList b2) + where + comparePairs (k1,v1) (k2,v2) = + compareComposeMany + [ + compare k1 k2, + compareVals v1 v2 + ] + + adjust = IMap.adjust + insert = IMap.insert + insertWith = IMap.insertWith + delete = IMap.delete + member = IMap.member + notMember = IMap.notMember + union = IMap.union + unionWith = IMap.unionWith + elems = IMap.elems + keys = IMap.keys + filter = IMap.filter + fold = IMap.fold + foldWithKey = IMap.foldWithKey + zipWith f b1 b2 = + applyF (IMap.toAscList b1) (IMap.toAscList b2) + where + applyF [] _ = [] + applyF _ [] = [] + applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) + | k1 == k2 = + (k1, f v1 v2) : (applyF rest1 rest2) + | k1 < k2 = applyF rest1 bl2 + | otherwise = applyF bl1 rest2 + zipWithDefault defaultValue f b1 b2 = + applyF (IMap.toAscList b1) (IMap.toAscList b2) + where + applyF [] [] = [] + applyF bl1@((k1,v1):rest1) [] = + (k1, f v1 defaultValue) : (applyF rest1 []) + applyF [] bl2@((k2,v2):rest2) = + (k2, f defaultValue v2) : (applyF [] rest2) + applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) + | k1 == k2 = + (k1, f v1 v2) : (applyF rest1 rest2) + | k1 < k2 = + (k1, f v1 defaultValue) : (applyF rest1 bl2) + | otherwise = + (k2, f defaultValue v2) : (applyF bl1 rest2) + zipWithDefaultSecond defaultValue f b1 b2 = + applyF (IMap.toAscList b1) (IMap.toAscList b2) + where + applyF [] _ = [] + applyF bl1@((k1,v1):rest1) [] = + (k1, f v1 defaultValue) : (applyF rest1 []) + applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) + | k1 == k2 = + (k1, f v1 v2) : (applyF rest1 rest2) + | k1 < k2 = + (k1, f v1 defaultValue) : (applyF rest1 bl2) + | otherwise = + applyF bl1 rest2 + findWithDefault = IMap.findWithDefault + lookup locspec var dom = + IMap.findWithDefault err var dom + where + err = + error $ + locspec ++ "DomainBox.IntMap lookup: domain box " ++ show dom + ++ " ignores variable " ++ show var + +instance (Show val1, Show val2) => + (DomainBoxMappable (Box val1) (Box val2) VarID val1 val2) + where + map = IMap.map + mapWithKey = IMap.mapWithKey + intersectionWith = IMap.intersectionWith + difference = IMap.difference + +instance (RA.ERIntApprox ira) => DomainIntBox (Box ira) VarID ira + where + compatible dom1 dom2 = + foldl (&&) True $ map snd $ + DBox.zipWith RA.equalIntervals dom1 dom2 + unify locspec dom1 dom2 + | DBox.compatible dom1 dom2 = + IMap.union dom1 dom2 + | otherwise = + error $ + locspec ++ "incompatible domains " ++ show dom1 ++ " and " ++ show dom2 + bestSplit domB = + (var, (varDom, pt)) + where + pt = + RA.defaultBisectPt varDom + (_, (varDom, var)) = + foldl findWidestVar (0, err) $ IMap.toList domB + err = + error $ "DomainBox: bestSplit: failed to find a split for " ++ show domB + findWidestVar (prevWidth, prevRes) (v, d) + | currWidth `RA.leqSingletons` prevWidth = (prevWidth, prevRes) + | otherwise = (currWidth, (d, v)) + where + currWidth = snd $ RA.bounds $ domHI - domLO + (domLO, domHI) = RA.bounds d + split domB var maybePt = + (IMap.insert var varDomL domB, + IMap.insert var varDomR domB) + where + varDomL = varDomLO RA.\/ pt + varDomR = pt RA.\/ varDomHI + pt = + case maybePt of + Nothing -> varDomMid + Just pt | pt `RA.refines` varDom -> pt + Just pt -> + error $ + "ER.DomainBox.IntMap: split given an invalid split point " + ++ show pt ++ " for the domain box " ++ show domB + ++ " and split variable " ++ show var + (varDomLO, varDomMid, varDomHI, _) = RA.exactMiddle varDom + varDom = DBox.lookup "DomainBox.IntMap: split: " var domB + classifyPosition dom sdom = + (away, touch, intersect, inside) + where + (away, touch, inside, intersect) = + foldl addDimension (True, True, True, False) awayTouchInsides + addDimension + (prevAway, prevTouch, prevInside, prevIntersect) + (thisAway, thisTouch, thisInside, thisIntersect) = + (prevAway && thisAway, + (prevTouch || prevAway) && (thisTouch || thisAway) && (prevTouch || thisTouch), + prevInside && thisInside, + prevIntersect || thisIntersect) + awayTouchInsides = + map snd $ + DBox.zipWith classifyRA dom sdom + classifyRA d sd = + (outsideNoTouch, outsideTouch, inside, + not (outsideNoTouch || outsideTouch || inside)) + where + outsideNoTouch = sdR < dL || dR < sdL + outsideTouch = sdR == dL || dR == sdL + inside = sdL =< dL && dR =< sdR + (==) = RA.eqSingletons + (<) = RA.ltSingletons + (=<) = RA.leqSingletons + (dL, dR) = RA.bounds d + (sdL, sdR) = RA.bounds sd + + + diff --git a/src/Data/Number/ER/BasicTypes/ExtendedInteger.hs b/src/Data/Number/ER/BasicTypes/ExtendedInteger.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes/ExtendedInteger.hs @@ -0,0 +1,125 @@ +{-| + Module : Data.Number.ER.BasicTypes.ExtendedInteger + Description : integer with infinities + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + An arbitrary sized integer type with additional +infinity and -infinity. + + To be imported qualified, usually with prefix EI. +-} +module Data.Number.ER.BasicTypes.ExtendedInteger +( + ExtendedInteger(..), + isInfinite, binaryLog, take +) +where + +import Prelude hiding (isInfinite, take) +import qualified Prelude + +data ExtendedInteger + = MinusInfinity | Finite Integer | PlusInfinity + deriving (Eq) + +isInfinite :: ExtendedInteger -> Bool +isInfinite MinusInfinity = True +isInfinite PlusInfinity = True +isInfinite _ = False + +{-| + the smallest integer i for which 2^i <= abs n +-} +binaryLog :: ExtendedInteger -> ExtendedInteger +binaryLog PlusInfinity = PlusInfinity +binaryLog MinusInfinity = PlusInfinity +binaryLog (Finite n) + | n < 0 = binaryLog (Finite (- n)) + | n == 0 = MinusInfinity + | otherwise = -- (n > 0) + -- how to do this fast? + intBinaryLog n + +intBinaryLog n + | n > 1 = 1 + (intBinaryLog (n `div` 2)) + | n == 1 = 0 + +instance Show ExtendedInteger where + show MinusInfinity = "-InfInt" + show PlusInfinity = "+InfInt" + show (Finite i) = show i + +take :: ExtendedInteger -> [a] -> [a] +take MinusInfinity _ = error "takeEI called with MinusInfinity" +take PlusInfinity list = list +take (Finite n) list = Prelude.take (fromInteger n) list + +instance Ord ExtendedInteger where + compare MinusInfinity MinusInfinity = EQ + compare MinusInfinity _ = LT + compare _ MinusInfinity = GT + compare PlusInfinity PlusInfinity = EQ + compare PlusInfinity _ = GT + compare _ PlusInfinity = LT + compare (Finite i1) (Finite i2) = + compare i1 i2 + +instance Num ExtendedInteger where + fromInteger i = Finite i + {- abs -} + abs MinusInfinity = PlusInfinity + abs PlusInfinity = PlusInfinity + abs (Finite i) = Finite $ abs i + {- signum -} + signum ei + | ei < 0 = -1 + | ei > 0 = 1 + | otherwise = 0 + {- negate -} + negate (Finite i) = Finite (-i) + negate MinusInfinity = PlusInfinity + negate PlusInfinity = MinusInfinity + {- addition -} + PlusInfinity + MinusInfinity = + error "cannot add PlusInfinity and MinusInfinity" + MinusInfinity + PlusInfinity = + error "cannot add PlusInfinity and MinusInfinity" + PlusInfinity + ei = PlusInfinity + ei + PlusInfinity = PlusInfinity + MinusInfinity + ei = MinusInfinity + ei + MinusInfinity = MinusInfinity + (Finite i1) + (Finite i2) = Finite $ i1 + i2 + {- multiplication -} + ei1 * ei2 | ei1 > ei2 = ei2 * ei1 + MinusInfinity * ei + | ei < 0 = PlusInfinity + | ei > 0 = MinusInfinity + | otherwise = error "cannot multiply MinusInfinity and 0" + ei * PlusInfinity + | ei < 0 = MinusInfinity + | ei > 0 = PlusInfinity + | otherwise = error "cannot multiply PlusInfinity and 0" + (Finite i1) * (Finite i2) = Finite $ i1 * i2 + +instance Enum ExtendedInteger where + toEnum i = Finite $ toInteger i + fromEnum (Finite i) = fromInteger i + fromEnum _ = error "infinite integers cannot be enumerated" + +instance Real ExtendedInteger where + toRational (Finite i) = toRational i + toRational _ = error "infinite integers cannot be converted to rational" + +instance Integral ExtendedInteger where + quotRem (Finite i) (Finite m) = + (Finite a, Finite b) + where + (a,b) = quotRem i m + quotRem _ _ = error "cannot make a quotient involving an infinite integer" + toInteger (Finite i) = i + toInteger _ = error "infinite integers cannot be converted to Integer" + diff --git a/src/Data/Number/ER/BasicTypes/PlusMinus.hs b/src/Data/Number/ER/BasicTypes/PlusMinus.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes/PlusMinus.hs @@ -0,0 +1,47 @@ +{-# LANGUAGE DeriveDataTypeable #-} +{-| + Module : Data.Number.ER.BasicTypes.PlusMinus + Description : mini sign datatype + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + A mini enumeration to represent the sign of different numbers and approximations. +-} +module Data.Number.ER.BasicTypes.PlusMinus where + +import Data.Typeable +import Data.Generics.Basics +import Data.Binary +--import BinaryDerive + +data PlusMinus = Minus | Plus + deriving (Eq, Ord, Typeable, Data) + +instance Show PlusMinus where + show Plus = "+" + show Minus = "-" + +{- the following has been generated by BinaryDerive -} +instance Binary PlusMinus where + put Minus = putWord8 0 + put Plus = putWord8 1 + get = do + tag_ <- getWord8 + case tag_ of + 0 -> return Minus + 1 -> return Plus + _ -> fail "no parse" +{- the above has been generated by BinaryDerive -} + +signNeg Plus = Minus +signNeg Minus = Plus + +signMult Plus s = s +signMult Minus s = signNeg s + +signToNum Plus = 1 +signToNum Minus = -1 diff --git a/src/Data/Number/ER/BasicTypes/Tests/Generate.hs b/src/Data/Number/ER/BasicTypes/Tests/Generate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/BasicTypes/Tests/Generate.hs @@ -0,0 +1,92 @@ +{-| + Module : Data.Number.ER.BasicTypes.Tests.Generate + Description : (testing) generating values for tests + Copyright : (c) 2007-2008 Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Instances of Arbitrary so that values + can be randomly generated for use in QuickCheck tests. +-} + +module Data.Number.ER.BasicTypes.Tests.Generate +where + +import Data.Number.ER.BasicTypes +import Data.Number.ER.BasicTypes.ExtendedInteger +import Data.Number.ER.BasicTypes.PlusMinus + +import Test.QuickCheck hiding (two, three) + +data Nat10 = Nat10 Int deriving (Show) +data Nat100 = Nat100 Int deriving (Show) + +data Ix10 = Ix10 EffortIndex deriving (Show) +data Ix20 = Ix20 EffortIndex deriving (Show) + +data Gran100 = Gran100 Granularity deriving (Show) +data Gran1000 = Gran1000 Granularity deriving (Show) + +data SmallRatio = SmallRatio Int Int deriving (Show) + +instance (Arbitrary Nat10) + where + arbitrary = + do + ix <- choose (0,10) + return $ Nat10 ix + coarbitrary (Nat10 ix) = + error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Nat10" + +instance (Arbitrary Nat100) + where + arbitrary = + do + ix <- choose (0,100) + return $ Nat100 ix + coarbitrary (Nat100 ix) = + error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Nat100" + +instance (Arbitrary Ix20) + where + arbitrary = + do + ix <- choose (2,20) + return $ Ix20 ix + coarbitrary (Ix20 ix) = + error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Ix20" + +instance (Arbitrary Ix10) + where + arbitrary = + do + ix <- choose (1,10) + return $ Ix10 ix + coarbitrary (Ix10 ix) = + error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for Ix10" + +instance (Arbitrary PlusMinus) + where + arbitrary = + do + isPlus <- arbitrary + case isPlus of + True -> return Plus + False -> return Minus + coarbitrary pm = + error "ER.BasicTypes.Tests.Generate: coarbitrary not implemented for PlusMinus" + +instance (Arbitrary SmallRatio) + where + arbitrary = + do + num <- choose (-1000000,1000000) + denom <- choose (1,1000000) + return $ SmallRatio num denom + coarbitrary (SmallRatio num denom) = + error "ERChebPoly: Generate: Arbitrary: coarbitrary not implemented for SmallRatio" + + diff --git a/src/Data/Number/ER/ExtendedInteger.hs b/src/Data/Number/ER/ExtendedInteger.hs deleted file mode 100644 --- a/src/Data/Number/ER/ExtendedInteger.hs +++ /dev/null @@ -1,125 +0,0 @@ -{-| - Module : Data.Number.ER.ExtendedInteger - Description : integer with infinities - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mik@konecny.aow.cz - Stability : experimental - Portability : portable - - An arbitrary sized integer type with additional +infinity and -infinity. - - To be imported qualified, usually with prefix EI. --} -module Data.Number.ER.ExtendedInteger -( - ExtendedInteger(..), - isInfinite, binaryLog, take -) -where - -import Prelude hiding (isInfinite, take) -import qualified Prelude - -data ExtendedInteger - = MinusInfinity | Finite Integer | PlusInfinity - deriving (Eq) - -isInfinite :: ExtendedInteger -> Bool -isInfinite MinusInfinity = True -isInfinite PlusInfinity = True -isInfinite _ = False - -{-| - the smallest integer i for which 2^i <= abs n --} -binaryLog :: ExtendedInteger -> ExtendedInteger -binaryLog PlusInfinity = PlusInfinity -binaryLog MinusInfinity = PlusInfinity -binaryLog (Finite n) - | n < 0 = binaryLog (Finite (- n)) - | n == 0 = MinusInfinity - | otherwise = -- (n > 0) - -- how to do this fast? - intBinaryLog n - -intBinaryLog n - | n > 1 = 1 + (intBinaryLog (n `div` 2)) - | n == 1 = 0 - -instance Show ExtendedInteger where - show MinusInfinity = "-InfInt" - show PlusInfinity = "+InfInt" - show (Finite i) = show i - -take :: ExtendedInteger -> [a] -> [a] -take MinusInfinity _ = error "takeEI called with MinusInfinity" -take PlusInfinity list = list -take (Finite n) list = Prelude.take (fromInteger n) list - -instance Ord ExtendedInteger where - compare MinusInfinity MinusInfinity = EQ - compare MinusInfinity _ = LT - compare _ MinusInfinity = GT - compare PlusInfinity PlusInfinity = EQ - compare PlusInfinity _ = GT - compare _ PlusInfinity = LT - compare (Finite i1) (Finite i2) = - compare i1 i2 - -instance Num ExtendedInteger where - fromInteger i = Finite i - {- abs -} - abs MinusInfinity = PlusInfinity - abs PlusInfinity = PlusInfinity - abs (Finite i) = Finite $ abs i - {- signum -} - signum ei - | ei < 0 = -1 - | ei > 0 = 1 - | otherwise = 0 - {- negate -} - negate (Finite i) = Finite (-i) - negate MinusInfinity = PlusInfinity - negate PlusInfinity = MinusInfinity - {- addition -} - PlusInfinity + MinusInfinity = - error "cannot add PlusInfinity and MinusInfinity" - MinusInfinity + PlusInfinity = - error "cannot add PlusInfinity and MinusInfinity" - PlusInfinity + ei = PlusInfinity - ei + PlusInfinity = PlusInfinity - MinusInfinity + ei = MinusInfinity - ei + MinusInfinity = MinusInfinity - (Finite i1) + (Finite i2) = Finite $ i1 + i2 - {- multiplication -} - ei1 * ei2 | ei1 > ei2 = ei2 * ei1 - MinusInfinity * ei - | ei < 0 = PlusInfinity - | ei > 0 = MinusInfinity - | otherwise = error "cannot multiply MinusInfinity and 0" - ei * PlusInfinity - | ei < 0 = MinusInfinity - | ei > 0 = PlusInfinity - | otherwise = error "cannot multiply PlusInfinity and 0" - (Finite i1) * (Finite i2) = Finite $ i1 * i2 - -instance Enum ExtendedInteger where - toEnum i = Finite $ toInteger i - fromEnum (Finite i) = fromInteger i - fromEnum _ = error "infinite integers cannot be enumerated" - -instance Real ExtendedInteger where - toRational (Finite i) = toRational i - toRational _ = error "infinite integers cannot be converted to rational" - -instance Integral ExtendedInteger where - quotRem (Finite i) (Finite m) = - (Finite a, Finite b) - where - (a,b) = quotRem i m - quotRem _ _ = error "cannot make a quotient involving an infinite integer" - toInteger (Finite i) = i - toInteger _ = error "infinite integers cannot be converted to Integer" - diff --git a/src/Data/Number/ER/Misc/STM.hs b/src/Data/Number/ER/Misc/STM.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Misc/STM.hs @@ -0,0 +1,42 @@ +{-| + Module : Data.Number.ER.Misc.STM + Description : some STM extras + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Miscelaneous utilities related to concurrency. +-} +module Data.Number.ER.Misc.STM where + +import Control.Concurrent as Concurrent +import Control.Concurrent.STM as STM + +modifyTVar tv update = + do + value <- readTVar tv + let newValue = update value + writeTVar tv newValue + return newValue + +modifyTVarGetOldVal tv update = + do + value <- readTVar tv + writeTVar tv $ update value + return value + +modifyTVarHasChanged tv update = + do + value <- readTVar tv + let newValue = update value + if value == newValue + then return False + else + do + writeTVar tv $ update value + return True + + diff --git a/src/Data/Number/ER/Misc/Tests.hs b/src/Data/Number/ER/Misc/Tests.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Misc/Tests.hs @@ -0,0 +1,54 @@ +{-| + Module : Data.Number.ER.Misc.Tests + Description : some QuickCheck extras + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Miscelaneous utilities related to testing. +-} +module Data.Number.ER.Misc.Tests + +where + +import Data.Number.ER.Misc + +import Test.QuickCheck +import Test.QuickCheck.Batch + +import System.IO + +erRunTests testsetName options initialise tests = + do + mapM (mkRunTest $ length tests) $ zip [1..] tests + return () + where + mkRunTest testCount (n, (testName, test)) = + do + initialise + putStr testDescr + result <- test options + putStrLn $ " result: " ++ show result +-- runTests testDescr options [test] + hFlush stdout + where + testDescr = + "(" ++ show n ++ "/" ++ show testCount ++ ") " ++ testsetName ++ ": " ++ testName ++ "\n" + +instance Show TestResult + where + show result = + case result of + TestOk msg ntest stamps -> + msg ++ " " ++ show ntest ++ " " -- ++ show stamps + TestExausted msg ntest stamps -> + msg ++ " " ++ show ntest ++ " " -- ++ show stamps + TestAborted exception -> + "aborted: " ++ show exception + TestFailed args ntest -> + "failed after " ++ show ntest ++ " tests" + ++ "\n args = " ++ show args + diff --git a/src/Data/Number/ER/MiscSTM.hs b/src/Data/Number/ER/MiscSTM.hs deleted file mode 100644 --- a/src/Data/Number/ER/MiscSTM.hs +++ /dev/null @@ -1,31 +0,0 @@ -{-| - Utilities related to concurrency. --} -module Data.Number.ER.MiscSTM where - -import Control.Concurrent as Concurrent -import Control.Concurrent.STM as STM - -modifyTVar tv update = - do - value <- readTVar tv - writeTVar tv $ update value - -modifyTVarGetOldVal tv update = - do - value <- readTVar tv - writeTVar tv $ update value - return value - -modifyTVarHasChanged tv update = - do - value <- readTVar tv - let newValue = update value - if value == newValue - then return False - else - do - writeTVar tv $ update value - return True - - diff --git a/src/Data/Number/ER/PlusMinus.hs b/src/Data/Number/ER/PlusMinus.hs deleted file mode 100644 --- a/src/Data/Number/ER/PlusMinus.hs +++ /dev/null @@ -1,47 +0,0 @@ -{-# LANGUAGE DeriveDataTypeable #-} -{-| - Module : Data.Number.ER.PlusMinus - Description : mini sign datatype - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mik@konecny.aow.cz - Stability : experimental - Portability : portable - - A mini enumeration to represent the sign of different numbers and approximations. --} -module Data.Number.ER.PlusMinus where - -import Data.Typeable -import Data.Generics.Basics -import Data.Binary ---import BinaryDerive - -data PlusMinus = Minus | Plus - deriving (Eq, Ord, Typeable, Data) - -instance Show PlusMinus where - show Plus = "+" - show Minus = "-" - -{- the following has been generated by BinaryDerive -} -instance Binary PlusMinus where - put Minus = putWord8 0 - put Plus = putWord8 1 - get = do - tag_ <- getWord8 - case tag_ of - 0 -> return Minus - 1 -> return Plus - _ -> fail "no parse" -{- the above has been generated by BinaryDerive -} - -signNeg Plus = Minus -signNeg Minus = Plus - -signMult Plus s = s -signMult Minus s = signNeg s - -signToNum Plus = 1 -signToNum Minus = -1 diff --git a/src/Data/Number/ER/Real/Approx.hs b/src/Data/Number/ER/Real/Approx.hs --- a/src/Data/Number/ER/Real/Approx.hs +++ b/src/Data/Number/ER/Real/Approx.hs @@ -16,47 +16,70 @@ approximations: * 'ERApprox' = - a *set* of approximated numbers whose size is - measured using some fixed measure - - * 'ERIntApprox' = - an *interval* of real numbers with finitely - representable endpoints + Approximating a real number by a *set* of real numbers + that includes the approximated number. + Precision is measured using some fixed measure on the sets. + Operations are "safe" wrt inclusion. + The sets can sometimes be "anti-consistent" - being smaller than + the empty set in the inclusion order. + + * 'ERInnerOuterApprox' = + Like 'ERApprox' with the addition of operations that are "inner rounded" + in the sense that each element of the rounded result set can + be obtained by the same operation performed on some elements of the arument set(s). + + * 'ERIntApprox' = + Like ERApprox but assuming that the sets are + *intervals* of real numbers with finitely + representable endpoints. To be imported qualified, usually with the synonym RA. -} module Data.Number.ER.Real.Approx ( ERApprox(..), - ERIntApprox(..), - effIx2ra, - splitIRA, --- checkShrinking, eqSingletons, leqSingletons, ltSingletons, + effIx2ra, + ERInnerOuterApprox(..), + ERIntApprox(..), + splitIRA, equalIntervals, exactMiddle, maxExtensionR2R, - maxExtensionInnerR2R + maxExtensionInnerR2R, + ERApproxApprox(..) ) where import Data.Number.ER.BasicTypes -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Typeable {-| - A type whose elements represent sets that can be used - to approximate a single extended real number with arbitrary precision. + A type whose elements represent sets that can be used + to approximate a single extended real number with arbitrary precision. + + Operations are "safe" with respect to inclusion, which means that + for any numbers admitted by the operand approximations the result + of the operation is admitted by the result approximation. + + The sets can sometimes be "anti-consistent" - being smaller than + the empty set in the inclusion order. + This can be understood as indicating that not only there is no correct real number + approximated here, but some numbers (ie those in interior of the set) + are excluded more strongly than the others. + Prime examples of such sets are directed "inverted" intervals such as [2,1]. + Such sets arise naturally from "inner rounded" operations - see 'ERInnerOuterApprox'. -} class (Fractional ra) => ERApprox ra where initialiseBaseArithmetic :: ra -> IO () getPrecision :: ra -> Precision {-^ - Precision is a measure of the set size. + Precision is a measure of the set size. It can be infinite. The default interpretation: @@ -65,53 +88,69 @@ -} getGranularity :: ra -> Granularity -- ^ the lower the granularity the bigger the rounding errors - setGranularity :: Granularity -> ra -> ra + setGranularityOuter :: Granularity -> ra -> ra -- ^ increase or safely decrease granularity - setMinGranularity :: Granularity -> ra -> ra + setMinGranularityOuter :: Granularity -> ra -> ra -- ^ ensure granularity is not below the first arg - isEmpty :: ra -> Bool - -- ^ true if this represents a computational error isBottom :: ra -> Bool - -- ^ true if this holds no information + -- ^ true if this approximation holds no information, ie it admits any real number + bottomApprox :: ra + -- ^ the bottom approximation - it admits any real number isExact :: ra -> Bool - -- ^ true if this is a singleton + -- ^ true if this approximation admits only one real number + isConsistent :: ra -> Bool + {- ^ true iff this approximation admits at least one real number -} + isAnticonsistent :: ra -> Bool + {- ^ true if this approximation is anti-consistent, which is a computational error + unless we used inner rounded operations -} + toggleConsistency :: ra -> ra + {- ^ + Toggle consistency - anti-consistency of the approximation. + Top is toggled with bottom. + Exact approximations are the only fixed points for this operation. + -} + isTop :: ra -> Bool + -- ^ true if this approximation is the most anti-consistent one + topApprox :: ra + -- ^ the top approximation - strongly rejects all real numbers isDisjoint :: ra -> ra -> Bool - isDisjoint a b = isEmpty $ a /\ b + isDisjoint a b = not $ isConsistent $ a /\ b isInteriorDisjoint :: ra -> ra -> Bool - isInteriorDisjoint a b = - isEmpty isect || isExact isect - where - isect = a /\ b - isBounded :: ra -> Bool - -- ^ true if the approximation excludes infinity - bottomApprox :: ra - -- ^ the bottom element - any number - emptyApprox :: ra - -- ^ the top element - error + isInteriorDisjoint a b = isAnticonsistent $ a /\ b + isBounded :: ra -> Bool + {- ^ + True iff the approximation excludes infinity + and, if anti-consistent, does not strongly exclude infinity. + -} + plusInfinity :: ra + -- ^ an exact approximation admitting only the positive infinity refines :: ra -> ra -> Bool -- ^ first arg is a subset of the second arg + maybeRefines :: ra -> ra -> Maybe Bool + -- ^ like 'refines' but usable for types where 'refines' is only partially decidable (/\) :: ra -> ra -> ra - -- ^ join; combining two approximations of the same number + -- ^ join; combining the information in two approximations of the same number intersectMeasureImprovement :: EffortIndex -> ra -> ra -> (ra, ra) {-^ - Like intersection but the second component: + First component of result is the intersection and the second component: - * measures improvement of the intersection relative to the first of the two approximations + * measures precision improvement of the intersection relative to the first argument * is a positive number: 1 means no improvement, 2 means doubled precision, etc. -} - equalReals :: ra -> ra -> Maybe Bool - -- ^ nothing if overlapping and not singletons + equalReals :: ra -> ra -> Maybe Bool + -- ^ semantic semi-decidable equality test compareReals :: ra -> ra -> Maybe Ordering - -- ^ nothing if overlapping and not singletons + -- ^ semantic semi-decidable comparison leqReals :: ra -> ra -> Maybe Bool - -- ^ nothing if overlapping on interior or by a wrong endpoint + -- ^ semantic semi-decidable less-than-or-equal comparison equalApprox :: ra -> ra -> Bool - -- ^ syntactic comparison + -- ^ syntactic equality test compareApprox :: ra -> ra -> Ordering -- ^ syntactic linear ordering double2ra :: Double -> ra + -- ^ safe approximate conversion showApprox :: Int {-^ number of relevant decimals to show -} -> Bool {-^ should show granularity -} -> @@ -120,26 +159,99 @@ String {-| - For a finite sequence of real approximations, determine - whether it is a shrinking sequence. --} -checkShrinking :: - (ERApprox ra) => - [ra] -> - Maybe (ra, ra) -checkShrinking [] = Nothing -checkShrinking [_] = Nothing -checkShrinking (a : b : rest) - | b `refines` a = checkShrinking (b : rest) - | otherwise = Just (a,b) + Assuming the arguments are singletons, equality is decidable. +-} +eqSingletons :: (ERApprox ra) => ra -> ra -> Bool +eqSingletons s1 s2 = + case equalReals s1 s2 of + Just b -> b + _ -> False +{-| + Assuming the arguments are singletons, @<=@ is decidable. +-} +leqSingletons :: (ERApprox ra) => ra -> ra -> Bool +leqSingletons s1 s2 = + case compareReals s1 s2 of + Just EQ -> True + Just LT -> True + _ -> False {-| + Assuming the arguments are singletons, @<@ is decidable. +-} +ltSingletons :: (ERApprox ra) => ra -> ra -> Bool +ltSingletons s1 s2 = + case compareReals s1 s2 of + Just LT -> True + _ -> False + +{-| + This function converts + an effort index to a real number approximation. + + Useful when an effort index is used in a formula + mixed with real approximations. +-} +effIx2ra :: + (ERApprox ra) => + EffortIndex -> ra +effIx2ra = fromInteger . toInteger + +{-| + A type whose elements represent some kind of nominal sets of real numbers + over which one can perform two kinds of arithmetic: + + * "outer rounded": arithmetic that approximates maximal extensions from outside (ie the 'ERApprox' arithmetic) + + * "inner rounded": arithmetic that approximates maximal extensions from inside, potentially leading to + anti-consistent set specifications (eg intervals whose endpoints are not in the usual order) + + Another explanation of the difference: + + * `outer': the approximation contains all the number(s) of interest + * `inner': all numbers eligible for the approximation are numbers of interest + + Ie inner rounded operations have the property that each real number admitted by the result can + be obtained as the exact result of the same operation performed on some real numbers admitted + by the operand approximations. + + While in "outer rounded" operations it is desirable to make the result set as small as + possible in order to reduce the amount of bogus result numbers, + in "inner rounded" operations it is desirable to make the result set as large as possible + to lose less of the genuinely feasible result numbers. + + Inner rounded arithmetic is useful eg for proving/disproving inclusions "f(x) subset g(x)" + where f and g are expressions using arithmetic extended to sets. + For proving the inclusion, we need an inner rounded approximation of g(x) + and for disproving the inclusion we need an inner rounded approximation of f(x). + + This is an abstraction of Kaucher's extended interval arithmetic + [Kaucher, E.: Interval Analysis in the Extended Interval Space IR, + Computing, Suppl. 2, 1980, pp. 33-49]. +-} +class (ERApprox xra) => ERInnerOuterApprox xra + where + (+:) :: xra -> xra -> xra + -- ^ inner rounded addition + (-:) :: xra -> xra -> xra + -- ^ inner rounded subtraction + a -: b = a +: (negate b) + (*:) :: xra -> xra -> xra + -- ^ inner rounded multiplication + (/:) :: xra -> xra -> xra + -- ^ inner rounded division + setGranularityInner :: Granularity -> xra -> xra + -- ^ increase or safely decrease granularity + setMinGranularityInner :: Granularity -> xra -> xra + -- ^ ensure granularity is not below the first arg + +{-| A type whose elements represent sets that can be used to approximate a recursive set of closed extended real number intervals with arbitrary precision. -} ---class (ERApprox sra) => SetOfRealsApprox sra where +--class (ERApprox sra) => ERSetApprox sra where -- (\/) :: sra -> sra -> sra -- ^ union; either approximation could be correct {-| @@ -174,34 +286,6 @@ (\/) :: ira -> ira -> ira {-| - Assuming the arguments are singletons, equality is decidable. --} -eqSingletons :: (ERApprox ra) => ra -> ra -> Bool -eqSingletons s1 s2 = - case equalReals s1 s2 of - Just b -> b - _ -> False - -{-| - Assuming the arguments are singletons, @<=@ is decidable. --} -leqSingletons :: (ERApprox ra) => ra -> ra -> Bool -leqSingletons s1 s2 = - case compareReals s1 s2 of - Just EQ -> True - Just LT -> True - _ -> False - -{-| - Assuming the arguments are singletons, @<@ is decidable. --} -ltSingletons :: (ERApprox ra) => ra -> ra -> Bool -ltSingletons s1 s2 = - case compareReals s1 s2 of - Just LT -> True - _ -> False - -{-| Return true if and only if the two intervals have equal endpoints. -} equalIntervals :: @@ -214,18 +298,6 @@ (d2L, d2U) = bounds d2 -{-| - This function converts - an effort index to a real number approximation. - - Useful when an effort index is used in a formula - mixed with real approximations. --} -effIx2ra :: - (ERApprox ra) => - EffortIndex -> ra -effIx2ra = fromInteger . toInteger - {-| Split an interval to a sequence of intervals whose union is the original interval using a given sequence of cut points. @@ -268,8 +340,8 @@ gran = max (getGranularity domL) (getGranularity domR) domM = (domL + domR) / 2 higherGran = gran + 1 - domLhg = setMinGranularity higherGran domL - domRhg = setMinGranularity higherGran domR + domLhg = setMinGranularityOuter higherGran domL + domRhg = setMinGranularityOuter higherGran domR domMhg = (domLhg + domRhg) / 2 @@ -278,16 +350,21 @@ given function. A maximal extension function has the property: f(I) = { f(x) | x in I }. Here we get this property only for the limit function for its 'EffortIndex' tending to infinity. + For finite effor indices the function may add *outer* rounding + but it should be reasonably small. -} maxExtensionR2R :: (ERIntApprox ira) => (EffortIndex -> ira -> [ira]) - {-^ returns a safe approximation of all extrema within the interval -} -> + {-^ returns an *outer* approximation of all extrema within the interval -} -> (EffortIndex -> ira -> ira) - {-^ a function behaving well on sequences that intersect to a point -} -> + {-^ an *outer* rounding function behaving well on sequences that intersect to a point -} -> (EffortIndex -> ira -> ira) - {- ^ a function behaving well on sequences that intersect to a non-empty interval -} + {- ^ an outer rounding function behaving well on sequences that intersect to a non-empty interval -} maxExtensionR2R getExtremes f ix x + | not $ isConsistent x = + toggleConsistency $ + maxExtensionInnerR2R getExtremes f ix $ toggleConsistency x | getPrecision x < effIx2prec ix = foldl1 (\/) $ [f ix xL, f ix xR] ++ (getExtremes ix x) -- x is thin enough (?), don't bother evaluating by endpoints and extrema: @@ -301,47 +378,44 @@ given function. A maximal extension function has the property: f(I) = { f(x) | x in I }. Here we get this property only for the limit function for its 'EffortIndex' tending to infinity. + For finite effor indices the function may include *inner* rounding + but it should be reasonably small. -} maxExtensionInnerR2R :: (ERIntApprox ira) => - (EffortIndex -> ira -> ([ira], (Maybe Bool, Maybe Bool))) - {-^ returns a safe approximation of all extrema within the interval - and an indication whether the function is increasing or decreasing - at the endpoints of the queried real approximation -} -> + (EffortIndex -> ira -> [ira]) + {-^ returns an *outer* approximation of all extrema within the interval -} -> (EffortIndex -> ira -> ira) - {-^ a function behaving well on sequences that intersect to a point -} -> - (EffortIndex -> ira -> ira) - {- ^ a function behaving well on sequences that intersect to a non-empty interval -} -maxExtensionInnerR2R getExtremesAndDirections f ix x = - case (isIncreasing, isDecreasing, compareReals leftVal rightVal) of --- (True, _, Just GT) -> emptyApprox --- (True, _, Nothing) -> emptyApprox --- (_, True, Just LT) -> emptyApprox --- (_, True, Nothing) -> emptyApprox - (True, _, _) -> fromBounds (leftVal, rightVal) - (_, True, _) -> fromBounds (rightVal, leftVal) - _ -> - (/\) ((-1) \/ 1) $ - foldl1 (\/) $ [leftVal, rightVal] ++ extremes + {-^ an *outer* rounding function behaving well on sequences that intersect to a point -} -> + (EffortIndex -> ira -> ira) + {- ^ an inner rounding function behaving well on sequences that intersect to a non-empty interval -} +maxExtensionInnerR2R getExtremes f ix x + | not $ isConsistent x = + toggleConsistency $ + maxExtensionR2R getExtremes f ix $ toggleConsistency x + | otherwise = + foldl1 (\/) $ map toggleConsistency $ [f ix xL, f ix xR] ++ (getExtremes ix x) where - (extremes, (maybeLowIncreasing, maybeHighIncreasing)) = - getExtremesAndDirections ix x - (Just lowIsIncreasing) = maybeLowIncreasing - (isIncreasing, isDecreasing) - | null extremes = (lowIsIncreasing, not lowIsIncreasing) - | otherwise = (False, False) - leftVal = - case maybeLowIncreasing of - Just True -> snd $ bounds $ f ix xL - Just False -> fst $ bounds $ f ix xL - Nothing -> emptyApprox - rightVal = - case maybeHighIncreasing of - Just True -> fst $ bounds $ f ix xR - Just False -> snd $ bounds $ f ix xR - Nothing -> emptyApprox (xL, xR) = bounds x - - - +{-| + A type whose elements are thought of as sets of approximations of real numbers. + + Eg intervals of intervals, eg [[0,3],[1,2]] containing all intervals + whose left endpoint is between 0 and 1 and the right endpoint is between 2 and 3. + The upper bound interval can sometimes be anti-consistent, + eg [[0,3],[2,1]] containing all intervals (consistent as well as anti-consistent) + with a left endpoint between [0,2] and the right endpoint between [1,3]. +-} +class ERApproxApprox xra + where + safeIncludes :: xra -> xra -> Bool + -- ^ safe inclusion of approximations + safeNotIncludes :: xra -> xra -> Bool + -- ^ safe negation of inclusion of approximations + includes :: xra -> xra -> Maybe Bool + -- ^ like 'safeIncludes' but usable for types where 'safeIncludes' is only partially decidable + includes aa1 aa2 + | safeIncludes aa1 aa2 = Just True + | safeNotIncludes aa1 aa2 = Just False + | otherwise = Nothing diff --git a/src/Data/Number/ER/Real/Approx/Elementary.hs b/src/Data/Number/ER/Real/Approx/Elementary.hs --- a/src/Data/Number/ER/Real/Approx/Elementary.hs +++ b/src/Data/Number/ER/Real/Approx/Elementary.hs @@ -12,13 +12,15 @@ -} module Data.Number.ER.Real.Approx.Elementary ( - ERApproxElementary(..) + ERApproxElementary(..), + ERInnerOuterApproxElementary(..) ) where import Prelude hiding (exp, log, sin, cos) import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:)) import Data.Number.ER.BasicTypes import Data.Number.ER.Real.Arithmetic.Elementary @@ -38,6 +40,8 @@ min ix = Prelude.min max :: EffortIndex -> ra -> ra -> ra max ix = Prelude.max + sqrt :: EffortIndex -> ra -> ra + sqrt = erSqrt_IR exp :: EffortIndex -> ra -> ra exp = erExp_IR log :: EffortIndex -> ra -> ra @@ -54,6 +58,38 @@ tan ix r = (sin ix r) / (cos ix r) atan :: EffortIndex -> ra -> ra atan = erATan_IR + +{-| + A class defining various common real number operations + in a approximation-aware fashion, ie introducing effort indices. + + All operations here have default implementations based on + "Data.Number.ER.Real.Arithmetic.Elementary". +-} +class (RA.ERIntApprox ra, RA.ERInnerOuterApprox ra, Ord ra) => (ERInnerOuterApproxElementary ra) + where + absInner :: EffortIndex -> ra -> ra + absInner ix = Prelude.abs + minInner :: EffortIndex -> ra -> ra -> ra + minInner ix = Prelude.min + maxInner :: EffortIndex -> ra -> ra -> ra + maxInner ix = Prelude.max + sqrtInner :: EffortIndex -> ra -> ra + sqrtInner = erSqrt_IR_Inner + expInner :: EffortIndex -> ra -> ra + expInner = erExp_IR_Inner + logInner :: EffortIndex -> ra -> ra + logInner = erLog_IR_Inner + (**:) :: EffortIndex -> ra -> ra -> ra + (**:) ix b e = expInner ix $ e *: (logInner ix b) + sinInner :: EffortIndex -> ra -> ra + sinInner = erSine_IR_Inner + cosInner :: EffortIndex -> ra -> ra + cosInner = erCosine_IR_Inner + tanInner :: EffortIndex -> ra -> ra + tanInner ix r = (sinInner ix r) /: (cosInner ix r) + atanInner :: EffortIndex -> ra -> ra + atanInner = erATan_IR_Inner diff --git a/src/Data/Number/ER/Real/Approx/Interval.hs b/src/Data/Number/ER/Real/Approx/Interval.hs --- a/src/Data/Number/ER/Real/Approx/Interval.hs +++ b/src/Data/Number/ER/Real/Approx/Interval.hs @@ -1,5 +1,6 @@ {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} {-| Module : Data.Number.ER.Real.Approx.Interval Description : safe interval arithmetic @@ -16,17 +17,16 @@ module Data.Number.ER.Real.Approx.Interval ( ERInterval(..), - normaliseERInterval, - intervalTimesInner, - intervalPlusInner, - intervalDivideInner + normaliseERIntervalOuter, + normaliseERIntervalInner ) where import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:)) import qualified Data.Number.ER.Real.Approx.Elementary as RAEL import qualified Data.Number.ER.Real.Base as B -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.BasicTypes import Data.Number.ER.Misc @@ -44,26 +44,22 @@ Type for arbitrary precision interval arithmetic. -} data ERInterval base = - ERIntervalEmpty -- ^ usually represents computation error (top element in the interval domain) - | ERIntervalAny -- ^ represents no knowledge of result (bottom element in the interval domain) - | ERInterval +-- ERIntervalEmpty -- ^ usually represents computation error (top element in the interval domain) +-- | ERIntervalAny -- ^ represents no knowledge of result (bottom element in the interval domain) + ERInterval { - erintv_left :: base, - erintv_right :: base + erintv_left :: !base, + erintv_right :: !base } deriving (Typeable, Data) {- the following has been generated by BinaryDerive -} instance (Binary a) => Binary (ERInterval a) where - put ERIntervalEmpty = putWord8 0 - put ERIntervalAny = putWord8 1 - put (ERInterval a b) = putWord8 2 >> put a >> put b + put (ERInterval a b) = putWord8 0 >> put a >> put b get = do tag_ <- getWord8 case tag_ of - 0 -> return ERIntervalEmpty - 1 -> return ERIntervalAny - 2 -> get >>= \a -> get >>= \b -> return (ERInterval a b) + 0 -> get >>= \a -> get >>= \b -> return (ERInterval a b) _ -> fail "no parse" {- the above has been generated by BinaryDerive -} @@ -73,52 +69,66 @@ * no NaNs as endpoints - * @l <= r@ - - * no (-Infty, +Infty) + Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic. + This version interprets NaN's as bottomApprox. -} -normaliseERInterval :: +normaliseERIntervalOuter :: (B.ERRealBase b) => ERInterval b -> ERInterval b -normaliseERInterval (ERInterval minusInfty plusInfty) - | B.isPlusInfinity plusInfty && B.isPlusInfinity (- minusInfty) = - ERIntervalAny -normaliseERInterval (ERInterval nan1 nan2) +normaliseERIntervalOuter (ERInterval nan1 nan2) | B.isERNaN nan1 && B.isERNaN nan2 = - ERIntervalAny -normaliseERInterval (ERInterval nan r) + RA.bottomApprox +normaliseERIntervalOuter (ERInterval nan r) | B.isERNaN nan = ERInterval (- B.plusInfinity) r -normaliseERInterval (ERInterval l nan) +normaliseERIntervalOuter (ERInterval l nan) | B.isERNaN nan = ERInterval l (B.plusInfinity) -normaliseERInterval (ERInterval l r) - | l > r = ERIntervalEmpty -normaliseERInterval i = i +normaliseERIntervalOuter i = i {-| + convert to a normal form, ie: + + * no NaNs as endpoints + + Note that inverted intervals are fully supported using Warmus-Kaucher arithmetic. + This version interprets NaN's as topApprox. +-} +normaliseERIntervalInner :: + (B.ERRealBase b) => + ERInterval b -> ERInterval b +normaliseERIntervalInner (ERInterval nan1 nan2) + | B.isERNaN nan1 && B.isERNaN nan2 = + RA.topApprox +normaliseERIntervalInner (ERInterval nan r) + | B.isERNaN nan = + ERInterval (B.plusInfinity) r +normaliseERIntervalInner (ERInterval l nan) + | B.isERNaN nan = + ERInterval l (- B.plusInfinity) +normaliseERIntervalInner i = i + +{-| erintvPrecision returns an approximation of the number of bits required to represent the mantissa of a normalised size of the interval: - > - log_2 ((r - l) / (1 + abs(r) + abs(l))) - Notice that this is +Infty for singleton and empty intervals - and -Infty for the whole real line. + Notice that this is +Infty for singleton and anti-consistent intervals + and -Infty for unbounded intervals. -} erintvPrecision :: (B.ERRealBase b) => ERInterval b -> EI.ExtendedInteger -erintvPrecision (ERInterval l r) = - -1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l)) -erintvPrecision ERIntervalEmpty = EI.PlusInfinity -erintvPrecision ERIntervalAny = EI.MinusInfinity +erintvPrecision i@(ERInterval l r) + | not $ RA.isConsistent i = EI.PlusInfinity + | not $ RA.isBounded i = EI.MinusInfinity + | otherwise = + -1 - (B.getApproxBinaryLog $ (r - l)) -- /(1 + abs r + abs l)) erintvGranularity :: (B.ERRealBase b) => ERInterval b -> Int -erintvGranularity ERIntervalAny = 0 -erintvGranularity ERIntervalEmpty = 0 erintvGranularity (ERInterval l r) = min (B.getGranularity l) (B.getGranularity r) @@ -132,9 +142,6 @@ ERInterval b -> ERInterval b -> Bool erintvEqualApprox (ERInterval l1 r1) (ERInterval l2 r2) = l1 == l2 && r1 == r2 -erintvEqualApprox ERIntervalEmpty ERIntervalEmpty = True -erintvEqualApprox ERIntervalAny ERIntervalAny = True -erintvEqualApprox _ _ = False {-| a syntactic linear order @@ -142,12 +149,6 @@ erintvCompareApprox :: (B.ERRealBase b) => ERInterval b -> ERInterval b -> Ordering -erintvCompareApprox ERIntervalEmpty ERIntervalEmpty = EQ -erintvCompareApprox ERIntervalEmpty _ = LT -erintvCompareApprox _ ERIntervalEmpty = GT -erintvCompareApprox ERIntervalAny ERIntervalAny = EQ -erintvCompareApprox ERIntervalAny _ = LT -erintvCompareApprox _ ERIntervalAny = GT erintvCompareApprox (ERInterval l1 r1) (ERInterval l2 r2) = case compare l1 l2 of EQ -> compare r1 r2 @@ -165,10 +166,6 @@ ERInterval b -> ERInterval b -> Maybe Bool -erintvEqualReals ERIntervalEmpty _ = Nothing -erintvEqualReals _ ERIntervalEmpty = Nothing -erintvEqualReals ERIntervalAny _ = Nothing -erintvEqualReals _ ERIntervalAny = Nothing erintvEqualReals (ERInterval l1 r1) (ERInterval l2 r2) | l1 == r1 && l2 == r2 && l1 == l2 = Just True | r1 < l2 || l1 > r2 = Just False @@ -184,10 +181,6 @@ ERInterval b -> ERInterval b -> Maybe Ordering -erintvCompareReals ERIntervalEmpty _ = Nothing -erintvCompareReals _ ERIntervalEmpty = Nothing -erintvCompareReals ERIntervalAny _ = Nothing -erintvCompareReals _ ERIntervalAny = Nothing erintvCompareReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) | r1 < l2 = Just LT | l1 > r2 = Just GT @@ -204,10 +197,6 @@ ERInterval b -> ERInterval b -> Maybe Bool -erintvLeqReals ERIntervalEmpty _ = Nothing -erintvLeqReals _ ERIntervalEmpty = Nothing -erintvLeqReals ERIntervalAny _ = Nothing -erintvLeqReals _ ERIntervalAny = Nothing erintvLeqReals i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) | r1 <= l2 = Just True | l1 > r2 = Just False @@ -237,34 +226,35 @@ Granularity -> (ERInterval b) -> (ERInterval b) -erintvDefaultBisectPt gran ERIntervalAny = 0 -erintvDefaultBisectPt gran ERIntervalEmpty = ERIntervalEmpty -erintvDefaultBisectPt gran (ERInterval l r) = +erintvDefaultBisectPt gran (ERInterval l r) = ERInterval m m where - m - | B.isPlusInfinity r = - if l < 0 - then 0 - else 2 * (B.setMinGranularity gran l) + 1 - | B.isPlusInfinity (-l) = - if r > 0 - then 0 - else 2 * (B.setMinGranularity gran r) - 1 - | otherwise = - ((B.setMinGranularity gran l) + r)/2 + m = + case (B.isMinusInfinity l, B.isPlusInfinity r, B.isPlusInfinity l, B.isMinusInfinity r) of + (True, True, _, _) -> 0 -- [-oo,+oo] + (True, _,_,True) -> B.minusInfinity -- [-oo,-oo] + (_, True,True,_) -> B.plusInfinity -- [+oo,+oo] + (True, _,_,_) | r > 0 -> 0 + (True, _,_,_) -> 2 * (B.setMinGranularity gran r) - 1 + (_,True,_,_) | l < 0 -> 0 + (_,True,_,_) -> 2 * (B.setMinGranularity gran l) + 1 + (_,_,True,_) | r < 0 -> 0 + (_,_,True,_) -> 2 * (B.setMinGranularity gran r) + 1 + (_,_,_,True) | l > 0 -> 0 + (_,_,_,True) -> 2 * (B.setMinGranularity gran l) - 1 + _ -> ((B.setMinGranularity gran l) + r)/2 -- no infinities erintvBisect :: - (B.ERRealBase b, RealFrac b) => + (B.ERRealBase b) => Granularity -> (Maybe (ERInterval b)) -> (ERInterval b) -> (ERInterval b, ERInterval b) -erintvBisect gran maybePt i = - (l RA.\/ m, m RA.\/ r) +erintvBisect gran maybePt i@(ERInterval l r) = + (ERInterval l mR, ERInterval mL r) where - (l,r) = RA.bounds i + ERInterval mL mR = m m = case maybePt of Just m -> m @@ -294,20 +284,10 @@ eg compare is not defined for overlapping intervals.) -} max i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) = - normaliseERInterval $ ERInterval (max l1 l2) (max r1 r2) - max ERIntervalEmpty _ = ERIntervalEmpty - max _ ERIntervalEmpty = ERIntervalEmpty - max ERIntervalAny ERIntervalAny = ERIntervalAny - max ERIntervalAny (ERInterval l r) = ERInterval l B.plusInfinity - max (ERInterval l r) ERIntervalAny = ERInterval l B.plusInfinity + ERInterval (max l1 l2) (max r1 r2) {- min: -} min i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) = - normaliseERInterval $ ERInterval (min l1 l2) (min r1 r2) - min ERIntervalEmpty _ = ERIntervalEmpty - min _ ERIntervalEmpty = ERIntervalEmpty - min ERIntervalAny ERIntervalAny = ERIntervalAny - min ERIntervalAny (ERInterval l r) = ERInterval (- B.plusInfinity) r - min (ERInterval l r) ERIntervalAny = ERInterval (- B.plusInfinity) r + ERInterval (min l1 l2) (min r1 r2) instance (B.ERRealBase b) => Show (ERInterval b) where @@ -316,10 +296,12 @@ erintvShow numDigits showGran showComponents interval = showERI interval where - showERI ERIntervalEmpty = "[NONE]" - showERI ERIntervalAny = "[ANY]" - showERI (ERInterval l r) + showERI (ERInterval l r) + | (B.isMinusInfinity r) && (B.isPlusInfinity r) = + "[ANY]" | l == r = "<" ++ showBase l ++ ">" + | l > r = + "[!" ++ showBase l ++ "," ++ showBase r ++ "!]" | otherwise = "[" ++ showBase l ++ "," ++ showBase r ++ "]" showBase = B.showDiGrCmp numDigits showGran showComponents @@ -331,195 +313,201 @@ H.toHtml $ show l | otherwise = H.simpleTable [] [] [[H.toHtml l],[H.toHtml r]] - toHtml i = H.toHtml $ show i instance (B.ERRealBase b) => Num (ERInterval b) where fromInteger n = - normaliseERInterval $ ERInterval (fromInteger n) (fromInteger n) + ERInterval (B.fromIntegerDown n) (B.fromIntegerUp n) {- abs -} abs (ERInterval l r) - | l < 0 && r > 0 = ERInterval 0 (max (-l) r) + | l <= 0 && r >= 0 = ERInterval 0 (max (-l) r) + | l >= 0 && r <= 0 = ERInterval (max l (-r)) 0 | r <= 0 = ERInterval (-r) (-l) | otherwise = ERInterval l r - abs ERIntervalAny = ERInterval 0 B.plusInfinity - abs ERIntervalEmpty = ERIntervalEmpty {- signum -} - signum i@(ERInterval l r) - | l < 0 && r > 0 = ERInterval (-1) 1 -- need many-valuedness via sequences of intervals - | r < 0 = ERInterval (-1) (-1) - | l > 0 = ERInterval 1 1 - | l == 0 && r == 0 = i - | l == 0 = ERInterval 0 1 - | r == 0 = ERInterval (-1) 0 - signum ERIntervalAny = ERInterval (-1) 1 - signum ERIntervalEmpty = ERIntervalEmpty + signum i@(ERInterval l r) = + error "ER.Real.Approx.Interval: signum not implemented for ERInterval" +-- | l < 0 && r > 0 = ERInterval (-1) 1 -- need many-valuedness via sequences of intervals +-- | r < 0 = ERInterval (-1) (-1) +-- | l > 0 = ERInterval 1 1 +-- | l == 0 && r == 0 = i +-- | l == 0 = ERInterval 0 1 +-- | r == 0 = ERInterval (-1) 0 {- negate -} negate (ERInterval l r) = (ERInterval (-r) (-l)) - negate ERIntervalEmpty = ERIntervalEmpty - negate ERIntervalAny = ERIntervalAny {- addition -} - i1 + i2 = fst $ intervalPlusOuterInner i1 i2 + i1@(ERInterval l1 r1) + i2@(ERInterval l2 r2) = + normaliseERIntervalOuter $ + ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2) {- multiplication -} - i1 * i2 = fst $ intervalTimesOuterInner i1 i2 - -{-| - Add two real approximations, assuming the approximations are `inner' - as opposed to `outer': - - * `outer': the approximation contains all the number(s) of interest - * `inner': all numbers eligible for the approximation are numbers of interest --} -intervalPlusInner :: - (B.ERRealBase b) => - (ERInterval b) -> - (ERInterval b) -> - (ERInterval b) -intervalPlusInner i1 i2 = snd $ intervalPlusOuterInner i1 i2 - -{-| - Multiply two real approximations, assuming the approximations are `inner' - as opposed to `outer': - - * `outer': the approximation contains all the number(s) of interest - * `inner': all numbers eligible for the approximation are numbers of interest --} -intervalTimesInner :: - (B.ERRealBase b) => - (ERInterval b) -> - (ERInterval b) -> - (ERInterval b) -intervalTimesInner i1 i2 = snd $ intervalTimesOuterInner i1 i2 - -intervalPlusOuterInner (ERInterval l1 r1) (ERInterval l2 r2) = - (normaliseERInterval $ - ERInterval (l1 `plusDown` l2) (r1 `plusUp` r2), - ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2)) -intervalPlusOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny) -intervalPlusOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny) -intervalPlusOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty) -intervalPlusOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty) - -intervalTimesOuterInner (ERInterval l1 r1) (ERInterval l2 r2) - | haveNan = (ERIntervalAny, ERIntervalAny) - | otherwise = - (normaliseERInterval $ - ERInterval minProdOuter maxProdOuter, - ERInterval minProdInner maxProdInner) - where - haveNan = or $ map B.isERNaN (prodsUp ++ prodsDown) - minProdOuter = foldl1 min prodsDown - maxProdOuter = foldl1 max prodsUp - minProdInner = foldl1 min prodsUp - maxProdInner = foldl1 max prodsDown - prodsDown = [l1 `timesDown` l2, l1 `timesDown` r2, r1 `timesDown` l2, r1 `timesDown` r2] - prodsUp = [l1 `timesUp` l2, l1 `timesUp` r2, r1 `timesUp` l2, r1 `timesUp` r2] -intervalTimesOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny) -intervalTimesOuterInner l1 ERIntervalAny = (ERIntervalAny, ERIntervalAny) -intervalTimesOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty) -intervalTimesOuterInner l1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty) + i1@(ERInterval l1 r1) * i2@(ERInterval l2 r2) = + normaliseERIntervalOuter $ + intervalTimes timesDown timesUp i1 i2 instance (B.ERRealBase b) => Fractional (ERInterval b) where fromRational rat = (fromInteger $ numerator rat) / (fromInteger $ denominator rat) {- division -} - i1 / i2 = - fst $ intervalDivideOuterInner i1 i2 - -intervalDivideInner :: - (B.ERRealBase b) => - (ERInterval b) -> - (ERInterval b) -> - (ERInterval b) -intervalDivideInner i1 i2 = snd $ intervalDivideOuterInner i1 i2 + recip i@(ERInterval l r) + | not $ RA.isConsistent i = + RA.toggleConsistency $ + 1 /: (RA.toggleConsistency i) + | 0 < l || r < 0 = + normaliseERIntervalOuter $ + ERInterval (1 `divideDown` r) (1 `divideUp` l) + | otherwise = + RA.bottomApprox -intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 r2) - | l2 < 0 && r2 > 0 = (ERIntervalAny, ERIntervalAny) - | haveNan = --- unsafePrint "ERInterval: /: haveNan" $ - (ERIntervalAny, ERIntervalAny) - | l2 == 0 && r2 > 0 && 1/l2 < 0 = -- minus 0 - intervalDivideOuterInner (ERInterval l1 r1) (ERInterval (-l2) r2) -- correct it to +0 - | r2 == 0 && l2 < 0 && 1/r2 > 0 = -- plus 0 - intervalDivideOuterInner (ERInterval l1 r1) (ERInterval l2 (-r2)) -- correct it to -0 + +instance (B.ERRealBase b) => RA.ERInnerOuterApprox (ERInterval b) + where + {- addition -} + i1@(ERInterval l1 r1) +: i2@(ERInterval l2 r2) = + normaliseERIntervalInner $ + ERInterval (l1 `plusUp` l2) (r1 `plusDown` r2) + {- multiplication -} + i1@(ERInterval l1 r1) *: i2@(ERInterval l2 r2) = + normaliseERIntervalInner $ + intervalTimes timesUp timesDown i1 i2 + {- division -} + i1@(ERInterval l1 r1) /: i2@(ERInterval l2 r2) + | not $ RA.isConsistent i2 = + (*:) i1 $ + RA.toggleConsistency $ + 1 / (RA.toggleConsistency i2) + | 0 < l2 || r2 < 0 = + (*:) i1 $ + normaliseERIntervalInner $ + ERInterval (1 `divideDown` r2) (1 `divideUp` l2) | otherwise = --- unsafePrintReturn --- ( --- "intervalDivideOuterInner: " --- ++ "\n divsUp = " ++ show divsUp --- ++ "\n divsDown = " ++ show divsDown --- ++ "\n result = " --- ) - ( - normaliseERInterval $ - ERInterval minDivOuter maxDivOuter - , - ERInterval minDivInner maxDivInner - ) - where - haveNan = or $ map B.isERNaN (divsUp ++ divsDown) - minDivOuter = foldl1 min divsDown - maxDivOuter = foldl1 max divsUp - minDivInner = foldl1 min divsUp - maxDivInner = foldl1 max divsDown - divsDown = [l1 `divideDown` l2, l1 `divideDown` r2, r1 `divideDown` l2, r1 `divideDown` r2] - divsUp = [l1 `divideUp` l2, l1 `divideUp` r2, r1 `divideUp` l2, r1 `divideUp` r2] -intervalDivideOuterInner ERIntervalAny i2 = (ERIntervalAny, ERIntervalAny) -intervalDivideOuterInner i1 ERIntervalAny = (ERIntervalAny, ERIntervalAny) -intervalDivideOuterInner ERIntervalEmpty i2 = (ERIntervalEmpty, ERIntervalEmpty) -intervalDivideOuterInner i1 ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty) + RA.bottomApprox + {- setMinGranularityInner -} + setMinGranularityInner gr (ERInterval l r) = + normaliseERIntervalInner $ + (ERInterval (B.setMinGranularity gr l) (negate $ B.setMinGranularity gr (-r))) + {- setGranularityInner -} + setGranularityInner gr (ERInterval l r) = + normaliseERIntervalInner $ + (ERInterval (B.setGranularity gr l) (negate $ B.setGranularity gr (- r))) +intervalTimes timesL timesR i1@(ERInterval l1 r1) i2@(ERInterval l2 r2) = + ERInterval l r + where + (l,r) = + case (compare l1 0, compare r1 0, l1 <= r1, compare l2 0, compare r2 0, l2 <= r2) of + -- i1 negative, i2 positive + (LT, LT, _, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` l2) + -- i1 negative, i2 negative + (LT, LT, _, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` l2) + -- i1 negative, i2 consistent and containing zero + (LT, LT, _, _, _, True) -> (l1 `timesL` r2, l1 `timesR` l2) + -- i1 negative, i2 inconsistent and anti-containing zero + (LT, LT, _, _, _, False) -> (r1 `timesL` r2, r1 `timesR` l2) + + -- i1 positive, i2 positive + (GT, GT, _, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` r2) + -- i1 positive, i2 negative + (GT, GT, _, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` r2) + -- i1 positive, i2 consistent and containing zero + (GT, GT, _, _, _, True) -> (r1 `timesL` l2, r1 `timesR` r2) + -- i1 positive, i2 inconsistent and anti-containing zero + (GT, GT, _, _, _, False) -> (l1 `timesL` l2, l1 `timesR` r2) + + -- i1 consistent and containing zero, i2 positive + (_, _, True, GT, GT, _) -> (l1 `timesL` r2, r1 `timesR` r2) + -- i1 consistent and containing zero, i2 negative + (_, _, True, LT, LT, _) -> (r1 `timesL` l2, l1 `timesR` l2) + -- i1 consistent and containing zero, i2 consistent and containing zero + (_, _, True, _, _, True) -> + (l,r) + where + l | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity + | otherwise = min l1r2 r1l2 + where + l1r2 = l1 `timesL` r2 + r1l2 = r1 `timesL` l2 + r | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity + | otherwise = max l1l2 r1r2 + where + l1l2 = l1 `timesR` l2 + r1r2 = r1 `timesR` r2 + -- i1 consistent and containing zero, i2 inconsistent and anti-containing zero + (_, _, True, _, _, False) -> (0, 0) + + -- i1 inconsistent and anti-containing zero, i2 positive + (_, _, False, GT, GT, _) -> (l1 `timesL` l2, r1 `timesR` l2) + -- i1 inconsistent and anti-containing zero, i2 negative + (_, _, False, LT, LT, _) -> (r1 `timesL` r2, l1 `timesR` r2) + -- i1 inconsistent and anti-containing zero, i2 consistent and containing zero + (_, _, False, _, _, True) -> (0, 0) + -- i1 inconsistent and anti-containing zero, i2 the same + (_, _, False, _, _, False) -> + (l,r) + where + l | B.isERNaN l1l2 || B.isERNaN r1r2 = B.plusInfinity + | otherwise = max l1l2 r1r2 + where + l1l2 = l1 `timesL` l2 + r1r2 = r1 `timesL` r2 + r | B.isERNaN l1r2 || B.isERNaN r1l2 = B.minusInfinity + | otherwise = min l1r2 r1l2 + where + l1r2 = l1 `timesR` r2 + r1l2 = r1 `timesR` l2 + + -instance (B.ERRealBase b, RealFrac b) => RA.ERApprox (ERInterval b) where +instance (B.ERRealBase b) => RA.ERApprox (ERInterval b) where initialiseBaseArithmetic _ = B.initialiseBaseArithmetic (0 :: b) getPrecision i = erintvPrecision i getGranularity i = erintvGranularity i {- setMinGranularity -} - setMinGranularity gr (ERInterval l r) = - normaliseERInterval $ + setMinGranularityOuter gr (ERInterval l r) = + normaliseERIntervalOuter $ (ERInterval (- (B.setMinGranularity gr (-l))) (B.setMinGranularity gr r)) - setMinGranularity _ i = i {- setGranularity -} - setGranularity gr (ERInterval l r) = - normaliseERInterval $ + setGranularityOuter gr (ERInterval l r) = + normaliseERIntervalOuter $ (ERInterval (- (B.setGranularity gr (-l))) (B.setGranularity gr r)) - setGranularity _ i = i - {- bottomApprox -} - bottomApprox = ERIntervalAny - {- emptyApprox -} - emptyApprox = ERIntervalEmpty - {- isEmpty -} - isEmpty ERIntervalEmpty = True - isEmpty _ = False {- isBottom -} - isBottom ERIntervalAny = True isBottom (ERInterval l r) = - B.isPlusInfinity r && B.isPlusInfinity (-l) - isBottom _ = False + B.isMinusInfinity l && B.isPlusInfinity r + {- bottomApprox -} + bottomApprox = + ERInterval B.minusInfinity B.plusInfinity {- isExact -} - isExact ERIntervalEmpty = False - isExact ERIntervalAny = False isExact (ERInterval l r) = l == r + {- isConsistent -} + isConsistent (ERInterval l r) = l <= r + {- isAnticonsistent -} + isAnticonsistent (ERInterval l r) = l >= r + {- toggleConsistency -} + toggleConsistency (ERInterval l r) = (ERInterval r l) + {- isTop -} + isTop (ERInterval l r) = + B.isPlusInfinity l && B.isMinusInfinity r + {- topApprox -} + topApprox = + ERInterval B.plusInfinity B.minusInfinity {- isBounded -} - isBounded ERIntervalEmpty = True - isBounded ERIntervalAny = False isBounded (ERInterval l r) = - (- B.plusInfinity) < l && r < B.plusInfinity + (- B.plusInfinity) < l && l < B.plusInfinity + && + (- B.plusInfinity) < r && r < B.plusInfinity + {- plusInfinity -} + plusInfinity = ERInterval B.plusInfinity B.plusInfinity + {- refines -} + refines (ERInterval l1 r1) (ERInterval l2 r2) = + l2 <= l1 && r1 <= r2 + {- maybeRefines -} + maybeRefines i1 i2 = Just $ RA.refines i1 i2 + {- intersection -} - ERIntervalEmpty /\ i = ERIntervalEmpty - i /\ ERIntervalEmpty = ERIntervalEmpty - ERIntervalAny /\ i = i - i /\ ERIntervalAny = i (ERInterval l1 r1) /\ (ERInterval l2 r2) = - normaliseERInterval $ ERInterval (max l1 l2) (min r1 r2) {- intersectMeasureImprovement -} - intersectMeasureImprovement _ ERIntervalEmpty i = (ERIntervalEmpty, 1) - intersectMeasureImprovement _ i ERIntervalEmpty = (ERIntervalEmpty, 1) - intersectMeasureImprovement _ ERIntervalAny i = (i, 1) - intersectMeasureImprovement _ i ERIntervalAny = (i, 1) intersectMeasureImprovement ix i1 i2 = (isec, impr) where @@ -529,18 +517,10 @@ | otherwise = i1Width / isecWidth i1Width = i1H - i1L isecWidth = isecH - isecL - (isecL, isecH) = RA.bounds $ RA.setMinGranularity gran isec - (i1L, i1H) = RA.bounds $ RA.setMinGranularity gran i1 - gran = effIx2gran ix - {- refines -} - refines _ ERIntervalAny = True - refines ERIntervalEmpty _ = True - refines ERIntervalAny (ERInterval l r) - | B.isPlusInfinity r && B.isPlusInfinity (-l) = True - refines ERIntervalAny _ = False - refines _ ERIntervalEmpty = False - refines (ERInterval l1 r1) (ERInterval l2 r2) = - l2 <= l1 && r1 <= r2 + (isecL, isecH) = RA.bounds $ RA.setMinGranularityOuter gran isec + (i1L, i1H) = RA.bounds $ RA.setMinGranularityOuter gran i1 + gran = effIx2gran ix + {- semantic comparisons -} equalReals = erintvEqualReals compareReals = erintvCompareReals @@ -556,50 +536,27 @@ {- formatting -} showApprox = erintvShow -instance (B.ERRealBase b, RealFrac b) => RA.ERIntApprox (ERInterval b) +instance (B.ERRealBase b) => RA.ERIntApprox (ERInterval b) where - doubleBounds ERIntervalAny = (- infinity, infinity) - where - infinity = 1/0 - doubleBounds ERIntervalEmpty = - error "ERInterval: doubleBounds: empty interval" doubleBounds (ERInterval l r) = - (B.toDouble l, B.toDouble r) - floatBounds ERIntervalAny = (- infinity, infinity) - where - infinity = 1/0 - floatBounds ERIntervalEmpty = - error "ERInterval: floatBounds: empty interval" + (negate $ B.toDouble (-l), B.toDouble r) floatBounds (ERInterval l r) = - (B.toFloat l, B.toFloat r) - integerBounds ERIntervalAny = (- infinity, infinity) - where - infinity = EI.PlusInfinity - integerBounds ERIntervalEmpty = - error "ERInterval: integerBounds: empty interval" + (negate $ B.toFloat (-l), B.toFloat r) integerBounds (ERInterval l r) = - (- (mkEI (- l)), mkEI r) + (negate $ mkEI (- l), mkEI r) where mkEI f | B.isPlusInfinity f = EI.PlusInfinity - | B.isPlusInfinity (-f) = EI.MinusInfinity + | B.isMinusInfinity f = EI.MinusInfinity | otherwise = ceiling f - defaultBisectPt dom = erintvDefaultBisectPt (RA.getGranularity dom + 1) dom + defaultBisectPt dom = + erintvDefaultBisectPt (RA.getGranularity dom + 1) dom bisectDomain maybePt dom = erintvBisect (RA.getGranularity dom + 1) maybePt dom {- \/ -} - ERIntervalEmpty \/ i = i - i \/ ERIntervalEmpty = i - ERIntervalAny \/ _ = ERIntervalAny - _ \/ ERIntervalAny = ERIntervalAny (ERInterval l1 r1) \/ (ERInterval l2 r2) = - normaliseERInterval $ ERInterval (min l1 l2) (max r1 r2) {- RA.bounds -} - bounds ERIntervalAny = - (ERInterval (-B.plusInfinity) (-B.plusInfinity), - ERInterval B.plusInfinity B.plusInfinity) - bounds ERIntervalEmpty = (ERIntervalEmpty, ERIntervalEmpty) bounds (ERInterval l r) = (ERInterval l l, ERInterval r r) {- RA.fromBounds -} @@ -607,9 +564,11 @@ | l1 == r1 && l2 == r2 = ERInterval l1 l2 fromBounds i1i2 = error $ - "ER.Real.Approx.Interval: fromBounds: bounds not thin: " + "ER.Real.Approx.Interval: fromBounds: bounds not exact: " ++ show i1i2 -instance (B.ERRealBase b, RealFrac b) => RAEL.ERApproxElementary (ERInterval b) +instance (B.ERRealBase b) => RAEL.ERApproxElementary (ERInterval b) +instance (B.ERRealBase b) => RAEL.ERInnerOuterApproxElementary (ERInterval b) -- all operations here have appropriate default implementations - + + diff --git a/src/Data/Number/ER/Real/Approx/OI.hs b/src/Data/Number/ER/Real/Approx/OI.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/OI.hs @@ -0,0 +1,56 @@ +{-# OPTIONS_GHC -fno-warn-missing-methods #-} +{-| + Module : Data.Number.ER.Real.Approx.OI + Description : outer and inner approximations of approximations + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + This module offers a transformation of a safely rounded real approximation type into + a type that approximates these approximations from outside as well as *inside*. +-} +module Data.Number.ER.Real.Approx.OI where + +import qualified Data.Number.ER.Real.Approx as RA + +{-| + A pair of approximations that form an "interval" in the lattice of + approximations. + + Eg outer = [1,4] inner = [3,2] can be thought of as the set of all + generalised intervals where the left endpoint is between 1 and 3 + and the right endpoint is between 2 and 4 (eg [1,4], [3,4], + [3,2], [3,3]). +-} +data ERApproxOI ra = + ERApproxOI + { + eroiOuter :: ra, + eroiInner :: ra + } + deriving (Eq, Ord) + +instance (RA.ERApprox ra) => (Show (ERApproxOI ra)) + where + show (ERApproxOI oi ii) = + "{ outer = " ++ show oi ++ "; inner = " ++ show ii ++ "}" + +instance (RA.ERApprox ra) => RA.ERApproxApprox (ERApproxOI ra) + where + safeIncludes (ERApproxOI oi1 ii1) (ERApproxOI oi2 ii2) = + oi2 `RA.refines` ii1 + safeNotIncludes (ERApproxOI oi1 ii1) (ERApproxOI oi2 ii2) = + not $ ii2 `RA.refines` oi1 + +{- TODO when required: -} +instance (RA.ERApprox ra) => (Num (ERApproxOI ra)) +instance (RA.ERApprox ra) => (Fractional (ERApproxOI ra)) +instance (RA.ERApprox ra) => RA.ERApprox (ERApproxOI ra) + where + (ERApproxOI oi1 ii1) `leqReals` (ERApproxOI oi2 ii2) = + oi1 `RA.leqReals` oi2 + + diff --git a/src/Data/Number/ER/Real/Approx/Sequence.hs b/src/Data/Number/ER/Real/Approx/Sequence.hs --- a/src/Data/Number/ER/Real/Approx/Sequence.hs +++ b/src/Data/Number/ER/Real/Approx/Sequence.hs @@ -174,7 +174,7 @@ fromInteger n = ConvergRealSeq sq where sq ix = - RA.setMinGranularity (effIx2gran ix) $ fromInteger n + RA.setMinGranularityOuter (effIx2gran ix) $ fromInteger n abs = pointwiseConvergRealSeq1 $ abs signum = pointwiseConvergRealSeq1 $ signum negate = pointwiseConvergRealSeq1 $ negate @@ -188,7 +188,7 @@ fromRational q = ConvergRealSeq sq where sq ix = - (RA.setMinGranularity (effIx2gran ix) num) / denom + (RA.setMinGranularityOuter (effIx2gran ix) num) / denom num = fromInteger $ numerator q denom = fromInteger $ denominator q recip = pointwiseConvergRealSeq1 $ recip diff --git a/src/Data/Number/ER/Real/Approx/Tests/Generate.hs b/src/Data/Number/ER/Real/Approx/Tests/Generate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Tests/Generate.hs @@ -0,0 +1,177 @@ +{-| + Module : Data.Number.ER.Real.Approx.Tests.Generate + Description : (testing) generating real approximations + Copyright : (c) 2009 Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Generic instances of 'Arbitrary' class for generating (almost) random instances. +-} + +module Data.Number.ER.Real.Approx.Tests.Generate where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.BasicTypes + +import Test.QuickCheck + +import qualified Data.List as List + +newtype RAThin ira = RAThin ira deriving (Show) +newtype RAConsistent ira = RAConsistent ira deriving (Show) +newtype RADirected ira = RADirected ira deriving (Show) + + +instance (RA.ERIntApprox ira) => Arbitrary (RAThin ira) + where + arbitrary = + sized arbitrarySized + where + arbitrarySized n + | n < 28 = + do + gran <- choose (8,20) + (f1,f2,f3) <- arbitrary + isInfty <- choose (-inftyChance,inftyChance) + pow <- choose (-10,10) + return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow + | n <= 68 = + do + gran <- choose (30,100) + (f1,f2,f3) <- arbitrary + isInfty <- choose (-inftyChance,inftyChance) + pow <- choose (-100,100) + return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow + | otherwise = + do + gran <- choose (400,1000) + (f1,f2,f3) <- arbitrary + isInfty <- choose (-inftyChance,inftyChance) + pow <- choose (-10000,10000) + return $ RAThin $ constructThinRA isInfty gran (f1,f2,f3) pow + coarbitrary _ = + error "ER.Real.Approx: Tests: coarbitrary not implemented" + +inftyChance = 15 + +constructThinRA :: + (RA.ERIntApprox ra) => + Granularity -> + Int -> + (Double, Double, Double) -> + Int -> + ra +constructThinRA gran isInfty (f1,f2,f3) pow + | isInfty == inftyChance = + RA.setGranularityOuter gran $ RA.plusInfinity + | isInfty == - inftyChance = + RA.setGranularityOuter gran $ negate $ RA.plusInfinity + | isInfty == 0 = + RA.setGranularityOuter gran 0 + | otherwise = + fst $ RA.bounds $ -- ensure thinness + (\ (Just a) -> a) $ List.find RA.isBounded results + where + results = [result1, result2, result3, result4, result5, 0] + result1 = (b1/b2) ^^ pow + b3 + result2 = b1 * b2 + b3 + result3 = b1 ^^ pow - b2 + result4 = b1 - b2 + result5 = b1 + [b1,b2,b3] = map cvt [f1,f2,f3] + cvt f = RA.setGranularityOuter gran $ RA.double2ra f + +instance (RA.ERIntApprox ira) => Arbitrary (RAConsistent ira) + where + arbitrary = + sized arbitrarySized + where + arbitrarySized n + | n < 28 = + do + gran <- choose (8,20) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-10,10) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-10,10) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RAConsistent $ t1 RA.\/ t2 + | n <= 68 = + do + gran <- choose (30,100) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-100,100) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-100,100) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RAConsistent $ t1 RA.\/ t2 + | otherwise = + do + gran <- choose (400,1000) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-10000,10000) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-10000,10000) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RAConsistent $ t1 RA.\/ t2 + coarbitrary _ = + error "ER.Real.Approx: Tests: coarbitrary not implemented" + +instance (RA.ERIntApprox ira) => Arbitrary (RADirected ira) + where + arbitrary = + sized arbitrarySized + where + arbitrarySized n + | n < 28 = + do + gran <- choose (8,20) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-10,10) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-10,10) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RADirected $ RA.fromBounds (t1, t2) + | n <= 68 = + do + gran <- choose (30,100) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-100,100) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-100,100) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RADirected $ RA.fromBounds (t1, t2) + | otherwise = + do + gran <- choose (400,1000) + (f11,f12,f13) <- arbitrary + isInfty1 <- choose (-inftyChance,inftyChance) + pow1 <- choose (-10000,10000) + (f21,f22,f23) <- arbitrary + isInfty2 <- choose (-inftyChance,inftyChance) + pow2 <- choose (-10000,10000) + let t1 = constructThinRA isInfty1 gran (f11,f12,f13) pow1 + let t2 = constructThinRA isInfty2 gran (f21,f22,f23) pow2 + return $ RADirected $ RA.fromBounds (t1, t2) + coarbitrary _ = + error "ER.Real.Approx: Tests: coarbitrary not implemented" + + diff --git a/src/Data/Number/ER/Real/Approx/Tests/Properties.hs b/src/Data/Number/ER/Real/Approx/Tests/Properties.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Tests/Properties.hs @@ -0,0 +1,266 @@ +{-| + Module : Data.Number.ER.Real.Base.Tests.Properties + Description : (testing) properties to check for real approximations + Copyright : (c) 2009 Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Properties of real approximations we want to check in tests. +-} + +module Data.Number.ER.Real.Approx.Tests.Properties +where + +import Data.Number.ER.Real.Approx.Tests.Generate +import Data.Number.ER.Real.Approx.Tests.Reporting +import Data.Number.ER.BasicTypes.Tests.Generate + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Real.Approx ((+:),(-:),(*:),(/:)) +import qualified Data.Number.ER.Real.Approx.Elementary as RAEL + +import Data.Number.ER.BasicTypes + +import Data.Number.ER.Misc + +import Test.QuickCheck + +type RAPropTupleUnary ira = + ira -> + String -> + ((Ix20, RAThin ira) -> Bool, + (Ix20, RAConsistent ira) -> Bool, + (Ix20, RAThin ira) -> Bool, + (Ix20, RAConsistent ira) -> Bool, + (Ix20, RADirected ira) -> Bool) + +props_ra_AMinusA_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_AMinusA_eq_oi = + props_ra_eq_oi_unary 0 (\ix a -> a - a) (\ix a -> a -: a) + +props_ra_ADivA_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_ADivA_eq_oi = + props_ra_eq_oi_unary 1 (\ix a -> a / a) (\ix a -> a /: a) + +props_ra_AddCommut_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_AddCommut_eq_oi = + props_ra_eq_oi_unary 0 commutDiff commutDiffInner + where + commutDiff ix a = + (a + b) - (b + a) + where + b = 1 / (a + 1) + commutDiffInner ix a = + (a +: b) -: (b +: a) + where + b = 1 / (a + 1) + +props_ra_MultCommut_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_MultCommut_eq_oi = + props_ra_eq_oi_unary 0 commutDiff commutDiffInner + where + commutDiff ix a = + (a * b) - (b * a) + where + b = 1 / (a + 1) + commutDiffInner ix a = + (a *: b) -: (b *: a) + where + b = 1 / (a + 1) + +props_ra_AddAssoc_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_AddAssoc_eq_oi = + props_ra_eq_oi_unary 0 assocDiff assocDiffInner + where + assocDiff ix a = + ((a + b) + c) - (a + (b + c)) + where + b = 1 / (a + 1) + c = (a - 1) + assocDiffInner ix a = + ((a +: b) +: c) -: (a +: (b +: c)) + where + b = 1 / (a + 1) + c = (a - 1) + +props_ra_MultAssoc_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_MultAssoc_eq_oi = + props_ra_eq_oi_unary 0 assocDiff assocDiffInner + where + assocDiff ix a = + ((a * b) * c) - (a * (b * c)) + where + b = 1 / (a + 1) + c = (a - 1) + assocDiffInner ix a = + ((a *: b) *: c) -: (a *: (b *: c)) + where + b = 1 / (a + 1) + c = (a - 1) + +props_ra_Distr_eq_oi :: + (RA.ERIntApprox ira, RA.ERInnerOuterApprox ira) => RAPropTupleUnary ira +props_ra_Distr_eq_oi = + props_ra_eq_oi_unary 0 distrDiff distrDiffInner + where + distrDiff ix a = + (a * (b + c)) - (a * b + a * c) + where + b = 1 / (a + 1) + c = (a - 1) + distrDiffInner ix a = + (a *: (b +: c)) -: ((a *: b) +: (a *: c)) + where + b = 1 / (a + 1) + c = (a - 1) + +props_ra_SinCos_eq_oi :: + (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira +props_ra_SinCos_eq_oi = + props_ra_eq_oi_unary 1 sincos sincosInner + where + sincos ix a = + (RAEL.sin ix a)^2 + (RAEL.cos ix a)^2 + sincosInner ix a = + (sina *: sina) +: (cosa *: cosa) + where + sina = RAEL.sinInner ix a + cosa = RAEL.cosInner ix a + +props_ra_TanATan_eq_oi :: + (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira +props_ra_TanATan_eq_oi = + props_ra_eq_oi_unary 0 tanAtan tanAtanInner + where + tanAtan ixP a = +-- unsafePrint +-- ( +-- "tanAtan: " +-- ++ "\n ix = " ++ show ix +-- ++ "\n a = " ++ show a +-- ++ "\n atan ix a = " ++ show tana +-- ++ "\n tan ix (atan ix a) = " ++ show tanatana +-- ) $ + tanatana - a + where + tanatana = RAEL.tan ix tana + tana = RAEL.atan ix a + ix = min 10 ixP + tanAtanInner ixP a = + (RAEL.tanInner ix $ RAEL.atanInner ix a) -: a + where + ix = min 10 ixP + +props_ra_LogExp_eq_oi :: + (RAEL.ERInnerOuterApproxElementary ira, RAEL.ERApproxElementary ira) => RAPropTupleUnary ira +props_ra_LogExp_eq_oi = + props_ra_eq_oi_unary 0 logExp logExpInner + where + logExp ixP a = +-- unsafePrint +-- ( +-- "logExp: " +-- ++ "\n ix = " ++ show ix +-- ++ "\n a = " ++ show a +-- ++ "\n exp ix a = " ++ show expa +-- ++ "\n log ix (exp ix a) = " ++ show logexpa +-- ) $ + logexpa - a + where + logexpa = RAEL.log ix expa + expa = RAEL.exp ix a + ix = min 10 ixP + logExpInner ixP a = + logexpa -: a + where + logexpa = RAEL.logInner ix expa + expa = RAEL.expInner ix a + ix = min 10 ixP + + +{------------------ auxiliary functions ------------------------} + +props_ra_eq_oi_unary constRes opOuter opInner sampleRA reportFileName = + (prop_Eq_Thin, prop_Eq_Consistent, + prop_OI_Thin, prop_OI_Consistent, prop_OI_Directed) + where + prop_Eq_Thin (Ix20 ix, RAThin a) = + raConsistentWithPrecise sampleRA (reportFileName ++ "_Eq_Thin") (ix,aId) 0 constRes resOuter + where + resOuter = opOuter ix a + aId = RA.showApprox 10 True True a + prop_Eq_Consistent (Ix20 ix, RAConsistent a) = + raConsistentWithPrecise sampleRA (reportFileName ++ "_Eq_Consistent") (ix,aId) 0 constRes resOuter + where + resOuter = opOuter ix a + aId = RA.showApprox 10 True True a + prop_OI_Thin (Ix20 ix, RAThin a) = + raIncludedIn sampleRA (reportFileName ++ "_OI_Thin") (ix, aId) 0 resInner resOuter + where + resOuter = opOuter ix a + resInner = opInner ix a + aId = RA.showApprox 10 True True a + prop_OI_Consistent (Ix20 ix, RAConsistent a) = + raIncludedIn sampleRA (reportFileName ++ "_OI_Consistent") (ix,aId) 0 resInner resOuter + where + resOuter = opOuter ix a + resInner = opInner ix a + aId = RA.showApprox 10 True True a + prop_OI_Directed (Ix20 ix, RADirected a) = + raIncludedIn sampleRA (reportFileName ++ "_OI_Directed") (ix, aId) 0 resInner resOuter + where + resOuter = opOuter ix a + resInner = opInner ix a + aId = RA.showApprox 10 True True a + +raConsistentWithPrecise sampleRA reportFileName caseId subId preciseVal approxVal + | result = + unsafeERTestReport reportFileName + (caseId, subId, preciseVal, approxVal) $ + result + | otherwise = + unsafePrint + ( + "raAntiIncludes failed" + ++ "\n caseId = " ++ show caseId + ++ "\n subId = " ++ show subId + ++ "\n preciseVal = " ++ show preciseVal + ++ "\n approxVal = " ++ show approxVal + ) $ + result + where + result = + (approxVal `RA.refines` preciseVal) + || + (preciseVal `RA.refines` approxVal) + _ = [sampleRA, approxVal] + +raIncludedIn sampleRA reportFileName caseId subId innerVal outerVal + | result = + unsafeERTestReport reportFileName + (caseId, subId, innerVal, outerVal) $ + result + | otherwise = + unsafePrint + ( + "raIncludes failed" + ++ "\n caseId = " ++ show caseId + ++ "\n subId = " ++ show subId + ++ "\n innerVal = " ++ show innerVal + ++ "\n outerVal = " ++ show outerVal + ) $ + result + where + result = innerVal `RA.refines` outerVal + _ = [sampleRA, innerVal] + + diff --git a/src/Data/Number/ER/Real/Approx/Tests/Reporting.hs b/src/Data/Number/ER/Real/Approx/Tests/Reporting.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Tests/Reporting.hs @@ -0,0 +1,167 @@ + +module Data.Number.ER.Real.Approx.Tests.Reporting + +where + +import qualified Data.Number.ER.Real.Approx as RA + +import Data.Number.ER.Misc + +import qualified Data.List as List +import Text.Regex +import System.IO + + +unsafeERTestReport :: + (Show tId, Show sId, RA.ERIntApprox ira) => + String -> + (tId, sId, ira, ira) -> + a -> a +unsafeERTestReport reportFileName (testId, subId, almostPreciseVal, approxVal) = + unsafeReport reportFileName $ + stdRepLine (testId, subId) (overestimation, detail) + where + overestimation = fst $ getOverestimation almostPreciseVal approxVal + detail = (almostPreciseVal, approxVal) + +stdRepLine (testId, subId) (overestimation, detail) = + "case=" ++ show testId + ++ ";pt=" ++ show subId + ++ ";ovest=" ++ show overestimation + ++ ";detail=" ++ show detail + +getOverestimation :: + (RA.ERIntApprox ira) => + ira -> ira -> (Double, (ira, ira)) +getOverestimation model res = + ((abs $ wMod - wRes) / (1 + (max 0 (wMod))), (model, res)) + where + wMod = hMod - lMod + wRes = hRes - lRes + (lMod, hMod) = RA.doubleBounds model + (lRes, hRes) = RA.doubleBounds res + +produceSummary :: String -> IO () +produceSummary filepath = + do + casesInfo <- parseReport filepath + writeFile summaryFilepath $ formatSummary casesInfo + return () + where + summaryFilepath = filepath ++ "-summary" + formatSummary casesInfo = + "all " ++ show casesCount ++ " cases:" + ++ "\n approx. average time per case: " ++ show timeInSeconds ++ " seconds" + ++ "\n approx. average per case average overestimation: " ++ show avgOverestimation + ++ "\n approx. average per case maximum overestimation: " ++ show maxOverestimation + ++ "\n\n removing the worst 5% of the cases (for each measure separately):" + ++ "\n approx. average time per case: " ++ show timeInSeconds95 ++ " seconds" + ++ "\n approx. average per case average overestimation: " ++ show avgOverestimation95 + ++ "\n approx. average per case maximum overestimation: " ++ show maxOverestimation95 + ++ "\n\n considering only the worst 50% but not the worst 5% of the cases (for each measure separately):" + ++ "\n approx. average time per case: " ++ show timeInSeconds45 ++ " seconds" + ++ "\n approx. average per case average overestimation: " ++ show avgOverestimation45 + ++ "\n approx. average per case maximum overestimation: " ++ show maxOverestimation45 + ++ "\n\n considering only the best 50% of the cases (for each measure separately):" + ++ "\n approx. average time per case: " ++ show timeInSeconds50 ++ " seconds" + ++ "\n approx. average per case average overestimation: " ++ show avgOverestimation50 + ++ "\n approx. average per case maximum overestimation: " ++ show maxOverestimation50 + ++ "\n\n" ++ (unlines $ map formatSummaryCase casesInfo) + where + (allTimes, (allAvgOvers, allMaxOvers)) = + mapSnd unzip $ unzip $ snd $ unzip casesInfo + timeInSeconds = (sum allTimes) / casesCountF + avgOverestimation = (sum allAvgOvers) / casesCountF + maxOverestimation = (sum allMaxOvers) / casesCountF + casesCount = length casesInfo + casesCountF :: Double + casesCountF = fromInteger $ toInteger casesCount + + timeInSeconds95 = (sum allTimes95) / casesCount95F + avgOverestimation95 = (sum allAvgOvers95) / casesCount95F + maxOverestimation95 = (sum allMaxOvers95) / casesCount95F + allTimes95 = drop fivePerCent $ reverse $ List.sort allTimes + allAvgOvers95 = drop fivePerCent $ reverse $ List.sort allAvgOvers + allMaxOvers95 = drop fivePerCent $ reverse $ List.sort allMaxOvers + casesCount95F = fromInteger $ toInteger $ casesCount - fivePerCent + fivePerCent = max 1 $ (5 * casesCount) `div` 100 + + timeInSeconds50 = (sum allTimes50) / casesCount50F + avgOverestimation50 = (sum allAvgOvers50) / casesCount50F + maxOverestimation50 = (sum allMaxOvers50) / casesCount50F + allTimes50 = drop fiftyPerCent $ reverse $ List.sort allTimes + allAvgOvers50 = drop fiftyPerCent $ reverse $ List.sort allAvgOvers + allMaxOvers50 = drop fiftyPerCent $ reverse $ List.sort allMaxOvers + casesCount50F = fromInteger $ toInteger $ casesCount - fiftyPerCent + fiftyPerCent = casesCount `div` 2 + + timeInSeconds45 = (sum allTimes45) / casesCount45F + avgOverestimation45 = (sum allAvgOvers45) / casesCount45F + maxOverestimation45 = (sum allMaxOvers45) / casesCount45F + allTimes45 = drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allTimes + allAvgOvers45 = drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allAvgOvers + allMaxOvers45 = drop fivePerCent $ reverse $ drop fiftyPerCent $ List.sort allMaxOvers + casesCount45F = fromInteger $ toInteger $ casesCount - fiftyPerCent - fivePerCent + formatSummaryCase (caseId, (timeInSeconds, (avgOverestimation, maxOverestimation))) = + "case " ++ caseId ++ ":" + ++ "\n approximate time = " ++ show timeInSeconds ++ " seconds" + ++ "\n average sampled overestimation = " ++ show avgOverestimation + ++ "\n maximal sampled overestimation = " ++ show maxOverestimation + parseReport :: String -> IO [(String, (Double, (Double, Double)))] + parseReport filepath = + withFile filepath ReadMode readFirstAndOtherLines + where + readFirstAndOtherLines h = + do + startLine <- hGetLine h + firstLine <- hGetLine h + readCases (firstLine, (getTime firstLine) - (getTime startLine)) h + readCases (currentLine, caseCompTime) h = + do + (caseOverestimations, maybeNextLineAndTime) <- readCase [] 0 currentLine + let caseInfo = (caseId, (caseCompTime, avgAndMax caseOverestimations)) + case maybeNextLineAndTime of + Nothing -> return [caseInfo] + Just (nextLine, nextCaseTime) -> + do + otherCases <- readCases (nextLine, nextCaseTime) h + return $ caseInfo : otherCases + where + avgAndMax ns = + (sum ns / (fromInteger $ toInteger $ length ns), foldl1 max ns) + caseId = getCaseId currentLine + readCase overestimationsSoFar currentTimeStep currentLine + | currentCaseId /= caseId = + return (overestimationsSoFar, Just (currentLine, currentTimeStep)) + | otherwise = + do + finished <- hIsEOF h + case finished of + True -> return (currentOverestimations, Nothing) + False -> + do + nextLine <- hGetLine h + let nextTimeStep = (getTime nextLine) - (getTime currentLine) + readCase currentOverestimations nextTimeStep nextLine + where + currentCaseId = getCaseId currentLine + currentOverestimations = + currentOverestimation : overestimationsSoFar + currentOverestimation = getOverestimation currentLine + getTime :: String -> Double + getTime line = + case reads line of + [(time,'s':_)] -> time + getCaseId :: String -> String + getCaseId line = + case matchRegex idRegex line of + Just [caseId] -> caseId + where + idRegex = mkRegex "case=([^;]*);" + getOverestimation :: String -> Double + getOverestimation line = + case matchRegex ovestRegex line of + Just [ovestS] -> read ovestS + where + ovestRegex = mkRegex "ovest=([^;]*);" + diff --git a/src/Data/Number/ER/Real/Approx/Tests/Run.hs b/src/Data/Number/ER/Real/Approx/Tests/Run.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Approx/Tests/Run.hs @@ -0,0 +1,100 @@ +{-| + Module : Data.Number.ER.Real.Approx.Tests.Run + Description : (testing) running all function enclosure base tests in a batch + Copyright : (c) 2007-2008 Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Support for running all function enclosure base tests in a batch. +-} +module Data.Number.ER.Real.Approx.Tests.Run +where + +import Data.Number.ER.Real.Approx.Tests.Generate +import Data.Number.ER.Real.Approx.Tests.Properties +import Data.Number.ER.Real.Approx.Tests.Reporting + +import qualified Data.Number.ER.Real.Approx as RA +import qualified Data.Number.ER.Real.Approx.Elementary as RAEL + +import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) + +import Data.Number.ER.Misc.Tests +import Data.Number.ER.Misc + +import Test.QuickCheck +import Test.QuickCheck.Batch + +import System.Directory +import qualified System.FilePath as FP +import Data.Time.Clock +import Data.Time.Calendar + +runRATests :: + (RAEL.ERApproxElementary ra, + RAEL.ERInnerOuterApproxElementary ra, + Ord ra) => + String -> ra -> IO () -> IO () +runRATests title sampleRA initialise = + do + (UTCTime (ModifiedJulianDay days) secs) <- getCurrentTime + let folder = "tests-" ++ title ++ "-" ++ (show days) ++ "-" ++ (show $ floor $ toRational secs) + createDirectory folder + erRunTests (title ++ " real approx tests") raTestOptions initialise (raTests sampleRA folder) + +raTestOptions = + TestOptions + { +-- no_of_tests = 10 +-- no_of_tests = 50 +-- no_of_tests = 200 + no_of_tests = 500 + , + length_of_tests = 240 * 3600 -- ie 4h time limit + , + debug_tests = False + } + +raTests sampleRA folder = + (propTuple "a-a=0" "AMinusAIsZero" props_ra_AMinusA_eq_oi) + ++ + (propTuple "a/a=1" "ADivAIsOne" props_ra_ADivA_eq_oi) + ++ + (propTuple "a+b=b+a" "AddCommut" props_ra_AddCommut_eq_oi) + ++ + (propTuple "a*b=b*a" "MultCommut" props_ra_MultCommut_eq_oi) + ++ + (propTuple "(a+b)+c=a+(b+c)" "AddAssoc" props_ra_AddAssoc_eq_oi) + ++ + (propTuple "(a*b)*c=a*(b*c)" "MultAssoc" props_ra_MultAssoc_eq_oi) + ++ + (propTuple "a*(b+c)=a*b+a*c" "Distr" props_ra_Distr_eq_oi) + ++ + (propTuple "log(exp(a))=a" "LogExp" props_ra_LogExp_eq_oi) + ++ + (propTuple "(sin a)^2+(cos a)^2=1" "SinCos" props_ra_SinCos_eq_oi) + ++ + (propTuple "tan(atan(a))=a" "TanATan" props_ra_TanATan_eq_oi) + where + propTuple testName testFileName propGen = + [ + (testName ++ ", equality, thin intervals", runR prop_eq_Thin $ filepath ++ "_Eq_Thin"), + (testName ++ ", equality, consistent intervals", runR prop_eq_Consistent $ filepath ++ "_Eq_Consistent"), + (testName ++ ", inner in outer, thin intervals", runR prop_oi_Thin $ filepath ++ "_OI_Thin"), + (testName ++ ", inner in outer, consistent intervals", runR prop_oi_Consistent $ filepath ++ "_OI_Consistent"), + (testName ++ ", inner in outer, directed intervals", runR prop_oi_Directed $ filepath ++ "_OI_Directed") + ] + where + (prop_eq_Thin, prop_eq_Consistent, prop_oi_Thin, prop_oi_Consistent, prop_oi_Directed) = + propGen sampleRA filepath + filepath = FP.combine folder testFileName + runR test filepath opts = + unsafeReport filepath "started" $ + do + result <- run test opts + produceSummary filepath + return result + diff --git a/src/Data/Number/ER/Real/Arithmetic/Elementary.hs b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs --- a/src/Data/Number/ER/Real/Arithmetic/Elementary.hs +++ b/src/Data/Number/ER/Real/Arithmetic/Elementary.hs @@ -47,6 +47,7 @@ import qualified Data.Number.ER.Real.Approx as RA import Data.Number.ER.BasicTypes +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Real.Arithmetic.Taylor -- import Data.Number.ER.Real.Arithmetic.Newton @@ -74,13 +75,10 @@ (RA.ERIntApprox ira, Ord ira) => EffortIndex -> ira -> ira -erSqr_R ix a - | RA.isEmpty a = - RA.emptyApprox - | otherwise = - max 0 $ a' * a' +erSqr_R ix a = + max 0 $ a' * a' where - a' = RA.setMinGranularity gran a + a' = RA.setMinGranularityOuter gran a gran = effIx2gran ix {- @@ -109,8 +107,6 @@ Integer -> ira -> ira erPow_R ix p a - | RA.isEmpty a = - RA.emptyApprox | p < 0 = 1 / erPow_R ix (-p) a | p == 0 = @@ -120,7 +116,7 @@ | otherwise = a' * (erPow_R ix (div (p - 1) 2) (erSqr_R ix a')) where - a' = RA.setMinGranularity gran a + a' = RA.setMinGranularityOuter gran a gran = effIx2gran ix {- @@ -145,7 +141,7 @@ EffortIndex -> ira -> ira erSqrt_IR_Inner = RA.maxExtensionInnerR2R - sqrtExtremaAndDirections + sqrtExtrema (\ ix x -> erSqrt_R ix x) sqrtExtrema ix x = fst $ sqrtExtremaAndDirections ix x @@ -163,11 +159,11 @@ EffortIndex -> ira -> ira erSqrtContFr_R ix a | aR == 0 = 0 - | aL == 1/0 = 1/0 - | aR `RA.ltSingletons` 0 = RA.emptyApprox + | aL == RA.plusInfinity = RA.plusInfinity + | aR `RA.ltSingletons` 0 = RA.topApprox | otherwise = contFrIter (ix + 3) $ - RA.setMinGranularity gran $ 0 RA.\/ aR -- assuming aR >= 0 + RA.setMinGranularityOuter gran $ 0 RA.\/ aR -- assuming aR >= 0 where gran = effIx2gran ix (aL, aR) = RA.bounds a @@ -185,10 +181,9 @@ (RA.ERIntApprox ira) => EffortIndex -> ira -> ira erSqrtNewton_R ix a - | RA.isEmpty a = RA.emptyApprox | aR == 0 = 0 - | aL == 1/0 = 1/0 - | aR `RA.ltSingletons` 0 = RA.emptyApprox + | aL == RA.plusInfinity = RA.plusInfinity + | aR `RA.ltSingletons` 0 = RA.topApprox | otherwise = x_i RA.\/ (a/x_i) where @@ -198,7 +193,7 @@ x_i = newtonIter ((ix `div` 10) + 5) $ - RA.setMinGranularity gran aR -- assuming aR >= 0 + RA.setMinGranularityOuter gran aR -- assuming aR >= 0 newtonIter i x_i | i == 0 = x_i | otherwise = @@ -230,7 +225,7 @@ EffortIndex -> Integer -> ira -> ira erRoot_IR_Inner ix p = RA.maxExtensionInnerR2R - (rootExtremaAndDirections p) + (rootExtrema p) (\ ix x -> erRoot_R ix p x) $ ix rootExtrema p ix x = fst $ rootExtremaAndDirections p ix x @@ -247,10 +242,9 @@ (RA.ERIntApprox ira, Ord ira) => EffortIndex -> Integer -> ira -> ira erRootNewton_R ix p a - | RA.isEmpty a = RA.emptyApprox | aR == 0 = 0 - | aL == 1/0 = 1/0 - | aR < 0 && even p = RA.emptyApprox + | aL == RA.plusInfinity = RA.plusInfinity + | aR < 0 && even p = RA.topApprox | aR < 0 = - erRootNewton_R ix p (-a) | p > 0 = x_i RA.\/ (a/x_i_pow_p_minus_1) @@ -265,7 +259,7 @@ (x_i, x_i_pow_p_minus_1) = newtonIter (ix + 5) $ - RA.setMinGranularity gran $ max 0 aR + RA.setMinGranularityOuter gran $ max 0 aR newtonIter i x_0 | i == 0 = @@ -294,10 +288,19 @@ erExp_R ix x | RA.isBounded x = +-- unsafePrintReturn +-- ( +-- "erExp_R: " +-- ++ "\n x = " ++ show x +-- ++ "\n xNear0 = " ++ show xNear0 +-- ++ "\n n = " ++ show n +-- ++ "\n erExp_Tay_Opt_R ix xNear0 = " ++ (show $ erExp_Tay_Opt_R ix xNear0) +-- ++ "\n result = " +-- ) $ erPow_IR ix n $ erExp_Tay_Opt_R ix xNear0 - | x `RA.refines` (-1/0) = 0 - | (-1/0) `RA.refines` x = + | x `RA.refines` (-RA.plusInfinity) = 0 + | (-RA.plusInfinity) `RA.refines` x = 0 RA.\/ (erExp_R ix (snd $ RA.bounds x)) | otherwise = RA.bottomApprox where @@ -316,7 +319,7 @@ erExp_IR = RA.maxExtensionR2R - (\ ix x -> []) + noExtrema erExp_R erExp_IR_Inner :: @@ -324,9 +327,11 @@ EffortIndex -> ira -> ira erExp_IR_Inner = RA.maxExtensionInnerR2R - (\ ix x -> ([], (Just True, Just True))) + noExtrema erExp_R +noExtrema ix x = [] + {- Log using Newton -} erLog_R :: @@ -352,7 +357,7 @@ erLog_IR_Inner = RA.maxExtensionInnerR2R - logExtremaAndDirections + logExtrema (\ ix x -> logDivSeries_R ix x) logExtrema ix x = fst $ logExtremaAndDirections ix x @@ -361,24 +366,40 @@ case RA.compareReals 0 x of Just LT -> ([], (Just True, Just True)) Just GT -> ([], (Nothing, Nothing)) - _ -> ([-1/0], (Nothing, Just True)) + _ -> ([-RA.plusInfinity], (Nothing, Just True)) {-| log using a fast converging series, designed to be used with singletons -} logDivSeries_R :: (RA.ERIntApprox ira, Ord ira) => EffortIndex -> ira -> ira logDivSeries_R ix x - | RA.isEmpty posx = RA.emptyApprox - | posx `RA.refines` 0 = -1/0 - | posx `RA.refines` (1/0) = 1/0 + | posx `RA.refines` 0 = -RA.plusInfinity + | 0 `RA.refines` posx = RA.bottomApprox + | posx `RA.refines` (RA.plusInfinity) = RA.plusInfinity | otherwise = case RA.compareReals posx 1 of Just LT -> - - (logDivSeries_R ix (recip posx)) +-- unsafePrint +-- ( +-- "logDivSeries_R: recursion via recip" +-- ) $ + negate $ + (logDivSeries_R ix posxRecipL) + RA.\/ + (logDivSeries_R ix posxRecipR) _ -> - nearLogx + 2 * t * (series ix (RA.setMinGranularity gran 1)) +-- unsafePrint +-- ( +-- "logDivSeries_R: using series" +-- ++ "\n posx = " ++ show posx +-- ++ "\n nearLogx = " ++ show nearLogx +-- ++ "\n remNearLogx = " ++ show remNearLogx +-- ++ "\n t = " ++ show t +-- ) $ + nearLogx + 2 * t * (series ix (RA.setMinGranularityOuter gran 1)) where gran = effIx2gran ix - posx = (RA.setMinGranularity gran x) RA./\ (0 RA.\/ (1/0)) + posx = (RA.setMinGranularityOuter gran x) RA./\ (0 RA.\/ (RA.plusInfinity)) + (posxRecipL, posxRecipR) = RA.bounds $ recip posx nearLogx = 0.69314718055994530941 * (fromInteger $ intLogUp 2 $ xCeiling) remNearLogx = @@ -411,14 +432,14 @@ -- _ -> -- erNewton_FullArgs -- ( \ i y -> (erExp_RA i y) - posx, erExp_RA) --- (RA.setMinGranularity gran nearLogx) --- (RA.setMinGranularity gran 1) +-- (RA.setMinGranularityOuter gran nearLogx) +-- (RA.setMinGranularityOuter gran 1) -- (fromInteger $ toInteger i) -- i -- where -- gran = effIx2gran i -- posx = --- RA.setMinGranularity gran x /\ (ira2ra $ 0 RA.\/ (1/0)) +-- RA.setMinGranularityOuter gran x /\ (ira2ra $ 0 RA.\/ (RA.plusInfinity)) -- nearLogx = -- 0.69314718055994530941 * (fromInteger $ intLog 2 $ xCeiling) -- xCeiling @@ -490,7 +511,7 @@ EffortIndex -> ira -> ira erSine_Tay_R ix x - | (1/0) `RA.refines` x || (-1/0) `RA.refines` x = + | (RA.plusInfinity) `RA.refines` x || (-RA.plusInfinity) `RA.refines` x = (-1) RA.\/ 1 | otherwise = erTaylor_R ix sine_coefSeq sine_error 0 x @@ -528,14 +549,14 @@ EffortIndex -> ira -> ira erSine_IR_Inner = - RA.maxExtensionInnerR2R sineExtremesAndDirections erSine_R + RA.maxExtensionInnerR2R sineExtremes erSine_R erCosine_IR_Inner :: (RA.ERIntApprox ira) => EffortIndex -> ira -> ira erCosine_IR_Inner = - RA.maxExtensionInnerR2R cosineExtremesAndDirections erCosine_R + RA.maxExtensionInnerR2R cosineExtremes erCosine_R sineExtremes ix x = fst $ sineExtremesAndDirections ix x cosineExtremes ix x = fst $ cosineExtremesAndDirections ix x @@ -592,18 +613,14 @@ EffortIndex -> ira -> ira erATan_IR = - RA.maxExtensionR2R atanExtremes erATan_R - -atanExtremes ix x = [] + RA.maxExtensionR2R noExtrema erATan_R erATan_IR_Inner :: (RA.ERIntApprox ira) => EffortIndex -> ira -> ira erATan_IR_Inner = - RA.maxExtensionInnerR2R atanExtremesAndDirections erATan_R - -atanExtremesAndDirections ix x = ([], (Just True, Just True)) + RA.maxExtensionInnerR2R noExtrema erATan_R {- atan using Euler's series: (x / (1 + x^2)) * (1 + t*2*1/(2*1 + 1)*(1 + t*2*2/(2*2 + 1)*(1 + ... (1 + t*2*n/(2*n+1)*(1 + ...))))) @@ -619,11 +636,13 @@ EffortIndex -> ira -> ira atanEuler_R ix x - | RA.isEmpty x = RA.emptyApprox + | x `RA.refines` RA.plusInfinity = RA.plusInfinity + | x `RA.refines` (- RA.plusInfinity) = - RA.plusInfinity + | not $ RA.isBounded x = RA.bottomApprox | x `RA.refines` ((-1.5) RA.\/ 1.5) = - (x / xSquarePlus1) * (series ix (RA.setMinGranularity gran 2)) + (x / xSquarePlus1) * (series ix (RA.setMinGranularityOuter gran 2)) | otherwise = -- too far from 0, needs atan(x) = 2*atan(x/(1+sqrt(1+x^2))) - 2 * (atanEuler_R ix $ x / (1 + (erSqrt_R ix $ 1 + x * x))) + 2 * (atanEuler_R ix $ x / (1 + sqrtXQuarePlus1)) where gran = effIx2gran ix series termsCount coeffBase @@ -636,6 +655,16 @@ xSquare = abs $ x * x xSquarePlus1 = xSquare + 1 xSquareOverXSquarePlus1 = xSquare / xSquarePlus1 + sqrtXQuarePlus1 = + iterateIx 10 EI.MinusInfinity + where + iterateIx ix prevPrec + | prevPrec == currentPrec = result + | otherwise = + iterateIx (ix * 2) currentPrec + where + result = erSqrt_R ix xSquarePlus1 + currentPrec = RA.getPrecision result --{- atan using Newton -} -- @@ -646,8 +675,8 @@ --atanNewton_RA i x = -- erNewton_FullArgs -- ( \ i y -> (erTan_RA i y) - x, erTanDeriv_RA) --- (RA.setMinGranularity (effIx2gran i) (x)) --- (RA.setMinGranularity (effIx2gran i) 1) +-- (RA.setMinGranularityOuter (effIx2gran i) (x)) +-- (RA.setMinGranularityOuter (effIx2gran i) 1) -- (fromInteger $ toInteger i) -- i @@ -718,9 +747,9 @@ -} where gran = max 0 (effIx2gran ix) + 10 - r1over64 = (RA.setMinGranularity gran 1) / 64 - r1over1024 = (RA.setMinGranularity gran 1) / 1024 - z = RA.setMinGranularity gran 0 + r1over64 = (RA.setMinGranularityOuter gran 1) / 64 + r1over1024 = (RA.setMinGranularityOuter gran 1) / 1024 + z = RA.setMinGranularityOuter gran 0 bellardTerms n nMax (mult, r4n, r10n) | n >= nMax = [] | otherwise = diff --git a/src/Data/Number/ER/Real/Arithmetic/LinearSolver.hs b/src/Data/Number/ER/Real/Arithmetic/LinearSolver.hs --- a/src/Data/Number/ER/Real/Arithmetic/LinearSolver.hs +++ b/src/Data/Number/ER/Real/Arithmetic/LinearSolver.hs @@ -19,8 +19,8 @@ where import qualified Data.Number.ER.Real.Approx as RA -import qualified Data.Number.ER.Real.DomainBox as DBox -import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) +import qualified Data.Number.ER.BasicTypes.DomainBox as DBox +import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) import Data.Number.ER.BasicTypes import Data.List @@ -104,10 +104,10 @@ i = DBox.lookup "ER: LinearSolver: splitBox: " k box k = widestVar box incrementGranularity x = - RA.setMinGranularity (RA.getGranularity x + 1) x + RA.setMinGranularityOuter (RA.getGranularity x + 1) x widestVar box = - fst $ DBox.bestSplit box + fst $ DBox.bestSplit box width i = snd $ RA.bounds (iR-iL) diff --git a/src/Data/Number/ER/Real/Arithmetic/Taylor.hs b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs --- a/src/Data/Number/ER/Real/Arithmetic/Taylor.hs +++ b/src/Data/Number/ER/Real/Arithmetic/Taylor.hs @@ -13,8 +13,9 @@ module Data.Number.ER.Real.Arithmetic.Taylor where import qualified Data.Number.ER.Real.Approx as RA -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.BasicTypes +import Data.Number.ER.Misc erTaylor_R @@ -41,7 +42,7 @@ -> ira -> ira erTaylor_R_FullArgs coefSeq derivBounds n a gran x = - rec_apTaylor (RA.setMinGranularity gran 0) 0 + rec_apTaylor (RA.setMinGranularityOuter gran 0) 0 where rec_apTaylor i j | n > j = (coefSeq(j)) + @@ -59,9 +60,24 @@ -> ira -> ira erExp_Tay_Opt_R ix x - | RA.isEmpty x = RA.emptyApprox - | x `RA.refines` (-1/0) = 0 -- -infty is not handled well by the Taylor formula - | otherwise = 1 + (te ix x (RA.setMinGranularity gran 1)) + | x `RA.refines` (-RA.plusInfinity) = +-- unsafePrintReturn +-- ( +-- "erExp_Tay_Opt_R (x `RA.refines` (-RA.plusInfinity)): " +-- ++ "\n x = " ++ show x +-- ++ "\n ix = " ++ show ix +-- ++ "\n result = " +-- ) $ + 0 -- -infty is not handled well by the Taylor formula + | otherwise = +-- unsafePrintReturn +-- ( +-- "erExp_Tay_Opt_R: " +-- ++ "\n x = " ++ show x +-- ++ "\n ix = " ++ show ix +-- ++ "\n result = " +-- ) $ + 1 + (te ix x (RA.setMinGranularityOuter gran 1)) where gran = effIx2gran ix te steps x i @@ -87,13 +103,13 @@ | xFloor == EI.MinusInfinity = 0 | otherwise = - ((RA.setMinGranularity gran 26)/10) ^^ xFloor + ((RA.setMinGranularityOuter gran 26)/10) ^^ xFloor -- lower estimate of e^x pow28xCeiling | xCeiling == EI.PlusInfinity = - (1/0) + (RA.plusInfinity) | otherwise = - ((RA.setMinGranularity gran 28)/10) ^^ xCeiling + ((RA.setMinGranularityOuter gran 28)/10) ^^ xCeiling -- upper estimate of e^x {- @@ -109,7 +125,7 @@ -> ira -> ira erSine_Tay_Opt_R ix x = - taylor_seg ix x (RA.setMinGranularity gran 1) + taylor_seg ix x (RA.setMinGranularityOuter gran 1) where gran = effIx2gran ix taylor_seg i x n -- 'i' for iterator @@ -127,7 +143,7 @@ => EffortIndex -> ira -> ira -erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularity gran 1) +erCosine_Tay_Opt_R ix x = taylor_seg ix x (RA.setMinGranularityOuter gran 1) where gran = effIx2gran ix taylor_seg i x n -- 'i' for iterator diff --git a/src/Data/Number/ER/Real/Base.hs b/src/Data/Number/ER/Real/Base.hs --- a/src/Data/Number/ER/Real/Base.hs +++ b/src/Data/Number/ER/Real/Base.hs @@ -22,15 +22,15 @@ where import Data.Number.ER.BasicTypes -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Typeable {-| This class is an abstraction of a subset of real numbers - with upwards rounded operations. + with *upwards rounded* operations. -} -class (Fractional rb, Ord rb) => ERRealBase rb +class (RealFrac rb, Ord rb) => ERRealBase rb where typeName :: rb -> String initialiseBaseArithmetic :: rb -> IO () @@ -49,9 +49,14 @@ isERNaN :: rb -> Bool erNaN :: rb isPlusInfinity :: rb -> Bool + isMinusInfinity :: rb -> Bool + isMinusInfinity = isPlusInfinity . negate plusInfinity :: rb minusInfinity :: rb minusInfinity = - plusInfinity + fromIntegerUp :: Integer -> rb + fromIntegerDown :: Integer -> rb + fromIntegerDown i = negate $ fromIntegerUp $ - i fromDouble :: Double -> rb toDouble :: rb -> Double fromFloat :: Float -> rb diff --git a/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs b/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs --- a/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs +++ b/src/Data/Number/ER/Real/Base/CombinedMachineAP.hs @@ -23,7 +23,7 @@ where import qualified Data.Number.ER.Real.Base as B -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Real.Base.MachineDouble import Data.Number.ER.Real.Base.Float import Data.Number.ER.BasicTypes @@ -232,6 +232,7 @@ isPlusInfinity = op1ERMachineAP (== 1/0) B.isPlusInfinity plusInfinity = B.fromDouble $ 1/0 + fromIntegerUp = fromInteger fromDouble d = ERMachineAPMachineDouble (B.defaultGranularity (0 :: b)) d toDouble = op1ERMachineAP id B.toDouble diff --git a/src/Data/Number/ER/Real/Base/Float.hs b/src/Data/Number/ER/Real/Base/Float.hs --- a/src/Data/Number/ER/Real/Base/Float.hs +++ b/src/Data/Number/ER/Real/Base/Float.hs @@ -19,8 +19,8 @@ ) where -import qualified Data.Number.ER.ExtendedInteger as EI -import Data.Number.ER.PlusMinus +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI +import Data.Number.ER.BasicTypes.PlusMinus import Data.Number.ER.Misc import Data.Number.ER.BasicTypes import qualified Data.Number.ER.Real.Base as B @@ -494,7 +494,8 @@ erNaN = ERFloatNaN (B.defaultGranularity zero) isPlusInfinity (ERFloatInfty _ Plus) = True isPlusInfinity _ = False - plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus + plusInfinity = ERFloatInfty (B.defaultGranularity zero) Plus + fromIntegerUp i = fromInteger i fromDouble d | isNaN d = ERFloatNaN (B.defaultGranularity zero) | otherwise = (fromRational . toRational) d diff --git a/src/Data/Number/ER/Real/Base/MPFR.hs b/src/Data/Number/ER/Real/Base/MPFR.hs --- a/src/Data/Number/ER/Real/Base/MPFR.hs +++ b/src/Data/Number/ER/Real/Base/MPFR.hs @@ -24,7 +24,7 @@ where import qualified Data.Number.ER.Real.Base as B -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Misc import Data.Binary @@ -65,6 +65,7 @@ erNaN = 0/0 isPlusInfinity x = M.isInfinite x && x > 0 plusInfinity = 1/0 + fromIntegerUp = fromInteger fromDouble = M.fromDouble M.Up 53 toDouble = M.toDouble M.Up fromFloat = B.fromDouble . fromRational . toRational diff --git a/src/Data/Number/ER/Real/Base/MachineDouble.hs b/src/Data/Number/ER/Real/Base/MachineDouble.hs --- a/src/Data/Number/ER/Real/Base/MachineDouble.hs +++ b/src/Data/Number/ER/Real/Base/MachineDouble.hs @@ -12,14 +12,14 @@ Make 'Double' an instance of 'B.ERRealBase' as much as possible. -} -module Data.Number.ER.Real.Base.MachineDouble +module Data.Number.ER.Real.Base.MachineDouble ( initMachineDouble ) where import qualified Data.Number.ER.Real.Base as B -import qualified Data.Number.ER.ExtendedInteger as EI +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI import Data.Number.ER.Misc import Foreign.C @@ -78,6 +78,8 @@ fromInteger $ intLogUp 2 $ ceiling d | d < 1 = negate $ fromInteger $ intLogUp 2 $ ceiling $ recip d + | otherwise = + error $ "ER.Real.Base.MachineDouble: getApproxBinaryLog: illegal argument " ++ show d getGranularity _ = 53 setMinGranularity _ = id setGranularity _ = id @@ -86,6 +88,14 @@ erNaN = 0/0 isPlusInfinity f = isInfinite f && f > 0 plusInfinity = 1/0 + fromIntegerUp i + | i <= floor nearest = nearest + | otherwise = nearestIncreased + where + nearestCeil = ceiling nearest + nearest = fromInteger i + nearestIncreased = encodeFloat (s+1) e + (s,e) = decodeFloat nearest fromDouble = fromRational . toRational toDouble = fromRational . toRational fromFloat = fromRational . toRational diff --git a/src/Data/Number/ER/Real/Base/Rational.hs b/src/Data/Number/ER/Real/Base/Rational.hs --- a/src/Data/Number/ER/Real/Base/Rational.hs +++ b/src/Data/Number/ER/Real/Base/Rational.hs @@ -24,8 +24,8 @@ import Prelude hiding (isNaN) import qualified Data.Number.ER.Real.Base as B -import qualified Data.Number.ER.ExtendedInteger as EI -import Data.Number.ER.PlusMinus +import qualified Data.Number.ER.BasicTypes.ExtendedInteger as EI +import Data.Number.ER.BasicTypes.PlusMinus import Data.Number.ER.Misc import Data.Ratio @@ -226,6 +226,7 @@ isPlusInfinity (Infinity Plus) = True isPlusInfinity _ = False plusInfinity = Infinity Plus + fromIntegerUp = fromInteger fromDouble = fromRational . toRational toDouble (Infinity Plus) = 1/0 toDouble (Infinity Minus) = -1/0 diff --git a/src/Data/Number/ER/Real/Base/Tests/Generate.hs b/src/Data/Number/ER/Real/Base/Tests/Generate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/ER/Real/Base/Tests/Generate.hs @@ -0,0 +1,90 @@ +{-| + Module : Data.Number.ER.Real.Base.Tests.Generate + Description : (testing) generating base real numbers + Copyright : (c) 2009 Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + Generic instances of 'Arbitrary' class for generating (almost) random instances + according to different distributions. +-} +module Data.Number.ER.Real.Base.Tests.Generate where + +import qualified Data.Number.ER.Real.Base as B +import Data.Number.ER.BasicTypes + +import Test.QuickCheck + +newtype BGran20 b = BGran20 b deriving Show +newtype BGran100 b = BGran100 b deriving Show +newtype BGran1000 b = BGran1000 b deriving Show + +instance (B.ERRealBase b) => Arbitrary (BGran20 b) + where + arbitrary = + do + gran <- choose (8,20) + (f1,f2,f3) <- arbitrary + pow <- choose (-10,10) + return $ BGran20 $ constructB gran (f1,f2,f3) pow + coarbitrary _ = + error "ER.Real.Base: Tests: coarbitrary not implemented" + +constructB :: + (B.ERRealBase b) => + Granularity -> + (Double, Double, Double) -> + Int -> + b +constructB gran (f1,f2,f3) pow = + (b1/b2) ^^ pow + b3 + where + [b1,b2,b3] = map cvt [f1,f2,f3] + cvt f = B.setGranularity gran $ B.fromDouble f + +instance (B.ERRealBase b) => Arbitrary (BGran100 b) + where + arbitrary = + sized arbitrarySized + where + arbitrarySized n + | n <= 28 = + do + (BGran20 b) <- arbitrary + return (BGran100 b) + | otherwise = + do + gran <- choose (30,100) + (f1,f2,f3) <- arbitrary + pow <- choose (-100,100) + return $ BGran100 $ constructB gran (f1,f2,f3) pow + coarbitrary _ = + error "ER.Real.Base: Tests: coarbitrary not implemented" + +instance (B.ERRealBase b) => Arbitrary (BGran1000 b) + where + arbitrary = + sized arbitrarySized + where + arbitrarySized n + | n <= 28 = + do + (BGran20 b) <- arbitrary + return (BGran1000 b) + | n <= 68 = + do + (BGran100 b) <- arbitrary + return (BGran1000 b) + | otherwise = + do + gran <- choose (400,1000) + (f1,f2,f3) <- arbitrary + pow <- choose (-10000,10000) + return $ BGran1000 $ constructB gran (f1,f2,f3) pow + coarbitrary _ = + error "ER.Real.Base: Tests: coarbitrary not implemented" + + diff --git a/src/Data/Number/ER/Real/DomainBox.hs b/src/Data/Number/ER/Real/DomainBox.hs deleted file mode 100644 --- a/src/Data/Number/ER/Real/DomainBox.hs +++ /dev/null @@ -1,181 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-| - Module : Data.Number.ER.Real.DomainBox - Description : portions of many-dimensional domains - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mik@konecny.aow.cz - Stability : experimental - Portability : portable - - Abstractions of the 'Box' datatype, often used to represent - sections of multi-dimensional function domains. - - To be imported qualified, usually with prefix DBox. - - VariableID(..) and DomainBox - are usually imported separately and not qualified. --} -module Data.Number.ER.Real.DomainBox -( - VariableID(..), - DomainBox(..), - DomainBoxMappable(..), - DomainIntBox(..) -) -where - -import qualified Data.Number.ER.Real.Approx as RA - -import Data.Number.ER.BasicTypes - -import qualified Data.Set as Set -import qualified Data.Map as Map - -import Prelude hiding (lookup) - - -{-| - A class abstracting a type of variable identifiers - for axes in function domains, polynomials etc. --} -class (Ord varid) => VariableID varid - where - newVarID :: Set.Set varid -> varid - defaultVar :: varid - defaultVar = newVarID Set.empty - showVar :: varid -> String - -{-| - A class abstracting a type of many-dimensional points, intervals - or anything indexed by a subset of dimensions. - - More generally, this class abstracts most of 'Data.Map.Map'. --} -class (VariableID varid) => DomainBox box varid val - | box -> varid val, varid val -> box - where - noinfo :: box - isNoinfo :: box -> Bool - size :: box -> Int - {-| constructor using 'defaultVar' -} - unary :: val -> box - singleton :: varid -> val -> box - toList :: box -> [(varid, val)] - fromList :: [(varid, val)] -> box - toAscList :: box -> [(varid, val)] - fromAscList :: [(varid, val)] -> box --- toMap :: box -> Map.Map varid val --- fromMap :: Map.Map varid val -> box - compare :: (val -> val -> Ordering) -> box -> box -> Ordering - insert :: varid -> val -> box -> box - insertWith :: (val -> val -> val) -> varid -> val -> box -> box - delete :: varid -> box -> box - member :: varid -> box -> Bool - notMember :: varid -> box -> Bool - union :: box -> box -> box - unionWith :: (val -> val -> val) -> box -> box -> box - keys :: box -> [varid] - elems :: box -> [val] - filter :: (val -> Bool) -> box -> box - fold :: (val -> a -> a) -> a -> box -> a - foldWithKey :: (varid -> val -> a -> a) -> a -> box -> a - {-| - for all variables that appear in both boxes, - apply the function and add the result to the list - -} - zipWith :: (val -> val -> a) -> box -> box -> [(varid, a)] - {-| - For all variables that appear in either of the two boxes, - apply the function and add the result to the list. - - Supply the default value when the variable is missing from either box. - -} - zipWithDefault :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] - {-| - For all variables that appear in the first box, - apply the function and add the result to the list. - - Supply the default value when the variable is missing from the second box. - -} - zipWithDefaultSecond :: val -> (val -> val -> a) -> box -> box -> [(varid, a)] - findWithDefault :: val -> varid -> box -> val - {-| - Pick the extents of a single variable in a domain box. - If there is no information for this variable, assume the - variable ranges over the whole real line. - -} - lookup :: - String {-^ identification of caller location to use in error messages -} -> - varid -> - box -> - val - -{-| - A class linking two domain box types that share the - index type so that boxes of the two types can be - converted etc. --} -class (DomainBox box1 varid val1, DomainBox box2 varid val2) => - DomainBoxMappable box1 box2 varid val1 val2 - where - map :: (val1 -> val2) -> box1 -> box2 - mapWithKey :: (varid -> val1 -> val2) -> box1 -> box2 - intersectionWith :: (val1 -> val2 -> val1) -> box1 -> box2 -> box1 - difference :: box1 -> box2 -> box1 - -{-| - A class abstracting a type of many-dimensional intervals. --} -class (DomainBox box varid ira) => DomainIntBox box varid ira - | box -> varid ira, varid ira -> box - where - {-| - Check whether the two domains specify the same - interval for each variable that they share. - -} - compatible :: - box -> - box -> - Bool - {-| - Assuming that two domains are compatible, take the - most information from both of the domains about the - ranges of variables. - -} - unify :: - String {-^ identification of caller location to use in error messages -} -> - box -> - box -> - box - {-| - Find the variable with the largest interval - and return it together with the default splitting point - in its domain. - -} - bestSplit :: - box -> - (varid, ira) - split :: - box -> - varid -> - ira -> - (box, box) - classifyPosition :: - box {-^ domain @d1@ -} -> - box {-^ domain @d2@ -} -> - (Bool, Bool, Bool, Bool) - {-^ - Answers to these (mutually exclusive) questions: - - * is @d1@ outside and /not/ touching @d2@? - - * is @d1@ outside and touching @d2@? - - * is @d1@ intersecting and not inside @d2@? - - * is @d1@ inside @d2@? - -} - diff --git a/src/Data/Number/ER/Real/DomainBox/IntMap.hs b/src/Data/Number/ER/Real/DomainBox/IntMap.hs deleted file mode 100644 --- a/src/Data/Number/ER/Real/DomainBox/IntMap.hs +++ /dev/null @@ -1,196 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-| - Module : Data.Number.ER.Real.DomainBox.IntMap - Description : implementation of DomainBox based on Data.IntMap - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mik@konecny.aow.cz - Stability : experimental - Portability : portable - - A simple implementation of the 'VariableID' and 'DomainBox' classes. --} -module Data.Number.ER.Real.DomainBox.IntMap -( - VarID, Box -) -where - -import qualified Data.Number.ER.Real.Approx as RA -import qualified Data.Number.ER.Real.DomainBox as DBox -import Data.Number.ER.Real.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) - -import Data.Number.ER.Misc - -import qualified Data.IntMap as IMap -import qualified Data.Set as Set - -type VarID = Int -type Box ira = IMap.IntMap ira - -instance VariableID VarID - where - newVarID prevVars - | Set.null prevVars = 0 - | otherwise = - 1 + (Set.findMax prevVars) - showVar v - | v == 0 = "x" - | otherwise = "x" ++ show v - -instance (Show val) => (DomainBox (Box val) VarID val) - where - noinfo = IMap.empty - isNoinfo = IMap.null - size = IMap.size - unary r = IMap.singleton defaultVar r - singleton = IMap.singleton - toList = IMap.toList - fromList = IMap.fromList - toAscList = IMap.toAscList - fromAscList = IMap.fromAscList --- toMap = id --- fromMap = id - compare compareVals b1 b2 = - compareListsWith comparePairs (IMap.toList b1) (IMap.toList b2) - where - comparePairs (k1,v1) (k2,v2) = - compareComposeMany - [ - compare k1 k2, - compareVals v1 v2 - ] - - insert = IMap.insert - insertWith = IMap.insertWith - delete = IMap.delete - member = IMap.member - notMember = IMap.notMember - union = IMap.union - unionWith = IMap.unionWith - elems = IMap.elems - keys = IMap.keys - filter = IMap.filter - fold = IMap.fold - foldWithKey = IMap.foldWithKey - zipWith f b1 b2 = - applyF (IMap.toAscList b1) (IMap.toAscList b2) - where - applyF [] _ = [] - applyF _ [] = [] - applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) - | k1 == k2 = - (k1, f v1 v2) : (applyF rest1 rest2) - | k1 < k2 = applyF rest1 bl2 - | otherwise = applyF bl1 rest2 - zipWithDefault defaultValue f b1 b2 = - applyF (IMap.toAscList b1) (IMap.toAscList b2) - where - applyF [] [] = [] - applyF bl1@((k1,v1):rest1) [] = - (k1, f v1 defaultValue) : (applyF rest1 []) - applyF [] bl2@((k2,v2):rest2) = - (k2, f defaultValue v2) : (applyF [] rest2) - applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) - | k1 == k2 = - (k1, f v1 v2) : (applyF rest1 rest2) - | k1 < k2 = - (k1, f v1 defaultValue) : (applyF rest1 bl2) - | otherwise = - (k2, f defaultValue v2) : (applyF bl1 rest2) - zipWithDefaultSecond defaultValue f b1 b2 = - applyF (IMap.toAscList b1) (IMap.toAscList b2) - where - applyF [] _ = [] - applyF bl1@((k1,v1):rest1) [] = - (k1, f v1 defaultValue) : (applyF rest1 []) - applyF bl1@((k1,v1):rest1) bl2@((k2,v2):rest2) - | k1 == k2 = - (k1, f v1 v2) : (applyF rest1 rest2) - | k1 < k2 = - (k1, f v1 defaultValue) : (applyF rest1 bl2) - | otherwise = - applyF bl1 rest2 - findWithDefault = IMap.findWithDefault - lookup locspec var dom = - IMap.findWithDefault err var dom - where - err = - error $ - locspec ++ "DomainBox.IntMap lookup: domain box " ++ show dom - ++ " ignores variable " ++ show var - -instance (Show val1, Show val2) => - (DomainBoxMappable (Box val1) (Box val2) VarID val1 val2) - where - map = IMap.map - mapWithKey = IMap.mapWithKey - intersectionWith = IMap.intersectionWith - difference = IMap.difference - -instance (RA.ERIntApprox ira) => DomainIntBox (Box ira) VarID ira - where - compatible dom1 dom2 = - foldl (&&) True $ map snd $ - DBox.zipWith RA.equalIntervals dom1 dom2 - unify locspec dom1 dom2 - | DBox.compatible dom1 dom2 = - IMap.union dom1 dom2 - | otherwise = - error $ - locspec ++ "incompatible domains " ++ show dom1 ++ " and " ++ show dom2 - bestSplit dom = - (var, pt) - where - pt = - RA.defaultBisectPt varDom - (_, (varDom, var)) = - foldl findWidestVar (0, err) $ IMap.toList dom - err = - error $ "DomainBox: bestSplit: failed to find a split for " ++ show dom - findWidestVar (prevWidth, prevRes) (v, d) - | currWidth `RA.leqSingletons` prevWidth = (prevWidth, prevRes) - | otherwise = (currWidth, (d, v)) - where - currWidth = snd $ RA.bounds $ domHI - domLO - (domLO, domHI) = RA.bounds d - split dom var pt = - (IMap.insert var varDomL dom, IMap.insert var varDomR dom) - where - varDomL = varDomLO RA.\/ varDomMid - varDomR = varDomMid RA.\/ varDomHI - (varDomLO, varDomMid, varDomHI, _) = RA.exactMiddle varDom - varDom = DBox.lookup "DomainBox.IntMap: split: " var dom - classifyPosition dom sdom = - (away, touch, intersect, inside) - where - (away, touch, inside, intersect) = - foldl addDimension (True, True, True, False) awayTouchInsides - addDimension - (prevAway, prevTouch, prevInside, prevIntersect) - (thisAway, thisTouch, thisInside, thisIntersect) = - (prevAway && thisAway, - (prevTouch || prevAway) && (thisTouch || thisAway) && (prevTouch || thisTouch), - prevInside && thisInside, - prevIntersect || thisIntersect) - awayTouchInsides = - map snd $ - DBox.zipWith classifyRA dom sdom - classifyRA d sd = - (outsideNoTouch, outsideTouch, inside, - not (outsideNoTouch || outsideTouch || inside)) - where - outsideNoTouch = sdR < dL || dR < sdL - outsideTouch = sdR == dL || dR == sdL - inside = sdL =< dL && dR =< sdR - (==) = RA.eqSingletons - (<) = RA.ltSingletons - (=<) = RA.leqSingletons - (dL, dR) = RA.bounds d - (sdL, sdR) = RA.bounds sd - - - diff --git a/tests/Demo.hs b/tests/Demo.hs deleted file mode 100644 --- a/tests/Demo.hs +++ /dev/null @@ -1,149 +0,0 @@ -{-# LANGUAGE CPP #-} -{-| - Module : Main - Description : simple examples of using AERN-Real - Copyright : (c) Michal Konecny - License : BSD3 - - Maintainer : mik@konecny.aow.cz - Stability : experimental - Portability : portable - - Simple examples of using AERN-Real --} -module Main where - -import qualified Data.Number.ER.Real as AERN -import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq) - -#ifdef USE_MPFR ---type B = AERN.BAP -- use pure Haskell floats -type B = AERN.BMAP -- use combination of double and pure Haskell floats ---type B = AERN.BMPFR -- use MPFR floats -#else ---type B = AERN.BAP -- use pure Haskell floats -type B = AERN.BMAP -- use combination of double and pure Haskell floats -#endif -type RA = AERN.RA B -type IRA = AERN.IRA B -type R = ConvergRealSeq IRA - -one :: R -one = 1 - -two :: R -two = 2 - -piSeq :: R -piSeq = ConvergRealSeq $ AERN.pi - -seqExp = convertFuncRA2Seq $ AERN.exp -seqSine = convertFuncRA2Seq $ AERN.sin -seqCosine = convertFuncRA2Seq $ AERN.cos - -main = - do - AERN.initialiseBaseArithmetic (0 :: RA) - putStrLn "****************************" - putStrLn "Testing interval arithmetic:" - putStrLn "****************************" - putStrLn "**** Fractions:" - putStrLn $ - "(default granularity, show internals) 1/3 =\n " ++ - AERN.showApprox 30 True True (1/3 :: RA) - putStrLn $ - "(granularity 50, show internals) 1/3 =\n " ++ - AERN.showApprox 30 True True ((AERN.setGranularity 50 1/3) :: RA) - putStrLn $ - "(granularity 100, show internals) 1/3 =\n " ++ - AERN.showApprox 40 True True ((AERN.setGranularity 100 1/3) :: RA) - putStrLn $ - "(granularity 100, do not show internals) 1/3 =\n " ++ - AERN.showApprox 40 True False ((AERN.setGranularity 100 1/3) :: RA) - putStrLn $ - "(granularity 100, default show) 1/3 =\n " ++ - show ((AERN.setGranularity 100 1/3) :: RA) - putStrLn "**** Exp:" - putStrLn $ - "(effort 5, granularity 50) exp 1 =\n " ++ - (show $ AERN.exp 5 (AERN.setGranularity 50 (1::RA))) - putStrLn $ - "(effort 10, granularity 50) exp 1 =\n " ++ - (show $ AERN.exp 10 (AERN.setGranularity 50 (1::RA))) - putStrLn $ - "(effort 10, granularity 100) exp 1 =\n " ++ - (show $ AERN.exp 10 (AERN.setGranularity 100 (1::RA))) - putStrLn $ - "(effort 20, granularity 50) exp 1 =\n " ++ - (show $ AERN.exp 20 (AERN.setGranularity 50 (1::RA))) - putStrLn $ - "(effort 20, granularity 100) exp 1 =\n " ++ - (show $ AERN.exp 20 (AERN.setGranularity 100 (1::RA))) - putStrLn "**** Pi:" - putStrLn $ - "(effort 10) pi =\n " ++ - (show $ (AERN.pi 10 :: RA)) - putStrLn $ - "(effort 50) pi =\n " ++ - (AERN.showApprox 20 True False $ (AERN.pi 50 :: RA)) - putStrLn $ - "(effort 100) pi =\n " ++ - (AERN.showApprox 35 True False $ (AERN.pi 100 :: RA)) - putStrLn $ - "(effort 200) pi =\n " ++ - (AERN.showApprox 65 True False $ (AERN.pi 200 :: RA)) - putStrLn $ - "(effort 400) pi =\n " ++ - (AERN.showApprox 125 True False $ (AERN.pi 400 :: RA)) - putStrLn "**** Sine:" - putStrLn $ - "(effort 10, granularity 50) sin 1 =\n " ++ - (show $ AERN.sin 10 (AERN.setGranularity 50 (1::RA))) - putStrLn $ - "(effort 10, granularity 100) sin 1 =\n " ++ - (show $ AERN.sin 10 (AERN.setGranularity 100 (1::RA))) - putStrLn "**** Integration:" - putStrLn $ - "(effort 10, granularity 50) integrate exp 0 1 =\n " ++ - (show $ AERN.integrateContAdapt_R AERN.exp 10 0 (AERN.setGranularity 50 (1::RA))) - putStrLn $ - "(effort 20, granularity 50) integrate exp 0 1 =\n " ++ - (show $ AERN.integrateContAdapt_R AERN.exp 20 0 (AERN.setGranularity 50 (1::RA))) --- putStrLn $ --- "(effort 30, granularity 50) integrate exp 0 1 =\n " ++ --- (show $ AERN.integrateContAdapt_R AERN.exp 30 0 (AERN.setGranularity 50 (1::RA))) - putStrLn "*****************************" - putStrLn "Testing convergent sequences:" - putStrLn "*****************************" --- putStrLn $ "1 =\n " ++ show one --- putStrLn $ "1 + 2 =\n " ++ (show $ one + two) - putStrLn "**** Fractions:" - putStrLn $ - "(precision 20) 1/3 =\n " ++ - (AERN.showConvergRealSeqAuto 20 $ one / 3) - putStrLn $ - "(precision 20) 100000000001/300000000000 =\n " ++ - (AERN.showConvergRealSeqAuto 20 $ (one + 100000000000)/300000000000 ) - putStrLn $ - "100000000001/300000000000 =? 1/3:\n " ++ - (show $ one/3 == 100000000001/300000000000) --- putStrLn $ "abs -1 = " ++ (show $ abs (- one)) --- putStrLn $ "neg 2 = " ++ (show $ negate two) --- putStrLn $ "1 + 2 = " ++ (show $ one + 2) - putStrLn "**** Elementary:" - putStrLn $ - "(precision 30) exp 1 =\n " ++ - (AERN.showConvergRealSeqAuto 30 $ seqExp one) - putStrLn $ - "(precision 500) pi =\n " ++ - (AERN.showConvergRealSeqAuto 500 $ piSeq) - putStrLn $ - "(precision 30) cosine(1) =\n " ++ - (AERN.showConvergRealSeqAuto 30 $ seqCosine one) - putStrLn $ - "(precision 30) sine(1) =\n " ++ - (AERN.showConvergRealSeqAuto 30 $ seqSine one) - putStrLn "**** Integration:" - putStrLn $ -- very slow for precision > 4 - "(precision 3) integrate exp 0 1 =\n " ++ - (AERN.showConvergRealSeqAuto 3 $ AERN.integrateCont AERN.exp 0 one) diff --git a/tests/Matrix.hs b/tests/Matrix.hs deleted file mode 100644 --- a/tests/Matrix.hs +++ /dev/null @@ -1,384 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE DeriveDataTypeable #-} -module Main - -where - -import qualified Data.Number.ER.Real as AERN -import Data.Number.ER.BasicTypes -import Data.Number.ER.Misc - -import Data.Maybe -import qualified Data.List as List -import qualified Data.Map as Map - -import qualified Data.Array.IArray as IAr -import qualified Data.Array.MArray as MAr -import qualified Data.Array.ST as STAr -import qualified Data.Ix as Ix -import qualified Data.Array.Base as BAr - -import Control.Monad.ST -import GHC.Arr - -#ifdef USE_MPFR -type B = AERN.BMPFR -- use MPFR floats -#else -type B = AERN.BAP -- use pure Haskell floats -#endif -type RA = AERN.RA B -type IRA = AERN.IRA B - -testMatrixN = 100 -incrementGran = (+) 50 - --- Hilbert 100x100 matrix: -addOneDiag = False -targetPrec = 167 -- approx 50 decimal digits after the point -initialGran = 2050 -- 100x100 ---initialGran = 2388 -- 100x100 Norbert's ---initialGran = 750 -- 50x50 ---initialGran = 300 -- 10x10 - ---targetPrec = 34 -- approx 10 decimal digits after the point ---initialGran = 1350 ---initialGran = 50 -- 50x50 - --- Hilbert matrix + 1: ---addOneDiag = True ---targetPrec = 167 -- approx 50 decimal digits after the point ---initialGran = 200 - ---targetPrec = 34 -- approx 10 decimal digits after the point ---initialGran = 50 - -main = - do - AERN.initialiseBaseArithmetic (0 :: RA) - putStrLn $ - "Inverting the " ++ show n ++ "x" ++ show n ++ " Hilbert matrix " - ++ "with target binary precision " ++ show targetPrec ++ "..." --- putStrLn $ --- "sorted matrix elements = \n" ++ (unlines $ map show elemsSortedByPrec) - putStrLn $ - "sum of all elements in inverted matrix = " ++ show (sum elems) --- putStrLn $ show (Matrix n n rarr) - where - n = testMatrixN - elems = IAr.elems rarr - elemsSortedByPrec = - List.sortBy comparePrec elems - where - comparePrec a b = - compare aPrecLO bPrecLO - where - aPrecLO = fst $ AERN.bounds $ aHI - aLO - (aLO, aHI) = AERN.bounds a - bPrecLO = fst $ AERN.bounds $ bHI - bLO - (bLO, bHI) = AERN.bounds b - rarr = - STAr.runSTArray $ - do - mInv@(Matrix _ _ rowsInv) <- - invert testMatrix --- m <- testMatrix initialGran --- mUnit@(Matrix _ _ rowsUnit) <- multM m mInv - return rowsInv - - -testMatrix :: - Granularity -> - ST s (STMatrix s IRA) -testMatrix gran = - do - marr <- MAr.newArray ((1,1),(n,n)) 0 - mapM (updateCell marr) assocsGran - return $ Matrix n n marr - where - assocsGran = map (mapSnd $ AERN.setMinGranularity gran) assocs - assocs = --- assocsMini - assocsHilbert gran n - assocsMini = - [((1,1),1), - ((1,2),3), - ((2,1),2), - ((2,2),0) - ] - n = testMatrixN - updateCell marr (ix, el) = - do - unsafeMatrixWrite marr n ix el - -assocsHilbert gran n = - [((i,j), coeff i j)| i <- [1..n], j <- [1..n]] - where - coeff i j - | addOneDiag && i == j = - 1 + oneOverIplusJ - | otherwise = - oneOverIplusJ - where - oneOverIplusJ = - recip $ (AERN.setMinGranularity gran $ iRA + jRA + 1) - iRA = fromInteger $ toInteger i - jRA = fromInteger $ toInteger j - - ---invert :: --- Precision -> --- () -> -invert getMatrix = - do - gaussElim getMatrixI - where - n = testMatrixN - getMatrixI gran = - do - m <- getMatrix gran - mI <- addIdentity m - return mI - -gaussElim getMatrix = - elimWithMinGran initialGran - where - elimWithMinGran workingGran = - do - mI@(Matrix colN rowN _) <- getMatrix workingGran - idPerm <- MAr.newListArray (1,rowN) [1..rowN] - elimAtRow mI 1 idPerm - where - elimAtRow mI@(Matrix colN rowN mIarr) i perm = - do - success <- ensureNonZeroDiag -- make sure (i,i) is non-zero by permuting - case success of - False -> -- failed - all elements contain 0 -> try larger granularity - unsafePrint ("failed to divide at granularity " ++ show workingGran) $ - elimWithMinGran (incrementGran workingGran) - True -> - do - normaliseRow - eliminateColumn - case i == rowN of - True -> - do - mInv <- permuteRowsDropCols perm testMatrixN mI - mPrec <- getMatrixPrecision mInv - case mPrec >= targetPrec of - False -> -- resulting precision insufficient - unsafePrint - ("insufficient precision " ++ show mPrec ++ - " at granularity " ++ show workingGran) $ - elimWithMinGran (incrementGran workingGran) - True -> - unsafePrint - ("precision " ++ show mPrec ++ - " succeeded at granularity " ++ show workingGran) - return mInv - False -> elimAtRow mI (i+1) perm - where - ensureNonZeroDiag = - do - maybeNonZeroIx <- findNonZeroRow - case maybeNonZeroIx of - Nothing -> - return False - Just ii -> - do - case ii > 0 of - True -> swap i (i + ii) perm - False -> return () - return True - findNonZeroRow = - do - elems <- mapM getElemPerm [(i,rowIx) | rowIx <- [i..rowN]] - return $ List.findIndex (\e -> not $ 0 `AERN.refines` e) elems - getElemPerm (colIx,rowIx) = - do - rowIxPerm <- unsafePermRead perm rowIx - unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) - - normaliseRow = - do - rowIxPerm <- unsafePermRead perm i - e <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) - unsafeMatrixWrite mIarr rowN (i, rowIxPerm) 1 - mapM (divideCellBy e rowIxPerm) [(i+1)..colN] - divideCellBy e rowIxPerm colIx = - do - e2 <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) - unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (e2/e) - - eliminateColumn = - do - iRowPerm <- unsafePermRead perm i - mapM (eliminateColumnRow iRowPerm) $ [1..(i-1)] ++ [(i+1)..rowN] - eliminateColumnRow iRowPerm rowIx = - do - rowIxPerm <- unsafePermRead perm rowIx - c <- unsafeMatrixRead mIarr rowN (i, rowIxPerm) -- remember old element for scaling i'th row - unsafeMatrixWrite mIarr rowN (i,rowIxPerm) 0 -- at column i we set 0 - mapM (eliminateColumnRowColumn iRowPerm rowIxPerm c) [(i+1)..colN] - eliminateColumnRowColumn iRowPerm rowIxPerm c colIx = - do - ei <- unsafeMatrixRead mIarr rowN (colIx, iRowPerm) -- at i'th row - er <- unsafeMatrixRead mIarr rowN (colIx, rowIxPerm) -- at current row - unsafeMatrixWrite mIarr rowN (colIx, rowIxPerm) (er - c * ei) -- eliminate by i'th row - - -swap :: - Int -> - Int -> - (STAr.STUArray s Int Int) -> - ST s () -swap i1 i2 perm = - do - a1 <- unsafePermRead perm i1 - a2 <- unsafePermRead perm i2 - unsafePermWrite perm i1 a2 - unsafePermWrite perm i2 a1 - - -unsafePermWrite permArr i e = - do - BAr.unsafeWrite permArr (i - 1) e - -unsafePermRead permArr i = - do - BAr.unsafeRead permArr (i - 1) - - -addIdentity :: - (STMatrix s IRA) -> - ST s (STMatrix s IRA) -addIdentity (Matrix colN rowN marr) = - do --- (_, (colN,rowN)) <- MAr.getBounds marr - mElems <- MAr.getElems marr - mIarr <- MAr.newListArray ((1,1),(colN+rowN,rowN)) $ mElems ++ (idElems rowN) - return $ Matrix (colN + rowN) rowN mIarr - where - idElems m = - 1 : (concat $ replicate (m-1) $ (replicate m 0) ++ [1]) - - -data Matrix marr el = - Matrix - { - mxRowN :: Int, - mxColN :: Int, - mxRows :: marr (ColIx,RowIx) el - } - -type ColIx = Int -type RowIx = Int - -type IMatrix el = - Matrix Array el - -type STMatrix s el = - Matrix (STArray s) el - -instance - (IAr.IArray marr el,-- IAr.IArray marr (marr Int el), - Show el) => - Show (Matrix marr el) - where - show (Matrix colN rowN rows) = - "\nMatrix:\n" ++ - (concat $ map showCol [1..colN]) - where --- (_,(colN,rowN)) = IAr.bounds rows - showCol colIx = - unlines $ - map showCell [(colIx, rowIx) | rowIx <- [1..rowN]] - showCell ix@(colIx, rowIx) = - (show ix) ++ - (replicate colIx '.') ++ - (show $ (IAr.!) rows ix) - -getMatrixPrecision (Matrix _ _ marr) = - do - elems <- MAr.getElems marr - return $ foldl1 min $ map AERN.getPrecision elems - -unsafeMatrixWrite marr rowN (i,j) e = - do - BAr.unsafeWrite marr (rowN*(i-1) + j-1) e --- MAr.writeArray marr (i,j) e - -unsafeMatrixRead marr rowN (i,j) = - do - BAr.unsafeRead marr (rowN*(i-1) + j-1) --- MAr.readArray marr (i,j) - -permuteRowsDropCols :: - (STAr.STUArray s Int Int) -> - Int {-^ drop this many first columns -} -> - (STMatrix s IRA) -> - ST s (STMatrix s IRA) -permuteRowsDropCols perm dropN (Matrix colN rowN marr) = - do --- (_, (colN,rowN)) <- MAr.getBounds marr - (_, permN) <- MAr.getBounds perm - rarr <- MAr.newArray ((1,1),(colN - dropN, permN)) 0 - mapM (copyElem marr rarr rowN) [(colIx, rowIx) | colIx <- [1..colN - dropN], rowIx <- [1..permN]] - return (Matrix (colN - dropN) permN rarr) - where - copyElem marr rarr rowN (colIx, rowIx) = - do - permRowIx <- unsafePermRead perm rowIx - e <- unsafeMatrixRead marr rowN (colIx + dropN, permRowIx) - unsafeMatrixWrite rarr rowN (colIx, rowIx) e - - -addM m1 m2 - | mxColN m1 == mxColN m2 && mxRowN m1 == mxRowN m2 = - do - marr <- MAr.newArray ((1,1),(colN, rowN)) 0 - mapM (addCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]] - return (Matrix colN rowN marr) - | otherwise = - error "Matrix: addM mismatch" - where - colN = mxColN m1 - rowN = mxRowN m1 - marr1 = mxRows m1 - marr2 = mxRows m2 - addCell marr (colIx, rowIx) = - do - elem1 <- unsafeMatrixRead marr1 rowN (colIx, rowIx) - elem2 <- unsafeMatrixRead marr2 rowN (colIx, rowIx) - unsafeMatrixWrite marr rowN (colIx, rowIx) (elem1 + elem2) - -multM m1 m2 - | colN1 == rowN2 = - do - marr <- MAr.newArray ((1,1),(colN, rowN)) 0 - mapM (multCell marr) [(c,r) | c <- [1..colN], r <- [1..rowN]] - return (Matrix colN rowN marr) - | otherwise = - error "Matrix: multM mismatch" - where - colN1 = mxColN m1 - rowN1 = mxRowN m1 - colN2 = mxColN m2 - rowN2 = mxRowN m2 - colN = colN2 - rowN = rowN1 - marr1 = mxRows m1 - marr2 = mxRows m2 - multCell marr (colIx, rowIx) = - do - elems1 <- mapM (getCell1 rowIx) [1..colN1] - elems2 <- mapM (getCell2 colIx) [1..rowN2] - unsafeMatrixWrite marr rowN (colIx, rowIx) (sum $ zipWith (*) elems1 elems2) - getCell1 rowIx colIx = - do - unsafeMatrixRead marr1 rowN1 (colIx, rowIx) - getCell2 rowIx colIx = - do - unsafeMatrixRead marr2 rowN2 (colIx, rowIx) - diff --git a/tests/Pi.hs b/tests/Pi.hs deleted file mode 100644 --- a/tests/Pi.hs +++ /dev/null @@ -1,43 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE UndecidableInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE DeriveDataTypeable #-} -module Main - -where - -import qualified Data.Number.ER.Real as AERN -import Data.Number.ER.Real (ConvergRealSeq(..), convertFuncRA2Seq) -import Data.Number.ER.BasicTypes -import Data.Number.ER.Misc - -import Data.Maybe - -#ifdef USE_MPFR -type B = AERN.BMPFR -- use MPFR floats -#else -type B = AERN.BAP -- use pure Haskell floats ---type B = AERN.BMAP -- use combination of double and pure Haskell floats -#endif -type RA = AERN.RA B -type IRA = AERN.IRA B - - -decimalPrec = 100 ---decimalPrec = 1000 -binaryPrec = - fromInteger $ toInteger $ - snd $ AERN.integerBounds $ - (fromInteger decimalPrec :: RA) * (AERN.log 100 10)/(AERN.log 100 2) - -main = - do - AERN.initialiseBaseArithmetic (0 :: RA) - putStrLn $ - show decimalPrec - ++ " decimal digits of pi = \n" - ++ (AERN.showConvergRealSeqAuto binaryPrec pi) - where - pi :: ConvergRealSeq RA - pi = ConvergRealSeq AERN.pi - diff --git a/tests/RunERIntervalTests.hs b/tests/RunERIntervalTests.hs new file mode 100644 --- /dev/null +++ b/tests/RunERIntervalTests.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE CPP #-} +{-| + Module : Main + Description : laucher for approximated exact real arithmetic tests + Copyright : (c) Michal Konecny + License : BSD3 + + Maintainer : mik@konecny.aow.cz + Stability : experimental + Portability : portable + + An executable for easy automated launch of tests + of approximated exact real arithmetic. +-} +module Main where + +import qualified Data.Number.ER.Real.Approx as RA +import Data.Number.ER.Real.Approx.Tests.Run +import Data.Number.ER.Real.DefaultRepr + +--import Data.Number.ER.Real.Approx.Tests.Properties +--import Data.Number.ER.Real.Approx.Tests.Generate + +main = + do + runRATests "interval-double" sampleRABM (RA.initialiseBaseArithmetic sampleRABM) + runRATests "interval-haskell" sampleRABM (RA.initialiseBaseArithmetic sampleRABAP) + runRATests "interval-haskell-double" sampleRABM (RA.initialiseBaseArithmetic sampleRABMAP) +#ifdef USE_MPFR + runRATests "interval-mpfr" sampleRABM (RA.initialiseBaseArithmetic sampleRABMPFR) +#endif + +sampleRABM :: RA BM +sampleRABAP :: RA BAP +sampleRABMAP :: RA BMAP +sampleRABM = 0 +sampleRABAP = 0 +sampleRABMAP = 0 + +#ifdef USE_MPFR +sampleRABMPFR :: RA BMPFR +sampleRABMPFR = 0 +#endif