diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown --- a/CHANGELOG.markdown +++ b/CHANGELOG.markdown @@ -1,3 +1,7 @@ +0.1.1 +----- +* Ported `Data.Approximate.Numerics` from [analytics](http://github.com/analytics) + 0.1 --- -* Ported from [analytics](http://github.com/analytics) +* Ported `Data.Approximate.Type` and `Data.Approximate.Mass` from [analytics](http://github.com/analytics) diff --git a/approximate.cabal b/approximate.cabal --- a/approximate.cabal +++ b/approximate.cabal @@ -1,6 +1,6 @@ name: approximate category: Numeric -version: 0.1 +version: 0.1.1 license: BSD3 cabal-version: >= 1.8 license-file: LICENSE @@ -62,12 +62,16 @@ vector >= 0.9 && < 0.11 exposed-modules: + Data.Approximate Data.Approximate.Type Data.Approximate.Mass + Data.Approximate.Numerics if flag(lib-Werror) ghc-options: -Werror + + c-sources: cbits/fast.c ghc-options: -Wall -fwarn-tabs -O2 hs-source-dirs: src diff --git a/cbits/fast.c b/cbits/fast.c new file mode 100644 --- /dev/null +++ b/cbits/fast.c @@ -0,0 +1,357 @@ +/* + * See http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/ + * + * All of these rely on being on a little endian machine, such as an Intel box. + * + * These can be _quite_ inaccurate. ~20% in many cases, but being much faster (~7x) may + * permit more loop iterations of tuning algorithms that only need approximate powers. + * + * This version of Ankerl's algorithm has been extended to provide optionally conservative (lower) bounds + * and also to generate a full linear interpolation across the entire significand rather than 'stair-step' + * at the expense of performing a 64 bit operation rather than a 32 bit one. This is cheap these days. + * + * 'exp' is further improved by using a suggestion by Nic Schraudolph: + * + * "You can get a much better approximation (piecewise rational instead of linear) at + * the cost of a single floating-point division by using better_exp(x) = exp(x/2)/exp(-x/2), + * where exp() is my published approximation but you don't need the additive constant anymore, + * you can use c=0. On machines with hardware division this is very attractive." -- Nic Schraudolph + * + * --Edward Kmett + * + * TODO: Incorporate the techniques from https://code.google.com/p/fastapprox/ to enable us + * to calculate more interesting approximate functions. They might need to be generalized to work on + * Double values where appropriate I suppose. + * + * Magic numbers: + * float /int : round(1<<23/log(2)) = 12102203, 127<<23 = 1065353216 + * double/int : round(1<<20/log(2)) = 1512775, 1023<<20 = 1072693248 + * double/long long: round(1<<52/log(2)) = 6497320848556798, 1023<<52 = 4607182418800017408 + * + * The fudge factors such that exp y <= exp_fast y: + * >>> ceiling (2^23 * (1 - (log (log 2) + 1)/log 2)) + * 722019 + * >>> ceiling (2^20 * (1 - (log (log 2) + 1)/log 2)) + * 90253 + * >>> ceiling (2^52 * (1 - (log (log 2) + 1)/log 2)) + * 387630818974388 + * + * The fudge factor such that exp_fast y <= exp y is uniformly -1 + * + * TODO: perform exponential doubling for pow based on better_exp_fast instead for better accuracy. + */ + +/* Schraudolph's published algorithm extended into the least significant bits to avoid the stair step. + double long long approximation: round 1<<52/log(2) 6497320848556798, + mask = 0x3ff0000000000000LL = 4607182418800017408LL + double approximation: round(1<<20/log(2)) = 1512775, 1023<<20 = 1072693248 +*/ + +/* 4607182418800017408 - 387630818974388 = 4606794787981043020 + +Exponent mask adapted to full 64 bit precision: +>>> 1023 * 2^52 +4607182418800017408 + +The fudge factor for conservative lower bound adapted to full 64 bit precision: +>>> round (2^52 * (1 - (log (log 2) + 1)/log 2)) +387630818974388 + +As a lower bound this is suitable for use when generating Mass and Precision estimates. +*/ +double exp_fast_lb(double a) { + union { double d; long long x; } u; + u.x = (long long)(6497320848556798LL * a + 4606794787981043020); + return u.d; +} + +/* 4607182418800017408 + 1 */ +double exp_fast_ub(double a) { + union { double d; long long x; } u; + u.x = (long long)(6497320848556798LL * a + 4607182418800017409); + return u.d; +} + +double exp_fast(double a) { + union { double d; long long x; } u; + u.x = (long long)(6497320848556798LL * a + 0x3fef127e83d16f12LL); + return u.d; +} + +double better_exp_fast(double a) { + union { double d; long long x; } u, v; + u.x = (long long)(3248660424278399LL * a + 0x3fdf127e83d16f12LL); + v.x = (long long)(0x3fdf127e83d16f12LL - 3248660424278399LL * a); + return u.d / v.d; +} + +/* Schraudolph's published algorithm */ +double exp_fast_schraudolph(double a) { + union { double d; int x[2]; } u; + u.x[1] = (int) (1512775 * a + 1072632447); + u.x[0] = 0; + return u.d; +} + +/* 1065353216 + 1 */ +float expf_fast_ub(float a) { + union { float f; int x; } u; + u.x = (int) (12102203 * a + 1065353217); + return u.f; +} + +/* Schraudolph's published algorithm with John's constants */ +/* 1065353216 - 486411 = 1064866805 */ +float expf_fast(float a) { + union { float f; int x; } u; + u.x = (int) (12102203 * a + 1064866805); + return u.f; +} + +// 1056478197 +double better_expf_fast(float a) { + union { float f; int x; } u, v; + u.x = (long long)(6051102 * a + 1056478197); + v.x = (long long)(1056478197 - 6051102 * a); + return u.f / v.f; +} + +/* 1065353216 - 722019 */ +float expf_fast_lb(float a) { + union { float f; int x; } u; + u.x = (int) (12102203 * a + 1064631197); + return u.f; +} + +/* Ankerl's inversion of Schraudolph's published algorithm, converted to explicit multiplication */ +double log_fast_ankerl(double a) { + union { double d; int x[2]; } u = { a }; + return (u.x[1] - 1072632447) * 6.610368362777016e-7; /* 1 / 1512775.0; */ +} + +double log_fast_ub(double a) { + union { double d; long long x; } u = { a }; + return (u.x - 4606794787981043020) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */ +} + +/* Ankerl's inversion of Schraudolph's published algorithm with my constants */ +double log_fast(double a) { + union { double d; long long x; } u = { a }; + return (u.x - 4606921278410026770) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */ +} + +double log_fast_lb(double a) { + union { double d; long long x; } u = { a }; + return (u.x - 4607182418800017409) * 1.539095918623324e-16; /* 1 / 6497320848556798.0; */ +} + + +/* 1065353216 - 722019 */ +float logf_fast_ub(float a) { + union { float f; int x; } u = { a }; + return (u.x - 1064631197) * 8.262958405176314e-8f; /* 1 / 12102203.0; */ +} + +/* Ankerl's adaptation of Schraudolph's published algorithm with John's constants */ +/* 1065353216 - 486411 = 1064866805 */ +float logf_fast(float a) { + union { float f; int x; } u = { a }; + return (u.x - 1064866805) * 8.262958405176314e-8f; /* 1 / 12102203.0; */ +} + +/* 1065353216 + 1 */ +float logf_fast_lb(float a) { + union { float f; int x; } u = { a }; + return (u.x - 1065353217) * 8.262958405176314e-8f; /* 1 / 12102203.0 */ +} + +/* Ankerl's version of Schraudolph's approximation. */ +double pow_fast_ankerl(double a, double b) { + union { double d; int x[2]; } u = { a }; + u.x[1] = (int)(b * (u.x[1] - 1072632447) + 1072632447); + u.x[0] = 0; + return u.d; +} + +/* + These constants are based loosely on the following comment off of Ankerl's blog: + + "I have used the same trick for float, not double, with some slight modification to the constants to suite IEEE754 float format. The first constant for float is 1<<23/log(2) and the second is 127<<23 (for double they are 1<<20/log(2) and 1023<<20)." -- John +*/ + +/* 1065353216 + 1 = 1065353217 ub */ +/* 1065353216 - 486411 = 1064866805 min RMSE */ +/* 1065353216 - 722019 = 1064631197 lb */ +float powf_fast(float a, float b) { + union { float d; int x; } u = { a }; + u.x = (int)(b * (u.x - 1064866805) + 1064866805); + return u.d; +} + +float powf_fast_lb(float a, float b) { + union { float d; int x; } u = { a }; + u.x = (int)(b * (u.x - 1065353217) + 1064631197); + return u.d; +} + +float powf_fast_ub(float a, float b) { + union { float d; int x; } u = { a }; + u.x = (int)(b * (u.x - 1064631197) + 1065353217); + return u.d; +} + +/* + Now that 64 bit arithmetic is cheap we can (try to) improve on Ankerl's algorithm. + + double long long approximation: round 1<<52/log(2) 6497320848556798, + mask = 0x3ff0000000000000LL = 4607182418800017408LL + +>>> round (2**52 * log (3 / (8 * log 2) + 1/2) / log 2 - 1/2) +261140389990638 +>>> 0x3ff0000000000000 - round (2**52 * log (3 / (8 * log 2) + 1/2) / log 2 - 1/2) +4606921278410026770 + +*/ + +double pow_fast_ub(double a, double b) { + union { double d; long long x; } u = { a }; + u.x = (long long)(b * (u.x - 4606794787981043020LL) + 4607182418800017409LL); + return u.d; +} + +double pow_fast(double a, double b) { + union { double d; long long x; } u = { a }; + u.x = (long long)(b * (u.x - 4606921278410026770LL) + 4606921278410026770LL); + return u.d; +} + +double pow_fast_lb(double a, double b) { + union { double d; long long x; } u = { a }; + u.x = (long long)(b * (u.x - 4607182418800017409LL) + 4606794787981043020LL); + return u.d; +} + +/* should be much more precise with large b, still ~3.3x faster. */ +double pow_fast_precise_ankerl(double a, double b) { + int flipped = 0; + if (b < 0) { + flipped = 1; + b = -b; + } + + /* calculate approximation with fraction of the exponent */ + int e = (int) b; + union { double d; int x[2]; } u = { a }; + u.x[1] = (int)((b - e) * (u.x[1] - 1072632447) + 1072632447); + u.x[0] = 0; + + double r = 1.0; + while (e) { + if (e & 1) { + r *= a; + } + a *= a; + e >>= 1; + } + + r *= u.d; + return flipped ? 1.0/r : r; +} + +/* should be much more precise with large b, still ~3.3x faster. */ +double pow_fast_precise(double a, double b) { + int flipped = 0; + if (b < 0) { + flipped = 1; + b = -b; + } + + /* calculate approximation with fraction of the exponent */ + int e = (int) b; + double d = exp_fast(b - e); + + double r = 1.0; + while (e) { + if (e & 1) r *= a; + a *= a; + e >>= 1; + } + + r *= d; + return flipped ? 1.0/r : r; +} + +double better_pow_fast_precise(double a, double b) { + int flipped = 0; + if (b < 0) { + flipped = 1; + b = -b; + } + + /* calculate approximation with fraction of the exponent */ + int e = (int) b; + double d = better_exp_fast(b - e); + + double r = 1.0; + while (e) { + if (e & 1) r *= a; + a *= a; + e >>= 1; + } + + r *= d; + return flipped ? 1.0/r : r; +} + + +/* should be much more precise with large b */ +float powf_fast_precise(float a, float b) { + int flipped = 0; + if (b < 0) { + flipped = 1; + b = -b; + } + + /* calculate approximation with fraction of the exponent */ + int e = (int) b; + union { float f; int x; } u = { a }; + u.x = (int)((b - e) * (u.x - 1065353216) + 1065353216); + + float r = 1.0f; + while (e) { + if (e & 1) { + r *= a; + } + a *= a; + e >>= 1; + } + + r *= u.f; + return flipped ? 1.0f/r : r; +} + +/* should be much more precise with large b */ +float better_powf_fast_precise(float a, float b) { + int flipped = 0; + if (b < 0) { + flipped = 1; + b = -b; + } + + /* calculate approximation with fraction of the exponent */ + int e = (int) b; + float f = better_expf_fast(b - e); + + float r = 1.0f; + while (e) { + if (e & 1) { + r *= a; + } + a *= a; + e >>= 1; + } + + r *= f; + return flipped ? 1.0f/r : r; +} + diff --git a/src/Data/Approximate.hs b/src/Data/Approximate.hs new file mode 100644 --- /dev/null +++ b/src/Data/Approximate.hs @@ -0,0 +1,18 @@ +-------------------------------------------------------------------- +-- | +-- Copyright : (c) Edward Kmett 2013 +-- License : BSD3 +-- Maintainer: Edward Kmett +-- Stability : experimental +-- Portability: non-portable +-- +-------------------------------------------------------------------- +module Data.Approximate + ( module Data.Approximate.Mass + , module Data.Approximate.Numerics + , module Data.Approximate.Type + ) where + +import Data.Approximate.Mass +import Data.Approximate.Numerics +import Data.Approximate.Type diff --git a/src/Data/Approximate/Numerics.hs b/src/Data/Approximate/Numerics.hs new file mode 100644 --- /dev/null +++ b/src/Data/Approximate/Numerics.hs @@ -0,0 +1,86 @@ +{-# LANGUAGE ForeignFunctionInterface #-} +-------------------------------------------------------------------- +-- | +-- Copyright : (c) Edward Kmett 2013 +-- License : BSD3 +-- Maintainer: Edward Kmett +-- Stability : experimental +-- Portability: non-portable +-- +-- These functions provide wildly inaccurate but very fast +-- approximations to common transcendental functions. +-- +-- The algorithms here are based on Martin Ankerl's optimized 'pow', +-- +-- which is in turn based on +-- +-------------------------------------------------------------------- +module Data.Approximate.Numerics + ( Fast(..) + , blog + ) where + +class Floating a => Fast a where + -- | Calculate an approximate log. + flog :: a -> a + flog_lb :: a -> a + flog_ub :: a -> a + -- | Calculate an approximate exp. + fexp :: a -> a + fexp_lb :: a -> a + fexp_ub :: a -> a + -- | Calculate an approximate pow. + fpow :: a -> a -> a + fpow_lb :: a -> a -> a + fpow_ub :: a -> a -> a + +instance Fast Float where + flog = logf_fast + flog_lb = logf_fast_lb + flog_ub = logf_fast_ub + fexp = better_expf_fast + fexp_lb = expf_fast_lb + fexp_ub = expf_fast_ub + fpow = better_powf_fast_precise + fpow_lb = powf_fast_lb + fpow_ub = powf_fast_ub + +instance Fast Double where + flog = log_fast + flog_lb = log_fast_lb + flog_ub = log_fast_ub + fexp = better_exp_fast + fexp_lb = exp_fast_lb + fexp_ub = exp_fast_ub + fpow = better_pow_fast_precise + fpow_lb = pow_fast_lb + fpow_ub = pow_fast_ub + +-- | Borchardt’s Algorithm from “Dead Reckoning: Calculating without instruments”. +-- +-- This is a remarkably bad approximate logarithm. +-- +-- 'flog' had better outperform it! It is provided merely for comparison. +blog :: Floating a => a -> a +blog x = 6 * (x - 1) / (x + 1 + 4 * sqrt x); + +foreign import ccall unsafe pow_fast_lb :: Double -> Double -> Double +foreign import ccall unsafe pow_fast_ub :: Double -> Double -> Double +foreign import ccall unsafe better_pow_fast_precise :: Double -> Double -> Double +foreign import ccall unsafe powf_fast_lb :: Float -> Float -> Float +foreign import ccall unsafe powf_fast_ub :: Float -> Float -> Float +foreign import ccall unsafe better_powf_fast_precise :: Float -> Float -> Float + +foreign import ccall unsafe better_exp_fast :: Double -> Double +foreign import ccall unsafe exp_fast_lb :: Double -> Double +foreign import ccall unsafe exp_fast_ub :: Double -> Double +foreign import ccall unsafe better_expf_fast :: Float -> Float +foreign import ccall unsafe expf_fast_lb :: Float -> Float +foreign import ccall unsafe expf_fast_ub :: Float -> Float + +foreign import ccall unsafe log_fast :: Double -> Double +foreign import ccall unsafe log_fast_lb :: Double -> Double +foreign import ccall unsafe log_fast_ub :: Double -> Double +foreign import ccall unsafe logf_fast :: Float -> Float +foreign import ccall unsafe logf_fast_lb :: Float -> Float +foreign import ccall unsafe logf_fast_ub :: Float -> Float