From 9e2b867b3f2e9afc9e9f9178788ae07f6be1f3c0 Mon Sep 17 00:00:00 2001 From: Brian Paul Date: Fri, 22 Aug 2008 15:13:47 -0600 Subject: gallium: new u_math.[ch] files for math functions So far, optimized/low-precision versions of exp(), exp2(), log2(), pow(). --- src/gallium/auxiliary/util/u_math.c | 60 +++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 src/gallium/auxiliary/util/u_math.c (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c new file mode 100644 index 0000000000..8bf6551b6e --- /dev/null +++ b/src/gallium/auxiliary/util/u_math.c @@ -0,0 +1,60 @@ +/************************************************************************** + * + * Copyright 2008 Tungsten Graphics, Inc., Cedar Park, Texas. + * All Rights Reserved. + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sub license, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject to + * the following conditions: + * + * The above copyright notice and this permission notice (including the + * next paragraph) shall be included in all copies or substantial portions + * of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NON-INFRINGEMENT. + * IN NO EVENT SHALL TUNGSTEN GRAPHICS AND/OR ITS SUPPLIERS BE LIABLE FOR + * ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, + * TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE + * SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + * + **************************************************************************/ + + + +#include "util/u_math.h" + + + +float pow2_table[POW2_TABLE_SIZE]; + + +static void +init_pow2_table(void) +{ + int i; + for (i = 0; i < POW2_TABLE_SIZE; i++) { + pow2_table[i] = pow(2.0, i / POW2_TABLE_SCALE); + } +} + + +/** + * One time init for math utilities. + */ +void +util_init_math(void) +{ + static boolean initialized = FALSE; + if (!initialized) { + init_pow2_table(); + initialized = TRUE; + } +} + + -- cgit v1.2.3 From e7ff7f78be6c14b6c48e451d6d1f597af379f8f8 Mon Sep 17 00:00:00 2001 From: Michal Krol Date: Sat, 23 Aug 2008 12:31:16 +0200 Subject: util: Silence compiler warnings on Windows. --- src/gallium/auxiliary/util/u_math.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c index 8bf6551b6e..0729114d6a 100644 --- a/src/gallium/auxiliary/util/u_math.c +++ b/src/gallium/auxiliary/util/u_math.c @@ -39,7 +39,7 @@ init_pow2_table(void) { int i; for (i = 0; i < POW2_TABLE_SIZE; i++) { - pow2_table[i] = pow(2.0, i / POW2_TABLE_SCALE); + pow2_table[i] = (float) pow(2.0, i / POW2_TABLE_SCALE); } } -- cgit v1.2.3 From 325cbeb29a63e3d71da00baeab864970fe3aa595 Mon Sep 17 00:00:00 2001 From: Brian Date: Sun, 9 Nov 2008 10:15:32 -0700 Subject: util: Fix util_fast_pow/exp2/log2. - Use a lookup table for log2. - Compute (float) (1 << ipart) by tweaking with the exponent directly to avoid integer overflow and float conversion. - Also table negative exponents to avoid float division and branching. - Implement util_fast_exp as function of util_fast_exp2. -------- Cherry-picked from gallium-0.2: 8415d06d90a197e16554dab98d160334fd9f9f93 This fixes some pow() glitches seen in fslight.c, spectex.c, etc. Conflicts: src/gallium/auxiliary/util/u_math.h --- src/gallium/auxiliary/util/u_math.c | 21 +++++-- src/gallium/auxiliary/util/u_math.h | 112 +++++++++++++++--------------------- 2 files changed, 64 insertions(+), 69 deletions(-) (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c index 0729114d6a..5b3cab4642 100644 --- a/src/gallium/auxiliary/util/u_math.c +++ b/src/gallium/auxiliary/util/u_math.c @@ -30,7 +30,7 @@ #include "util/u_math.h" - +/** 2^x, for x in [-1.0, 1.0[ */ float pow2_table[POW2_TABLE_SIZE]; @@ -38,9 +38,21 @@ static void init_pow2_table(void) { int i; - for (i = 0; i < POW2_TABLE_SIZE; i++) { - pow2_table[i] = (float) pow(2.0, i / POW2_TABLE_SCALE); - } + for (i = 0; i < POW2_TABLE_SIZE; i++) + pow2_table[i] = (float) pow(2.0, (i - POW2_TABLE_OFFSET) / POW2_TABLE_SCALE); +} + + +/** log2(x), for x in [1.0, 2.0[ */ +float log2_table[LOG2_TABLE_SIZE]; + + +static void +init_log2_table(void) +{ + unsigned i; + for (i = 0; i < LOG2_TABLE_SIZE; i++) + log2_table[i] = (float) log2(1.0 + i * (1.0 / LOG2_TABLE_SIZE)); } @@ -53,6 +65,7 @@ util_init_math(void) static boolean initialized = FALSE; if (!initialized) { init_pow2_table(); + init_log2_table(); initialized = TRUE; } } diff --git a/src/gallium/auxiliary/util/u_math.h b/src/gallium/auxiliary/util/u_math.h index 196aeb28fa..be7303e550 100644 --- a/src/gallium/auxiliary/util/u_math.h +++ b/src/gallium/auxiliary/util/u_math.h @@ -174,8 +174,10 @@ static INLINE float logf( float f ) -#define POW2_TABLE_SIZE 256 -#define POW2_TABLE_SCALE ((float) (POW2_TABLE_SIZE-1)) +#define POW2_TABLE_SIZE_LOG2 9 +#define POW2_TABLE_SIZE (1 << POW2_TABLE_SIZE_LOG2) +#define POW2_TABLE_OFFSET (POW2_TABLE_SIZE/2) +#define POW2_TABLE_SCALE ((float)(POW2_TABLE_SIZE/2)) extern float pow2_table[POW2_TABLE_SIZE]; @@ -186,98 +188,78 @@ util_init_math(void); union fi { float f; - int i; - unsigned ui; + int32_t i; + uint32_t ui; }; /** - * Fast approximation to exp(x). - * Compute with base 2 exponents: exp(x) = exp2(log2(e) * x) - * Note: log2(e) is a constant, k = 1.44269 - * So, exp(x) = exp2(k * x); + * Fast version of 2^x * Identity: exp2(a + b) = exp2(a) * exp2(b) - * Let ipart = int(k*x) - * Let fpart = k*x - ipart; - * So, exp2(k*x) = exp2(ipart) * exp2(fpart) + * Let ipart = int(x) + * Let fpart = x - ipart; + * So, exp2(x) = exp2(ipart) * exp2(fpart) * Compute exp2(ipart) with i << ipart * Compute exp2(fpart) with lookup table. */ static INLINE float -util_fast_exp(float x) +util_fast_exp2(float x) { - if (x >= 0.0f) { - float k = 1.44269f; /* = log2(e) */ - float kx = k * x; - int ipart = (int) kx; - float fpart = kx - (float) ipart; - float y = (float) (1 << ipart) - * pow2_table[(int) (fpart * POW2_TABLE_SCALE)]; - return y; - } - else { - /* exp(-x) = 1.0 / exp(x) */ - float k = -1.44269f; - float kx = k * x; - int ipart = (int) kx; - float fpart = kx - (float) ipart; - float y = (float) (1 << ipart) - * pow2_table[(int) (fpart * POW2_TABLE_SCALE)]; - return 1.0f / y; - } + int32_t ipart; + float fpart, mpart; + union fi epart; + + if(x > 129.00000f) + return 3.402823466e+38f; + + if(x < -126.99999f) + return 0.0f; + + ipart = (int32_t) x; + fpart = x - (float) ipart; + + /* same as + * epart.f = (float) (1 << ipart) + * but faster and without integer overflow for ipart > 31 */ + epart.i = (ipart + 127 ) << 23; + + mpart = pow2_table[POW2_TABLE_OFFSET + (int)(fpart * POW2_TABLE_SCALE)]; + + return epart.f * mpart; } /** - * Fast version of 2^x - * XXX the above function could be implemented in terms of this one. + * Fast approximation to exp(x). */ static INLINE float -util_fast_exp2(float x) +util_fast_exp(float x) { - if (x >= 0.0f) { - int ipart = (int) x; - float fpart = x - (float) ipart; - float y = (float) (1 << ipart) - * pow2_table[(int) (fpart * POW2_TABLE_SCALE)]; - return y; - } - else { - /* exp(-x) = 1.0 / exp(x) */ - int ipart = (int) -x; - float fpart = -x - (float) ipart; - float y = (float) (1 << ipart) - * pow2_table[(int) (fpart * POW2_TABLE_SCALE)]; - return 1.0f / y; - } + const float k = 1.44269f; /* = log2(e) */ + return util_fast_exp2(k * x); } -/** - * Based on code from http://www.flipcode.com/totd/ - */ +#define LOG2_TABLE_SIZE_LOG2 8 +#define LOG2_TABLE_SIZE (1 << LOG2_TABLE_SIZE_LOG2) +extern float log2_table[LOG2_TABLE_SIZE]; + + static INLINE float -util_fast_log2(float val) +util_fast_log2(float x) { union fi num; - int log_2; - num.f = val; - log_2 = ((num.i >> 23) & 255) - 128; - num.i &= ~(255 << 23); - num.i += 127 << 23; - num.f = ((-1.0f/3) * num.f + 2) * num.f - 2.0f/3; - return num.f + log_2; + float epart, mpart; + num.f = x; + epart = (float)(((num.i & 0x7f800000) >> 23) - 127); + mpart = log2_table[(num.i & 0x007fffff) >> (23 - LOG2_TABLE_SIZE_LOG2)]; + return epart + mpart; } static INLINE float util_fast_pow(float x, float y) { - /* XXX these tests may need adjustment */ - if (y >= 3.0f && (-0.02f <= x && x <= 0.02f)) - return 0.0f; - if (y >= 50.0f && (-0.9f <= x && x <= 0.9f)) - return 0.0f; return util_fast_exp2(util_fast_log2(x) * y); } -- cgit v1.2.3 From 2276dcf05f7e0ae13ba434615cf7f34dc06b2afe Mon Sep 17 00:00:00 2001 From: Brian Paul Date: Mon, 10 Nov 2008 08:24:45 -0700 Subject: gallium: fix typos in comments --- src/gallium/auxiliary/util/u_math.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c index 5b3cab4642..9c5f616ceb 100644 --- a/src/gallium/auxiliary/util/u_math.c +++ b/src/gallium/auxiliary/util/u_math.c @@ -30,7 +30,7 @@ #include "util/u_math.h" -/** 2^x, for x in [-1.0, 1.0[ */ +/** 2^x, for x in [-1.0, 1.0] */ float pow2_table[POW2_TABLE_SIZE]; @@ -43,7 +43,7 @@ init_pow2_table(void) } -/** log2(x), for x in [1.0, 2.0[ */ +/** log2(x), for x in [1.0, 2.0] */ float log2_table[LOG2_TABLE_SIZE]; -- cgit v1.2.3 From ff42991c720bc1cfbf72194447fde0bebbd65b85 Mon Sep 17 00:00:00 2001 From: Brian Date: Mon, 10 Nov 2008 20:22:36 -0700 Subject: gallium: fix comment again. A half-closed interval was intended. Never saw the [a,b[ notation before. --- src/gallium/auxiliary/util/u_math.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c index 9c5f616ceb..d1571cd1fc 100644 --- a/src/gallium/auxiliary/util/u_math.c +++ b/src/gallium/auxiliary/util/u_math.c @@ -30,7 +30,7 @@ #include "util/u_math.h" -/** 2^x, for x in [-1.0, 1.0] */ +/** 2^x, for x in [-1.0, 1.0) */ float pow2_table[POW2_TABLE_SIZE]; @@ -43,7 +43,7 @@ init_pow2_table(void) } -/** log2(x), for x in [1.0, 2.0] */ +/** log2(x), for x in [1.0, 2.0) */ float log2_table[LOG2_TABLE_SIZE]; -- cgit v1.2.3 From 7e584a70c492698be18bf4d6372b50d1a1c38385 Mon Sep 17 00:00:00 2001 From: Brian Paul Date: Fri, 14 Nov 2008 12:55:05 -0700 Subject: gallium: increase table size for fast log/pow functions The various conformance tests pass now. --- src/gallium/auxiliary/util/u_math.c | 2 +- src/gallium/auxiliary/util/u_math.h | 8 +++++--- 2 files changed, 6 insertions(+), 4 deletions(-) (limited to 'src/gallium/auxiliary/util/u_math.c') diff --git a/src/gallium/auxiliary/util/u_math.c b/src/gallium/auxiliary/util/u_math.c index d1571cd1fc..2811475fa0 100644 --- a/src/gallium/auxiliary/util/u_math.c +++ b/src/gallium/auxiliary/util/u_math.c @@ -52,7 +52,7 @@ init_log2_table(void) { unsigned i; for (i = 0; i < LOG2_TABLE_SIZE; i++) - log2_table[i] = (float) log2(1.0 + i * (1.0 / LOG2_TABLE_SIZE)); + log2_table[i] = (float) log2(1.0 + i * (1.0 / LOG2_TABLE_SCALE)); } diff --git a/src/gallium/auxiliary/util/u_math.h b/src/gallium/auxiliary/util/u_math.h index aee69ab7ba..ac11d7001b 100644 --- a/src/gallium/auxiliary/util/u_math.h +++ b/src/gallium/auxiliary/util/u_math.h @@ -246,8 +246,9 @@ util_fast_exp(float x) } -#define LOG2_TABLE_SIZE_LOG2 8 -#define LOG2_TABLE_SIZE (1 << LOG2_TABLE_SIZE_LOG2) +#define LOG2_TABLE_SIZE_LOG2 16 +#define LOG2_TABLE_SCALE (1 << LOG2_TABLE_SIZE_LOG2) +#define LOG2_TABLE_SIZE (LOG2_TABLE_SCALE + 1) extern float log2_table[LOG2_TABLE_SIZE]; @@ -258,7 +259,8 @@ util_fast_log2(float x) float epart, mpart; num.f = x; epart = (float)(((num.i & 0x7f800000) >> 23) - 127); - mpart = log2_table[(num.i & 0x007fffff) >> (23 - LOG2_TABLE_SIZE_LOG2)]; + /* mpart = log2_table[mantissa*LOG2_TABLE_SCALE + 0.5] */ + mpart = log2_table[((num.i & 0x007fffff) + (1 << (22 - LOG2_TABLE_SIZE_LOG2))) >> (23 - LOG2_TABLE_SIZE_LOG2)]; return epart + mpart; } -- cgit v1.2.3