logfloat 0.12.1 → 0.14.0.2
raw patch · 10 files changed
Files
- CHANGELOG +113/−0
- LICENSE +1/−1
- README.md +238/−0
- logfloat.cabal +61/−32
- src/Data/Number/LogFloat.hs +321/−204
- src/Data/Number/LogFloat/Raw.hs +405/−0
- src/Data/Number/PartialOrd.hs +33/−17
- src/Data/Number/RealToFrac.hs +49/−16
- src/Data/Number/Transfinite.hs +41/−32
- src/Hugs/RealFloat.hs +10/−13
+ CHANGELOG view
@@ -0,0 +1,113 @@+0.14.0.2 (2026-02-28):+ - Nudged version bounds for GHC 9.14 (didn't actually need to+ change the upper bounds on 'base' and 'array', because they+ were already lenient)+ - Silenced GHC >= 9.10 warnings about importing 'Data.List'+0.14.0.1 (2024-08-29):+ - Nudged version bounds for ghc 9.10+0.14.0 (2021-03-19):+ - Separating out "Data.Number.LogFloat.Raw"+ - Added `Tested-With: GHC == 9.2.4, 9.4.4, 9.6.1` (didn't actually need to+ nudge the upper bound on 'base', because it's already lenient)+ - Silencing warnings about rewrite rules, as a hack for working around:+ <https://gitlab.haskell.org/ghc/ghc/-/issues/10595>+0.13.4 (2017-06-18):+ - Added Read LogFloat instance. (h/t Rob Zinkov)+ - Corrected a bug where `product` returns NaN on `logFloat 0`.+0.13.3.3 (2015-10-02):+ - ???+0.13.3.2 (2015-08-06):+ - Fixed the buggy Show LogFloat instance+0.13.3.1 (2015-05-30):+ - Moved VERSION to CHANGELOG+0.13.3 (2015-03-29):+ - Got clean compile on GHC 7.10.1+0.13.2 (2015-03-23):+ - Preliminary fixes for dealing with type-roles in GHC 7.10+0.13.1 (2015-03-10):+ - Fixed a major bug in sum+0.13 (2015-02-17):+ - monomorphized logFloat, logToLogFloat, fromLogFloat, logFromLogFloat+ - added sum, product, pow+ - cleaned up some extraneous rewrite rules+ - removed support for older versions of GHC, in order to clean things up++0.12.1 (2010-03-19):+ - Fixed some NaN injection bugs in the Num instance. These would+ only have been triggered by non-probabilistic uses of LogFloat+ (i.e., if you used @logFloat infinity@)+0.12.0.3 (2009-04-03):+ - Added notes to INSTALL about Windows compatibility+0.12.0.2 (2009-04-03):+ - Updated TODO and VERSION+0.12.0.1 (2009-03-25):+ - Disabled -fvia-C so the FFI can be used in GHC 6.10+0.12 (2009-03-10):+ - Added log1p.+ - Added Storable instance (for GHC only)+ - removed orphaned toRational/fromRational rules. The+ RealToFrac module obviates those optimizations.+ - Adjusted Real LogFloat instance to throw errors on transfinite values++0.11.2 (2009-03-09):+ - Moved log/exp rules from LogFloat into Transfinite+0.11.1.1 (2009-03-09):+ - Documentation updates+0.11.1 (2009-03-08):+ - Added IArray UArray LogFloat instance (thanks to Felipe Lessa).+0.11 (2009-01-29):+ - Moved the RealToFrac class from Data.Number.Transfinite to+ Data.Number.RealToFrac. This breaks backwards compatibility if+ you depended directly on that module for the class. If you+ imported the class from Data.Number.LogFloat then you're still+ fine.+ - Cabal build fail due to missing new module (Fixed in 0.11.0.1).++0.10 (2009-01-29):+ - Hugs' Prelude definitions of isInfinite and isNaN for Float+ and Double are buggy. The new Hugs.RealFloat module provides+ correct implementations (defaulting to the Prelude definition+ for non-Hugs compilers).+ - The PartialOrd and Transfinite instances for Double and Float+ have been updated to use these corrected functions.+ - Added maxPO and minPO to PartialOrd and added comparingPO to+ Data.Number.PartialOrd.+ - Minor changes to circumvent bugs in Hackage's new version of+ Haddock. You should be able to see all the documentation now.++0.9.1 (2008-08-30):+ - Fixed some PartialOrd stuff and sanitized documentation+ - Build fail on GHC 6.10 due to GHC.Prim (Fixed in 0.9.1.1).+0.9.0 (2008-08-30):+ - s/toFractional/realToFrac/g.+ - Also moved realToFrac and log to Transfinite+ - Build fail on GHC 6.10 due to GHC.Prim++0.8.6 (2008-08-17):+ - Removed buggy RULES+ - Build fail on GHC 6.10 due to -Werror and removal of -fno-warn-orphans+0.8.5 (2008-08-17):+ - Gave up and converted from lhs to hs so Hackage docs work+0.8.4 (2008-08-17):+ - Broke out Transfinite+0.8.3 (2008-08-16):+ - Documentation updates+0.8.2 (2008-08-14):+ - Announced release+0.8:+ - Did a bunch of tweaking. Things should be decent now++0.7:+ - Haddockified+0.6:+ - Fixed monomorphism of LogFloat constructors/destructors.+0.5:+ - Added optimization rules.+0.4:+ - Translated from Perl to Haskell at revision 2007.12.20.+0.3:+ - Converted extensive comments to POD format.+0.2:+ - Did a bunch of profiling, optimizing, and debugging.+0.1:+ - Initial version (in Perl) created for hw5 for NLP with Jason Eisner.
LICENSE view
@@ -1,4 +1,4 @@-Copyright (c) 2007--2010, wren ng thornton.+Copyright (c) 2007--2015, wren gayle romano. ALL RIGHTS RESERVED. Redistribution and use in source and binary forms, with or without
+ README.md view
@@ -0,0 +1,238 @@+logfloat+========+[](https://github.com/wrengr/logfloat/actions?query=workflow%3Aci+-event%3Apull_request)+[](https://hackage.haskell.org/package/logfloat)+[](https://stackage.org/lts/package/logfloat)+[](https://stackage.org/nightly/package/logfloat)++This package provides a type for storing numbers in the log-domain,+primarily useful for preventing underflow when multiplying many+probabilities as in HMMs and other probabilistic models. The package+also provides modules for dealing with floating numbers correctly.++As of version 0.13.3, we've dropped support for Hugs and GHC < 7.8.+Nothing major has changed, so they should still work; it's just+that they're no longer officially supported. (In practice, GHC 7.6+still works just fine, and is still being tested against. It's just+that this support will be droped in the future, once it becomes a+burden.) Thus, this version of the library provides a transitional+point between backwards compatability and adding new features (see+below).++## Install++In general, this is a simple package and should be easy to install.+The specifics are a bit murky however, since we use CPP and the FFI+and want to retain portability, and due to the rapid evolution of+Cabal and other development tools. However, you should be able to+use the standard:++ $> cabal install logfloat+++### FFI Problems++The logfloat package uses the FFI to access functions in libm to+improve accuracy. These functions are part of the ISO C 1999 standard+and are supported on most systems, however they're not part of the+ANSI C standard nor the System V standard and so they might be+unavailable on some systems. If you encounter errors during build+about not being able to find `<math.h>`, not having prototypes for+functions, or linking errors, then you are probably on such a system.++In order to use pure Haskell and disable the FFI, pass -f-useFFI+during the configure phase to disable the useFFI flag. This will+degrade the accuracy of certain operations, but should still compile+cleanly.++If you run into errors like "error: 'log1p' redeclared as different+kind of symbol" this is due to a bug in GHC 6.10.1 where you can't+use both -fvia-C and the FFI. This bug should be circumvented for+the moment, but resolving it is still a work in progress. See the+bug report (which is resolved in GHC 6.10.2):++ http://hackage.haskell.org/trac/ghc/ticket/3117+++### Windows FFI++The logfloat package builds and installs cleanly with GHC 8.0 and+later on Windows (without needing Cygwin nor Mingw/Msys), and the+installed package can be used in compiled programs without any+issues.++Prior to GHC 8.0, the package does not work from the GHCi debugger+by default and gives errors like "can't load .so/.DLL for: m (addDLL:+could not load DLL)". This is a long-standing issue with GHCi having+to do with the fact that, on Windows, the so-called "libm" file does+not actually contain anything (as it does on POSIX) and the C+functions we use are instead placed in `libmingwex.a` (which comes+bundled with GHC). The problem is that `ghc` (the compiler) knows+to pull in `libmingwex.a`, whereas `ghci` (the interactive debugger)+does not. All of this is true at least as far back as Windows XP+and GHC 6.10.1. Some more information can be found in the ticket:++ https://ghc.haskell.org/trac/ghc/ticket/3242++The most reliable workaround at this point, alas, is to compile the+library with FFI support disabled in order to be able to use it in+GHCi, and then to recompile with FFI enabled whenever you need to+ship a compiled program.++A less-reliable workaround (i.e., needs to be tailored for your+system based on your paths, and you'll need to have `gcc` installed)+which avoids the need for recompilation is to generate the DLL+yourself by running the following two commands and then placing the+resulting `m.dll` into your path.++ ar -x libmingwex.a+ gcc -shared *.o -o m.dll+++### Testing++If you want to run the test suite, use the following standard method+(with `runhaskell Setup.hs` in lieu of `cabal`, if necessary):++ $> cd logfloat+ $> cabal configure --enable-tests --enable-coverage+ $> cabal build+ $> cabal test --keep-tix-files++The results of the code coverage are in+`./dist/hpc/vanilla/html/logfloat-$VERSION/hpc_index.html`. If+you're not interested in the coverage of the test suite, then you+needn't pass the `--enable-coverage` nor `--keep-tix-files` flags.+Note that older versions of cabal used the flag name+`--enable-library-coverage` instead of `--enable-coverage`. And+IIRC hpc integration in cabal was broken for ghc-7.6.+++### Haddock Problems++In Cabal 1.2 there is a bug in the handling of building Haddock+documentation when CPP is involved. These issues have been fixed+in Cabal 1.6, but here are the instructions if you're on older+systems.++In Cabal 1.2 the cpp-options field is not passed to Haddock, and+therefore any macros defined there are not seen, which can cause+Haddock to fail. The old workaround was to define CPP macros in the+ghc-options field which *does* get passed to Haddock. This is now+considered bad style and is forbidden by Hackage. It appears that+passing --haddock-option=... or --ghc-option=... flags during+configure does not have the same effect as defining the field.++Therefore, in order to properly compile Haddock documentation on+Cabal 1.2, you should go into logfloat.cabal and uncomment the+ghc-options fields which declare CPP macros -D__USE_FFI__ and+-D__HUGS__=200609. If you know of a better workaround for this+configuration, contact the maintainer.+++### Building for Hugs (September 2006)++If you're feeling adventurous and want to try and get the library+to compile under Hugs, here's how you used to be able to compile+things:++ runhaskell Setup.hs configure --hugs \+ --with-cpphs="`which cpphs-hugs`" \+ --ffihugs-options="-98 +o" \+ --ffihugs-option=-F'cpp -P -traditional -D__HUGS__=200609 -D__USE_FFI__'+ runhaskell Setup.hs build+ runhaskell Setup.hs copy+ runhaskell Setup.hs register++If you need to disable the FFI due to issues with not being able+to find <math.h>, not having prototypes for functions, or linking+errors, be sure *not* to pass -D__USE_FFI__ to the cpp filter for+ffihugs.++Notably, Hugs installs cpphs under the name "cpphs-hugs" by default.+The `which` command will try to resolve the location, assuming it's+on your $PATH. If it isn't, then change the --with-cpphs= flag to+point to where cpphs is installed on your system.++Some additional details about difficulties with building for Hugs+can be found in this blog post and the bug reports:++ http://winterkoninkje.livejournal.com/60707.html+ http://hackage.haskell.org/trac/hackage/ticket/526+ http://hackage.haskell.org/trac/hackage/ticket/527++And while previous versions of Cabal could build this package (with+the above commandline to work around those bugs), Hugs support is+broken in Cabal 1.8 (including at least 1.8.0.2 through 1.8.0.6):++ http://hackage.haskell.org/trac/hackage/ticket/633+++## Changes: Version 0.13.3+ (2015-03-29) vs 0.12.1 (2010-03-19)++* Monomorphized `logFloat`, `logToLogFloat`, `fromLogFloat`, and+`logFromLogFloat`: that is, they all take/return `Double` now. The+change was made to help reduce the need for explicit type signatures.+It shouldn't really affect most users, since it seems noone was+really making use of the polymorphism provided by previous versions.+To get the previous behavior back, just explicitly add calls to+`realToFrac` wherever necessary.++* Fixed some instances to get them to compile under the new role-based+type system of GHC 7.10++* Cleaned up various extraneous rewrite rules, specializations, etc++* Added the functions `sum`, `product`, and `pow`. Both sum and+product preserve more precision than the fold-based definitions in+the Prelude. Moreover, sum is *much* faster than the Prelude version,+since it only requires crossing the log/exp boundary n+1 times,+instead of 2\*(n-1) times. The only downside is that sum requires+two passes over the input and thus is not amenable to list fusion.++* (Version 0.13.3.2; 2015-08-06) Fixed the `Show LogFloat` instance+to produce parentheses in the right place.+++## Upcoming changes (0.14+)++* Since the `Data.Number.RealToFrac` module is no longer required+by any of the others, it will probably be forked off to a separate+package in order to improve portability of the rest of the package+by removing the need for MPTCs.++* There's long been clamoring for adding a+vector:`Data.Vector.Unboxed.Unbox` instance. I've been reluctant+to add such an instance due to wanting to retain backwards compatibility+and portability. Having dropped support for Hugs and older versions+of GHC, I'm now willing to add them in.++The logfloat library is conceptually quite simple, and thus to+whatever extent possible I'd still like to retain portability to+non-GHC compilers. So if you are interested in using logfloat with+another compiler/interpreter but run into problems (e.g., due to+the type families required by the vector library), please get in+touch and I'll try to get things to work.+++## Compatibility / Portability++The package is compatible with GHC 7.8.3 and 7.10.1. It may still+compile with older versions of GHC (or even Hugs!), however they+are no longer officially supported.++The package is not compatible with nhc98 and Yhc because+`Data.Number.RealToFrac` uses MPTCs. However, that module is no+longer required by any others, and all the other modules should be+compatible with these compilers. Thus, it should be fairly easy to+port. If you do so, please let me know and I'll try to incorporate+support for them.+++## Links++* [Website](http://wrengr.org/)+* [Blog](http://winterkoninkje.dreamwidth.org/)+* [Twitter](https://twitter.com/wrengr)+* [Hackage](http://hackage.haskell.org/package/logfloat)+* [GitHub](https://github.com/wrengr/logfloat)
logfloat.cabal view
@@ -1,62 +1,91 @@+Cabal-Version: 2.2+-- Cabal >=2.2 is required for:+-- <https://cabal.readthedocs.io/en/latest/cabal-package.html#common-stanzas>+-- Since 2.1, the Cabal-Version must be the absolutely first thing+-- in the file, even before comments. Also, no longer uses ">=".+-- <https://github.com/haskell/cabal/issues/4899>+ ------------------------------------------------------------------- wren ng thornton <wren@community.haskell.org> ~ 2010.03.19+-- wren gayle romano <wren@cpan.org> ~ 2026-02-26 ---------------------------------------------------------------- Name: logfloat-Version: 0.12.1-Cabal-Version: >= 1.2+Version: 0.14.0.2 Build-Type: Simple-Stability: experimental-Copyright: Copyright (c) 2007--2010 wren ng thornton-License: BSD3+Stability: provisional+Homepage: https://wrengr.org/software/hackage.html+Bug-Reports: https://github.com/wrengr/logfloat/issues+Author: wren gayle romano+Maintainer: wren@cpan.org+Copyright: 2007–2026 wren romano+-- Cabal-2.2 requires us to say "BSD-3-Clause" not "BSD3"+License: BSD-3-Clause License-File: LICENSE-Author: wren ng thornton-Maintainer: wren@community.haskell.org-Homepage: http://code.haskell.org/~wren/+ Category: Data, Math, Natural Language Processing, Statistics-Synopsis: Log-domain floating point numbers-Description: This module presents a type for storing numbers- in the log-domain. The main reason for doing this- is to prevent underflow when multiplying many- probabilities as is done in Hidden Markov Models.- It is also helpful for preventing overflow.+Synopsis:+ Log-domain floating point numbers+Description:+ This module presents a type for storing numbers in the log-domain.+ The main reason for doing this is to prevent underflow when+ multiplying many probabilities as is done in Hidden Markov+ Models. It is also helpful for preventing overflow. +Extra-source-files:+ README.md, CHANGELOG -Flag splitBase- Description: base-3.0 broke out array and other packages- Default: False+-- This should work as far back as GHC 7.6.1, but we don't verify that by CI.+-- <https://github.com/wrengr/logfloat/actions?query=workflow%3Aci>+Tested-With:+ GHC ==8.0.2,+ GHC ==8.2.2,+ GHC ==8.4.4,+ GHC ==8.6.5,+ GHC ==8.8.4,+ GHC ==8.10.3,+ GHC ==9.0.1,+ GHC ==9.2.4,+ GHC ==9.4.8,+ GHC ==9.6.5,+ GHC ==9.8.2,+ GHC ==9.10.1,+ GHC ==9.12.1,+ GHC ==9.14.1 +----------------------------------------------------------------+Source-Repository head+ Type: git+ Location: https://github.com/wrengr/logfloat.git++---------------------------------------------------------------- Flag useFFI- Description: Use FFI to link C's log1p. Improves accuracy and- same performance. Disable only if your compiler- doesn't support FFI. Default: True+ Description:+ Use FFI to link C's log1p. Improves accuracy while maintaining+ the same performance. Disable only if your compiler doesn't+ support FFI. +---------------------------------------------------------------- Library+ Default-Language: Haskell2010 Hs-Source-Dirs: src Exposed-Modules: Data.Number.LogFloat+ , Data.Number.LogFloat.Raw , Data.Number.RealToFrac , Data.Number.Transfinite , Data.Number.PartialOrd , Hugs.RealFloat- if flag(splitBase)- Build-depends: base >= 3.0, array- else- Build-depends: base < 3.0 + Build-Depends: base >= 4.6 && < 5+ -- TODO: do we need any version constraints on array?+ , array+ if flag(useFFI) -- BUG (Cabal 1.2 + Haddock): See the INSTALL file. --GHC-Options: -D__USE_FFI__ CPP-Options: -D__USE_FFI__ includes: math.h extra-libraries: m- - -- BUG (GHC 6.10.1): Can't mix FFI and -fvia-C- -- <http://hackage.haskell.org/trac/ghc/ticket/3117>- -- This is a momentary hack to keep it for 6.8- -- TODO: see if -fasm suffices for the benefits- if impl(ghc < 6.10)- GHC-Options: -fvia-C -optc-O3 -- BUG (Cabal <= 1.6 + Hugs): See the INSTALL file. Hugs-Options: -98 +o
src/Data/Number/LogFloat.hs view
@@ -1,32 +1,25 @@--- FlexibleContexts needed by our RealToFrac contexts--- CPP needed for IArray UArray instance+-- CPP and GeneralizedNewtypeDeriving are needed for IArray UArray instance -- FFI is for log1p--- N.B. can't mix FFI and -fvia-C under ghc==6.10.1--- <http://hackage.haskell.org/trac/ghc/ticket/3117>--- TODO: figure out how to conditionally use "-fvia-C -optc-O3"--- for ghc6.8 when not using Cabal to compile.-{-# LANGUAGE FlexibleContexts- , CPP- , ForeignFunctionInterface- #-}+{-# LANGUAGE CPP, ForeignFunctionInterface, MultiParamTypeClasses #-}+-- We don't put these in LANGUAGE, because it's CPP guarded for GHC only+-- HACK: ScopedTypeVariables and InstanceSigs are for GHC 7.10 only...+{-# OPTIONS_GHC+ -XGeneralizedNewtypeDeriving+ -XScopedTypeVariables+ -XInstanceSigs+ #-} --- Removed -Wall because -fno-warn-orphans was removed in GHC 6.10 {-# OPTIONS_GHC -Wall -fwarn-tabs #-} --- Unfortunately we need -fglasgow-exts in order to actually pick--- up on the rules (see -ddump-rules). The -frewrite-rules flag--- doesn't do what you want.--- cf <http://hackage.haskell.org/trac/ghc/ticket/2213>--- cf <http://www.mail-archive.com/glasgow-haskell-users@haskell.org/msg14313.html>-{-# OPTIONS_GHC -O2 -fexcess-precision -fglasgow-exts #-}+{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-} ------------------------------------------------------------------- ~ 2010.03.19+-- ~ 2021.10.17 -- | -- Module : Data.Number.LogFloat--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable (with CPP, FFI) --@@ -51,26 +44,28 @@ ( -- * Exceptional numeric values module Data.Number.Transfinite- , module Data.Number.RealToFrac- + -- * @LogFloat@ data type- , LogFloat+ , LogFloat() -- ** Isomorphism to normal-domain , logFloat , fromLogFloat -- ** Isomorphism to log-domain , logToLogFloat , logFromLogFloat- + -- ** Additional operations+ , sum, product+ , pow+ -- * Accurate versions of logarithm\/exponentiation , log1p, expm1 ) where -import Prelude hiding (log, realToFrac, isInfinite, isNaN)+import Prelude hiding (log, sum, product, isInfinite, isNaN) -import Data.Number.RealToFrac import Data.Number.Transfinite import Data.Number.PartialOrd+import Data.Number.LogFloat.Raw -- GHC can derive (IArray UArray LogFloat), but Hugs needs to coerce@@ -81,9 +76,15 @@ -- Hugs (Sept 2006) doesn't use the generic wrapper in base:Unsafe.Coerce -- so we'll just have to go back to the original source. #ifdef __HUGS__-import Hugs.IOExts (unsafeCoerce)+import Hugs.IOExts (unsafeCoerce) #elif __NHC__ import NonStdUnsafeCoerce (unsafeCoerce)+#elif __GLASGOW_HASKELL__ >= 710+-- For when the *heap* representations are the same+--import Data.Coerce (coerce)+-- For when the *unboxed array* storage representations are the same+import Unsafe.Coerce (unsafeCoerce)+import Data.Ix (Ix) #endif #ifdef __GLASGOW_HASKELL__@@ -95,9 +96,27 @@ -- | A @LogFloat@ is just a 'Double' with a special interpretation. -- The 'logFloat' function is presented instead of the constructor, -- in order to ensure semantic conversion. At present the 'Show'--- instance will convert back to the normal-domain, and so will+-- instance will convert back to the normal-domain, and hence will -- underflow at that point. This behavior may change in the future.+-- At present, the 'Read' instance parses things in the normal-domain+-- and then converts them to the log-domain. Again, this behavior+-- may change in the future. --+-- Because 'logFloat' performs the semantic conversion, we can use+-- operators which say what we /mean/ rather than saying what we're+-- actually doing to the underlying representation. That is,+-- equivalences like the following are true[1] thanks to type-class+-- overloading:+--+-- > logFloat (p + q) == logFloat p + logFloat q+-- > logFloat (p * q) == logFloat p * logFloat q+--+-- (Do note, however, that subtraction can and negation will throw+-- errors: since @LogFloat@ can only represent the positive half of+-- 'Double'. 'Num' is the wrong abstraction to put at the bottom+-- of the numeric type-class hierarchy; but alas, we're stuck with+-- it.)+-- -- Performing operations in the log-domain is cheap, prevents -- underflow, and is otherwise very nice for dealing with miniscule -- probabilities. However, crossing into and out of the log-domain@@ -105,56 +124,128 @@ -- particular, if you're doing a series of multiplications as in -- @lp * logFloat q * logFloat r@ it's faster to do @lp * logFloat -- (q * r)@ if you're reasonably sure the normal-domain multiplication--- won't underflow, because that way you enter the log-domain only--- once, instead of twice.+-- won't underflow; because that way you enter the log-domain only+-- once, instead of twice. Also note that, for precision, if you're+-- doing more than a few multiplications in the log-domain, you+-- should use 'product' rather than using '(*)' repeatedly. -- -- Even more particularly, you should /avoid addition/ whenever--- possible. Addition is provided because it's necessary at times--- and the proper implementation is not immediately transparent.--- However, between two @LogFloat@s addition requires crossing the--- exp\/log boundary twice; with a @LogFloat@ and a regular number--- it's three times since the regular number needs to enter the--- log-domain first. This makes addition incredibly slow. Again,--- if you can parenthesize to do plain operations first, do it!+-- possible. Addition is provided because sometimes we need it, and+-- the proper implementation is not immediately apparent. However,+-- between two @LogFloat@s addition requires crossing the exp\/log+-- boundary twice; with a @LogFloat@ and a 'Double' it's three+-- times, since the regular number needs to enter the log-domain+-- first. This makes addition incredibly slow. Again, if you can+-- parenthesize to do normal-domain operations first, do it!+--+-- [1] That is, true up-to underflow and floating point fuzziness.+-- Which is, of course, the whole point of this module. newtype LogFloat = LogFloat Double deriving ( Eq , Ord -- Should we really perpetuate the Ord lie? #ifdef __GLASGOW_HASKELL__- -- At least GHC 6.8.2 can derive these (without- -- GeneralizedNewtypeDeriving). The H98 Report doesn't include- -- them among the options for automatic derivation though.+ -- The H98 Report doesn't include these among the options for+ -- automatic derivation. But GHC >= 6.8.2 (at least) can derive+ -- them (even without GeneralizedNewtypeDeriving). However,+ -- GHC >= 7.10 can't derive @IArray UArray@ thanks to the new+ -- type-role stuff! since 'UArray' is declared to be nominal+ -- in both arguments... and that seems to be necessary:+ -- cf., <https://ghc.haskell.org/trac/ghc/ticket/9220>+#if __GLASGOW_HASKELL__ < 710 , IArray UArray+#endif , Storable #endif ) +#if __GLASGOW_HASKELL__ >= 710+-- TODO: this version should also work for NHC and Hugs, I think...+-- HACK: we should be able to just unsafeCoerce the functions+-- themselves, instead of coercing the inputs and the outputs; but,+-- GHC 7.10 seems to get confused about trying to coerce the index+-- types too... To fix this we give explicit signatures, as below,+-- but this requires both ScopedTypeVariables and InstanceSigs; and+-- I'm not sure when InstanceSigs was introduced. -#if __HUGS__ || __NHC__+instance IArray UArray LogFloat where+ {-# INLINE bounds #-}+ bounds :: forall i. Ix i => UArray i LogFloat -> (i, i)+ bounds = unsafeCoerce (bounds :: UArray i Double -> (i, i))++ {-# INLINE numElements #-}+ numElements :: forall i. Ix i => UArray i LogFloat -> Int+ numElements = unsafeCoerce (numElements :: UArray i Double -> Int)++ {-# INLINE unsafeArray #-}+ unsafeArray+ :: forall i. Ix i => (i,i) -> [(Int,LogFloat)] -> UArray i LogFloat+ unsafeArray = unsafeCoerce (unsafeArray+ :: (i,i) -> [(Int,Double)] -> UArray i Double)++ {-# INLINE unsafeAt #-}+ unsafeAt :: forall i. Ix i => UArray i LogFloat -> Int -> LogFloat+ unsafeAt = unsafeCoerce (unsafeAt :: UArray i Double -> Int -> Double)++ {-# INLINE unsafeReplace #-}+ unsafeReplace+ :: forall i. Ix i+ => UArray i LogFloat -> [(Int,LogFloat)] -> UArray i LogFloat+ unsafeReplace = unsafeCoerce (unsafeReplace+ :: UArray i Double -> [(Int,Double)] -> UArray i Double)++ {-# INLINE unsafeAccum #-}+ unsafeAccum+ :: forall i e. Ix i+ => (LogFloat -> e -> LogFloat)+ -> UArray i LogFloat -> [(Int,e)] -> UArray i LogFloat+ unsafeAccum = unsafeCoerce (unsafeAccum+ :: (Double -> e -> Double)+ -> UArray i Double -> [(Int,e)] -> UArray i Double)++ {-# INLINE unsafeAccumArray #-}+ unsafeAccumArray+ :: forall i e. Ix i+ => (LogFloat -> e -> LogFloat)+ -> LogFloat -> (i,i) -> [(Int,e)] -> UArray i LogFloat+ unsafeAccumArray = unsafeCoerce (unsafeAccumArray+ :: (Double -> e -> Double)+ -> Double -> (i,i) -> [(Int,e)] -> UArray i Double)++#elif __HUGS__ || __NHC__ -- TODO: Storable instance. Though Foreign.Storable isn't in Hugs(Sept06) +-- TODO: depend on my @pointless-fun@ package rather than repeating things here... -- These two operators make it much easier to read the instance. -- Hopefully inlining everything will get rid of the eta overhead. -- <http://matt.immute.net/content/pointless-fun>+(~>) :: (a -> b) -> (d -> c) -> (b -> d) -> a -> c {-# INLINE (~>) #-} infixr 2 ~> f ~> g = (. f) . (g .) -{-# INLINE ($.) #-}-infixl 1 $.-($.) = flip ($)+($::) :: a -> (a -> b) -> b+{-# INLINE ($::) #-}+infixl 1 $::+($::) = flip ($) {-# INLINE logFromLFAssocs #-} logFromLFAssocs :: [(Int, LogFloat)] -> [(Int, Double)]+#if __GLASGOW_HASKELL__ >= 710+logFromLFAssocs = coerce+#else logFromLFAssocs = unsafeCoerce+#endif +-- HACK: can't coerce, cf: <https://ghc.haskell.org/trac/ghc/ticket/9220> {-# INLINE logFromLFUArray #-} logFromLFUArray :: UArray a LogFloat -> UArray a Double logFromLFUArray = unsafeCoerce -- Named unsafe because it could allow injecting NaN if misused+-- HACK: can't coerce, cf: <https://ghc.haskell.org/trac/ghc/ticket/9220> {-# INLINE unsafeLogToLFUArray #-} unsafeLogToLFUArray :: UArray a Double -> UArray a LogFloat unsafeLogToLFUArray = unsafeCoerce@@ -162,9 +253,9 @@ -- Named unsafe because it could allow injecting NaN if misused {-# INLINE unsafeLogToLFFunc #-} unsafeLogToLFFunc :: (LogFloat -> a -> LogFloat) -> (Double -> a -> Double)-unsafeLogToLFFunc = ($. unsafeLogToLogFloat ~> id ~> logFromLogFloat)+unsafeLogToLFFunc = ($:: unsafeLogToLogFloat ~> id ~> logFromLogFloat) --- | Remove the extranious 'isNaN' test of 'logToLogFloat', when+-- | Remove the extraneous 'isNaN' test of 'logToLogFloat', when -- we know it's safe. {-# INLINE unsafeLogToLogFloat #-} unsafeLogToLogFloat :: Double -> LogFloat@@ -174,124 +265,118 @@ instance IArray UArray LogFloat where {-# INLINE bounds #-} bounds = bounds . logFromLFUArray- + -- Apparently this method was added in base-2.0/GHC-6.6 but Hugs -- (Sept 2006) doesn't have it. Not sure about NHC's base-#if __HUGS__ > 200609+#if (!(defined(__HUGS__))) || (__HUGS__ > 200609) {-# INLINE numElements #-} numElements = numElements . logFromLFUArray #endif- + {-# INLINE unsafeArray #-}- unsafeArray =- unsafeArray $. id ~> logFromLFAssocs ~> unsafeLogToLFUArray- + unsafeArray = unsafeArray $:: id ~> logFromLFAssocs ~> unsafeLogToLFUArray+ {-# INLINE unsafeAt #-}- unsafeAt =- unsafeAt $. logFromLFUArray ~> id ~> unsafeLogToLogFloat- + unsafeAt = unsafeAt $:: logFromLFUArray ~> id ~> unsafeLogToLogFloat+ {-# INLINE unsafeReplace #-}- unsafeReplace =- unsafeReplace $. logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray- + unsafeReplace = unsafeReplace+ $:: logFromLFUArray ~> logFromLFAssocs ~> unsafeLogToLFUArray+ {-# INLINE unsafeAccum #-}- unsafeAccum =- unsafeAccum $. unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray- + unsafeAccum = unsafeAccum+ $:: unsafeLogToLFFunc ~> logFromLFUArray ~> id ~> unsafeLogToLFUArray+ {-# INLINE unsafeAccumArray #-}- unsafeAccumArray =- unsafeAccumArray $. unsafeLogToLFFunc ~> logFromLogFloat ~> id ~> id ~> unsafeLogToLFUArray+ unsafeAccumArray = unsafeAccumArray+ $:: unsafeLogToLFFunc ~> logFromLogFloat ~> id ~> id ~> unsafeLogToLFUArray #endif -- TODO: the Nothing branch should never be reachable. Once we get--- a test suite up and going to *verify* the never-NaN invariant,+-- a test suite up and going to /verify/ the never-NaN invariant, -- we should be able to eliminate the branch and the isNaN checks. instance PartialOrd LogFloat where- cmp (LogFloat x) (LogFloat y) + cmp (LogFloat x) (LogFloat y) | isNaN x || isNaN y = Nothing | otherwise = Just $! x `compare` y +instance Read LogFloat where+ readsPrec p s =+ [(LogFloat (log x), r) | (x, r) <- readsPrec p s, not (isNaN x), x >= 0] ---------------------------------------------------------------- -- | Reduce the number of constant string literals we need to store. errorOutOfRange :: String -> a {-# NOINLINE errorOutOfRange #-}-errorOutOfRange fun = error $! "Data.Number.LogFloat."++fun- ++ ": argument out of range"+errorOutOfRange fun =+ error $! "Data.Number.LogFloat."++fun++ ": argument out of range" -- Both guards are redundant due to the subsequent call to -- 'Data.Number.Transfinite.log' at all use sites. However we use -- this function to give local error messages. Perhaps we should -- catch the exception and throw the new message instead? Portability? -guardNonNegative :: String -> Double -> Double-guardNonNegative fun x | isNaN x || x < 0 = errorOutOfRange fun- | otherwise = x+guardNonNegative :: String -> Double -> Double+guardNonNegative fun x+ | isNaN x || x < 0 = errorOutOfRange fun+ | otherwise = x -guardIsANumber :: String -> Double -> Double-guardIsANumber fun x | isNaN x = errorOutOfRange fun- | otherwise = x+guardIsANumber :: String -> Double -> Double+guardIsANumber fun x+ | isNaN x = errorOutOfRange fun+ | otherwise = x ---------------------------------------------------------------- -- | Constructor which does semantic conversion from normal-domain--- to log-domain. Throws errors on negative input.-logFloat :: (Real a, RealToFrac a Double) => a -> LogFloat-{-# SPECIALIZE logFloat :: Double -> LogFloat #-}-logFloat = LogFloat . log . guardNonNegative "logFloat" . realToFrac+-- to log-domain. Throws errors on negative and NaN inputs. If @p@+-- is non-negative, then following equivalence holds:+--+-- > logFloat p == logToLogFloat (log p)+--+-- If @p@ is NaN or negative, then the two sides differ only in+-- which error is thrown.+logFloat :: Double -> LogFloat+{-# INLINE [0] logFloat #-}+-- TODO: should we use NOINLINE or [~0] to avoid the possibility of code bloat?+logFloat = LogFloat . log . guardNonNegative "logFloat" --- This is simply a polymorphic version of the 'LogFloat' data--- constructor. We present it mainly because we hide the constructor--- in order to make the type a bit more opaque. If the polymorphism--- turns out to be a performance liability because the rewrite rules--- can't remove it, then we need to rethink all four--- constructors\/destructors.--- -- | Constructor which assumes the argument is already in the--- log-domain. Throws errors on @notANumber@ input.-logToLogFloat :: (Real a, RealToFrac a Double) => a -> LogFloat-{-# SPECIALIZE logToLogFloat :: Double -> LogFloat #-}-logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" . realToFrac+-- log-domain. Throws errors on @notANumber@ inputs.+logToLogFloat :: Double -> LogFloat+logToLogFloat = LogFloat . guardIsANumber "logToLogFloat" --- | Return our log-domain value back into normal-domain. Beware--- of overflow\/underflow.-fromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a)- => LogFloat -> a-{-# SPECIALIZE fromLogFloat :: LogFloat -> Double #-}-fromLogFloat (LogFloat x) = realToFrac (exp x)+-- | Semantically convert our log-domain value back into the+-- normal-domain. Beware of overflow\/underflow. The following+-- equivalence holds (without qualification):+--+-- > fromLogFloat == exp . logFromLogFloat+--+fromLogFloat :: LogFloat -> Double+{-# INLINE [0] fromLogFloat #-}+-- TODO: should we use NOINLINE or [~0] to avoid the possibility of code bloat?+fromLogFloat (LogFloat x) = exp x -- | Return the log-domain value itself without conversion.-logFromLogFloat :: (Fractional a, Transfinite a, RealToFrac Double a)- => LogFloat -> a-{-# SPECIALIZE logFromLogFloat :: LogFloat -> Double #-}-logFromLogFloat (LogFloat x) = realToFrac x+logFromLogFloat :: LogFloat -> Double+logFromLogFloat (LogFloat x) = x -- These are our module-specific versions of "log\/exp" and "exp\/log"; -- They do the same things but also have a @LogFloat@ in between--- the logarithm and exponentiation.------ In order to ensure these rules fire we may need to delay inlining--- of the four con-\/destructors, like we do for 'realToFrac'.--- Unfortunately, I'm not entirely sure whether they will be inlined--- already or not (and whether they are may be fragile) and I don't--- want to inline them excessively and lead to code bloat in the--- off chance that we could prune some of it away.--- TODO: thoroughly investigate this.+-- the logarithm and exponentiation. In order to ensure these rules+-- fire, we have to delay the inlining on two of the four+-- con-\/destructors. {-# RULES -- Out of log-domain and back in "log/fromLogFloat" forall x. log (fromLogFloat x) = logFromLogFloat x-"log.fromLogFloat" log . fromLogFloat = logFromLogFloat- "logFloat/fromLogFloat" forall x. logFloat (fromLogFloat x) = x-"logFloat.fromLogFloat" logFloat . fromLogFloat = id -- Into log-domain and back out "fromLogFloat/logFloat" forall x. fromLogFloat (logFloat x) = x-"fromLogFloat.logFloat" fromLogFloat . logFloat = id #-} ----------------------------------------------------------------@@ -305,61 +390,43 @@ -- log-domain value instead. instance Show LogFloat where- show (LogFloat x) = let y = exp x- in y `seq` "LogFloat "++show y+ showsPrec p (LogFloat x) =+ let y = exp x in y `seq`+ showParen (p > 9)+ ( showString "logFloat "+ . showsPrec 11 y+ ) ------------------------------------------------------------------- Technically these should use 'Foreign.C.CDouble' however there's--- no isomorphism provided to normal 'Double'. The former is--- documented as being a newtype of the later, and so this should--- be safe.--#ifdef __USE_FFI__-#define LOG1P_WHICH_VERSION specialized version.-#else-#define LOG1P_WHICH_VERSION naive version! \- Contact the maintainer with any FFI difficulties.-#endif----- | Definition: @log1p == log . (1+)@. The C language provides a--- special definition for 'log1p' which is more accurate than doing--- the naive thing, especially for very small arguments. For example,--- the naive version underflows around @2 ** -53@, whereas the--- specialized version underflows around @2 ** -1074@. This function--- is used by ('+') and ('-') on @LogFloat@.+-- | A curried function for converting arbitrary pairs into ordered+-- pairs. The continuation recieves the minimum first and the maximum+-- second. ----- /This installation was compiled to use the LOG1P_WHICH_VERSION/--#ifdef __USE_FFI__-foreign import ccall unsafe "math.h log1p"- log1p :: Double -> Double-#else-log1p :: Double -> Double-{-# INLINE log1p #-}-log1p x = log (1 + x)-#endif+-- This combinator is primarily intended to reduce repetition in+-- the source code; but hopefully it should also help reduce bloat+-- in the compiled code, by sharing the continuation and just+-- swapping the variables in place. Of course, if the continuation+-- is very small, then requiring a join point after the conditional+-- swap may end up being more expensive than simply duplicating the+-- continuation. Also, given as we're inlining it, I'm not sure+-- whether GHC will decide to keep the sharing we introduced or+-- whether it'll end up duplicating the continuation into the two+-- call sites.+ordered :: Ord a => a -> a -> (a -> a -> b) -> b+ordered x y k+ | x <= y = k x y+ | otherwise = k y x+ -- N.B., the implementation of @(>=)@ in Hugs (Sept2006) will+ -- always returns True if either argument isNaN. This does not+ -- constitute a bug for us, since we maintain the invariant that+ -- values wrapped by 'LogFloat' are not NaN.+{-# INLINE ordered #-} --- | Definition: @expm1 == (subtract 1) . exp@. The C language--- provides a special definition for 'expm1' which is more accurate--- than doing the naive thing, especially for very small arguments.--- This function isn't needed internally, but is provided for--- symmetry with 'log1p'.+-- TODO: Do we need to add explicit INLINE pragmas here? Or will+-- GHC automatically see that they're small enough to want inlining? ----- /This installation was compiled to use the LOG1P_WHICH_VERSION/--#ifdef __USE_FFI__-foreign import ccall unsafe "math.h expm1"- expm1 :: Double -> Double-#else-expm1 :: Double -> Double-{-# INLINE expm1 #-}-expm1 x = exp x - 1-#endif------------------------------------------------------------------ -- These all work without causing underflow. However, do note that -- they tend to induce more of the floating-point fuzz than using -- regular floating numbers because @exp . log@ doesn't really equal@@ -367,33 +434,33 @@ -- multiplying many small numbers (and preventing overflow for -- multiplying many large numbers) so we're not too worried about -- +\/- 4e-16.--instance Num LogFloat where - -- N.B. In Hugs (Sept2006) the (>=) always returns True if- -- either isNaN. This does not constitute a bug since we- -- maintain the invariant that values wrapped by 'LogFloat'- -- are not NaN.- +instance Num LogFloat where (*) (LogFloat x) (LogFloat y)- | isInfinite x- && isInfinite y- && x == negate y = LogFloat negativeInfinity -- @0*infinity == 0@- | otherwise = LogFloat (x+y)- + | isInfinite x && isInfinite y && x == negate y =+ LogFloat negativeInfinity -- @0 * infinity == 0@+ | otherwise =+ -- This includes the @0 * 0 == 0@ and @infty * infty == infty@+ -- cases, since @(+)@ treats them appropriately.+ LogFloat (x + y)+ (+) (LogFloat x) (LogFloat y)- | x == y- && isInfinite x- && isInfinite y = LogFloat x -- @0+0 == 0@, @infinity+infinity == infinity@- | x >= y = LogFloat (x + log1p (exp (y - x)))- | otherwise = LogFloat (y + log1p (exp (x - y)))- + | isInfinite x && isInfinite y && x == y =+ LogFloat x -- @0 + 0 == 0@ and @infty + infty == infty@+ | otherwise =+ -- This includes the @0 + infinity == infinity@ case,+ -- since 'log1pexp' (and 'ordered') treats them appropriately.+ ordered x y $ \n m ->+ LogFloat (m + log1pexp (n - m))++ -- TODO: give a better error message in the (infinity,infinity) case.+ -- TODO: does 'log1mexp' handle the (+infty,-infty) cases correctly? (-) (LogFloat x) (LogFloat y)- | x == negativeInfinity- && y == negativeInfinity = LogFloat negativeInfinity -- @0-0 == 0@+ | x == negativeInfinity && y == negativeInfinity =+ LogFloat negativeInfinity -- @0 - 0 == 0@ | otherwise =- -- Will throw error if x < y or if (x,y) is (infinity,infinity)- LogFloat (guardIsANumber "(-)" (x + log1p (negate (exp (y - x)))))- + ordered x y $ \n m ->+ LogFloat (guardIsANumber "(-)" (m + log1mexp (n - m)))+ signum (LogFloat x) | x == negativeInfinity = 0 | x > negativeInfinity = 1@@ -402,29 +469,26 @@ -- broke the invariant. That shouldn't be possible and -- so noone else bothers to check, but we check here just -- in case.- + -- TODO: wouldn't @not (isNaN x)@ be a better guard to use?+ negate _ = errorOutOfRange "negate"- abs = id- - fromInteger = LogFloat . log- . guardNonNegative "fromInteger" . fromInteger+ fromInteger = LogFloat . log . guardNonNegative "fromInteger" . fromInteger instance Fractional LogFloat where- -- n/0 == infinity is handled seamlessly for us. We must catch 0/0 and infinity/infinity NaNs, and handle 0/infinity.+ -- @n / 0 == infinity@ is handled seamlessly for us. We must catch+ -- @0 / 0@ and @infinity / infinity@ NaNs, and handle @0 / infinity@. (/) (LogFloat x) (LogFloat y)- | x == y- && isInfinite x- && isInfinite y = errorOutOfRange "(/)"+ | isInfinite x && isInfinite y && x == y = errorOutOfRange "(/)" | x == negativeInfinity = LogFloat negativeInfinity -- @0/infinity == 0@- | otherwise = LogFloat (x-y)- + | otherwise = LogFloat (x - y)+ fromRational = LogFloat . log . guardNonNegative "fromRational" . fromRational --- Just for fun. The more coersion functions the better. Though+-- Just for fun. The more coercion functions the better. Though -- Rationals are very buggy when it comes to transfinite values instance Real LogFloat where toRational (LogFloat x)@@ -434,18 +498,71 @@ ex = exp x -{- -- Commented out because I'm not sure about requiring MPTCs. Of course, those are already required by "Data.Number.Transfinite" so it's pretty moot...+----------------------------------------------------------------+-- | /O(1)/. Compute powers in the log-domain; that is, the following+-- equivalence holds (modulo underflow and all that):+--+-- > logFloat (p ** m) == logFloat p `pow` m+--+-- /Since: 0.13/+pow :: LogFloat -> Double -> LogFloat+{-# INLINE pow #-}+infixr 8 `pow`+pow (LogFloat x) m+ | isNaN mx = LogFloat 0+ | otherwise = LogFloat mx+ where+ -- N.B., will be NaN when @x == negativeInfinity && m == 0@+ -- (which is true when m is -0 as well as +0). We check for NaN+ -- after multiplying, rather than checking this precondition+ -- before multiplying, in an attempt to simplify/optimize the+ -- generated code.+ -- TODO: benchmark.+ mx = m * x --- LogFloat->LogFloat is already given via generic (a->a)--- No LogFloat->Rational since LogFloat can have 'infinity'--- Can't have LogFloat->a using fromLogFloat because Hugs dislikes incoherence. Adding an explicit LogFloat->LogFloat instance doesn't help like it does for GHC. -instance RealToFrac LogFloat Double where- realToFrac = fromLogFloat- -instance RealToFrac LogFloat Float where- realToFrac = fromLogFloat--}+-- N.B., the default implementation of (**) for Complex is wrong.+-- It can be fixed by using the definition:+-- > x ** y = if x == 0 then 0 else exp (log x * y)+-- cf., <https://ghc.haskell.org/trac/ghc/ticket/8539>+-- TODO: Is this relevant to us at all?+++-- TODO: check out ekmett's compensated library.+++-- Some good test cases:+-- for @logsumexp == log . sum . map exp@:+-- logsumexp[0,1,0] should be about 1.55+-- for correctness of avoiding underflow:+-- logsumexp[1000,1001,1000] ~~ 1001.55 == 1000 + 1.55+-- logsumexp[-1000,-999,-1000] ~~ -998.45 == -1000 + 1.55+--+-- | /O(n)/. Compute the sum of a finite list of 'LogFloat's, being+-- careful to avoid underflow issues. That is, the following+-- equivalence holds (modulo underflow and all that):+--+-- > logFloat . sum == sum . map logFloat+--+-- /N.B./, this function requires two passes over the input. Thus,+-- it is not amenable to list fusion, and hence will use a lot of+-- memory when summing long lists.+--+-- /Since: 0.13/+sum :: [LogFloat] -> LogFloat+sum = LogFloat . logSumExp . fmap logFromLogFloat+++-- | /O(n)/. Compute the product of a finite list of 'LogFloat's,+-- being careful to avoid numerical error due to loss of precision.+-- That is, the following equivalence holds (modulo underflow and+-- all that):+--+-- > logFloat . product == product . map logFloat+--+-- /Since: 0.13/+product :: [LogFloat] -> LogFloat+product = LogFloat . kahanSum . fmap logFromLogFloat ---------------------------------------------------------------- ----------------------------------------------------------- fin.
+ src/Data/Number/LogFloat/Raw.hs view
@@ -0,0 +1,405 @@+{-# LANGUAGE CPP, ForeignFunctionInterface, BangPatterns #-}+{-# OPTIONS_GHC -Wall -fwarn-tabs #-}+{-# OPTIONS_GHC -O2 -fexcess-precision -fenable-rewrite-rules #-}++----------------------------------------------------------------+-- ~ 2026-02-28+-- |+-- Module : Data.Number.LogFloat.Raw+-- Copyright : 2007--2026 wren romano+-- License : BSD-3-Clause+-- Maintainer : wren@cpan.org+-- Stability : provisional+-- Portability : portable (with CPP, FFI)+--+-- This module provides implementations for computing various+-- logarithmic and exponential functions without losing precision+-- (as the naive implementations do). These are the \"raw\"+-- implementations; i.e., sans newtypes and other conveniences.+-- Since the lack of newtypes means we can't rely on types to clarify+-- things, we use the traditional baroque names for things. The+-- design considerations behind (most of) these implementations are+-- documented at:+-- <https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf>+--+-- In base-4.9.0.0 GHC added some of these to the 'Floating' class+-- exported from "Numeric". Alas, they provide default definitions+-- using the naive implementations, so one can't really rely on the+-- 'Floating' methods being precision preserving. Overall, the+-- specific instance for 'Double' looks fine (though they use+-- different cutoffs for 'log1pexp' for some reason); but it's easy+-- enough to reimplement here, to make absolutely sure we're getting+-- the right thing.+--+-- @since: 0.14.0+----------------------------------------------------------------+module Data.Number.LogFloat.Raw+ (+ -- * Logarithmic\/exponential basics+ expm1+ , log1p+ , log1mexp+ , log1pexp+ -- * Summation+ , logSumExp+ , kahanSum+ -- , neumaierSum+ -- * Softmax+ , logSoftmax+ , softmax+ -- * Sigmoid and related functions+ , sigmoid+ , logit+ , logitExp+ ) where++#if __GLASGOW_HASKELL__ < 910+import Data.List (foldl')+#endif+import Data.Number.Transfinite (negativeInfinity)++----------------------------------------------------------------+-- Technically these should use 'Foreign.C.CDouble' however there's+-- no isomorphism provided to normal 'Double'. The former is+-- documented as being a newtype of the later, and so this should+-- be safe.++#ifdef __USE_FFI__+#define LOG1P_WHICH_VERSION FFI version.+#else+#define LOG1P_WHICH_VERSION naive version! \+ Contact the maintainer with any FFI difficulties.+#endif+++-- | Compute @log (1 + x)@ without losing precision.+--+-- Standard C libraries provide a special definition for this+-- function, which is more accurate than doing the naive thing,+-- especially for very small arguments. For example, the naive+-- version underflows around @2 ** -53@, whereas the specialized+-- version underflows around @2 ** -1074@.+--+-- N.B. The @statistics:Statistics.Math@ module provides a pure+-- Haskell implementation of @log1p@ for those who are interested.+-- We do not copy it here because it relies on the @vector@ package+-- which is non-portable. If there is sufficient interest, a portable+-- variant of that implementation could be made. Contact the+-- maintainer if the FFI and naive implementations are insufficient+-- for your needs.+--+-- /This installation was compiled to use the LOG1P_WHICH_VERSION/++#ifdef __USE_FFI__+-- TODO: verify that the Haddock comes out as intended...+foreign import ccall unsafe "math.h log1p"+ log1p+ :: Double -- ^ N.B., only defined on the @[-1,infty]@ interval.+ -> Double+#else+-- See @statistics@:"Statistics.Math" for a more accurate Haskell+-- implementation.+log1p+ :: Double -- ^ N.B., only defined on the @[-1,infty]@ interval.+ -> Double+{-# INLINE [0] log1p #-}+log1p x = log (1 + x)+#endif+++-- | Compute @exp x - 1@ without losing precision.+--+-- Standard C libraries provide a special definition for 'expm1'+-- which is more accurate than doing the naive thing, especially+-- for very small arguments.+--+-- /This installation was compiled to use the LOG1P_WHICH_VERSION/++#ifdef __USE_FFI__+foreign import ccall unsafe "math.h expm1"+ expm1 :: Double -> Double+#else+expm1 :: Double -> Double+{-# INLINE [0] expm1 #-}+expm1 x = exp x - 1+#endif+++-- CPP guarded because they won't fire if we're using the FFI versions.+-- TODO: can we get them to fire if we to the standard thing about+-- naming the FFI version @c_foo@ and then defining a Haskell+-- function @foo = c_foo@?+#if !defined(__USE_FFI__)+{-# RULES+-- Into log-domain and back out+"expm1/log1p" forall x. expm1 (log1p x) = x++-- Out of log-domain and back in+"log1p/expm1" forall x. log1p (expm1 x) = x+ #-}+#endif+++-- | Compute @log (1 - exp x)@ without losing precision.+log1mexp+ :: Double -- ^ N.B., only defined on the @[-infty,0]@ interval.+ -> Double+log1mexp x+ | x <= log 2 = (log . negate . expm1) x+ | otherwise = (log1p . negate . exp) x+{-# INLINE log1mexp #-}+++-- | Compute @log (1 + exp x)@ without losing precision. Algebraically+-- this is @0 ⊔ x@, which is the log-domain's analogue of @1 + x@.+log1pexp :: Double -> Double+log1pexp x+ | x <= -37 = exp x+ | x <= 18 = log1p (exp x)+ | x <= 33.3 = x + exp (negate x)+ | otherwise = x+{-# INLINE log1pexp #-}+++-- TODO: bring back 'expm1c' and 'log1pc'+++----------------------------------------------------------------+-- | The logistic function; aka, the inverse of 'logit'.+-- > sigmoid x = 1 / (1 + exp (-x))+-- > sigmoid x = exp x / (exp x + 1)+-- > sigmoid x = (1 + tanh (x/2)) / 2+sigmoid :: Double -> Double+sigmoid x = (1 + tanh (x/2)) / 2+{-# INLINE sigmoid #-}+-- We prefer the 'tanh'-based definition because it's (exactly!)+-- symmetric about zero, whereas the naive version isn't (due to+-- floating-point fuzz).+-- TODO(b/68203642): Properly analyze the accuracy and precision+-- of the 'tanh' version.+++-- | The quantile function; aka, the inverse of 'sigmoid'.+-- > logit x = log (x / (1 - x))+-- > logit x = 2 * atanh (2*x - 1)+logit+ :: Double -- ^ N.B., only defined on the @[0,1]@ interval.+ -> Double+logit x = 2 * atanh (2*x - 1)+{-# INLINE logit #-}+-- TODO(b/68203642): properly analyze the precision of the 'atanh' version.+++-- | A variant of 'logit' for when the argument is already in the+-- log-domain; hence, @logitExp = logit . exp@+logitExp+ :: Double -- ^ N.B., only defined on the @[-infty,0]@ interval.+ -> Double+logitExp x = x - log1mexp x+{-# INLINE logitExp #-}+-- TODO(b/68203642): properly analyze the precision of this+-- implementation with respect to the 'logit' implementation.+++----------------------------------------------------------------+-- TODO: double check that everything inlines away, so this data+-- type doesn't introduce any slowdown.+--+-- | A helper type for 'logSumExp'. As a semigroup, this is isomorphic to:+-- @(WrappedMonoid (Sum Int), Max Double)@; however, we strictify and+-- flatten everything to improve performance.+data LSE = LSE+ {-# UNPACK #-}!Int -- The length, minus one.+ {-# UNPACK #-}!Double -- The maximum.++-- | Compute the length and maximum of a list. This is a semigroup+-- reduction. However we roll it ourselves rather than using the+-- semigroup class: since that would incur an otherwise unnecessary+-- dependency on @base >= 4.9.0.0@.+foldLSE :: Double -> [Double] -> LSE+foldLSE = foldl' step . LSE 0+ where+ step (LSE lm1 m) x = LSE (lm1 + 1) (m `max` x)+++-- TODO: expose a single-pass version for the special case where+-- the first element of the list is (promised to be) the maximum+-- element?+--+-- | /O(n)/. Log-domain summation, aka: @(log . sum . fmap exp)@.+-- Algebraically this is @⨆ xs@, which is the log-domain equivalent+-- of @∑ xs@.+--+-- /N.B./, this function requires two passes over the input. Thus,+-- it is not amenable to list fusion, and hence will use a lot of+-- memory when summing long lists.+logSumExp :: [Double] -> Double+logSumExp [] = (-1)/0+logSumExp xs0@(x:xs) =+ case foldLSE x xs of+ LSE lm1 m+ | isInfinite m -> m+ | otherwise ->+ -- TODO: push the addition of @lm1@ into the 'kahanSum',+ -- but making sure to add it in only at the very end.+ -- TODO: would using 'neumaierSum' be better? Should+ -- we factor the summation function out as an argument?+ -- TODO: is using 'log1p' here /really/ any better than+ -- just using 'log'?+ -- TODO: does that 'fmap' properly fuse into the+ -- 'kahanSum', or need we inline it ourselves?+ m + log1p (fromIntegral lm1 + kahanSum (fmap (expm1 . subtract m) xs0))++{-+-- TODO(wrengr): Compare precision of the following implementations.+-- We need to make sure to structure it in such a way that the @m@+-- doesn't obliterate the whole purpose of using @exp (x - m)@ in+-- the first place; but supposing we can do that, then it might+-- could help++sumExp = exp . logSumExp++sumExp [] = 0+sumExp xs0@(x:xs) =+ case foldLSE x xs of+ LSE lm1 m+ | isInfinite m -> m+ | otherwise ->+ exp m * kahanSum (fromIntegral lm1 : fmap (expm1 . subtract m) xs0)+-}+++----------------------------------------------------------------+-- | /O(n)/. Log-domain softmax, aka: @(fmap log . softmax)@.+--+-- /N.B./, this requires three passes over the data: two for the+-- 'logSumExp', and a third for the normalization itself. Thus,+-- it is not amenable to list fusion, and hence will use a lot of+-- memory when summing long lists.+logSoftmax :: [Double] -> [Double]+logSoftmax xs = let z = logSumExp xs in z `seq` fmap (subtract z) xs+-- TODO(wrengr): alternatively we could use a variant of 'logSumExp'+-- which doesn't add the maximum back in, and do the final rescaling+-- by subtracting both the maximum and the summation; that is, a more+-- efficient\/straightforward variant of:+-- > logSoftmax xs =+-- > subtract z <$> xs' -- aka @subtract (m + z) <$> xs@+-- > where+-- > m = maximum xs+-- > xs' = subtract m <$> xs+-- > z = logSumExp xs'+-- This works because for any constant @c@, @softmax xs == softmax ((+c)+-- <$> xs)@. Of course, I don't know that doing that would really help+-- precision by much (given the improved performance of using 'logSumExp'+-- in the first place), and saving a single add won't really matter+-- performance-wise. Perhaps if instead of the thing just proposed+-- about avoiding adding the max back in, what if instead we did things+-- exactly as written above: so we subtract off the maximum, but then+-- also do 'logSumExp' such that it subtracts off the maximum of those+-- differences. We could get the top-2 maxima in a single pass without+-- much extra work; but again, unclear whether it'd really help...+++-- | /O(n)/. Normal-domain softmax:+-- > softmax xs = [ exp x / sum [ exp y | y <- xs] | x <- xs ]+--+-- /N.B./, this requires three passes over the data: same as 'logSoftmax'.+softmax :: [Double] -> [Double]+softmax = fmap exp . logSoftmax+{-# INLINE softmax #-}+++----------------------------------------------------------------+-- TODO: double check that everything inlines away, so this data+-- type doesn't introduce any slowdown.+--+-- | A helper type for 'kahanSum'. As a data type, this is really+-- just so we can phrase things as using 'foldl''.+data Kahan = Kahan+ {-# UNPACK #-}!Double -- The total.+ {-# UNPACK #-}!Double -- The error correction.++kahanZero :: Kahan+kahanZero = Kahan 0 0+{-# INLINE kahanZero #-}++-- DONOTSUBMIT: if @x == negativeInfinity@ then our use case demands we return negativeInfinity (so that @0 * infinity == 0@ as desired). But moreover, we really want to short-circuit things to avoid even scanning the rest of the list. To do that, we need to re-inline everything and use recursion directly instead of using 'foldl''.+kahanPlus :: Kahan -> Double -> Kahan+kahanPlus (Kahan t c) x = Kahan t' c'+ where+ -- Beware this getting incorrectly optimized away by constant folding!+ x' = x - c+ t' = t + x'+ c' = (t' - t) - x'+{-# INLINE kahanPlus #-}++fromKahan :: Kahan -> Double+fromKahan (Kahan t _) = t+{-# INLINE fromKahan #-}++-- | /O(n)/. Floating-point summation, via Kahan's algorithm. This+-- is nominally equivalent to 'sum', but greatly mitigates the+-- problem of losing precision.+--+-- /N.B./, this only requires a single pass over the data; but we+-- use a strict left fold for performance, so it's still not amenable+-- to list fusion.+kahanSum :: [Double] -> Double+{-+-- Alas, this implementation loses the optimization below where we+-- avoid NaN and short-circuit to return @LogFloat -infty@ aka 0.+kahanSum = fromKahan . foldl' kahanPlus kahanZero+-}+kahanSum = go kahanZero+ where+ go tc _ | tc `seq` False = undefined+ go tc [] = fromKahan tc+ go tc (x:xs)+ -- Avoid NaN when there's a negInfty in the list. N.B.,+ -- this causes zero to annihilate infinity.+ | x == negativeInfinity = negativeInfinity+ | otherwise = go (kahanPlus tc x) xs+++-- TODO: bring back the 'neumaierSum'+++-- This version *completely* eliminates rounding errors and loss+-- of significance due to catastrophic cancellation during summation.+-- <http://code.activestate.com/recipes/393090/> Also see the other+-- implementations given there. For Python's actual C implementation,+-- see math_fsum in+-- <http://svn.python.org/view/python/trunk/Modules/mathmodule.c?view=markup>+--+-- For merely *mitigating* errors rather than completely eliminating+-- them, see <http://code.activestate.com/recipes/298339/>.+--+-- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@+{-+-- For proof of correctness, see+-- <www-2.cs.cmu.edu/afs/cs/project/quake/public/papers/robust-arithmetic.ps>+def msum(xs):+ partials = [] # sorted, non-overlapping partial sums+ # N.B., the actual C implementation uses a 32 array, doubling size as needed+ for x in xs:+ i = 0+ for y in partials: # for(i = j = 0; j < n; j++)+ if abs(x) < abs(y):+ x, y = y, x+ hi = x + y+ lo = y - (hi - x)+ if lo != 0.0:+ partials[i] = lo+ i += 1+ x = hi+ # does an append of x while dropping all the partials after+ # i. The C version does n=i; and leaves the garbage in place+ partials[i:] = [x]+ # BUG: this last step isn't entirely correct and can lose+ # precision <http://stackoverflow.com/a/2704565/358069>+ return sum(partials, 0.0)+-}+++----------------------------------------------------------------+----------------------------------------------------------- fin.
src/Data/Number/PartialOrd.hs view
@@ -1,24 +1,28 @@-{-# LANGUAGE OverlappingInstances+{-# LANGUAGE CPP , FlexibleInstances , UndecidableInstances #-} +#if __GLASGOW_HASKELL__ < 710+{-# LANGUAGE OverlappingInstances #-}+#endif+ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} ------------------------------------------------------------------- ~ 2009.01.29+-- ~ 2021.10.17 -- | -- Module : Data.Number.PartialOrd--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : semi-portable (OverlappingInstances,...)--- +-- -- The Prelude's 'Ord' class for dealing with ordered types is often -- onerous to use because it requires 'Eq' as well as a total -- ordering. While such total orderings are common, partial orderings--- are moreso. This module presents a class for partially ordered+-- are more so. This module presents a class for partially ordered -- types. ---------------------------------------------------------------- module Data.Number.PartialOrd@@ -45,49 +49,49 @@ class PartialOrd a where -- | like 'compare' cmp :: a -> a -> Maybe Ordering- + -- | like ('>') gt :: a -> a -> Maybe Bool gt x y = case x `cmp` y of Just GT -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like ('>=') ge :: a -> a -> Maybe Bool ge x y = case x `cmp` y of Just LT -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('==') eq :: a -> a -> Maybe Bool eq x y = case x `cmp` y of Just EQ -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like ('/=') ne :: a -> a -> Maybe Bool ne x y = case x `cmp` y of Just EQ -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('<=') le :: a -> a -> Maybe Bool le x y = case x `cmp` y of Just GT -> Just False Just _ -> Just True Nothing -> Nothing- + -- | like ('<') lt :: a -> a -> Maybe Bool lt x y = case x `cmp` y of Just LT -> Just True Just _ -> Just False Nothing -> Nothing- + -- | like 'max'. The default instance returns the left argument -- when they're equal. maxPO :: a -> a -> Maybe a@@ -96,7 +100,7 @@ GT -> Just x EQ -> Just x LT -> Just y- + -- | like 'min'. The default instance returns the left argument -- when they're equal. minPO :: a -> a -> Maybe a@@ -108,7 +112,11 @@ infix 4 `gt`, `ge`, `eq`, `ne`, `le`, `lt`, `maxPO`, `minPO` -instance (Ord a) => PartialOrd a where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPABLE #-}+#endif+ (Ord a) => PartialOrd a where cmp x y = Just $! x `compare` y gt x y = Just $! x > y ge x y = Just $! x >= y@@ -126,11 +134,19 @@ -- returns @Just Eq@ for @notANumber@ then CPP was run wrongly. -- -- The instances inherited from Ord are wrong. So we'll fix them.-instance PartialOrd Float where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ PartialOrd Float where cmp x y | isNaN x || isNaN y = Nothing | otherwise = Just $! x `compare` y -instance PartialOrd Double where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ PartialOrd Double where cmp x y | isNaN x || isNaN y = Nothing | otherwise = Just $! x `compare` y
src/Data/Number/RealToFrac.hs view
@@ -1,25 +1,30 @@ -- Needed to ensure correctness, and because we can't guarantee rules fire -- The MagicHash is for unboxed primitives (-fglasgow-exts also works) -- We only need MagicHash if on GHC, but we can't hide it in an #ifdef-{-# LANGUAGE MultiParamTypeClasses- , OverlappingInstances+{-# LANGUAGE CPP+ , MultiParamTypeClasses , FlexibleInstances- , CPP- , MagicHash #-} +#if __GLASGOW_HASKELL__ < 710+{-# LANGUAGE OverlappingInstances #-}+#endif++-- We don't put these in LANGUAGE, because it's CPP guarded for GHC only+{-# OPTIONS_GHC -XMagicHash #-}+ {-# OPTIONS_GHC -Wall -fwarn-tabs #-} ------------------------------------------------------------------- ~ 2009.01.29+-- ~ 2021.10.17 -- | -- Module : Data.Number.RealToFrac--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : semi-portable (CPP, MPTC, OverlappingInstances)--- +-- -- This module presents a type class for generic conversion between -- numeric types, generalizing @realToFrac@ in order to overcome -- problems with pivoting through 'Rational'@@ -75,7 +80,11 @@ instance (Real a, Fractional a) => RealToFrac a a where realToFrac = id -instance (Real a, Transfinite a, Fractional b, Transfinite b)+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPABLE #-}+#endif+ (Real a, Transfinite a, Fractional b, Transfinite b) => RealToFrac a b where realToFrac x@@ -86,31 +95,55 @@ #ifdef __GLASGOW_HASKELL__-instance RealToFrac Int Float where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Int Float where {-# INLINE realToFrac #-} realToFrac (I# i) = F# (int2Float# i) -instance RealToFrac Int Double where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Int Double where {-# INLINE realToFrac #-} realToFrac (I# i) = D# (int2Double# i) -instance RealToFrac Integer Float where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Integer Float where -- TODO: is there a more primitive way? {-# INLINE realToFrac #-} realToFrac j = Prelude.realToFrac j -instance RealToFrac Integer Double where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Integer Double where -- TODO: is there a more primitive way? {-# INLINE realToFrac #-} realToFrac j = Prelude.realToFrac j -instance RealToFrac Float Double where+instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Float Double where {-# INLINE realToFrac #-} realToFrac (F# f) = D# (float2Double# f)- -instance RealToFrac Double Float where++instance+#if __GLASGOW_HASKELL__ >= 710+ {-# OVERLAPPING #-}+#endif+ RealToFrac Double Float where {-# INLINE realToFrac #-} realToFrac (D# d) = F# (double2Float# d) #endif
src/Data/Number/Transfinite.hs view
@@ -1,29 +1,33 @@ {-# OPTIONS_GHC -Wall -fwarn-tabs #-}---- Unfortunately we need -fglasgow-exts in order to actually pick--- up on the rules (see -ddump-rules). The -frewrite-rules flag--- doesn't do what you want.--- <http://hackage.haskell.org/trac/ghc/ticket/2213>--- <http://www.mail-archive.com/glasgow-haskell-users@haskell.org/msg14313.html>-{-# OPTIONS_GHC -fglasgow-exts #-}+{-# OPTIONS_GHC -O2 -fenable-rewrite-rules #-} +-- FIXME(2023-03-19): Since recent versions of GHC the rewrite rules+-- in this file generate warnings that they may never fire because+-- the rule "Class op exp" may fire first. Although the warning+-- suggests adding a phase limit, that doesn't actually help because+-- the "Class op exp" rule is a built-in. So it's not actually+-- clear how to silence this warning:+-- <https://gitlab.haskell.org/ghc/ghc/-/issues/10595>+--+-- So for now we just silence the warnings.+{-# OPTIONS_GHC -fno-warn-inline-rule-shadowing #-} ------------------------------------------------------------------- ~ 2009.03.09+-- ~ 2021.10.17 -- | -- Module : Data.Number.Transfinite--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable--- +-- -- This module presents a type class for numbers which have -- representations for transfinite values. The idea originated from -- the IEEE-754 floating-point special values, used by -- "Data.Number.LogFloat". However not all 'Fractional' types -- necessarily support transfinite values. In particular, @Ratio@ -- types including 'Rational' do not have portable representations.--- +-- -- For the Glasgow compiler (GHC 6.8.2), "GHC.Real" defines @1%0@ -- and @0%0@ as representations for 'infinity' and 'notANumber', -- but most operations on them will raise exceptions. If 'toRational'@@ -36,7 +40,7 @@ -- * <http://www.haskell.org/pipermail/haskell-prime/2006-February/000791.html> -- -- * <http://www.haskell.org/pipermail/haskell-prime/2006-February/000821.html>--- +-- -- Hugs (September 2006) stays closer to the haskell98 spec and -- offers no way of constructing those values, raising arithmetic -- overflow errors if attempted.@@ -75,24 +79,24 @@ -- is compiled correctly. class (PartialOrd a) => Transfinite a where- + -- | A transfinite value which is greater than all finite values. -- Adding or subtracting any finite value is a no-op. As is -- multiplying by any non-zero positive value (including -- @infinity@), and dividing by any positive finite value. Also -- obeys the law @negate infinity = negativeInfinity@ with all -- appropriate ramifications.- + infinity :: a- - ++ -- | A transfinite value which is less than all finite values. -- Obeys all the same laws as @infinity@ with the appropriate -- changes for the sign difference.- + negativeInfinity :: a- - ++ -- | An exceptional transfinite value for dealing with undefined -- results when manipulating infinite values. The following -- operations must return @notANumber@, where @inf@ is any value@@ -112,11 +116,11 @@ -- -- * @inf \/ inf@ --- -- * @inf `div` inf@+ -- * @inf \`div\` inf@ -- -- * @0 \/ 0@ --- -- * @0 `div` 0@+ -- * @0 \`div\` 0@ -- -- Additionally, any mathematical operations on @notANumber@ -- must also return @notANumber@, and any equality or ordering@@ -125,14 +129,14 @@ -- for 'Eq'; thus, 'eq' and 'ne' are preferred over ('==') and -- ('/=')). Since it returns false for equality, there may be -- more than one machine representation of this `value'.- + notANumber :: a- - ++ -- | Return true for both @infinity@ and @negativeInfinity@, -- false for all other values. isInfinite :: a -> Bool- + -- | Return true only for @notANumber@. isNaN :: a -> Bool @@ -163,7 +167,7 @@ -- This function will raise an error when taking the log of negative -- numbers, rather than returning 'notANumber' as the newer GHC -- implementation does. The reason being that typically this is a--- logical error, and @notANumber@ allows the error to propegate+-- logical error, and @notANumber@ allows the error to propagate -- silently. -- -- In order to improve portability, the 'Transfinite' class is@@ -184,8 +188,8 @@ -- be noted in case your code depends on the implementation details. log :: (Floating a, Transfinite a) => a -> a-{-# SPECIALIZE log :: Double -> Double #-}-{-# SPECIALIZE log :: Float -> Float #-}+{-# INLINE [0] log #-}+-- TODO: should we use NOINLINE or [~0] to avoid the possibility of code bloat? log x = case x `cmp` 0 of Just GT -> Prelude.log x Just EQ -> negativeInfinity@@ -201,12 +205,17 @@ ---------------------------------------------------------------- -- These rules moved here from "LogFloat" in v0.11.2+--+-- FIXME(2023-03-19): Since recent versions of GHC these rules+-- generate warnings that they may never fire because the rule+-- "Class op exp" may fire first. Although the warning suggests+-- adding a phase limit, that doesn't actually help because the+-- "Class op exp" rule is a built-in. So it's not actually clear+-- how to silence this warning:+-- <https://gitlab.haskell.org/ghc/ghc/-/issues/10595> {-# RULES "log/exp" forall x. log (exp x) = x-"log.exp" log . exp = id- "exp/log" forall x. exp (log x) = x-"exp.log" exp . log = id #-} -- We'd like to be able to take advantage of general rule versions
src/Hugs/RealFloat.hs view
@@ -11,15 +11,15 @@ #define REALFLOAT_VERSION normal Prelude version. This could be buggy. #endif ------------------------------------------------------------------- ~ 2010.03.19+-- ~ 2021.10.17 -- | -- Module : Hugs.RealFloat--- Copyright : Copyright (c) 2007--2010 wren ng thornton+-- Copyright : Copyright (c) 2007--2021 wren gayle romano -- License : BSD3--- Maintainer : wren@community.haskell.org+-- Maintainer : wren@cpan.org -- Stability : stable -- Portability : portable (with CPP)--- +-- -- Hugs (September 2006) has buggy definitions for 'Prelude.isNaN' -- and 'Prelude.isInfinite' on Float and Double. If this module is -- run through CPP with the macro @__HUGS__@ set to a value no@@ -44,9 +44,7 @@ import qualified Prelude ---------------------------------------------------------------- -isInfinite :: (RealFloat a) => a -> Bool-{-# SPECIALIZE isInfinite :: Double -> Bool #-}-{-# SPECIALIZE isInfinite :: Float -> Bool #-}+isInfinite :: (RealFloat a) => a -> Bool {-# INLINE isInfinite #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isInfinite x = (1/0) == abs x@@ -56,8 +54,6 @@ isNaN :: (RealFloat a) => a -> Bool-{-# SPECIALIZE isNaN :: Double -> Bool #-}-{-# SPECIALIZE isNaN :: Float -> Bool #-} {-# INLINE isNaN #-} #if defined(__HUGS__) && (__HUGS__ <= 200609) isNaN x = compareEQ x 0 && compareEQ x 1@@ -65,10 +61,11 @@ -- | In Hugs (September 2006), 'compare' always returns @EQ@ if one -- of the arguments is not a number. Thus, if a number is @compareEQ@ -- against multiple different numbers, then it must be @isNaN@.-compareEQ :: (Ord a) => a -> a -> Bool-compareEQ x y = case compare x y of- EQ -> True- _ -> False+compareEQ :: (Ord a) => a -> a -> Bool+compareEQ x y =+ case compare x y of+ EQ -> True+ _ -> False #else isNaN = Prelude.isNaN #endif