Mesa (master): glsl 1.30: Fix numerical instabilities in asinh

[Paul Berry]

Module: Mesa
Branch: master
Commit: 9c7552729971463af015a767e547c43f94fe1212
URL:    http://cgit.freedesktop.org/mesa/mesa/commit/?id=9c7552729971463af015a767e547c43f94fe1212

Author: Paul Berry <stereotype441 at gmail.com>
Date:   Mon Sep 26 15:51:39 2011 -0700

glsl 1.30: Fix numerical instabilities in asinh

The formula we were previously using for asinh:

    asinh x = ln(x + sqrt(x * x + 1))

is numerically unstable: when x is a large negative value, the quantity

    x + sqrt(x * x + 1)

is a small positive value (on the order of 1/(2|x|)).  Since the
logarithm function is very sensitive in this range, any error in the
computation of the square root manifests as a large error in the
result.

This patch changes to the equivalent formula:

    asinh x = sign(x) * ln(abs(x) + sqrt(x * x + 1))

which is only slightly more expensive to compute, and is numerically
stable for all x.

Fixes piglit tests
spec/glsl-1.30/execution/built-in-functions/[fv]s-asinh-*.

Reviewed-by: Chad Versace <chad at chad-versace.us>
Acked-by: Kenneth Graunke <kenneth at whitecape.org>
Reviewed-by: Eric Anholt <eric at anholt.net>

---

 src/glsl/builtins/ir/asinh |   40 ++++++++++++++++++++++++++++++++++++----
 1 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/src/glsl/builtins/ir/asinh b/src/glsl/builtins/ir/asinh
index 93f73cc..d2dc710 100644
--- a/src/glsl/builtins/ir/asinh
+++ b/src/glsl/builtins/ir/asinh
@@ -2,20 +2,52 @@
    (signature float
      (parameters
        (declare (in) float x))
-     ((return (expression float log (expression float + (var_ref x) (expression float sqrt (expression float + (expression float * (var_ref x) (var_ref x)) (constant float (1)))))))))
+     ((return (expression float *
+               (expression float sign (var_ref x))
+               (expression float log
+                (expression float +
+                 (expression float abs (var_ref x))
+                 (expression float sqrt
+                  (expression float +
+                   (expression float * (var_ref x) (var_ref x))
+                   (constant float (1))))))))))
 
    (signature vec2
      (parameters
        (declare (in) vec2 x))
-     ((return (expression vec2 log (expression vec2 + (var_ref x) (expression vec2 sqrt (expression vec2 + (expression vec2 * (var_ref x) (var_ref x)) (constant float (1)))))))))
+     ((return (expression vec2 *
+               (expression vec2 sign (var_ref x))
+               (expression vec2 log
+                (expression vec2 +
+                 (expression vec2 abs (var_ref x))
+                 (expression vec2 sqrt
+                  (expression vec2 +
+                   (expression vec2 * (var_ref x) (var_ref x))
+                   (constant float (1))))))))))
 
    (signature vec3
      (parameters
        (declare (in) vec3 x))
-     ((return (expression vec3 log (expression vec3 + (var_ref x) (expression vec3 sqrt (expression vec3 + (expression vec3 * (var_ref x) (var_ref x)) (constant float (1)))))))))
+     ((return (expression vec3 *
+               (expression vec3 sign (var_ref x))
+               (expression vec3 log
+                (expression vec3 +
+                 (expression vec3 abs (var_ref x))
+                 (expression vec3 sqrt
+                  (expression vec3 +
+                   (expression vec3 * (var_ref x) (var_ref x))
+                   (constant float (1))))))))))
 
    (signature vec4
      (parameters
        (declare (in) vec4 x))
-     ((return (expression vec4 log (expression vec4 + (var_ref x) (expression vec4 sqrt (expression vec4 + (expression vec4 * (var_ref x) (var_ref x)) (constant float (1)))))))))
+     ((return (expression vec4 *
+               (expression vec4 sign (var_ref x))
+               (expression vec4 log
+                (expression vec4 +
+                 (expression vec4 abs (var_ref x))
+                 (expression vec4 sqrt
+                  (expression vec4 +
+                   (expression vec4 * (var_ref x) (var_ref x))
+                   (constant float (1))))))))))
 ))