approximate 0.1 → 0.1.1
raw patch · 5 files changed
+471/−2 lines, 5 files
Files
- CHANGELOG.markdown +5/−1
- approximate.cabal +5/−1
- cbits/fast.c +357/−0
- src/Data/Approximate.hs +18/−0
- src/Data/Approximate/Numerics.hs +86/−0
CHANGELOG.markdown view
@@ -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)
approximate.cabal view
@@ -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
+ cbits/fast.c view
@@ -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;+}+
+ src/Data/Approximate.hs view
@@ -0,0 +1,18 @@+--------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2013+-- License : BSD3+-- Maintainer: Edward Kmett <ekmett@gmail.com>+-- 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
+ src/Data/Approximate/Numerics.hs view
@@ -0,0 +1,86 @@+{-# LANGUAGE ForeignFunctionInterface #-}+--------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2013+-- License : BSD3+-- Maintainer: Edward Kmett <ekmett@gmail.com>+-- 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',+-- <http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/>+-- which is in turn based on+-- <http://nic.schraudolph.org/pubs/Schraudolph99.pdf>+--------------------------------------------------------------------+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