rounded (empty) → 0.1
raw patch · 18 files changed
+2885/−0 lines, 18 filesdep +basedep +ghc-primdep +hgmpsetup-changed
Dependencies added: base, ghc-prim, hgmp, long-double, reflection, rounded, singletons
Files
- CHANGELOG.markdown +9/−0
- LICENSE +165/−0
- README.markdown +67/−0
- Setup.hs +2/−0
- cbits/wrappers.c +79/−0
- rounded.cabal +108/−0
- src/Numeric/MPFR/Raw.hs +25/−0
- src/Numeric/MPFR/Raw/Safe.hs +116/−0
- src/Numeric/MPFR/Raw/Unsafe.hs +116/−0
- src/Numeric/MPFR/Types.hsc +127/−0
- src/Numeric/Rounded.hs +95/−0
- src/Numeric/Rounded/Internal.hs +610/−0
- src/Numeric/Rounded/Interval.hs +717/−0
- src/Numeric/Rounded/Precision.hs +84/−0
- src/Numeric/Rounded/Rounding.hs +93/−0
- src/Numeric/Rounded/Simple.hs +384/−0
- test.hs +58/−0
- test.txt +30/−0
+ CHANGELOG.markdown view
@@ -0,0 +1,9 @@+0.1+-----+* Bindings rewritten as a straight FFI binding without a patched `MPFR`.+* Added `Simple` module.++0.0.1+-----+* Repository initialized based on older work on a package called `precision`+* Daniel Peebles was able to get a patched `MPFR` to work
+ LICENSE view
@@ -0,0 +1,165 @@+ GNU LESSER GENERAL PUBLIC LICENSE+ Version 3, 29 June 2007++ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.+++ This version of the GNU Lesser General Public License incorporates+the terms and conditions of version 3 of the GNU General Public+License, supplemented by the additional permissions listed below.++ 0. Additional Definitions. ++ As used herein, "this License" refers to version 3 of the GNU Lesser+General Public License, and the "GNU GPL" refers to version 3 of the GNU+General Public License.++ "The Library" refers to a covered work governed by this License,+other than an Application or a Combined Work as defined below.++ An "Application" is any work that makes use of an interface provided+by the Library, but which is not otherwise based on the Library.+Defining a subclass of a class defined by the Library is deemed a mode+of using an interface provided by the Library.++ A "Combined Work" is a work produced by combining or linking an+Application with the Library. The particular version of the Library+with which the Combined Work was made is also called the "Linked+Version".++ The "Minimal Corresponding Source" for a Combined Work means the+Corresponding Source for the Combined Work, excluding any source code+for portions of the Combined Work that, considered in isolation, are+based on the Application, and not on the Linked Version.++ The "Corresponding Application Code" for a Combined Work means the+object code and/or source code for the Application, including any data+and utility programs needed for reproducing the Combined Work from the+Application, but excluding the System Libraries of the Combined Work.++ 1. Exception to Section 3 of the GNU GPL.++ You may convey a covered work under sections 3 and 4 of this License+without being bound by section 3 of the GNU GPL.++ 2. Conveying Modified Versions.++ If you modify a copy of the Library, and, in your modifications, a+facility refers to a function or data to be supplied by an Application+that uses the facility (other than as an argument passed when the+facility is invoked), then you may convey a copy of the modified+version:++ a) under this License, provided that you make a good faith effort to+ ensure that, in the event an Application does not supply the+ function or data, the facility still operates, and performs+ whatever part of its purpose remains meaningful, or++ b) under the GNU GPL, with none of the additional permissions of+ this License applicable to that copy.++ 3. Object Code Incorporating Material from Library Header Files.++ The object code form of an Application may incorporate material from+a header file that is part of the Library. You may convey such object+code under terms of your choice, provided that, if the incorporated+material is not limited to numerical parameters, data structure+layouts and accessors, or small macros, inline functions and templates+(ten or fewer lines in length), you do both of the following:++ a) Give prominent notice with each copy of the object code that the+ Library is used in it and that the Library and its use are+ covered by this License.++ b) Accompany the object code with a copy of the GNU GPL and this license+ document.++ 4. Combined Works.++ You may convey a Combined Work under terms of your choice that,+taken together, effectively do not restrict modification of the+portions of the Library contained in the Combined Work and reverse+engineering for debugging such modifications, if you also do each of+the following:++ a) Give prominent notice with each copy of the Combined Work that+ the Library is used in it and that the Library and its use are+ covered by this License.++ b) Accompany the Combined Work with a copy of the GNU GPL and this license+ document.++ c) For a Combined Work that displays copyright notices during+ execution, include the copyright notice for the Library among+ these notices, as well as a reference directing the user to the+ copies of the GNU GPL and this license document.++ d) Do one of the following:++ 0) Convey the Minimal Corresponding Source under the terms of this+ License, and the Corresponding Application Code in a form+ suitable for, and under terms that permit, the user to+ recombine or relink the Application with a modified version of+ the Linked Version to produce a modified Combined Work, in the+ manner specified by section 6 of the GNU GPL for conveying+ Corresponding Source.++ 1) Use a suitable shared library mechanism for linking with the+ Library. A suitable mechanism is one that (a) uses at run time+ a copy of the Library already present on the user's computer+ system, and (b) will operate properly with a modified version+ of the Library that is interface-compatible with the Linked+ Version. ++ e) Provide Installation Information, but only if you would otherwise+ be required to provide such information under section 6 of the+ GNU GPL, and only to the extent that such information is+ necessary to install and execute a modified version of the+ Combined Work produced by recombining or relinking the+ Application with a modified version of the Linked Version. (If+ you use option 4d0, the Installation Information must accompany+ the Minimal Corresponding Source and Corresponding Application+ Code. If you use option 4d1, you must provide the Installation+ Information in the manner specified by section 6 of the GNU GPL+ for conveying Corresponding Source.)++ 5. Combined Libraries.++ You may place library facilities that are a work based on the+Library side by side in a single library together with other library+facilities that are not Applications and are not covered by this+License, and convey such a combined library under terms of your+choice, if you do both of the following:++ a) Accompany the combined library with a copy of the same work based+ on the Library, uncombined with any other library facilities,+ conveyed under the terms of this License.++ b) Give prominent notice with the combined library that part of it+ is a work based on the Library, and explaining where to find the+ accompanying uncombined form of the same work.++ 6. Revised Versions of the GNU Lesser General Public License.++ The Free Software Foundation may publish revised and/or new versions+of the GNU Lesser General Public License from time to time. Such new+versions will be similar in spirit to the present version, but may+differ in detail to address new problems or concerns.++ Each version is given a distinguishing version number. If the+Library as you received it specifies that a certain numbered version+of the GNU Lesser General Public License "or any later version"+applies to it, you have the option of following the terms and+conditions either of that published version or of any later version+published by the Free Software Foundation. If the Library as you+received it does not specify a version number of the GNU Lesser+General Public License, you may choose any version of the GNU Lesser+General Public License ever published by the Free Software Foundation.++ If the Library as you received it specifies that a proxy can decide+whether future versions of the GNU Lesser General Public License shall+apply, that proxy's public statement of acceptance of any version is+permanent authorization for you to choose that version for the+Library.
+ README.markdown view
@@ -0,0 +1,67 @@+rounded+=======++[](http://travis-ci.org/ekmett/rounded)++This package provides properly rounded floating point numbers of arbitrary precision.+It does so by wrapping the GNU MPFR library.++Phantom types carry the information about the precision and rounding mode, letting you treat properly rounded floating+point numbers as instances of `Num` or `Floating`, like any other numeric type in Haskell.++Unlike other attempts to port MPFR to Haskell, this library does not require you to cripple `Integer` performance+or link your code in an unnatural way.++Usage+-----++```haskell+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+import Numeric.Rounded+import Data.Proxy+```++To use a 53 bit significand (the same size as used by a Double), and round down intermediate results:++```haskell+>>> pi :: Rounded TowardZero Double+3.141592653589793+```++We can also round away from zero, or use other rounding modes.++```haskell+>>> pi :: Rounded AwayFromZero Double+3.1415926535897936+```++We can specify the significand size directly using type literals in GHC:++```haskell+>>> kCatalan :: Rounded TowardZero 128+0.915965594177219015054603514932384110773+```++You can also specify a dynamic significand size at runtime:++```haskell+>>> reifyPrecision 512 (\(_ :: Proxy p) -> show (logBase 10 2 :: Rounded TowardNearest p))+"0.3010299956639811952137388947244930267681898814621085413104274611271081892744245094869272521181861720406844771914309953790947678811335235059996923337046956"+```++or a dynamic rounding mode:++```haskell+ghci> reifyRounding TowardZero (\(_ :: Proxy r) -> show (logBase 10 2 :: Rounded r 512))+"0.30102999566398119521373889472449302676818988146210854131042746112710818927442450948692725211818617204068447719143099537909476788113352350599969233370469556"+```++Contact Information+-------------------++Please, feel free to contact me with questions, concerns, or bug fixes.++I can be reached as ekmett via github or as edwardk on the #haskell IRC channel on irc.freenode.net.++-Edward Kmett
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ cbits/wrappers.c view
@@ -0,0 +1,79 @@+#include <mpfr.h>++static int get_flags(void)+{+ return+ ((!!mpfr_underflow_p ()) << 0) |+ ((!!mpfr_overflow_p ()) << 1) |+ ((!!mpfr_divby0_p ()) << 2) |+ ((!!mpfr_nanflag_p ()) << 3) |+ ((!!mpfr_inexflag_p ()) << 4) |+ ((!!mpfr_erangeflag_p()) << 5);+}++static void set_flags(int flags)+{+ if (flags & (1 << 0)) mpfr_set_underflow (); else mpfr_clear_underflow ();+ if (flags & (1 << 1)) mpfr_set_overflow (); else mpfr_clear_overflow ();+ if (flags & (1 << 2)) mpfr_set_divby0 (); else mpfr_clear_divby0 ();+ if (flags & (1 << 3)) mpfr_set_nanflag (); else mpfr_clear_nanflag ();+ if (flags & (1 << 4)) mpfr_set_inexflag (); else mpfr_clear_inexflag ();+ if (flags & (1 << 5)) mpfr_set_erangeflag(); else mpfr_clear_erangeflag();+}++int wrapped_mpfr_get_z(mpz_t rop, mpfr_t op, mpfr_rnd_t rnd, int *flags)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ int retval = mpfr_get_z(rop, op, rnd);+ *flags = get_flags();+ set_flags(flags0);+ return retval;+}++mpfr_exp_t wrapped_mpfr_get_z_2exp(mpz_t rop, mpfr_t op, int *flags)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ mpfr_exp_t retval = mpfr_get_z_2exp(rop, op);+ *flags = get_flags();+ set_flags(flags0);+ return retval;+}++int wrapped_mpfr_set_ld(mpfr_t rop, long double *op, mpfr_rnd_t rnd, int *flags)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ int retval = mpfr_set_ld(rop, *op, rnd);+ *flags = get_flags();+ set_flags(flags0);+ return retval;+}++void wrapped_mpfr_get_ld(long double *rop, mpfr_t op, mpfr_rnd_t rnd, int *flags)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ *rop = mpfr_get_ld(op, rnd);+ *flags = get_flags();+ set_flags(flags0);+}++void wrapped_mpfr_get_ld_2exp(long double *rop, long *e, mpfr_t op, mpfr_rnd_t rnd, int *flags)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ *rop = mpfr_get_ld_2exp(e, op, rnd);+ *flags = get_flags();+ set_flags(flags0);+}++int wrapped_mpfr_cmp_ld(mpfr_t rop, long double *op)+{+ int flags0 = get_flags();+ mpfr_clear_flags();+ int retval = mpfr_cmp_ld(rop, *op);+ set_flags(flags0);+ return retval;+}
+ rounded.cabal view
@@ -0,0 +1,108 @@+name: rounded+version: 0.1+synopsis: Correctly-rounded arbitrary-precision floating-point arithmetic+homepage: https://github.com/ekmett/rounded+bug-reports: https://github.com/ekmett/rounded/issues+license: LGPL+license-file: LICENSE+author: Edward A. Kmett, Daniel G. Peebles, Claude Heiland-Allen+maintainer: Claude Heiland-Allen <claude@mathr.co.uk>+copyright: Copyright (C) 2012-2014 Edward A. Kmett, Daniel G. Peebles;+ Copyright (C) 2013-2018 Claude Heiland-Allen+category: Numeric, Math+build-type: Simple+cabal-version: 1.22+tested-with: GHC == 7.10.3, GHC == 8.0.2, GHC == 8.2.2, GHC == 8.4.4, GHC == 8.6.1+description:+ This package provides numeric instances for MPFR that use+ \"Implicit Configurations\" from+ <http://www.cs.rutgers.edu/~ccshan/prepose/prepose.pdf>+ to choose a 'Rounding' and 'Precision'. For those that do not want to+ use reflection, explicit instances are provided for common precisions+ and for the built-in rounding modes.+ .+ This package should work correctly with GHC 7.10.1 or later.+ .+ >>> import Numeric.Rounded+ >>> :set -XDataKinds+ >>> exp pi :: Rounded TowardZero 512+ 23.140692632779269005729086367948547380266106242600211993445046409524342350690452783516971997067549219675952704801087773144428044414693835844717445879609842++extra-source-files: README.markdown CHANGELOG.markdown test.hs test.txt++source-repository head+ type: git+ location: git://github.com/ekmett/rounded.git++source-repository this+ type: git+ location: git://github.com/ekmett/rounded.git+ tag: rounded-0.1++library+ exposed-modules:+ Numeric.Rounded+ Numeric.Rounded.Interval+ Numeric.Rounded.Simple+ Numeric.MPFR.Types+ Numeric.MPFR.Raw+ Numeric.MPFR.Raw.Safe+ Numeric.MPFR.Raw.Unsafe++ other-modules:+ Numeric.Rounded.Internal+ Numeric.Rounded.Rounding+ Numeric.Rounded.Precision++ build-depends:+ base >= 4.8 && < 4.13,+ ghc-prim >= 0.4 && < 0.6,+ reflection >= 2.1.2 && < 2.2,+ singletons >= 2.1 && < 2.6,+ hgmp >= 0.1.1 && < 0.2,+ long-double >= 0.1 && < 0.2++ ghc-options: -Wall -fspec-constr -funbox-strict-fields -fno-warn-unticked-promoted-constructors+ build-tools: hsc2hs++ extra-libraries: mpfr gmp++ hs-source-dirs: src+ c-sources: cbits/wrappers.c++ default-language: Haskell2010+ other-extensions:+ CPP+ DataKinds+ DeriveDataTypeable+ DeriveGeneric+ EmptyDataDecls+ FlexibleContexts+ FlexibleInstances+ ForeignFunctionInterface+ GeneralizedNewtypeDeriving+ KindSignatures+ MagicHash+ MultiParamTypeClasses+ PatternSynonyms+ PolyKinds+ Rank2Types+ RankNTypes+ RoleAnnotations+ ScopedTypeVariables+ TypeFamilies+ TypeSynonymInstances+ UnboxedTuples+ UndecidableInstances++test-suite rounded-test+ type: exitcode-stdio-1.0+ main-is: test.hs+ build-depends:+ base,+ long-double,+ rounded+ default-language: Haskell2010+ other-extensions:+ DataKinds+ ScopedTypeVariables
+ src/Numeric/MPFR/Raw.hs view
@@ -0,0 +1,25 @@+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.MPFR.Raw+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2017 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+-- This module contains FFI imports. In a future release they may switch+-- at runtime between safe and unsafe calls depending on cost estimates, but+-- for now thie module simply re-exports 'Numeric.MPFR.Raw.Safe'.+--+-- Note: beware issues with the GHC threaded runtime if you need the MPFR+-- status flags (which use OS thread local storage aka TLS). Even if you use+-- 'Control.Concurrent.forkOS' to create a bound thread whose FFI will all be+-- done by the same OS thread, sparks created with 'Control.Parallel.par'+-- may run on a different OS thread and lead to unpredictable behaviour.++module Numeric.MPFR.Raw+ ( module Numeric.MPFR.Raw.Safe+ ) where++import Numeric.MPFR.Raw.Safe
+ src/Numeric/MPFR/Raw/Safe.hs view
@@ -0,0 +1,116 @@+{-# LANGUAGE ForeignFunctionInterface #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.MPFR.Raw.Unsafe+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2017 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+-- This module contains FFI imports as safe ccalls.+----------------------------------------------------------------------------+module Numeric.MPFR.Raw.Safe where++import Foreign.C (CInt(..), CIntMax(..), CSize(..), CChar(..), CLong(..))++import GHC.Exts (Ptr)++import Numeric.LongDouble (LongDouble)+import Numeric.GMP.Types (MPZ, MPQ)++import Numeric.MPFR.Types++foreign import ccall safe "mpfr_init2" mpfr_init2 :: Ptr MPFR -> MPFRPrec -> IO ()+foreign import ccall safe "mpfr_clear" mpfr_clear :: Ptr MPFR -> IO ()++foreign import ccall safe "mpfr_sgn" mpfr_sgn :: Ptr MPFR -> IO CInt+foreign import ccall safe "mpfr_get_d" mpfr_get_d :: Ptr MPFR -> MPFRRnd -> IO Double+foreign import ccall safe "wrapped_mpfr_get_z" wrapped_mpfr_get_z :: Ptr MPZ -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt++foreign import ccall safe "mpfr_get_str" mpfr_get_str :: Ptr CChar -> Ptr MPFRExp -> Int -> CSize -> Ptr MPFR -> MPFRRnd -> IO (Ptr CChar)+foreign import ccall safe "mpfr_free_str" mpfr_free_str :: Ptr CChar -> IO ()++foreign import ccall safe "mpfr_set_z" mpfr_set_z :: Ptr MPFR -> Ptr MPZ -> MPFRRnd -> IO CInt+foreign import ccall safe "__gmpfr_set_sj" mpfr_set_sj :: Ptr MPFR -> CIntMax -> MPFRRnd -> IO CInt+foreign import ccall safe "mpfr_set_q" mpfr_set_q :: Ptr MPFR -> Ptr MPQ -> MPFRRnd -> IO CInt+foreign import ccall safe "mpfr_set_d" mpfr_set_d :: Ptr MPFR -> Double -> MPFRRnd -> IO CInt++type Test = Ptr MPFR -> IO CInt++foreign import ccall safe "mpfr_nan_p" mpfr_nan_p :: Test+foreign import ccall safe "mpfr_inf_p" mpfr_inf_p :: Test+foreign import ccall safe "mpfr_zero_p" mpfr_zero_p :: Test+foreign import ccall safe "mpfr_signbit" mpfr_signbit :: Test++foreign import ccall safe "wrapped_mpfr_get_z_2exp" wrapped_mpfr_get_z_2exp :: Ptr MPZ -> Ptr MPFR -> Ptr CInt -> IO MPFRExp+foreign import ccall safe "mpfr_set_z_2exp" mpfr_set_z_2exp :: Ptr MPFR -> Ptr MPZ -> MPFRExp -> MPFRRnd -> IO CInt++type Constant = Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall safe "mpfr_const_pi" mpfr_const_pi :: Constant+foreign import ccall safe "mpfr_const_log2" mpfr_const_log2 :: Constant+foreign import ccall safe "mpfr_const_euler" mpfr_const_euler :: Constant+foreign import ccall safe "mpfr_const_catalan" mpfr_const_catalan :: Constant++type Unary = Ptr MPFR -> Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall safe "mpfr_set" mpfr_set :: Unary++foreign import ccall safe "mpfr_abs" mpfr_abs :: Unary+foreign import ccall safe "mpfr_neg" mpfr_neg :: Unary+foreign import ccall safe "mpfr_log" mpfr_log :: Unary+foreign import ccall safe "mpfr_exp" mpfr_exp :: Unary+foreign import ccall safe "mpfr_sqrt" mpfr_sqrt :: Unary+foreign import ccall safe "mpfr_sin" mpfr_sin :: Unary+foreign import ccall safe "mpfr_cos" mpfr_cos :: Unary+foreign import ccall safe "mpfr_tan" mpfr_tan :: Unary+foreign import ccall safe "mpfr_asin" mpfr_asin :: Unary+foreign import ccall safe "mpfr_acos" mpfr_acos :: Unary+foreign import ccall safe "mpfr_atan" mpfr_atan :: Unary+foreign import ccall safe "mpfr_sinh" mpfr_sinh :: Unary+foreign import ccall safe "mpfr_cosh" mpfr_cosh :: Unary+foreign import ccall safe "mpfr_tanh" mpfr_tanh :: Unary+foreign import ccall safe "mpfr_asinh" mpfr_asinh :: Unary+foreign import ccall safe "mpfr_acosh" mpfr_acosh :: Unary+foreign import ccall safe "mpfr_atanh" mpfr_atanh :: Unary+foreign import ccall safe "mpfr_log1p" mpfr_log1p :: Unary+foreign import ccall safe "mpfr_expm1" mpfr_expm1 :: Unary+foreign import ccall safe "mpfr_rint" mpfr_rint :: Unary++type Unary' = Ptr MPFR -> Ptr MPFR -> IO CInt++foreign import ccall safe "mpfr_trunc" mpfr_trunc :: Unary'+foreign import ccall safe "mpfr_ceil" mpfr_ceil :: Unary'+foreign import ccall safe "mpfr_floor" mpfr_floor :: Unary'++type Binary = Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall safe "mpfr_add" mpfr_add :: Binary+foreign import ccall safe "mpfr_sub" mpfr_sub :: Binary+foreign import ccall safe "mpfr_mul" mpfr_mul :: Binary+foreign import ccall safe "mpfr_div" mpfr_div :: Binary+foreign import ccall safe "mpfr_min" mpfr_min :: Binary+foreign import ccall safe "mpfr_max" mpfr_max :: Binary+foreign import ccall safe "mpfr_atan2" mpfr_atan2 :: Binary++foreign import ccall safe "mpfr_modf" mpfr_modf :: Binary++type Comparison = Ptr MPFR -> Ptr MPFR -> IO CInt++foreign import ccall safe "mpfr_cmp" mpfr_cmp :: Comparison+foreign import ccall safe "mpfr_equal_p" mpfr_equal_p :: Comparison+foreign import ccall safe "mpfr_lessgreater_p" mpfr_lessgreater_p :: Comparison+foreign import ccall safe "mpfr_less_p" mpfr_less_p :: Comparison+foreign import ccall safe "mpfr_greater_p" mpfr_greater_p :: Comparison+foreign import ccall safe "mpfr_lessequal_p" mpfr_lessequal_p :: Comparison+foreign import ccall safe "mpfr_greaterequal_p" mpfr_greaterequal_p :: Comparison++foreign import ccall safe "mpfr_nextabove" mpfr_nextabove :: Ptr MPFR -> IO ()+foreign import ccall safe "mpfr_nextbelow" mpfr_nextbelow :: Ptr MPFR -> IO ()++foreign import ccall safe "wrapped_mpfr_get_ld" wrapped_mpfr_get_ld :: Ptr LongDouble -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt+foreign import ccall safe "wrapped_mpfr_get_ld_2exp" wrapped_mpfr_get_ld_2exp :: Ptr LongDouble -> Ptr CLong -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt+foreign import ccall safe "wrapped_mpfr_set_ld" wrapped_mpfr_set_ld :: Ptr MPFR -> Ptr LongDouble -> MPFRRnd -> IO CInt+foreign import ccall safe "wrapped_mpfr_cmp_ld" wrapped_mpfr_cmp_ld :: Ptr MPFR -> Ptr LongDouble -> IO CInt
+ src/Numeric/MPFR/Raw/Unsafe.hs view
@@ -0,0 +1,116 @@+{-# LANGUAGE ForeignFunctionInterface #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.MPFR.Raw.Unsafe+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2017 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+-- This module contains FFI imports as unsafe ccalls.+----------------------------------------------------------------------------+module Numeric.MPFR.Raw.Unsafe where++import Foreign.C (CInt(..), CIntMax(..), CSize(..), CChar(..), CLong(..))++import GHC.Exts (Ptr)++import Numeric.LongDouble (LongDouble)+import Numeric.GMP.Types (MPZ, MPQ)++import Numeric.MPFR.Types++foreign import ccall unsafe "mpfr_init2" mpfr_init2 :: Ptr MPFR -> MPFRPrec -> IO ()+foreign import ccall unsafe "mpfr_clear" mpfr_clear :: Ptr MPFR -> IO ()++foreign import ccall unsafe "mpfr_sgn" mpfr_sgn :: Ptr MPFR -> IO CInt+foreign import ccall unsafe "mpfr_get_d" mpfr_get_d :: Ptr MPFR -> MPFRRnd -> IO Double+foreign import ccall unsafe "wrapped_mpfr_get_z" wrapped_mpfr_get_z :: Ptr MPZ -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt++foreign import ccall unsafe "mpfr_get_str" mpfr_get_str :: Ptr CChar -> Ptr MPFRExp -> Int -> CSize -> Ptr MPFR -> MPFRRnd -> IO (Ptr CChar)+foreign import ccall unsafe "mpfr_free_str" mpfr_free_str :: Ptr CChar -> IO ()++foreign import ccall unsafe "mpfr_set_z" mpfr_set_z :: Ptr MPFR -> Ptr MPZ -> MPFRRnd -> IO CInt+foreign import ccall unsafe "__gmpfr_set_sj" mpfr_set_sj :: Ptr MPFR -> CIntMax -> MPFRRnd -> IO CInt+foreign import ccall unsafe "mpfr_set_q" mpfr_set_q :: Ptr MPFR -> Ptr MPQ -> MPFRRnd -> IO CInt+foreign import ccall unsafe "mpfr_set_d" mpfr_set_d :: Ptr MPFR -> Double -> MPFRRnd -> IO CInt++type Test = Ptr MPFR -> IO CInt++foreign import ccall unsafe "mpfr_nan_p" mpfr_nan_p :: Test+foreign import ccall unsafe "mpfr_inf_p" mpfr_inf_p :: Test+foreign import ccall unsafe "mpfr_zero_p" mpfr_zero_p :: Test+foreign import ccall unsafe "mpfr_signbit" mpfr_signbit :: Test++foreign import ccall unsafe "wrapped_mpfr_get_z_2exp" wrapped_mpfr_get_z_2exp :: Ptr MPZ -> Ptr MPFR -> Ptr CInt -> IO MPFRExp+foreign import ccall unsafe "mpfr_set_z_2exp" mpfr_set_z_2exp :: Ptr MPFR -> Ptr MPZ -> MPFRExp -> MPFRRnd -> IO CInt++type Constant = Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall unsafe "mpfr_const_pi" mpfr_const_pi :: Constant+foreign import ccall unsafe "mpfr_const_log2" mpfr_const_log2 :: Constant+foreign import ccall unsafe "mpfr_const_euler" mpfr_const_euler :: Constant+foreign import ccall unsafe "mpfr_const_catalan" mpfr_const_catalan :: Constant++type Unary = Ptr MPFR -> Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall unsafe "mpfr_set" mpfr_set :: Unary++foreign import ccall unsafe "mpfr_abs" mpfr_abs :: Unary+foreign import ccall unsafe "mpfr_neg" mpfr_neg :: Unary+foreign import ccall unsafe "mpfr_log" mpfr_log :: Unary+foreign import ccall unsafe "mpfr_exp" mpfr_exp :: Unary+foreign import ccall unsafe "mpfr_sqrt" mpfr_sqrt :: Unary+foreign import ccall unsafe "mpfr_sin" mpfr_sin :: Unary+foreign import ccall unsafe "mpfr_cos" mpfr_cos :: Unary+foreign import ccall unsafe "mpfr_tan" mpfr_tan :: Unary+foreign import ccall unsafe "mpfr_asin" mpfr_asin :: Unary+foreign import ccall unsafe "mpfr_acos" mpfr_acos :: Unary+foreign import ccall unsafe "mpfr_atan" mpfr_atan :: Unary+foreign import ccall unsafe "mpfr_sinh" mpfr_sinh :: Unary+foreign import ccall unsafe "mpfr_cosh" mpfr_cosh :: Unary+foreign import ccall unsafe "mpfr_tanh" mpfr_tanh :: Unary+foreign import ccall unsafe "mpfr_asinh" mpfr_asinh :: Unary+foreign import ccall unsafe "mpfr_acosh" mpfr_acosh :: Unary+foreign import ccall unsafe "mpfr_atanh" mpfr_atanh :: Unary+foreign import ccall unsafe "mpfr_log1p" mpfr_log1p :: Unary+foreign import ccall unsafe "mpfr_expm1" mpfr_expm1 :: Unary+foreign import ccall unsafe "mpfr_rint" mpfr_rint :: Unary++type Unary' = Ptr MPFR -> Ptr MPFR -> IO CInt++foreign import ccall unsafe "mpfr_trunc" mpfr_trunc :: Unary'+foreign import ccall unsafe "mpfr_ceil" mpfr_ceil :: Unary'+foreign import ccall unsafe "mpfr_floor" mpfr_floor :: Unary'++type Binary = Ptr MPFR -> Ptr MPFR -> Ptr MPFR -> MPFRRnd -> IO CInt++foreign import ccall unsafe "mpfr_add" mpfr_add :: Binary+foreign import ccall unsafe "mpfr_sub" mpfr_sub :: Binary+foreign import ccall unsafe "mpfr_mul" mpfr_mul :: Binary+foreign import ccall unsafe "mpfr_div" mpfr_div :: Binary+foreign import ccall unsafe "mpfr_min" mpfr_min :: Binary+foreign import ccall unsafe "mpfr_max" mpfr_max :: Binary+foreign import ccall unsafe "mpfr_atan2" mpfr_atan2 :: Binary++foreign import ccall unsafe "mpfr_modf" mpfr_modf :: Binary++type Comparison = Ptr MPFR -> Ptr MPFR -> IO CInt++foreign import ccall unsafe "mpfr_cmp" mpfr_cmp :: Comparison+foreign import ccall unsafe "mpfr_equal_p" mpfr_equal_p :: Comparison+foreign import ccall unsafe "mpfr_lessgreater_p" mpfr_lessgreater_p :: Comparison+foreign import ccall unsafe "mpfr_less_p" mpfr_less_p :: Comparison+foreign import ccall unsafe "mpfr_greater_p" mpfr_greater_p :: Comparison+foreign import ccall unsafe "mpfr_lessequal_p" mpfr_lessequal_p :: Comparison+foreign import ccall unsafe "mpfr_greaterequal_p" mpfr_greaterequal_p :: Comparison++foreign import ccall unsafe "mpfr_nextabove" mpfr_nextabove :: Ptr MPFR -> IO ()+foreign import ccall unsafe "mpfr_nextbelow" mpfr_nextbelow :: Ptr MPFR -> IO ()++foreign import ccall unsafe "wrapped_mpfr_get_ld" wrapped_mpfr_get_ld :: Ptr LongDouble -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt+foreign import ccall unsafe "wrapped_mpfr_get_ld_2exp" wrapped_mpfr_get_ld_2exp :: Ptr LongDouble -> Ptr CLong -> Ptr MPFR -> MPFRRnd -> Ptr CInt -> IO CInt+foreign import ccall unsafe "wrapped_mpfr_set_ld" wrapped_mpfr_set_ld :: Ptr MPFR -> Ptr LongDouble -> MPFRRnd -> IO CInt+foreign import ccall unsafe "wrapped_mpfr_cmp_ld" wrapped_mpfr_cmp_ld :: Ptr MPFR -> Ptr LongDouble -> IO CInt
+ src/Numeric/MPFR/Types.hsc view
@@ -0,0 +1,127 @@+#include <ghc-gmp.h>+#include <mpfr.h>++#if __GLASGOW_HASKELL__ < 800+#let alignment t = "%lu", (unsigned long)offsetof(struct {char x__; t (y__); }, y__)+#endif++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE PatternSynonyms #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.MPFR.Types+-- Copyright : (C) 2015-2018 Claude Heiland-Allen+-- License : LGPL+-- Stability : experimental+-- Portability : non-portable+--+-- MPFR types.+----------------------------------------------------------------------------+module Numeric.MPFR.Types where++import Data.Data+import Data.Typeable+import Data.Bits+import Data.Ix+import Data.Int+import Data.Word++import Foreign (Storable(..), Ptr)++import Numeric.GMP.Types (MPLimb(..))++-- | Version numbers.+pattern MPFR_VERSION_MAJOR = #const MPFR_VERSION_MAJOR+pattern MPFR_VERSION_MINOR = #const MPFR_VERSION_MINOR+pattern MPFR_VERSION_PATCHLEVEL = #const MPFR_VERSION_PATCHLEVEL+pattern MPFR_VERSION = #const MPFR_VERSION++-- | @mpfr_t@+data MPFR = MPFR+ { mpfrPrec :: !MPFRPrec+ , mpfrSign :: !MPFRSign+ , mpfrExp :: !MPFRExp+ , mpfrD :: !(Ptr MPLimb)+ }++instance Storable MPFR where+ sizeOf _ = (#size __mpfr_struct)+ alignment _ = (#alignment __mpfr_struct)+ peek ptr = do+ prec <- (#peek __mpfr_struct, _mpfr_prec) ptr+ sign <- (#peek __mpfr_struct, _mpfr_sign) ptr+ exp' <- (#peek __mpfr_struct, _mpfr_exp) ptr+ d <- (#peek __mpfr_struct, _mpfr_d) ptr+ return (MPFR{ mpfrPrec = prec, mpfrSign = sign, mpfrExp = exp', mpfrD = d })+ poke ptr (MPFR{ mpfrPrec = prec, mpfrSign = sign, mpfrExp = exp', mpfrD = d }) = do+ (#poke __mpfr_struct, _mpfr_prec) ptr prec+ (#poke __mpfr_struct, _mpfr_sign) ptr sign+ (#poke __mpfr_struct, _mpfr_exp) ptr exp'+ (#poke __mpfr_struct, _mpfr_d) ptr d++-- | @mpfr_int@+newtype MPFRInt = MPFRInt (#type mpfr_int)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+-- | @mpfr_uint@+newtype MPFRUInt = MPFRUInt (#type mpfr_uint)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)++-- | @mpfr_long@+newtype MPFRLong = MPFRLong (#type mpfr_long)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+-- | @mpfr_ulong@+newtype MPFRULong = MPFRULong (#type mpfr_ulong)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)++-- | @mpfr_size_t@+newtype MPFRSize = MPFRSize (#type mpfr_size_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)++-- | @mpfr_prec_t@+newtype MPFRPrec = MPFRPrec (#type mpfr_prec_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+-- | @mpfr_uprec_t@+newtype MPFRUPrec = MPFRUPrec (#type mpfr_uprec_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+pattern MPFR_PREC_MIN = (#const MPFR_PREC_MIN)+pattern MPFR_PREC_MAX = (#const MPFR_PREC_MAX)++-- | @mpfr_sign_t@+newtype MPFRSign = MPFRSign (#type mpfr_sign_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)++-- | @mpfr_exp_t@+newtype MPFRExp = MPFRExp (#type mpfr_exp_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+-- | @mpfr_uexp_t@+newtype MPFRUExp = MPFRUExp (#type mpfr_uexp_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+pattern MPFR_EMAX_DEFAULT = (#const MPFR_EMAX_DEFAULT)+pattern MPFR_EMIN_DEFAULT = (#const MPFR_EMIN_DEFAULT)++-- | @mpfr_rnd_t@+newtype MPFRRnd = MPFRRnd (#type mpfr_rnd_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+-- | round to nearest, with ties to even+pattern MPFR_RNDN = (#const MPFR_RNDN)+-- | round toward zero+pattern MPFR_RNDZ = (#const MPFR_RNDZ)+-- | round toward +Inf+pattern MPFR_RNDU = (#const MPFR_RNDU)+-- | round toward -Inf+pattern MPFR_RNDD = (#const MPFR_RNDD)+-- | round away from zero+pattern MPFR_RNDA = (#const MPFR_RNDA)+-- | faithful rounding (not implemented yet)+pattern MPFR_RNDF = (#const MPFR_RNDF)+-- | round to nearest, with ties away from zero (mpfr_round) (do not use)+pattern MPFR_RNDNA = (#const MPFR_RNDNA)++-- | @mpfr_kind_t@+newtype MPFRKind = MPFRKind (#type mpfr_kind_t)+ deriving (Eq, Ord, Read, Show, Enum, Bounded, Num, Integral, Real, Ix, Bits, FiniteBits, Data, Typeable, Storable)+pattern MPFR_NAN_KIND = (#const MPFR_NAN_KIND)+pattern MPFR_INF_KIND = (#const MPFR_INF_KIND)+pattern MPFR_ZERO_KIND = (#const MPFR_ZERO_KIND)+pattern MPFR_REGULAR_KIND = (#const MPFR_REGULAR_KIND)
+ src/Numeric/Rounded.hs view
@@ -0,0 +1,95 @@+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Rounded+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2018 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+----------------------------------------------------------------------------+module Numeric.Rounded+ (+ -- * Floating point numbers with a specified rounding mode and precision+ Rounded()+ , fromInt+ , fromDouble+ , fromLongDouble+ , toDouble+ , toLongDouble+ , toInteger'+ , precRound+ -- * Precision+ , Precision(precision)+ , Bytes+ , reifyPrecision+ -- * Rounding+ , Rounding(rounding)+ , RoundingMode(..)+ , reifyRounding+ -- * Useful Constants+ , kLog2+ , kEuler+ , kCatalan+ -- * Combinators that are oblivious to precision+ , (.+.)+ , (.-.)+ , (.*.)+ , abs'+ , negate'+ , decodeFloat'+ , succUlp+ , predUlp+ -- * Mixed-precision operations+ , (!+!)+ , (!-!)+ , (!*!)+ , (!/!)+ , abs_+ , negate_+ , compare_+ , (!==!)+ , (!/=!)+ , (!>=!)+ , (!<=!)+ , (!>!)+ , (!<!)+ , min_+ , max_+ , sqrt_+ , exp_+ , expm1_+ , log_+ , log1p_+ , sin_+ , cos_+ , tan_+ , asin_+ , acos_+ , atan_+ , atan2_+ , sinh_+ , cosh_+ , tanh_+ , asinh_+ , acosh_+ , atanh_+ , truncate_+ , round_+ , ceiling_+ , floor_+ -- * Foreign Function Interface+ , withInRounded+ , withInOutRounded+ , withInOutRounded_+ , withOutRounded+ , withOutRounded_+ , unsafeWithOutRounded+ , unsafeWithOutRounded_+ , peekRounded+ ) where++import Numeric.Rounded.Internal+import Numeric.Rounded.Precision+import Numeric.Rounded.Rounding
+ src/Numeric/Rounded/Internal.hs view
@@ -0,0 +1,610 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE RoleAnnotations #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples #-}+{-# OPTIONS_HADDOCK not-home #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Rounded.Internal+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2018 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+----------------------------------------------------------------------------+module Numeric.Rounded.Internal where++import Control.Exception (bracket, bracket_, throwIO, ArithException(Overflow))+import Data.Bits (shiftL, testBit)+import Data.Coerce (coerce)+import Data.Int (Int32)+import Data.Proxy (Proxy(..))+import Data.Ratio ((%))++import Foreign (with, alloca, allocaBytes, peek, sizeOf, nullPtr)+import Foreign.C (CInt(..), CIntMax(..))+import Foreign.C.String (peekCString)+import Numeric.LongDouble (LongDouble)++import System.IO.Unsafe (unsafePerformIO)++import GHC.Prim+ ( ByteArray#+ , sizeofByteArray#+ , copyByteArrayToAddr#+ , newByteArray#+ , copyAddrToByteArray#+ , unsafeFreezeByteArray#+ )+import GHC.Types (IO(..))+import GHC.Exts (Ptr(..), Int(..))++#if MIN_VERSION_base(4,9,0)+import Numeric (Floating(..))+#endif+import Numeric (readSigned, readFloat)++import Numeric.GMP.Utils (withInInteger, withOutInteger, withOutInteger_, withInRational)+import Numeric.GMP.Types (MPLimb)++import Numeric.MPFR.Types+import Numeric.MPFR.Raw++import Numeric.Rounded.Precision+import Numeric.Rounded.Rounding++type role Rounded phantom nominal++-- | A properly rounded floating-point number with a given rounding mode and precision.+--+-- You can 'Data.Coerce.coerce' to change rounding modes, but not precision.+data Rounded (r :: RoundingMode) p = Rounded+ { roundedPrec :: !MPFRPrec+ , roundedSign :: !MPFRSign+ , roundedExp :: !MPFRExp+ , roundedLimbs :: !ByteArray#+ }++-- | Round to 'Double' with the given rounding mode.+toDouble :: (Rounding r, Precision p) => Rounded r p -> Double+toDouble x = unsafePerformIO $ in_ x $ \xfr -> mpfr_get_d xfr (rnd x)+-- this syntax is strange, but it seems to be the way it works...+{-# RULES "realToFrac/toDouble" forall (x :: (Rounding r, Precision p) => Rounded r p) . realToFrac x = toDouble x #-}++-- | Round to 'LongDouble' with the given rounding mode.+toLongDouble :: (Rounding r, Precision p) => Rounded r p -> LongDouble+toLongDouble x = unsafePerformIO $ in_ x $ \xfr -> with 0 $ \yfr -> with 0 $ \ffr -> wrapped_mpfr_get_ld yfr xfr (rnd x) ffr >> peek yfr+-- this syntax is strange, but it seems to be the way it works...+{-# RULES "realToFrac/toLongDouble" forall (x :: (Rounding r, Precision p) => Rounded r p) . realToFrac x = toLongDouble x #-}++-- | Round to a different precision with the given rounding mode.+precRound :: (Rounding r, Precision p1, Precision p2) => Rounded r p1 -> Rounded r p2+precRound x = unsafePerformIO $ do+ (Just y, _) <- in_ x $ \xfr -> out_ $ \yfr ->+ mpfr_set yfr xfr (rnd x)+ return y+-- TODO figure out correct syntax (if even possible) to allow RULE+-- {-# RULES "realToFrac/precRound" realToFrac = precRound #-}++toString :: (Rounding r, Precision p) => Rounded r p -> String+-- FIXME: what do about unsightly 0.1 -> 0.1000...0002 or 9.999...9995e-2 issues+toString x = unsafePerformIO $ do+ (s, e) <- in_ x $ \xfr -> with 0 $ \eptr -> do+ s <- bracket (mpfr_get_str nullPtr eptr 10 0 xfr (fromIntegral (fromEnum TowardNearest))) mpfr_free_str peekCString+ e <- peek eptr+ return (s, fromIntegral e)+ return $ case () of+ _ | isNaN x -> "NaN"+ | isInfinite x && sgn' == GT -> "Infinity"+ | isInfinite x -> "-Infinity"+ | isNegativeZero x -> "-0.0"+ | sgn' == EQ -> "0.0"+ | e < 0 ||+ e >= threshold -> sign ++ take 1 digits ++ "." +++ dropTrailingZeroes (take (n - 1) (drop 1 digits0)) +++ "e" ++ show (e - 1)+ | e == 0 -> sign ++ "0." +++ dropTrailingZeroes digits+ | e < threshold -> sign ++ take e digits0 ++ "." +++ dropTrailingZeroes (take (n - e) (drop e digits0))+ where+ sgn' = sgn x+ sign = case sgn' of+ GT -> ""+ EQ -> ""+ LT -> "-"+ threshold = 8+ n = length digits+ digits = case take 1 s of+ "-" -> drop 1 s+ _ -> s+ digits0 = digits ++ repeat '0'+ dropTrailingZeroes a = case dropWhile ('0' ==) (reverse a) of+ "" -> "0"+ b -> reverse b++instance (Rounding r, Precision p) => Show (Rounded r p) where+ showsPrec p x = showParen (p >= 7 && take 1 s == "-") (s ++) -- FIXME: precedence issues?+ where s = toString x++instance (Rounding r, Precision p) => Read (Rounded r p) where+ -- apparently this handles parens without any extra fuss+ readsPrec _ = readSigned readFloat -- FIXME: precedence issues?++unary+ :: (Rounding r, Precision p1, Precision p2)+ => Unary -> Rounded r p1 -> Rounded r p2+unary f a = unsafePerformIO $ do+ (Just c, _) <- in_ a $ \afr ->+ out_ $ \cfr ->+ f cfr afr (rnd a)+ return c++unary' :: Rounding r => Unary -> Rounded r p -> Rounded r p+unary' f a = unsafePerformIO $ do+ (Just c, _) <- in_ a $ \afr ->+ out_' (roundedPrec a) $ \cfr ->+ f cfr afr (rnd a)+ return c++unary'' :: Unary -> Rounded r p -> Rounded r p+unary'' f a = unsafePerformIO $ do+ (Just c, _) <- in_ a $ \afr ->+ out_' (roundedPrec a) $ \cfr ->+ f cfr afr (fromIntegral (fromEnum TowardNearest))+ return c++abs' :: Rounded r p -> Rounded r p+abs' = unary'' mpfr_abs++negate' :: Rounded r p -> Rounded r p+negate' = unary'' mpfr_neg++(.-.), (.+.), (.*.) :: Rounding r => Rounded r p -> Rounded r p -> Rounded r p+(.-.) = binary' mpfr_sub+(.+.) = binary' mpfr_add+(.*.) = binary' mpfr_mul++infixl 6 .+., .-.+infixl 7 .*.+++abs_, negate_, log_, exp_, sqrt_,+ sin_, cos_, tan_, asin_, acos_, atan_,+ sinh_, cosh_, tanh_, asinh_, acosh_, atanh_,+ log1p_, expm1_+ :: (Rounding r, Precision p1, Precision p2)+ => Rounded r p1 -> Rounded r p2+abs_ = unary mpfr_abs+negate_ = unary mpfr_neg+log_ = unary mpfr_log+exp_ = unary mpfr_exp+sqrt_ = unary mpfr_sqrt+sin_ = unary mpfr_sin+cos_ = unary mpfr_cos+tan_ = unary mpfr_tan+asin_ = unary mpfr_asin+acos_ = unary mpfr_acos+atan_ = unary mpfr_atan+sinh_ = unary mpfr_sinh+cosh_ = unary mpfr_cosh+tanh_ = unary mpfr_tanh+asinh_ = unary mpfr_asinh+acosh_ = unary mpfr_acosh+atanh_ = unary mpfr_atanh+log1p_ = unary mpfr_log1p+expm1_ = unary mpfr_expm1++binary+ :: (Rounding r, Precision p1, Precision p2, Precision p3)+ => Binary -> Rounded r p1 -> Rounded r p2 -> Rounded r p3+binary f a b = unsafePerformIO $ do+ (Just c, _) <- in_ a $ \afr ->+ in_ b $ \bfr ->+ out_ $ \cfr ->+ f cfr afr bfr (rnd a)+ return c++min_, max_, (!+!), (!-!), (!*!), (!/!), atan2_+ :: (Rounding r, Precision p1, Precision p2, Precision p3)+ => Rounded r p1 -> Rounded r p2 -> Rounded r p3+min_ = binary mpfr_min+max_ = binary mpfr_max+(!+!) = binary mpfr_add+(!-!) = binary mpfr_sub+(!*!) = binary mpfr_mul+(!/!) = binary mpfr_div+atan2_ = binary mpfr_atan2++infixl 6 !+!, !-!+infixl 7 !*!, !/!++binary' :: Rounding r => Binary -> Rounded r p -> Rounded r p -> Rounded r p+binary' f a b = unsafePerformIO $ do+ (Just c, _) <- in_ a $ \afr ->+ in_ b $ \bfr ->+ out_' (roundedPrec a) $ \cfr ->+ f cfr afr bfr (rnd a)+ return c++cmp' :: Comparison -> Rounded r p1 -> Rounded r p2 -> CInt+cmp' f a b = unsafePerformIO $+ in_ a $ \afr ->+ in_ b $ \bfr -> do+ f afr bfr++cmp :: Comparison -> Rounded r p1 -> Rounded r p2 -> Bool+cmp f a b = cmp' f a b /= 0++(!==!), (!/=!), (!<=!), (!>=!), (!<!), (!>!)+ :: (Precision p1, Precision p2)+ => Rounded r p1 -> Rounded r p2 -> Bool+(!==!) = cmp mpfr_equal_p+(!/=!) = cmp mpfr_lessgreater_p+(!<=!) = cmp mpfr_lessequal_p+(!>=!) = cmp mpfr_greaterequal_p+(!<!) = cmp mpfr_less_p+(!>!) = cmp mpfr_greater_p++infix 4 !==!, !/=!, !<=!, !>=!, !<!, !>!++compare_ :: (Precision p1, Precision p2) => Rounded r p1 -> Rounded r p2 -> Ordering+compare_ a b = compare (cmp' mpfr_cmp a b) 0++instance Eq (Rounded r p) where+ (==) = cmp mpfr_equal_p+ (/=) = cmp mpfr_lessgreater_p++instance Rounding r => Ord (Rounded r p) where+ compare a b = compare (cmp' mpfr_cmp a b) 0+ (<=) = cmp mpfr_lessequal_p+ (>=) = cmp mpfr_greaterequal_p+ (<) = cmp mpfr_less_p+ (>) = cmp mpfr_greater_p+ min = binary' mpfr_min+ max = binary' mpfr_max++sgn :: (Rounding r, Precision p) => Rounded r p -> Ordering+sgn x = compare (unsafePerformIO $ in_ x mpfr_sgn) 0++instance (Rounding r, Precision p) => Num (Rounded r p) where+ (+) = (.+.)+ (-) = (.-.)+ (*) = (.*.)+ negate = negate'+ fromInteger j = r where+ r = unsafePerformIO $ do+ if toInteger (minBound :: CIntMax) <= j && j <= toInteger (maxBound :: CIntMax)+ then do+ (Just x, _) <- out_ $ \jfr -> mpfr_set_sj jfr (fromInteger j :: CIntMax) (rnd r)+ return x+ else do+ (Just x, _) <- withInInteger j $ \jz -> out_ $ \jfr -> mpfr_set_z jfr jz (rnd r)+ return x+ abs = abs'+ signum x = case sgn x of+ LT -> -1+ EQ -> 0+ GT -> 1++instance (Rounding r, Precision p) => Fractional (Rounded r p) where+ fromRational q = r where -- TODO small integer optimisation+ r = unsafePerformIO $ do+ (Just x, _) <- withInRational q $ \qq -> out_ $ \qfr -> mpfr_set_q qfr qq (rnd r)+ return x+ (/) = (!/!)++-- | Construct a properly rounded floating point number from an 'Int'.+fromInt :: (Rounding r, Precision p) => Int -> Rounded r p+fromInt i = r+ where+ r = unsafePerformIO $ do+ (Just x, _) <- out_ $ \xfr -> mpfr_set_sj xfr (fromIntegral i) (rnd r)+ return x+-- TODO figure out correct syntax (if even possible) to allow RULE+-- {-# RULES "fromIntegral/fromInt" fromIntegral = fromInt #-}++-- | Construct a rounded floating point number directly from a 'Double'.+fromDouble :: (Rounding r, Precision p) => Double -> Rounded r p+fromDouble d = r+ where+ r = unsafePerformIO $ do+ (Just x, _) <- out_ $ \xfr -> mpfr_set_d xfr d (rnd r)+ return x+-- TODO figure out correct syntax (if even possible) to allow RULE+-- {-# RULES "realToFrac/fromDouble" realToFrac = fromDouble #-}++-- | Construct a rounded floating point number directly from a 'LongDouble'.+fromLongDouble :: (Rounding r, Precision p) => LongDouble -> Rounded r p+fromLongDouble d = r+ where+ r = unsafePerformIO $ do+ (Just x, _) <- out_ $ \xfr -> with d $ \dp -> wrapped_mpfr_set_ld xfr dp (rnd r)+ return x+-- TODO figure out correct syntax (if even possible) to allow RULE+-- {-# RULES "realToFrac/fromLongDouble" realToFrac = fromLongDouble #-}+++inplace :: (Ptr MPFR -> IO ()) -> Rounded r p -> Rounded r p+inplace f y = unsafePerformIO $ do+ (Just x, _) <- out_' (roundedPrec y) $ \xfr -> in_ y $ \yfr -> do+ _ <- mpfr_set xfr yfr (fromIntegral (fromEnum TowardNearest))+ f xfr+ return x++succUlp, predUlp :: Rounded r p -> Rounded r p+succUlp = inplace mpfr_nextabove+predUlp = inplace mpfr_nextbelow++constant :: (Rounding r, Precision p) => Constant -> Rounded r p+constant k = r where+ r = unsafePerformIO $ do+ (Just x, _) <- out_ $ \xfr -> k xfr (rnd r)+ return x+++instance (Rounding r, Precision p) => Floating (Rounded r p) where+ pi = kPi+ exp = exp_+ sqrt = sqrt_+ log = log_+ sin = sin_+ tan = tan_+ cos = cos_+ asin = asin_+ atan = atan_+ acos = acos_+ sinh = sinh_+ tanh = tanh_+ cosh = cosh_+ asinh = asinh_+ atanh = atanh_+ acosh = acosh_+#if MIN_VERSION_base(4,9,0)+ log1p = log1p_+ expm1 = expm1_+#endif++toRational' :: Precision p => Rounded r p -> Rational+toRational' r+ | e > 0 = fromIntegral (s `shiftL` e)+ | otherwise = s % (1 `shiftL` negate e)+ where (s, e) = decodeFloat' r++instance (Rounding r, Precision p) => Real (Rounded r p) where+ toRational = toRational'++modf :: (Rounding r, Precision p) => Rounded r p -> (Rounded r p, Rounded r p)+modf x = unsafePerformIO $ do+ (Just y, (Just z, _)) <- in_ x $ \xfr ->+ out_ $ \yfr ->+ out_ $ \zfr ->+ mpfr_modf yfr zfr xfr (rnd x)+ return (y, z)++-- | Round to 'Integer' using the specified rounding mode. Throws 'Overflow' if+-- the result cannot be represented (for example, infinities or NaN).+toInteger' :: (Rounding r, Precision p) => Rounded r p -> Integer+toInteger' x = unsafePerformIO $+ withOutInteger_ $ \yz ->+ in_ x $ \xfr ->+ with 0 $ \flagsptr -> do+ e <- wrapped_mpfr_get_z yz xfr (rnd x) flagsptr+ flags <- peek flagsptr+ case testBit flags erangeBit of+ False -> return e+ True -> throwIO Overflow++instance (Rounding r, Precision p) => RealFrac (Rounded r p) where+ properFraction r = (fromInteger (toInteger' i), f) where (i, f) = modf r+ truncate = roundFunc truncate_+ round = roundFunc round_+ ceiling = roundFunc ceiling_+ floor = roundFunc floor_++roundFunc :: (Integral i, Precision p) => (Rounded TowardNearest p -> Rounded TowardNearest p) -> Rounded r p -> i+roundFunc f = fromInteger . toInteger' . f . coerce++unary_ :: (Precision p1, Precision p2) => (Ptr MPFR -> Ptr MPFR -> IO CInt) -> Rounded r p1 -> Rounded r p2+unary_ f x = unsafePerformIO $ do+ Just y <- withInRounded x $ \xp -> withOutRounded_ $ \yp -> f yp xp+ return y++truncate_, ceiling_, floor_, round_ :: (Precision p1, Precision p2) => Rounded r p1 -> Rounded r p2+truncate_ = unary_ mpfr_trunc+ceiling_ = unary_ mpfr_ceil+floor_ = unary_ mpfr_floor+round_ = unary_ (\yp xp -> mpfr_rint yp xp (fromIntegral (fromEnum TowardNearest)))++tst :: (Precision p) => Test -> Rounded r p -> Bool+tst f x = unsafePerformIO $ in_ x $ \xfr -> do+ t <- f xfr+ return (t /= 0)++decodeFloat' :: Rounded r p -> (Integer, Int)+decodeFloat' x = case (unsafePerformIO $ do+ in_ x $ \xfr -> withOutInteger $ \xz -> with 0 $ \flagsptr -> do+ e <- wrapped_mpfr_get_z_2exp xz xfr flagsptr+ flags <- peek flagsptr+ case testBit flags erangeBit of+ False -> return (fromIntegral e)+ True -> throwIO Overflow) of+ (0, _) -> (0, 0) -- mpfr_get_z_2exp returns emin instead of 0 for exponent+ me -> me++encodeFloat' :: (Rounding r, Precision p) => Integer -> Int -> Rounded r p+encodeFloat' j e = r where+ r = unsafePerformIO $ do+ (Just x, _) <- withInInteger j $ \jz -> out_ $ \xfr -> mpfr_set_z_2exp xfr jz (fromIntegral e) (rnd r)+ return x++instance (Rounding r, Precision p) => RealFloat (Rounded r p) where+ floatRadix _ = 2+ floatDigits = self where+ self _ = p+ p = precision (0 `asTypeIn` self)+ asTypeIn :: a -> (a -> b) -> a+ asTypeIn = const++ -- FIXME: this should do for now, but the real ones can change...+ -- FIXME: do these need to be offset to match Haskell conventions?+ floatRange _ = (MPFR_EMIN_DEFAULT, MPFR_EMAX_DEFAULT)++ decodeFloat = decodeFloat'+ encodeFloat = encodeFloat'+ isNaN = tst mpfr_nan_p+ isInfinite = tst mpfr_inf_p+ isDenormalized _ = False+ isNegativeZero r = tst mpfr_zero_p r && tst mpfr_signbit r+ isIEEE _ = True -- is this a lie? it mostly behaves like an IEEE float, despite being much bigger+ atan2 = atan2_++kPi :: (Rounding r, Precision p) => Rounded r p+kPi = constant mpfr_const_pi++-- | Natural logarithm of 2+kLog2 :: (Rounding r, Precision p) => Rounded r p+kLog2 = constant mpfr_const_log2++-- | 0.577...+kEuler :: (Rounding r, Precision p) => Rounded r p+kEuler = constant mpfr_const_euler++-- | 0.915...+kCatalan :: (Rounding r, Precision p) => Rounded r p+kCatalan = constant mpfr_const_catalan+++in_' :: Rounded r p -> (MPFR -> IO a) -> IO a+in_' (Rounded p s e l) f = withByteArray l $ \ptr _bytes -> f MPFR+ { mpfrPrec = p+ , mpfrSign = s+ , mpfrExp = e+ , mpfrD = ptr+ }++in_ :: Rounded r p -> (Ptr MPFR -> IO a) -> IO a+in_ x f = in_' x $ \y -> with y f+++out_' :: MPFRPrec -> (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p), a)+out_' p f = allocaBytes (precBytes p) $ \d -> with+ MPFR{ mpfrPrec = p, mpfrSign = 0, mpfrExp = 0, mpfrD = d } $ \ptr -> do+ a <- f ptr+ MPFR{ mpfrPrec = p', mpfrSign = s', mpfrExp = e', mpfrD = d' } <- peek ptr+ if p /= p' then return (Nothing, a) else+ asByteArray d' (precBytes p') $ \l' -> return (Just (Rounded p' s' e' l'), a)++out_ :: Precision p => (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p), a)+out_ f = r where+ r = out_' prec f+ prec = fromIntegral (precision (t r))+ t :: IO (Maybe t, a) -> t+ t _ = undefined+++-- | Use a value as a /constant/ @mpfr_t@ (attempts to modify it may explode,+-- changing the precision will explode).+withInRounded :: Rounded r p -> (Ptr MPFR -> IO a) -> IO a+withInRounded = in_++-- | Allocates and initializes a new @mpfr_t@, if the precision matches after+-- the action then it is peeked and returned. Otherwise you get 'Nothing'.+withOutRounded :: Precision p => (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p), a)+withOutRounded f = r where+ r = alloca $ \ptr -> bracket_ (mpfr_init2 ptr prec) (mpfr_clear ptr) $ do+ a <- f ptr+ MPFR{ mpfrPrec = prec', mpfrSign = s, mpfrExp = e, mpfrD = d } <- peek ptr+ if prec /= prec'+ then return (Nothing, a)+ else asByteArray d (precBytes prec) $ \l ->+ return (Just (Rounded prec s e l), a)+ prec = fromIntegral (precision (t r))+ t :: IO (Maybe b, a) -> b+ t _ = undefined++-- | Allocates and initializes a new @mpfr_t@, if the precision matches after+-- the action then it is peeked and returned. Otherwise you get 'Nothing'.+-- The result of the action is ignored.+withOutRounded_ :: Precision p => (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p))+withOutRounded_ = fmap fst . withOutRounded++-- | Like 'withOutRounded' but with the limbs allocated by GHC, which should be+-- slightly faster. However, it will crash if MPFR tries to reallocate the+-- limbs, so the action must not try to change the precision or clear it, etc.+unsafeWithOutRounded :: Precision p => (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p), a)+unsafeWithOutRounded = out_++-- | Like 'withOutRounded_' but with the limbs allocated by GHC, which should be+-- slightly faster. However, it will crash if MPFR tries to reallocate the+-- limbs, so the action must not try to change the precision or clear it, etc.+unsafeWithOutRounded_ :: Precision p => (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p))+unsafeWithOutRounded_ = fmap fst . out_++-- | Allocates and initializes a new @mpfr_t@ to the value. If the precision matches after+-- the action then it is peeked and returned. Otherwise you get 'Nothing'.+withInOutRounded :: Precision p => Rounded r p -> (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p), a)+-- FIXME: optimize to reduce copying+withInOutRounded i f =+ withOutRounded $ \ofr ->+ in_ i $ \ifr -> do+ _ <- mpfr_set ofr ifr (fromIntegral (fromEnum TowardNearest))+ f ofr++-- | Allocates and initializes a new @mpfr_t@ to the value. If the precision matches after+-- the action then it is peeked and returned. Otherwise you get 'Nothing'. The result+-- ot the action is ignored.+withInOutRounded_ :: Precision p => Rounded r p -> (Ptr MPFR -> IO a) -> IO (Maybe (Rounded r p))+withInOutRounded_ x = fmap fst . withInOutRounded x++-- | Peek an @mpfr_t@ at its actual precision, reified.+peekRounded :: Rounding r => Ptr MPFR -> (forall (p :: *) . Precision p => Rounded r p -> IO a) -> IO a+peekRounded ptr f = do+ MPFR{ mpfrPrec = p', mpfrSign = s', mpfrExp = e', mpfrD = d' } <- peek ptr+ asByteArray d' (precBytes p') $ \l' -> reifyPrecision (fromIntegral p') (wrap f (Rounded p' s' e' l'))+ where+ wrap :: forall (p :: *) (r :: RoundingMode) (a :: *) . (Rounding r, Precision p) => (forall (q :: *) . Precision q => Rounded r q -> IO a) -> Rounded r p -> Proxy p -> IO a+ wrap g r = \_proxy -> g r+++-- "The number of limbs in use is controlled by _mpfr_prec, namely ceil(_mpfr_prec/mp_bits_per_limb)."+-- <http://www.mpfr.org/mpfr-current/mpfr.html#Internals>+precBytes :: MPFRPrec -> Int+precBytes prec = bytesPerLimb * ((fromIntegral prec + bitsPerLimb1) `div` bitsPerLimb)+bytesPerLimb :: Int+bytesPerLimb = sizeOf (undefined :: MPLimb)+bitsPerLimb :: Int+bitsPerLimb = bytesPerLimb * 8+bitsPerLimb1 :: Int+bitsPerLimb1 = bitsPerLimb - 1+++erangeBit :: Int+erangeBit = 5 -- sync with cbits/wrappers.c++rnd :: Rounding r => Rounded r p -> MPFRRnd+rnd = fromIntegral . fromEnum . rounding . proxyRounding++proxyRounding :: Rounded r p -> Proxy r+proxyRounding _ = Proxy++withByteArray :: ByteArray# -> (Ptr a -> Int -> IO r) -> IO r+withByteArray ba# f = do+ let bytes = I# (sizeofByteArray# ba#)+ allocaBytes bytes $ \ptr@(Ptr addr#) -> do+ IO (\s -> (# copyByteArrayToAddr# ba# 0# addr# (sizeofByteArray# ba#) s, () #))+ f ptr bytes++asByteArray :: Ptr a -> Int -> (ByteArray# -> IO r) -> IO r+asByteArray (Ptr addr#) (I# bytes#) f = do+ IO $ \s# -> case newByteArray# bytes# s# of+ (# s'#, mba# #) ->+ case unsafeFreezeByteArray# mba# (copyAddrToByteArray# addr# mba# 0# bytes# s'#) of+ (# s''#, ba# #) -> case f ba# of IO r -> r s''#
+ src/Numeric/Rounded/Interval.hs view
@@ -0,0 +1,717 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE PolyKinds #-}+module Numeric.Rounded.Interval where++import Control.Applicative+import Numeric.Rounded+import Data.Coerce+import Data.Typeable+import GHC.Generics+import Prelude hiding (elem, notElem)++data Interval p+ = I (Rounded TowardNegInf p) (Rounded TowardInf p)+ | Empty+ deriving (Typeable, Generic)+++-- TODO: take from mpfr+fmod :: RealFrac a => a -> a -> a+fmod a b = a - q*b where+ q = realToFrac (truncate $ a / b :: Integer)+{-# INLINE fmod #-}+++instance Precision p => Num (Interval p) where+ I a b + I a' b' = I (a + a') (b + b')+ _ + _ = Empty+ I a b - I a' b' = I (a - coerce b') (b - coerce a')+ _ - _ = Empty+ negate (I a b) = I (coerce (negate b)) (coerce (negate a))+ negate Empty = Empty+ I a b * I a' b' =+ I (minimum [a * a', a * coerce b', coerce b * a', coerce b * coerce b'])+ (maximum [coerce a * coerce a', coerce a * b', b * coerce a', b * b'])+ _ * _ = Empty+ abs x@(I a b)+ | a >= 0 = x+ | b <= 0 = negate x+ | otherwise = I 0 (max (negate (coerce a)) b)+ abs Empty = Empty+ {-# INLINE abs #-}+ signum = increasing signum+ {-# INLINE signum #-}+ fromInteger = I <$> fromInteger <*> fromInteger++-- | lift a monotone increasing function over a given interval+increasing :: (forall r. Rounding r => Rounded r a -> Rounded r b) -> Interval a -> Interval b+increasing f (I a b) = I (f a) (f b)+increasing _ Empty = Empty++-- -- | lift a monotone decreasing function over a given interval+decreasing :: (forall r. Rounding r => Rounded r a -> Rounded r b) -> Interval a -> Interval b+decreasing f (I a b) = I (coerce (f b)) (coerce (f a))+decreasing _ Empty = Empty+--++(...) :: Rounded TowardNegInf p -> Rounded TowardInf p -> Interval p+a ... b+ | coerce a <= b = I a b+ | otherwise = Empty++infixl 6 +/-++(+/-) :: Rounded r p -> Rounded r' p -> Interval p+a +/- b = (coerce a .-. coerce b) ... (coerce a .+. coerce b)++negInfinity :: Fractional a => a+negInfinity = (-1)/0+{-# INLINE negInfinity #-}++posInfinity :: Fractional a => a+posInfinity = 1/0+{-# INLINE posInfinity #-}++-- | create a non-empty interval or fail+interval :: Rounded TowardNegInf p -> Rounded TowardInf p -> Maybe (Interval p)+interval a b+ | coerce a <= b = Just $ I a b+ | otherwise = Nothing+{-# INLINE interval #-}++-- | The whole real number line+--+-- >>> whole+-- -Infinity ... Infinity+whole :: Precision p => Interval p+whole = I negInfinity posInfinity+{-# INLINE whole #-}++-- | An empty interval+--+-- >>> empty+-- Empty+empty :: Interval p+empty = Empty+{-# INLINE empty #-}++-- | Check if an interval is empty+--+-- >>> null (1 ... 5)+-- False+--+-- >>> null (1 ... 1)+-- False+--+-- >>> null empty+-- True+null :: Interval p -> Bool+null Empty = True+null _ = False+{-# INLINE null #-}++-- | The infimum (lower bound) of an interval+--+-- >>> inf (1.0 ... 20.0)+-- 1.0+--+-- >>> inf empty+-- *** Exception: empty interval+inf :: Interval p -> Rounded TowardNegInf p+inf (I a _) = a+inf Empty = error "empty interval"+{-# INLINE inf #-}++-- | The supremum (upper bound) of an interval+--+-- >>> sup (1.0 ... 20.0)+-- 20.0+--+-- >>> sup empty+-- *** Exception: empty interval+sup :: Interval p -> Rounded TowardInf p+sup (I _ b) = b+sup Empty = error "empty interval"+{-# INLINE sup #-}++-- | Is the interval a singleton point?+-- N.B. This is fairly fragile and likely will not hold after+-- even a few operations that only involve singletons+--+-- >>> singular (singleton 1)+-- True+--+-- >>> singular (1.0 ... 20.0)+-- False+singular :: Interval p -> Bool+singular Empty = False+singular (I a b) = coerce a == b+{-# INLINE singular #-}++instance Eq (Interval p) where+ (==) = (==!)+ {-# INLINE (==) #-}++instance Precision p => Ord (Interval p) where+ compare Empty Empty = EQ+ compare Empty _ = LT+ compare _ Empty = GT+ compare (I ax bx) (I ay by)+ | coerce bx < ay = LT+ | coerce ax > by = GT+ | coerce bx == ay && coerce ax == by = EQ+ | otherwise = error "ambiguous comparison"+ {-# INLINE compare #-}++ max (I a b) (I a' b') = I (max a a') (max b b')+ max Empty i = i+ max i Empty = i+ {-# INLINE max #-}++ min (I a b) (I a' b') = I (min a a') (min b b')+ min Empty _ = Empty+ min _ Empty = Empty+ {-# INLINE min #-}++-- | 'realToFrac' will use the midpoint+instance Precision p => Real (Interval p) where+ toRational Empty = error "empty interval"+ toRational (I ra rb) = a + (b - a) / 2 where+ a = toRational ra+ b = toRational rb+ {-# INLINE toRational #-}++instance Precision p => Show (Interval p) where+ showsPrec _ Empty = showString "Empty"+ showsPrec n (I a b) =+ showParen (n > 3) $+ showsPrec 3 a .+ showString " ... " .+ showsPrec 3 b++-- | Calculate the width of an interval.+--+-- >>> width (1 ... 20)+-- 19 ... 19+--+-- >>> width (singleton 1)+-- 0 ... 0+--+-- >>> width empty+-- 0 ... 0+width :: Precision p => Interval p -> Rounded TowardInf p+width (I a b) = b - coerce a+width Empty = 0+{-# INLINE width #-}++-- | Magnitude+--+-- >>> magnitude (1 ... 20)+-- 20+--+-- >>> magnitude (-20 ... 10)+-- 20+--+-- >>> magnitude (singleton 5)+-- 5+--+-- throws 'EmptyInterval' if the interval is empty.+--+-- >>> magnitude empty+-- *** Exception: empty interval+magnitude :: Precision p => Interval p -> Rounded TowardInf p+magnitude = sup . abs+{-# INLINE magnitude #-}++-- | \"mignitude\"+--+-- >>> mignitude (1 ... 20)+-- 1+--+-- >>> mignitude (-20 ... 10)+-- 0+--+-- >>> mignitude (singleton 5)+-- 5+--+-- throws 'EmptyInterval' if the interval is empty.+--+-- >>> mignitude empty+-- *** Exception: empty interval+mignitude :: Precision p => Interval p -> Rounded TowardNegInf p -- TowardZero?+mignitude = inf . abs+{-# INLINE mignitude #-}++-- | Construct a symmetric interval.+--+-- >>> symmetric 3+-- -3 ... 3+symmetric :: Rounded TowardInf p -> Interval p+symmetric b = coerce (negate' b) ... b++-- | Hausdorff distance between intervals.+--+-- >>> distance (1 ... 7) (6 ... 10)+-- 0+--+-- >>> distance (1 ... 7) (15 ... 24)+-- 8+--+-- >>> distance (1 ... 7) (-10 ... -2)+-- 3+--+-- >>> distance Empty (1 ... 1)+-- *** Exception: empty interval+distance :: Precision p => Interval p -> Interval p -> Rounded TowardNegInf p -- TowardZero?+distance i1 i2 = mignitude (i1 - i2)++-- | Inflate an interval by enlarging it at both ends.+--+-- >>> inflate 3 (-1 ... 7)+-- -4 ... 10+--+-- >>> inflate (-2) (0 ... 4)+-- -2 ... 6+--+-- >>> inflate 1 empty+-- Empty+inflate :: Precision p => Rounded TowardInf p -> Interval p -> Interval p+inflate x y = symmetric x + y++{-++-- | Deflate an interval by shrinking it from both ends.+--+-- >>> deflate 3.0 (-4.0 ... 10.0)+-- -1.0 ... 7.0+--+-- >>> deflate 2.0 (-1.0 ... 1.0)+-- Empty+--+-- >>> deflate 1.0 empty+-- Empty+deflate :: => a -> Interval a -> Interval a+deflate _ Empty = Empty+deflate x (I a b) | a' <= b' = I a' b'+ | otherwise = Empty+ where+ a' = a + x+ b' = b - x++-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '<' y@+--+-- >>> (5 ... 10 :: Interval Double) <! (20 ... 30 :: Interval Double)+-- True+--+-- >>> (5 ... 10 :: Interval Double) <! (10 ... 30 :: Interval Double)+-- False+--+-- >>> (20 ... 30 :: Interval Double) <! (5 ... 10 :: Interval Double)+-- False+(<!) :: Precision p => Interval p -> Interval p -> Bool+I _ bx <! I ay _ = coerce bx < ay+_ <! _ = True+{-# INLINE (<!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '<=' y@+--+-- >>> (5 ... 10 :: Interval Double) <=! (20 ... 30 :: Interval Double)+-- True+--+-- >>> (5 ... 10 :: Interval Double) <=! (10 ... 30 :: Interval Double)+-- True+--+-- >>> (20 ... 30 :: Interval Double) <=! (5 ... 10 :: Interval Double)+-- False+(<=!) :: Precision p => Interval p -> Interval p -> Bool+I _ bx <=! I ay _ = coerce bx <= ay+_ <=! _ = True+{-# INLINE (<=!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '==' y@+--+-- Only singleton intervals or empty intervals can return true+--+-- >>> (singleton 5 :: Interval Double) ==! (singleton 5 :: Interval Double)+-- True+--+-- >>> (5 ... 10 :: Interval Double) ==! (5 ... 10 :: Interval Double)+-- False+(==!) :: Interval p -> Interval p -> Bool+I ax bx ==! I ay by = coerce bx == ay && coerce ax == by+_ ==! _ = True+{-# INLINE (==!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '/=' y@+--+-- >>> (5 ... 15 :: Interval Double) /=! (20 ... 40 :: Interval Double)+-- True+--+-- >>> (5 ... 15 :: Interval Double) /=! (15 ... 40 :: Interval Double)+-- False+(/=!) :: Interval p -> Interval p -> Bool+I ax bx /=! I ay by = bx < coerce ay || coerce ax > by+_ /=! _ = True+{-# INLINE (/=!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '>' y@+--+-- >>> (20 ... 40 :: Interval Double) >! (10 ... 19 :: Interval Double)+-- True+--+-- >>> (5 ... 20 :: Interval Double) >! (15 ... 40 :: Interval Double)+-- False+(>!) :: Precision p => Interval p -> Interval p -> Bool+I ax _ >! I _ by = ax > coerce by+_ >! _ = True+{-# INLINE (>!) #-}++-- | For all @x@ in @X@, @y@ in @Y@. @x '>=' y@+--+-- >>> (20 ... 40 :: Interval Double) >=! (10 ... 20 :: Interval Double)+-- True+--+-- >>> (5 ... 20 :: Interval Double) >=! (15 ... 40 :: Interval Double)+-- False+(>=!) :: Precision p => Interval p -> Interval p -> Bool+I ax _ >=! I _ by = coerce ax >= by+_ >=! _ = True++-- | Determine if a point is in the interval.+--+-- >>> elem 3.2 (1 ... 5)+-- True+--+-- >>> elem 5 (1 ... 5)+-- True+--+-- >>> elem 1 (1 ... 5)+-- True+--+-- >>> elem 8 (1 ... 5)+-- False+--+-- >>> elem 5 empty+-- False+--+elem :: Rounded TowardZero p -> Interval p -> Bool+elem x (I a b) = coerce x >= a && coerce x <= b+elem _ Empty = False+{-# INLINE elem #-}++-- | Determine if a point is not included in the interval+--+-- >>> notElem 8 (1.0 ... 5.0)+-- True+--+-- >>> notElem 1.4 (1.0 ... 5.0)+-- False+--+-- And of course, nothing is a member of the empty interval.+--+-- >>> notElem 5 empty+-- True+notElem :: Rounded TowardZero p -> Interval p -> Bool+notElem x xs = not (elem x xs)+{-# INLINE notElem #-}+++-- | For all @x@ in @X@, @y@ in @Y@. @x `op` y@+certainly :: Precision p => (forall b. Ord b => b -> b -> Bool) -> Interval p -> Interval p -> Bool+certainly cmp l r+ | lt && eq && gt = True+ | lt && eq = l <=! r+ | lt && gt = l /=! r+ | lt = l <! r+ | eq && gt = l >=! r+ | eq = l ==! r+ | gt = l >! r+ | otherwise = False+ where+ lt = cmp False True+ eq = cmp True True+ gt = cmp True False+{-# INLINE certainly #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<' y@?+(<?) :: Precision p => Interval p -> Interval p -> Bool+Empty <? _ = False+_ <? Empty = False+I ax _ <? I _ by = coerce ax < by+{-# INLINE (<?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '<=' y@?+(<=?) :: Precision p => Interval p -> Interval p -> Bool+Empty <=? _ = False+_ <=? Empty = False+I ax _ <=? I _ by = coerce ax <= by+{-# INLINE (<=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '==' y@?+(==?) :: Interval a -> Interval a -> Bool+I ax bx ==? I ay by = coerce ax <= by && coerce bx >= ay+_ ==? _ = False+{-# INLINE (==?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '/=' y@?+(/=?) :: Interval a -> Interval a -> Bool+I ax bx /=? I ay by = coerce ax /= by || coerce bx /= ay+_ /=? _ = False+{-# INLINE (/=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>' y@?+(>?) :: Precision p => Interval p -> Interval p -> Bool+I _ bx >? I ay _ = bx > coerce ay+_ >? _ = False+{-# INLINE (>?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x '>=' y@?+(>=?) :: Precision p => Interval p -> Interval p -> Bool+I _ bx >=? I ay _ = bx >= coerce ay+_ >=? _ = False+{-# INLINE (>=?) #-}++-- | Does there exist an @x@ in @X@, @y@ in @Y@ such that @x `op` y@?+possibly :: Precision p => (forall b. Ord b => b -> b -> Bool) -> Interval p -> Interval p -> Bool+possibly cmp l r+ | lt && eq && gt = True+ | lt && eq = l <=? r+ | lt && gt = l /=? r+ | lt = l <? r+ | eq && gt = l >=? r+ | eq = l ==? r+ | gt = l >? r+ | otherwise = False+ where+ lt = cmp LT EQ+ eq = cmp EQ EQ+ gt = cmp GT EQ+{-# INLINE possibly #-}++-- | Check if interval @X@ totally contains interval @Y@+--+-- >>> (20 ... 40 :: Interval Double) `contains` (25 ... 35 :: Interval Double)+-- True+--+-- >>> (20 ... 40 :: Interval Double) `contains` (15 ... 35 :: Interval Double)+-- False+contains :: Precision p => Interval p -> Interval p -> Bool+contains _ Empty = True+contains (I ax bx) (I ay by) = ax <= ay && by <= bx+contains Empty I{} = False+{-# INLINE contains #-}++-- | Flipped version of `contains`. Check if interval @X@ a subset of interval @Y@+--+-- >>> (25 ... 35 :: Interval Double) `isSubsetOf` (20 ... 40 :: Interval Double)+-- True+--+-- >>> (20 ... 40 :: Interval Double) `isSubsetOf` (15 ... 35 :: Interval Double)+-- False+isSubsetOf :: Precision p => Interval p -> Interval p -> Bool+isSubsetOf = flip contains+{-# INLINE isSubsetOf #-}++-- | Calculate the intersection of two intervals.+--+-- >>> intersection (1 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)+-- 5.0 ... 10.0+intersection :: Precision p => Interval p -> Interval p -> Interval p+intersection x@(I a b) y@(I a' b')+ | x /=! y = Empty+ | otherwise = I (max a a') (min b b')+intersection _ _ = Empty+{-# INLINE intersection #-}++-- | Calculate the convex hull of two intervals+--+-- >>> hull (0 ... 10 :: Interval Double) (5 ... 15 :: Interval Double)+-- 0.0 ... 15.0+--+-- >>> hull (15 ... 85 :: Interval Double) (0 ... 10 :: Interval Double)+-- 0.0 ... 85.0+hull :: Precision p => Interval p -> Interval p -> Interval p+hull (I a b) (I a' b') = I (min a a') (max b b')+hull Empty x = x+hull x Empty = x+{-# INLINE hull #-}++-- | Bisect an interval at its midpoint.+--+-- >>> bisect (10.0 ... 20.0)+-- (10.0 ... 15.0,15.0 ... 20.0)+--+-- >>> bisect (singleton 5.0)+-- (5.0 ... 5.0,5.0 ... 5.0)+--+-- >>> bisect Empty+-- (Empty,Empty)+bisect :: Precision p => Interval p -> (Interval p, Interval p)+bisect Empty = (Empty,Empty)+bisect (I a b) = (a...coerce m, succUlp m...b) where m = a + (coerce b - a) / 2+{-# INLINE bisect #-}++-- @'divNonZero' X Y@ assumes @0 `'notElem'` Y@+divNonZero :: Precision p => Interval p -> Interval p -> Interval p+divNonZero (I a b) (I a' b') =+ minimum [a / a', a / coerce b', coerce b / a', coerce b / coerce b']+ ...+ maximum [coerce a / coerce a', coerce a / b', b / coerce a', b / b']+divNonZero _ _ = Empty++-- @'divPositive' X y@ assumes y > 0, and divides @X@ by [0 ... y]+divPositive :: Precision p => Interval p -> Rounded TowardInf p -> Interval p+divPositive Empty _ = Empty+divPositive x@(I a b) y+ | a == 0 && b == 0 = x+ | b < 0 || isNegativeZero b = negInfinity ... (b / y)+ | a < 0 = whole+ | otherwise = (a / coerce y) ... posInfinity+{-# INLINE divPositive #-}++-- divNegative assumes y < 0 and divides the interval @X@ by [y ... 0]+divNegative :: Precision p => Interval p -> Rounded TowardNegInf p -> Interval p+divNegative Empty _ = Empty+divNegative x@(I a b) y+ | a == 0 && b == 0 = negate x -- flip negative zeros+ | b < 0 || isNegativeZero b = (coerce b / y) ... posInfinity+ | a < 0 = whole+ | otherwise = negInfinity ... (coerce a / coerce y)+{-# INLINE divNegative #-}++divZero :: Precision p => Interval p -> Interval p+divZero x@(I a b)+ | a == 0 && b == 0 = x+ | otherwise = whole+divZero Empty = Empty+{-# INLINE divZero #-}++instance Precision p => Fractional (Interval p) where+ -- TODO: check isNegativeZero properly?+ _ / Empty = Empty+ x / y@(I a b)+ | 0 `notElem` y = divNonZero x y+ | iz && sz = error "divide by zero"+ | iz = divPositive x b+ | sz = divNegative x a+ | otherwise = divZero x+ where+ iz = a == 0+ sz = b == 0+ recip Empty = Empty+ recip (I a b) = min (recip $ coerce a) (recip $ coerce b) ... max (recip $ coerce a) (recip $ coerce b)+ {-# INLINE recip #-}+ fromRational = I <$> fromRational <*> fromRational+ {-# INLINE fromRational #-}++midpoint :: Precision p => Interval p -> Rounded TowardNegInf p+midpoint (I a b) = a + (coerce b - a) / 2+midpoint _ = 0/0 -- TODO: use mpfr's nan++instance Precision p => RealFrac (Interval p) where+ properFraction x = (b, x - fromIntegral b)+ where b = truncate (midpoint x)+ {-# INLINE properFraction #-}+ ceiling x = ceiling (sup x)+ {-# INLINE ceiling #-}+ floor x = floor (inf x)+ {-# INLINE floor #-}+ round x = round (midpoint x)+ {-# INLINE round #-}+ truncate x = truncate (midpoint x)+ {-# INLINE truncate #-}++instance Precision p => Floating (Interval p) where+ pi = I pi pi+ {-# INLINE pi #-}++ exp = increasing exp+ {-# INLINE exp #-}++ log (I a b) = (if a > 0 then log a else negInfinity) ... log b+ log Empty = Empty+ {-# INLINE log #-}++ cos Empty = Empty+ cos x+ | width t >= pi = negate 1 ... 1+ | inf t >= pi = negate $ cos (t - pi)+ | sup t <= pi = decreasing cos t+ | sup t <= 2 * pi = negate 1 ... cos (((pi * 2 - sup t) `min` coerce (inf t)))+ | otherwise = negate 1 ... 1+ where+ t = fmod x (pi * 2)+ {-# INLINE cos #-}++ sin Empty = Empty+ sin x = cos (x - pi / 2)+ {-# INLINE sin #-}++ tan Empty = Empty+ tan x+ | inf t' <= negate pi / 2 || sup t' >= pi / 2 = whole+ | otherwise = increasing tan x+ where+ t = x `fmod` pi+ t' | t >= pi / 2 = t - pi+ | otherwise = t+ {-# INLINE tan #-}++ asin Empty = Empty+ asin (I a b)+ | b < -1 || a > 1 = Empty+ | otherwise =+ (if a <= -1 then negate pi / 2 else asin a)+ ...+ (if b >= 1 then pi / 2 else asin b)+ {-# INLINE asin #-}++ acos Empty = Empty+ acos (I a b)+ | b < -1 || a > 1 = Empty+ | otherwise =+ (if b >= 1 then 0 else acos (coerce b))+ ...+ (if a < -1 then pi else acos (coerce a))+ {-# INLINE acos #-}++ atan = increasing atan+ {-# INLINE atan #-}++ sinh = increasing sinh+ {-# INLINE sinh #-}++ cosh Empty = Empty+ cosh x@(I a b)+ | b < 0 = decreasing cosh x+ | a >= 0 = increasing cosh x+ | otherwise = I 0 $ cosh $ if negate a > coerce b then coerce a else b+ {-# INLINE cosh #-}++ tanh = increasing tanh+ {-# INLINE tanh #-}++ asinh = increasing asinh+ {-# INLINE asinh #-}++ acosh Empty = Empty+ acosh (I a b)+ | b < 1 = Empty+ | otherwise = I lo $ acosh b+ where lo | a <= 1 = 0+ | otherwise = acosh a+ {-# INLINE acosh #-}++ atanh Empty = Empty+ atanh (I a b)+ | b < -1 || a > 1 = Empty+ | otherwise =+ (if a <= - 1 then negInfinity else atanh a)+ ...+ (if b >= 1 then posInfinity else atanh b)+ {-# INLINE atanh #-}
+ src/Numeric/Rounded/Precision.hs view
@@ -0,0 +1,84 @@+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# OPTIONS_HADDOCK not-home #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Rounded.Precision+-- Copyright : (C) 2012 Edward Kmett+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+----------------------------------------------------------------------------+module Numeric.Rounded.Precision+ ( Precision(..)+ , reifyPrecision+ , Bytes+ ) where++import Data.Proxy+import Data.Reflection+import Foreign.C.Types+import GHC.TypeLits++import Numeric.LongDouble (LongDouble)+import Numeric.MPFR.Types++-- | This class is used to specify the number of bits of precision that are maintained in the+-- significand of a properly 'Numeric.Rounded.Rounded' floating point number.+class Precision p where+ precision :: proxy p -> Int++floatPrecision :: RealFloat a => p a -> Int+floatPrecision p = fromIntegral (floatDigits (proxyArg p)) where+ proxyArg :: p a -> a+ proxyArg _ = undefined+{-# INLINE floatPrecision #-}++instance Precision Float where+ precision = floatPrecision++instance Precision CFloat where+ precision = floatPrecision++instance Precision Double where+ precision = floatPrecision++instance Precision CDouble where+ precision = floatPrecision++instance Precision LongDouble where+ precision = floatPrecision++instance KnownNat n => Precision (n :: Nat) where+ precision p = max MPFR_PREC_MIN . min MPFR_PREC_MAX $ fromInteger (natVal p)++data Bytes (n :: Nat)++instance KnownNat n => Precision (Bytes n) where+ precision _ = max MPFR_PREC_MIN . min MPFR_PREC_MAX $ 8 * fromInteger (natVal (undefined :: Bytes n))++data ReifiedPrecision (s :: *)++retagReifiedPrecision :: (Proxy s -> a) -> proxy (ReifiedPrecision s) -> a+retagReifiedPrecision f _ = f Proxy+{-# INLINE retagReifiedPrecision #-}++instance Reifies s Int => Precision (ReifiedPrecision s) where+ precision = retagReifiedPrecision reflect++reifyPrecision :: Int -> (forall (p :: *). Precision p => Proxy p -> a) -> a+reifyPrecision m f = reify m (go f) where+ go :: Reifies p Int => (Proxy (ReifiedPrecision p) -> a) -> proxy p -> a+ go g _ = g Proxy+{-# INLINE reifyPrecision #-}
+ src/Numeric/Rounded/Rounding.hs view
@@ -0,0 +1,93 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# OPTIONS_HADDOCK not-home #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Rounded.Rounding+-- Copyright : (C) 2012-2014 Edward Kmett+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+-- Various rounding modes+----------------------------------------------------------------------------+module Numeric.Rounded.Rounding+ ( Rounding(..)+ , RoundingMode(..)+ , reifyRounding+ ) where++import Data.Data+import Data.Singletons++import Numeric.MPFR.Types++data RoundingMode+ = TowardNearestWithTiesAwayFromZero -- ^ currently unsupported placeholder+ | TowardNearest -- ^ roundTiesToEven in IEEE 754-2008+ | TowardZero -- ^ roundTowardZero in IEEE 754-2008+ | TowardInf -- ^ roundTowardPositive in IEEE 754-2008+ | TowardNegInf -- ^ roundTowardNegative in IEEE 754-2008+ | AwayFromZero -- ^ round away from zero+ | Faithfully -- ^ currently unsupported placeholder+ deriving (Eq,Ord,Show,Read,Data,Typeable)++class Rounding (r :: RoundingMode) where rounding :: Proxy r -> RoundingMode+instance Rounding TowardNearest where rounding _ = TowardNearest+instance Rounding TowardZero where rounding _ = TowardZero+instance Rounding TowardInf where rounding _ = TowardInf+instance Rounding TowardNegInf where rounding _ = TowardNegInf+instance Rounding AwayFromZero where rounding _ = AwayFromZero+instance Rounding Faithfully where rounding _ = Faithfully+instance Rounding TowardNearestWithTiesAwayFromZero where rounding _ = TowardNearestWithTiesAwayFromZero++instance Enum RoundingMode where+ toEnum MPFR_RNDNA = TowardNearestWithTiesAwayFromZero+ toEnum MPFR_RNDN = TowardNearest+ toEnum MPFR_RNDZ = TowardZero+ toEnum MPFR_RNDU = TowardInf+ toEnum MPFR_RNDD = TowardNegInf+ toEnum MPFR_RNDA = AwayFromZero+ toEnum MPFR_RNDF = Faithfully+ toEnum _ = error "out of range"++ fromEnum TowardNearestWithTiesAwayFromZero = MPFR_RNDNA+ fromEnum TowardNearest = MPFR_RNDN+ fromEnum TowardZero = MPFR_RNDZ+ fromEnum TowardInf = MPFR_RNDU+ fromEnum TowardNegInf = MPFR_RNDD+ fromEnum AwayFromZero = MPFR_RNDA+ fromEnum Faithfully = MPFR_RNDF++instance Bounded RoundingMode where+ minBound = TowardNearestWithTiesAwayFromZero+ maxBound = Faithfully++newtype instance Sing (m :: RoundingMode) = SRounding RoundingMode++instance SingI TowardNearestWithTiesAwayFromZero where sing = SRounding TowardNearestWithTiesAwayFromZero+instance SingI TowardNearest where sing = SRounding TowardNearest+instance SingI TowardZero where sing = SRounding TowardZero+instance SingI TowardInf where sing = SRounding TowardInf+instance SingI TowardNegInf where sing = SRounding TowardNegInf+instance SingI AwayFromZero where sing = SRounding AwayFromZero+instance SingI Faithfully where sing = SRounding Faithfully++reifyRounding :: RoundingMode -> (forall s. Rounding s => Proxy s -> r) -> r+reifyRounding TowardNearestWithTiesAwayFromZero f = f (Proxy :: Proxy TowardNearestWithTiesAwayFromZero)+reifyRounding TowardNearest f = f (Proxy :: Proxy TowardNearest)+reifyRounding TowardZero f = f (Proxy :: Proxy TowardZero)+reifyRounding TowardInf f = f (Proxy :: Proxy TowardInf)+reifyRounding TowardNegInf f = f (Proxy :: Proxy TowardNegInf)+reifyRounding AwayFromZero f = f (Proxy :: Proxy AwayFromZero)+reifyRounding Faithfully f = f (Proxy :: Proxy Faithfully)+{-# INLINE reifyRounding #-}
+ src/Numeric/Rounded/Simple.hs view
@@ -0,0 +1,384 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE RankNTypes #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.Rounded.Simple+-- Copyright : (C) 2012-2014 Edward Kmett, Daniel Peebles+-- (C) 2013-2018 Claude Heiland-Allen+-- License : LGPL+-- Maintainer : Claude Heiland-Allen <claude@mathr.co.uk>+-- Stability : experimental+-- Portability : non-portable+--+-- This module provides an interface without advanced type system features,+-- that may be more convenient if the precision is changed often.+----------------------------------------------------------------------------+module Numeric.Rounded.Simple+ (+ -- * Floating point numbers with a specified rounding mode and precision+ Rounded()+ , reifyRounded+ , simplify+ , fromInt+ , fromDouble+ , fromLongDouble+ , toDouble+ , toLongDouble+ , toInteger'+ , precRound+ -- * Precision+ , Precision+ , precision+ -- * Rounding+ , RoundingMode(..)+ -- * Useful Constants+ , kPi+ , kLog2+ , kEuler+ , kCatalan+ -- * Ord+ , min_+ , max_+ -- * Num+ , add_+ , sub_+ , mul_+ , negate_+ , abs_+-- , signum_+ , fromInteger'+ -- * Fractional+ , div_+-- , recip_+ , fromRational'+ -- * Real+ , toRational'+ -- * RealFrac+ , properFraction_+ , truncate_+ , round_+ , ceiling_+ , floor_+ -- * Floating+ , sqrt_+ , exp_+ , expm1_+ , log_+ , log1p_+ , sin_+ , cos_+ , tan_+ , asin_+ , acos_+ , atan_+ , sinh_+ , cosh_+ , tanh_+ , asinh_+ , acosh_+ , atanh_+ -- * RealFloat+ , atan2_+ , floatRadix'+ , floatDigits'+ , floatRange'+ , decodeFloat'+ , encodeFloat'+ , exponent'+ , significand'+ , scaleFloat'+ , isNaN'+ , isInfinite'+ , isDenormalized'+ , isNegativeZero'+ , isIEEE'+ -- * Show+ , show'+ -- * Read+ , read'+ -- * Foreign Function Interface+ , withInRounded+ , withInOutRounded+ , withInOutRounded_+ , withOutRounded+ , withOutRounded_+ , peekRounded+ ) where++import Control.Exception (bracket_)+import Foreign (Ptr(..), alloca)+import GHC.Prim ( ByteArray# )++import Numeric.LongDouble (LongDouble)++import Numeric.MPFR.Types+import Numeric.MPFR.Raw (mpfr_init2, mpfr_clear, mpfr_set)+import qualified Numeric.Rounded as R+import qualified Numeric.Rounded.Internal as R+import Numeric.Rounded.Rounding++type Precision = Int++-- | A properly rounded floating-point number with a given rounding mode and precision.+data Rounded = Rounded+ { roundedPrec :: !MPFRPrec+ , _roundedSign :: !MPFRSign+ , _roundedExp :: !MPFRExp+ , _roundedLimbs :: !ByteArray#+ }++precision :: Rounded -> Precision+precision = fromIntegral . roundedPrec++reifyRounded :: Rounded -> (forall p . R.Precision p => R.Rounded r p -> a) -> a+reifyRounded (Rounded p s e l) f = R.reifyPrecision (fromIntegral p) (\q -> f (g q (R.Rounded p s e l)))+ where+ g :: R.Precision q => proxy q -> R.Rounded s q -> R.Rounded s q+ g _ x = x++simplify :: R.Rounded r p -> Rounded+simplify (R.Rounded p s e l) = Rounded p s e l++constant :: (forall r p . (R.Rounding r, R.Precision p) => R.Rounded r p) -> RoundingMode -> Precision -> Rounded+constant f r q = R.reifyRounding r (\pr -> R.reifyPrecision q (\pq -> g pr pq f))+ where+ g :: (R.Rounding r, R.Precision p) => proxy1 r -> proxy2 p -> R.Rounded r p -> Rounded+ g _ _ b = simplify b++kPi, kLog2, kEuler, kCatalan :: RoundingMode -> Precision -> Rounded++kPi = constant pi+kLog2 = constant R.kLog2+kEuler = constant R.kEuler+kCatalan = constant R.kCatalan++unary :: (forall r p q . (R.Rounding r, R.Precision p, R.Precision q) => R.Rounded r p -> R.Rounded r q) -> RoundingMode -> Precision -> Rounded -> Rounded+unary f r q a = R.reifyRounding r (\pr -> R.reifyPrecision q (\pq -> reifyRounded a (\ra -> g pr pq f ra)))+ where+ g :: (R.Rounding r, R.Precision p, R.Precision q) => proxy1 r -> proxy2 q -> (R.Rounded r p -> R.Rounded r q) -> R.Rounded r p -> Rounded+ g _ _ h b = simplify (h b)++abs_, negate_, log_, exp_, sqrt_,+ sin_, cos_, tan_, asin_, acos_, atan_,+ sinh_, cosh_, tanh_, asinh_, acosh_, atanh_,+ log1p_, expm1_,+ precRound :: RoundingMode -> Precision -> Rounded -> Rounded++abs_ = unary R.abs_+negate_ = unary R.negate_+log_ = unary R.log_+exp_ = unary R.exp_+sqrt_ = unary R.sqrt_+sin_ = unary R.sin_+cos_ = unary R.cos_+tan_ = unary R.tan_+asin_ = unary R.asin_+acos_ = unary R.acos_+atan_ = unary R.atan_+sinh_ = unary R.sinh_+cosh_ = unary R.cosh_+tanh_ = unary R.tanh_+asinh_ = unary R.asinh_+acosh_ = unary R.acosh_+atanh_ = unary R.atanh_+log1p_ = unary R.log1p_+expm1_ = unary R.expm1_+precRound = unary R.precRound++fromInt :: RoundingMode -> Precision -> Int -> Rounded+fromInt = fromX R.fromInt++fromDouble :: RoundingMode -> Precision -> Double -> Rounded+fromDouble = fromX R.fromDouble++fromLongDouble :: RoundingMode -> Precision -> LongDouble -> Rounded+fromLongDouble = fromX R.fromLongDouble++fromInteger' :: RoundingMode -> Precision -> Integer -> Rounded+fromInteger' = fromX fromInteger++fromRational' :: RoundingMode -> Precision -> Rational -> Rounded+fromRational' = fromX fromRational++fromX :: (forall r p . (R.Rounding r, R.Precision p) => x -> R.Rounded r p) -> RoundingMode -> Precision -> x -> Rounded+fromX f r p x = R.reifyRounding r (\pr -> R.reifyPrecision p (\pp -> g pr pp (f x)))+ where+ g :: (R.Rounding r, R.Precision p) => proxy1 r -> proxy2 p -> R.Rounded r p -> Rounded+ g _ _ x = simplify x++binary :: (forall r p q pq . (R.Rounding r, R.Precision p, R.Precision q, R.Precision pq) => R.Rounded r p -> R.Rounded r q -> R.Rounded r pq) -> RoundingMode -> Precision -> Rounded -> Rounded -> Rounded+binary f r pq a b = R.reifyRounding r (\pr -> R.reifyPrecision pq (\ppq -> reifyRounded a (\ra -> reifyRounded b (\rb -> g pr ppq f ra rb))))+ where+ g :: (R.Rounding r, R.Precision p, R.Precision q, R.Precision pq) => proxy1 r -> proxy2 pq -> (R.Rounded r p -> R.Rounded r q -> R.Rounded r pq) -> R.Rounded r p -> R.Rounded r q -> Rounded+ g _ _ h x y = simplify (h x y)++binary' :: (forall r p q pq . (R.Rounding r, R.Precision p, R.Precision q, R.Precision pq) => R.Rounded r p -> R.Rounded r q -> R.Rounded r pq) -> Rounded -> Rounded -> Rounded+binary' f a b = binary f R.TowardNearest (precision a `max` precision b) a b++min_, max_, add_, sub_, mul_, div_, atan2_ :: RoundingMode -> Precision -> Rounded -> Rounded -> Rounded++min_ = binary R.min_+max_ = binary R.max_+add_ = binary (R.!+!)+sub_ = binary (R.!-!)+mul_ = binary (R.!*!)+div_ = binary (R.!/!)+atan2_ = binary R.atan2_++unary' :: (forall r p . (R.Rounding r, R.Precision p) => R.Rounded r p -> a) -> RoundingMode -> Rounded -> a+unary' f r a = R.reifyRounding r (\pr -> reifyRounded a (\ra -> g pr f ra))+ where+ g :: (R.Rounding r, R.Precision p) => proxy r -> (R.Rounded r p -> a) -> R.Rounded r p -> a+ g _ h x = h x++unary'' :: (forall r p . (R.Rounding r, R.Precision p) => R.Rounded r p -> a) -> Rounded -> a+unary'' f a = unary' f R.TowardNearest a++toDouble :: RoundingMode -> Rounded -> Double+toDouble = unary' R.toDouble++toLongDouble :: RoundingMode -> Rounded -> LongDouble+toLongDouble = unary' R.toLongDouble++toInteger' :: RoundingMode -> Rounded -> Integer+toInteger' = unary' R.toInteger'++-- Real++toRational' :: RoundingMode -> Rounded -> Rational+toRational' = unary' toRational++-- RealFloat++floatRadix' :: Rounded -> Integer+floatRadix' = unary'' floatRadix++floatDigits' :: Rounded -> Int+floatDigits' = unary'' floatDigits++floatRange' :: Rounded -> (Int, Int)+floatRange' = unary'' floatRange++decodeFloat' :: Rounded -> (Integer, Int)+decodeFloat' = unary'' decodeFloat++encodeFloat' :: RoundingMode -> Precision -> Integer -> Int -> Rounded+encodeFloat' r p m e = R.reifyRounding r (\rp -> R.reifyPrecision p (\pp -> g rp pp (encodeFloat m e)))+ where+ g :: R.Precision p => proxy1 r -> proxy2 p -> R.Rounded r p -> Rounded+ g _ _ x = simplify x++exponent' :: Rounded -> Int+exponent' = unary'' exponent++significand' :: Rounded -> Rounded+significand' = unary'' (\a -> simplify (significand a))++scaleFloat' :: Int -> Rounded -> Rounded+scaleFloat' n = unary'' (\a -> simplify (scaleFloat n a))++isNaN' :: Rounded -> Bool+isNaN' = unary'' isNaN++isInfinite' :: Rounded -> Bool+isInfinite' = unary'' isInfinite++isDenormalized' :: Rounded -> Bool+isDenormalized' = unary'' isDenormalized++isNegativeZero' :: Rounded -> Bool+isNegativeZero' = unary'' isNegativeZero++isIEEE' :: Rounded -> Bool+isIEEE' = unary'' isIEEE++-- RealFrac++properFraction_ :: Integral i => Rounded -> (i, Rounded)+properFraction_ a = reifyRounded a g+ where+ g :: (Integral j, R.Precision p) => R.Rounded R.TowardNearest p -> (j, Rounded)+ g ra = case properFraction ra of (i, b) -> (i, simplify b)++truncate_, ceiling_, floor_, round_ :: Precision -> Rounded -> Rounded+truncate_ = unary R.truncate_ TowardNearest+round_ = unary R.round_ TowardNearest+ceiling_ = unary R.ceiling_ TowardNearest+floor_ = unary R.floor_ TowardNearest++type Comparison = Rounded -> Rounded -> Bool++cmp :: (forall p q . (R.Precision p, R.Precision q) => R.Rounded R.TowardNearest p -> R.Rounded R.TowardNearest q -> Bool) -> Comparison+cmp f a b = reifyRounded a (\ra -> reifyRounded b (\rb -> f ra rb))++instance Eq Rounded where+ (==) = cmp (R.!==!)+ (/=) = cmp (R.!/=!)++instance Ord Rounded where+ compare a b = reifyRounded a (\ra -> reifyRounded b (\rb -> R.compare_ ra rb))+ (<) = cmp (R.!<!)+ (<=) = cmp (R.!<=!)+ (>) = cmp (R.!>!)+ (>=) = cmp (R.!>=!)+ max = binary' R.max_+ min = binary' R.min_++-- Show++show' :: Rounded -> String+show' = unary'' show++-- Read++read' :: RoundingMode -> Precision -> String -> Rounded+read' r p s = R.reifyRounding r (\pr -> R.reifyPrecision p (\pp -> g pr pp (read s)))+ where+ g :: (R.Rounding r, R.Precision p) => proxy1 r -> proxy2 p -> R.Rounded r p -> Rounded+ g _ _ x = simplify x++-- Foreign Function Interface++-- | Use a value as a /constant/ @mpfr_t@ (attempts to modify it may explode,+-- changing the precision will explode).+withInRounded :: Rounded -> (Ptr MPFR -> IO a) -> IO a+withInRounded a f = reifyRounded a (\ra -> R.withInRounded ra f)++-- | Allocates and initializes a new @mpfr_t@, after the action it is peeked+-- and returned.+withOutRounded :: Precision -> (Ptr MPFR -> IO a) -> IO (Rounded, a)+withOutRounded prec f = r where+ r = alloca $ \ptr -> bracket_ (mpfr_init2 ptr (fromIntegral prec)) (mpfr_clear ptr) $ do+ a <- f ptr+ m <- peekRounded ptr+ return (m, a)++-- | Allocates and initializes a new @mpfr_t@, after the action it is peeked+-- and returned.+-- The result of the action is ignored.+withOutRounded_ :: Precision -> (Ptr MPFR -> IO a) -> IO Rounded+withOutRounded_ p = fmap fst . withOutRounded p++-- | Allocates and initializes a new @mpfr_t@ to the value. After the action+-- it is peeked and returned.+withInOutRounded :: Rounded -> (Ptr MPFR -> IO a) -> IO (Rounded, a)+-- FIXME: optimize to reduce copying+withInOutRounded i f =+ withOutRounded (fromIntegral (roundedPrec i)) $ \ofr ->+ withInRounded i $ \ifr -> do+ _ <- mpfr_set ofr ifr (fromIntegral (fromEnum TowardNearest))+ f ofr++-- | Allocates and initializes a new @mpfr_t@ to the value. After the action+-- it is peeked and returned.+-- The result of the action is ignored.+withInOutRounded_ :: Rounded -> (Ptr MPFR -> IO a) -> IO Rounded+withInOutRounded_ x = fmap fst . withInOutRounded x++-- | Peek an @mpfr_t@ at its actual precision.+peekRounded :: Ptr MPFR -> IO Rounded+peekRounded ptr = R.peekRounded ptr f+ where+ f :: R.Precision p => R.Rounded TowardNearest p -> IO Rounded+ f mr = return (simplify mr)
+ test.hs view
@@ -0,0 +1,58 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}++import Control.Exception (catch, SomeException)+import Data.Proxy (Proxy(..))+import System.Exit (exitSuccess, exitFailure)++import Numeric.LongDouble (LongDouble)++import Numeric.Rounded++filename :: String+filename = "test.txt"++main :: IO ()+main = do+ golden <- readFile filename `Control.Exception.catch` update+ if golden == test then exitSuccess else exitFailure++update :: SomeException -> IO String+update _ = writeFile filename test >> readFile filename++pf :: RealFrac a => a -> (Integer, a)+pf = properFraction++test :: String+test = unlines+ [ show (exp pi :: Rounded TowardZero 512)+ , show (pi :: Rounded TowardZero Double)+ , show (pi :: Rounded AwayFromZero Double)+ , show (kCatalan :: Rounded TowardZero 128)+ , (reifyPrecision 512 (\(_ :: Proxy p) -> show (logBase 10 2 :: Rounded TowardNearest p)))+ , (reifyRounding TowardZero (\(_ :: Proxy r) -> show (logBase 10 2 :: Rounded r 512)))+ , show (fromDouble pi - pi :: Rounded TowardNearest 64)+ , show (fromInt 100000000 :: Rounded TowardNearest Float)+ , show (fromInt 123456789 :: Rounded TowardNearest Float)+ , show (realToFrac (pi :: Rounded TowardNearest 512) :: Double)+ , show . pf $ (-2.5 :: Rational)+ , show . pf $ (-1.5 :: Rational)+ , show . pf $ (-0.5 :: Rational)+ , show . pf $ ( 0.5 :: Rational)+ , show . pf $ ( 1.5 :: Rational)+ , show . pf $ ( 2.5 :: Rational)+ , show . pf $ (-2.5 :: Rounded TowardNearest Float)+ , show . pf $ (-1.5 :: Rounded TowardNearest Float)+ , show . pf $ (-0.5 :: Rounded TowardNearest Float)+ , show . pf $ ( 0.5 :: Rounded TowardNearest Float)+ , show . pf $ ( 1.5 :: Rounded TowardNearest Float)+ , show . pf $ ( 2.5 :: Rounded TowardNearest Float)+ , show . pf $ (-(2^23 + 0.5) :: Rounded TowardNearest Float)+ , show . pf $ (-(2^22 + 0.5) :: Rounded TowardNearest Float)+ , show . pf $ (-(2^21 + 0.5) :: Rounded TowardNearest Float)+ , show . pf $ ( (2^21 + 0.5) :: Rounded TowardNearest Float)+ , show . pf $ ( (2^22 + 0.5) :: Rounded TowardNearest Float)+ , show . pf $ ( (2^23 + 0.5) :: Rounded TowardNearest Float)+ , show (fromLongDouble pi == (pi :: Rounded TowardNearest LongDouble))+ , show (pi == toLongDouble (pi :: Rounded TowardNearest LongDouble))+ ]
+ test.txt view
@@ -0,0 +1,30 @@+23.1406926327792690057290863679485473802661062426002119934450464095243423506904527835169719970675492196759527048010877731444280444146938358447174458796098418+3.1415926535897931+3.1415926535897936+0.9159655941772190150546035149323841107734+0.301029995663981195213738894724493026768189881462108541310427461127108189274424509486927252118186172040684477191430995379094767881133523505999692333704695598+0.301029995663981195213738894724493026768189881462108541310427461127108189274424509486927252118186172040684477191430995379094767881133523505999692333704695561+-1.22514845490862001043e-16+1.0e8+1.23456792e8+3.141592653589793+(-2,(-1) % 2)+(-1,(-1) % 2)+(0,(-1) % 2)+(0,1 % 2)+(1,1 % 2)+(2,1 % 2)+(-2,-0.5)+(-1,-0.5)+(0,-0.5)+(0,0.5)+(1,0.5)+(2,0.5)+(-8388608,-0.0)+(-4194304,-0.5)+(-2097152,-0.5)+(2097152,0.5)+(4194304,0.5)+(8388608,0.0)+True+True