diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown new file mode 100644 --- /dev/null +++ b/CHANGELOG.markdown @@ -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 diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,165 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 3, 29 June 2007 + + Copyright (C) 2007 Free Software Foundation, Inc. + 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. diff --git a/README.markdown b/README.markdown new file mode 100644 --- /dev/null +++ b/README.markdown @@ -0,0 +1,67 @@ +rounded +======= + +[![Build Status](https://secure.travis-ci.org/ekmett/rounded.png?branch=master)](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 diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/cbits/wrappers.c b/cbits/wrappers.c new file mode 100644 --- /dev/null +++ b/cbits/wrappers.c @@ -0,0 +1,79 @@ +#include + +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; +} diff --git a/rounded.cabal b/rounded.cabal new file mode 100644 --- /dev/null +++ b/rounded.cabal @@ -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 +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 + + 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 diff --git a/src/Numeric/MPFR/Raw.hs b/src/Numeric/MPFR/Raw.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MPFR/Raw.hs @@ -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 +-- 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 diff --git a/src/Numeric/MPFR/Raw/Safe.hs b/src/Numeric/MPFR/Raw/Safe.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MPFR/Raw/Safe.hs @@ -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 +-- 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 diff --git a/src/Numeric/MPFR/Raw/Unsafe.hs b/src/Numeric/MPFR/Raw/Unsafe.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/MPFR/Raw/Unsafe.hs @@ -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 +-- 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 diff --git a/src/Numeric/MPFR/Types.hsc b/src/Numeric/MPFR/Types.hsc new file mode 100644 --- /dev/null +++ b/src/Numeric/MPFR/Types.hsc @@ -0,0 +1,127 @@ +#include +#include + +#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) diff --git a/src/Numeric/Rounded.hs b/src/Numeric/Rounded.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Rounded.hs @@ -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 +-- 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_ + , (!==!) + , (!/=!) + , (!>=!) + , (!<=!) + , (!>!) + , (! +-- 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_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)." +-- +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''# diff --git a/src/Numeric/Rounded/Interval.hs b/src/Numeric/Rounded/Interval.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Rounded/Interval.hs @@ -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) >> (5 ... 10 :: Interval Double) >> (20 ... 30 :: Interval Double) Interval p -> Interval p -> Bool +I _ bx >> (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 = 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@? +( Interval p -> Interval p -> Bool +Empty 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 = 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 #-} diff --git a/src/Numeric/Rounded/Precision.hs b/src/Numeric/Rounded/Precision.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Rounded/Precision.hs @@ -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 +-- 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 #-} diff --git a/src/Numeric/Rounded/Rounding.hs b/src/Numeric/Rounded/Rounding.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Rounded/Rounding.hs @@ -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 +-- 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 #-} diff --git a/src/Numeric/Rounded/Simple.hs b/src/Numeric/Rounded/Simple.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/Rounded/Simple.hs @@ -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 +-- 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.!>=!) + 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) diff --git a/test.hs b/test.hs new file mode 100644 --- /dev/null +++ b/test.hs @@ -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)) + ] diff --git a/test.txt b/test.txt new file mode 100644 --- /dev/null +++ b/test.txt @@ -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