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.10.0.1 +Version: 0.10.0.2 Cabal-Version: >= 1.2 Build-Type: Simple License: BSD3 @@ -26,13 +26,13 @@ . For an architectural overview, see module "Data.Number.ER.Real". . - Simple examples of usage can be found in folder @demos@. + Simple examples of usage can be found in folder @examples@. . 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 + examples/Demo.hs examples/Pi.hs examples/Matrix.hs tests/RunERIntervalTests.hs ChangeLog diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,4 @@ +0.10.0.2: 29 July 2009: renamed "demos" folder to "examples" 0.10.0.1: 28 July 2009: fixed a few bugs in meta-data 0.10.0: 28 July 2009 * switching to beta status diff --git a/demos/Demo.hs b/demos/Demo.hs deleted file mode 100644 --- a/demos/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.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 deleted file mode 100644 --- a/demos/Matrix.hs +++ /dev/null @@ -1,385 +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.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 deleted file mode 100644 --- a/demos/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 -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/examples/Demo.hs b/examples/Demo.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/Matrix.hs b/examples/Matrix.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/Pi.hs b/examples/Pi.hs new file mode 100644 --- /dev/null +++ b/examples/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 +