diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 --- /dev/null +++ b/CHANGELOG @@ -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: + +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. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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 diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,238 @@ +logfloat +======== +[![CI Status](https://github.com/wrengr/logfloat/actions/workflows/ci.yml/badge.svg)](https://github.com/wrengr/logfloat/actions?query=workflow%3Aci+-event%3Apull_request) +[![Hackage version](https://img.shields.io/hackage/v/logfloat.svg?style=flat)](https://hackage.haskell.org/package/logfloat) +[![Stackage LTS version](https://stackage.org/package/logfloat/badge/lts)](https://stackage.org/lts/package/logfloat) +[![Stackage Nightly version](https://stackage.org/package/logfloat/badge/nightly)](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 ``, 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 , 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) diff --git a/logfloat.cabal b/logfloat.cabal --- a/logfloat.cabal +++ b/logfloat.cabal @@ -1,62 +1,91 @@ +Cabal-Version: 2.2 +-- Cabal >=2.2 is required for: +-- +-- Since 2.1, the Cabal-Version must be the absolutely first thing +-- in the file, even before comments. Also, no longer uses ">=". +-- + ---------------------------------------------------------------- --- wren ng thornton ~ 2010.03.19 +-- wren gayle romano ~ 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. +-- +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 - -- - -- 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 diff --git a/src/Data/Number/LogFloat.hs b/src/Data/Number/LogFloat.hs --- a/src/Data/Number/LogFloat.hs +++ b/src/Data/Number/LogFloat.hs @@ -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 --- --- 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 --- cf -{-# 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., +#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. -- +(~>) :: (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: {-# 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: {-# 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., +-- 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. diff --git a/src/Data/Number/LogFloat/Raw.hs b/src/Data/Number/LogFloat/Raw.hs new file mode 100644 --- /dev/null +++ b/src/Data/Number/LogFloat/Raw.hs @@ -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: +-- +-- +-- 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. +-- Also see the other +-- implementations given there. For Python's actual C implementation, +-- see math_fsum in +-- +-- +-- For merely *mitigating* errors rather than completely eliminating +-- them, see . +-- +-- A good test case is @msum([1, 1e100, 1, -1e100] * 10000) == 20000.0@ +{- +-- For proof of correctness, see +-- +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 + return sum(partials, 0.0) +-} + + +---------------------------------------------------------------- +----------------------------------------------------------- fin. diff --git a/src/Data/Number/PartialOrd.hs b/src/Data/Number/PartialOrd.hs --- a/src/Data/Number/PartialOrd.hs +++ b/src/Data/Number/PartialOrd.hs @@ -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 diff --git a/src/Data/Number/RealToFrac.hs b/src/Data/Number/RealToFrac.hs --- a/src/Data/Number/RealToFrac.hs +++ b/src/Data/Number/RealToFrac.hs @@ -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 diff --git a/src/Data/Number/Transfinite.hs b/src/Data/Number/Transfinite.hs --- a/src/Data/Number/Transfinite.hs +++ b/src/Data/Number/Transfinite.hs @@ -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. --- --- -{-# 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: +-- +-- +-- 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 @@ -- * -- -- * --- +-- -- 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: +-- {-# 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 diff --git a/src/Hugs/RealFloat.hs b/src/Hugs/RealFloat.hs --- a/src/Hugs/RealFloat.hs +++ b/src/Hugs/RealFloat.hs @@ -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