[Mesa-dev] [PATCH] glsl: improve accuracy of atan()
Erik Faye-Lund
kusmabite at gmail.com
Wed Sep 24 04:35:34 PDT 2014
On Wed, Sep 24, 2014 at 4:10 AM, Ian Romanick <idr at freedesktop.org> wrote:
> On 09/23/2014 03:39 PM, Erik Faye-Lund wrote:
>> Our current atan()-approximation is pretty inaccurate at 1.0, so
>> let's try to improve the situation by doing a direct approximation
>> without going through atan.
>>
>> This new implementation uses an 11th degree polynomial to approximate
>> atan in the [-1..1] range, and the following identitiy to reduce the
>> entire range to [-1..1]:
>>
>> atan(x) = 0.5 * pi * sign(x) - atan(1.0 / x)
>>
>> This range-reduction idea is taken from the paper "Fast computation
>> of Arctangent Functions for Embedded Applications: A Comparative
>> Analysis" (Ukil et al. 2011).
>>
>> The polynomial that approximates atan(x) is:
>>
>> x * 0.9999793128310355 - x^3 * 0.3326756418091246 +
>> x^5 * 0.1938924977115610 - x^7 * 0.1173503194786851 +
>> x^9 * 0.0536813784310406 - x^11 * 0.0121323213173444
>>
>> This polynomial was found with the following GNU Octave script:
>>
>> x = linspace(0, 1);
>> y = atan(x);
>> n = [1, 3, 5, 7, 9, 11];
>> format long;
>> polyfitc(x, y, n)
>>
>> The polyfitc function is not built-in, but too long to include here.
>> It can be downloaded from the following URL:
>>
>> http://www.mathworks.com/matlabcentral/fileexchange/47851-constraint-polynomial-fit/content/polyfitc.m
>>
>> This fixes the following piglit test:
>> shaders/glsl-const-folding-01
>>
>> Signed-off-by: Erik Faye-Lund <kusmabite at gmail.com>
>
> Bugzilla: https://bugs.freedesktop.org/show_bug.cgi?id=49713
>
Ah! Thanks a lot for that link, very useful.
>> We could also go down that path, by introducing an ir_unop_atan that
>> gets lowered to a polynomial before code-generation. That would benefit
>> drivers for hardware with atan-support (at least Mali supports this,
>> according to http://limadriver.org/Lima+ISA/), but it worries me a bit to
>> do constant folding with different precision than execution. But perhaps
>> we already have that problem, only a bit more subtle?
>
> We used to implement constant folding for many built-in functions this
> way. Built-in functions like atan were detected in the constant folding
> process, and C library functions were used. Commit 363c14ae0 changed
> this to simplify the constant folding code. This test case began
> failing at that time, and Vinson submitted bug #49713.
Aha, thanks for clearing that up. However, just to be clear, my
alternative suggestion isn't to go back in that direction. It's to
make atan a middle-class citizen, like ir_unop_[exp|log]. These always
gets lowered to ir_unop_[exp|log]2, but in this case the lowering
would go to arithmetic.
Something like this (also no piglit regressions):
--8<--
diff --git a/src/glsl/builtin_functions.cpp b/src/glsl/builtin_functions.cpp
index 9be7f6d..cc6c8b3 100644
--- a/src/glsl/builtin_functions.cpp
+++ b/src/glsl/builtin_functions.cpp
@@ -442,6 +442,7 @@ private:
ir_swizzle *matrix_elt(ir_variable *var, int col, int row);
ir_expression *asin_expr(ir_variable *x);
+ void do_atan(ir_factory &body, const glsl_type *type, ir_variable
*res, operand y_over_x);
/**
* Call function \param f with parameters specified as the linked
@@ -2596,6 +2597,7 @@ builtin_builder::_degrees(const glsl_type *type)
UNOP(sin, ir_unop_sin, always_available)
UNOP(cos, ir_unop_cos, always_available)
+UNOP(atan, ir_unop_atan, always_available)
ir_function_signature *
builtin_builder::_tan(const glsl_type *type)
@@ -2684,11 +2686,7 @@ builtin_builder::_atan2(const glsl_type *type)
ir_factory outer_then(&outer_if->then_instructions, mem_ctx);
/* Then...call atan(y/x) */
- ir_variable *y_over_x =
outer_then.make_temp(glsl_type::float_type, "y_over_x");
- outer_then.emit(assign(y_over_x, div(y, x)));
- outer_then.emit(assign(r, mul(y_over_x, rsq(add(mul(y_over_x, y_over_x),
- imm(1.0f))))));
- outer_then.emit(assign(r, asin_expr(r)));
+ body.emit(assign(r, expr(ir_unop_atan, div(y, x))));
/* ...and fix it up: */
ir_if *inner_if = new(mem_ctx) ir_if(less(x, imm(0.0f)));
@@ -2712,21 +2710,6 @@ builtin_builder::_atan2(const glsl_type *type)
}
ir_function_signature *
-builtin_builder::_atan(const glsl_type *type)
-{
- ir_variable *y_over_x = in_var(type, "y_over_x");
- MAKE_SIG(type, always_available, 1, y_over_x);
-
- ir_variable *t = body.make_temp(type, "t");
- body.emit(assign(t, mul(y_over_x, rsq(add(mul(y_over_x, y_over_x),
- imm(1.0f))))));
-
- body.emit(ret(asin_expr(t)));
-
- return sig;
-}
-
-ir_function_signature *
builtin_builder::_sinh(const glsl_type *type)
{
ir_variable *x = in_var(type, "x");
diff --git a/src/glsl/ir.cpp b/src/glsl/ir.cpp
index 739a9f4..1ec7789 100644
--- a/src/glsl/ir.cpp
+++ b/src/glsl/ir.cpp
@@ -247,6 +247,7 @@ ir_expression::ir_expression(int op, ir_rvalue *op0)
case ir_unop_cos:
case ir_unop_sin_reduced:
case ir_unop_cos_reduced:
+ case ir_unop_atan:
case ir_unop_dFdx:
case ir_unop_dFdx_coarse:
case ir_unop_dFdx_fine:
@@ -513,6 +514,7 @@ static const char *const operator_strs[] = {
"cos",
"sin_reduced",
"cos_reduced",
+ "atan",
"dFdx",
"dFdxCoarse",
"dFdxFine",
diff --git a/src/glsl/ir.h b/src/glsl/ir.h
index 8003f88..c715203 100644
--- a/src/glsl/ir.h
+++ b/src/glsl/ir.h
@@ -1200,6 +1200,7 @@ enum ir_expression_operation {
ir_unop_cos,
ir_unop_sin_reduced, /**< Reduced range sin. [-pi, pi] */
ir_unop_cos_reduced, /**< Reduced range cos. [-pi, pi] */
+ ir_unop_atan,
/*@}*/
/**
diff --git a/src/glsl/ir_constant_expression.cpp
b/src/glsl/ir_constant_expression.cpp
index 1e8b3a3..ca3d4d7 100644
--- a/src/glsl/ir_constant_expression.cpp
+++ b/src/glsl/ir_constant_expression.cpp
@@ -730,6 +730,12 @@ ir_expression::constant_expression_value(struct
hash_table *variable_context)
}
break;
+ case ir_unop_atan:
+ for (unsigned c = 0; c < op[0]->type->components(); c++) {
+ data.f[c] = atanf(op[0]->value.f[c]);
+ }
+ break;
+
case ir_unop_neg:
for (unsigned c = 0; c < op[0]->type->components(); c++) {
switch (this->type->base_type) {
diff --git a/src/glsl/ir_validate.cpp b/src/glsl/ir_validate.cpp
index 97a581d..9cd04d4 100644
--- a/src/glsl/ir_validate.cpp
+++ b/src/glsl/ir_validate.cpp
@@ -317,6 +317,7 @@ ir_validate::visit_leave(ir_expression *ir)
case ir_unop_cos:
case ir_unop_sin_reduced:
case ir_unop_cos_reduced:
+ case ir_unop_atan:
case ir_unop_dFdx:
case ir_unop_dFdx_coarse:
case ir_unop_dFdx_fine:
diff --git a/src/glsl/lower_instructions.cpp b/src/glsl/lower_instructions.cpp
index 6842853..4879424 100644
--- a/src/glsl/lower_instructions.cpp
+++ b/src/glsl/lower_instructions.cpp
@@ -112,6 +112,7 @@
*/
#include "main/core.h" /* for M_LOG2E */
+#include "program/prog_instruction.h" /* for SWIZZLE_XYZW */
#include "glsl_types.h"
#include "ir.h"
#include "ir_builder.h"
@@ -145,6 +146,7 @@ private:
void carry_to_arith(ir_expression *);
void borrow_to_arith(ir_expression *);
void sat_to_clamp(ir_expression *);
+ void atan_to_arith(ir_expression *);
};
} /* anonymous namespace */
@@ -508,6 +510,75 @@ lower_instructions_visitor::sat_to_clamp(ir_expression *ir)
this->progress = true;
}
+void
+lower_instructions_visitor::atan_to_arith(ir_expression *ir)
+{
+ ir_variable *y_over_x = new(ir) ir_variable(ir->type, "y_over_x",
ir_var_temporary);
+ ir_variable *x = new(ir) ir_variable(ir->type, "x", ir_var_temporary);
+ ir_variable *xx = new(ir) ir_variable(ir->type, "xx", ir_var_temporary);
+ ir_variable *tmp = new(ir) ir_variable(ir->type, "tmp", ir_var_temporary);
+
+ ir_instruction &i = *base_ir;
+
+ /* Copy <y_over_x> argument */
+ i.insert_before(y_over_x);
+ i.insert_before(assign(y_over_x, ir->operands[0]));
+
+ /*
+ * range-reduction, first step:
+ *
+ * / y_over_x if |y_over_x| <= 1.0;
+ * x = <
+ * \ 1.0 / y_over_x otherwise
+ */
+ i.insert_before(x);
+ i.insert_before(assign(x, div(min2(abs(y_over_x),
+ new(ir) ir_constant(1.0f)),
+ max2(abs(y_over_x),
+ new(ir) ir_constant(1.0f)))));
+
+ /*
+ * approximate atan by evaluating polynomial:
+ *
+ * x * 0.9999793128310355 - x^3 * 0.3326756418091246 +
+ * x^5 * 0.1938924977115610 - x^7 * 0.1173503194786851 +
+ * x^9 * 0.0536813784310406 - x^11 * 0.0121323213173444
+ */
+ i.insert_before(xx);
+ i.insert_before(assign(xx, mul(x, x)));
+ ir_constant *term1 = new(ir) ir_constant(+0.9999793128310355f);
+ ir_constant *term3 = new(ir) ir_constant(-0.3326756418091246f);
+ ir_constant *term5 = new(ir) ir_constant(+0.1938924977115610f);
+ ir_constant *term7 = new(ir) ir_constant(-0.1173503194786851f);
+ ir_constant *term9 = new(ir) ir_constant(+0.0536813784310406f);
+ ir_constant *term11 = new(ir) ir_constant(-0.0121323213173444f);
+ i.insert_before(tmp);
+ i.insert_before(assign(tmp,
mul(add(mul(add(mul(add(mul(add(mul(add(mul(term11,
+ xx),
+ term9),
+ xx),
+ term7),
+ xx),
+ term5),
+ xx),
+ term3),
+ xx),
+ term1),
+ x)));
+
+ /* range-reduction and sign fixup */
+ ir->operation = ir_binop_mul;
+ ir->operands[0] = add(tmp,
+ csel(greater(abs(y_over_x),
+ new(ir) ir_constant(1.0f,
ir->type->components())),
+ add(mul(tmp,
+ new(ir) ir_constant(-2.0f)),
+ new(ir) ir_constant((float) M_PI_2)),
+ new(ir) ir_constant(0.0f,
ir->type->components())));
+ ir->operands[1] = sign(y_over_x);
+ this->progress = true;
+}
+
ir_visitor_status
lower_instructions_visitor::visit_leave(ir_expression *ir)
{
@@ -569,6 +640,10 @@ lower_instructions_visitor::visit_leave(ir_expression *ir)
sat_to_clamp(ir);
break;
+ case ir_unop_atan:
+ // TODO: make a flag for this also
+ atan_to_arith(ir);
+ break;
+
default:
return visit_continue;
}
diff --git a/src/mesa/drivers/dri/i965/brw_fs_channel_expressions.cpp
b/src/mesa/drivers/dri/i965/brw_fs_channel_expressions.cpp
index cb0a079..5b56f9a 100644
--- a/src/mesa/drivers/dri/i965/brw_fs_channel_expressions.cpp
+++ b/src/mesa/drivers/dri/i965/brw_fs_channel_expressions.cpp
@@ -431,6 +431,7 @@
ir_channel_expressions_visitor::visit_leave(ir_assignment *ir)
case ir_unop_unpack_unorm_2x16:
case ir_unop_unpack_unorm_4x8:
case ir_unop_unpack_half_2x16:
+ case ir_unop_atan:
case ir_binop_ldexp:
case ir_binop_vector_extract:
case ir_triop_vector_insert:
diff --git a/src/mesa/drivers/dri/i965/brw_fs_visitor.cpp
b/src/mesa/drivers/dri/i965/brw_fs_visitor.cpp
index 2d5318a..2a8a6fb 100644
--- a/src/mesa/drivers/dri/i965/brw_fs_visitor.cpp
+++ b/src/mesa/drivers/dri/i965/brw_fs_visitor.cpp
@@ -585,6 +586,8 @@ fs_visitor::visit(ir_expression *ir)
case ir_unop_cos_reduced:
emit_math(SHADER_OPCODE_COS, this->result, op[0]);
break;
+ case ir_unop_atan:
+ unreachable("not reached: should be handled by atan_to_arith");
case ir_unop_dFdx:
emit(FS_OPCODE_DDX, this->result, op[0], fs_reg(BRW_DERIVATIVE_BY_HINT));
diff --git a/src/mesa/drivers/dri/i965/brw_vec4_visitor.cpp
b/src/mesa/drivers/dri/i965/brw_vec4_visitor.cpp
index 473b3c5..8ed09b8 100644
--- a/src/mesa/drivers/dri/i965/brw_vec4_visitor.cpp
+++ b/src/mesa/drivers/dri/i965/brw_vec4_visitor.cpp
@@ -1365,6 +1366,8 @@ vec4_visitor::visit(ir_expression *ir)
case ir_unop_cos_reduced:
emit_math(SHADER_OPCODE_COS, result_dst, op[0]);
break;
+ case ir_unop_atan:
+ unreachable("not reached: should be handled by atan_to_arith");
case ir_unop_dFdx:
case ir_unop_dFdx_coarse:
diff --git a/src/mesa/program/ir_to_mesa.cpp b/src/mesa/program/ir_to_mesa.cpp
index 83cd1ae..991127d 100644
--- a/src/mesa/program/ir_to_mesa.cpp
+++ b/src/mesa/program/ir_to_mesa.cpp
@@ -1115,6 +1117,9 @@ ir_to_mesa_visitor::visit(ir_expression *ir)
case ir_unop_cos_reduced:
emit_scs(ir, OPCODE_COS, result_dst, op[0]);
break;
+ case ir_unop_atan:
+ assert(!"not reached: should be handled by atan_to_arith");
+ break;
case ir_unop_dFdx:
emit(ir, OPCODE_DDX, result_dst, op[0]);
diff --git a/src/mesa/state_tracker/st_glsl_to_tgsi.cpp
b/src/mesa/state_tracker/st_glsl_to_tgsi.cpp
index 7c2f0fc..ce3744b 100644
--- a/src/mesa/state_tracker/st_glsl_to_tgsi.cpp
+++ b/src/mesa/state_tracker/st_glsl_to_tgsi.cpp
@@ -1401,6 +1403,9 @@ glsl_to_tgsi_visitor::visit(ir_expression *ir)
case ir_unop_cos_reduced:
emit_scs(ir, TGSI_OPCODE_COS, result_dst, op[0]);
break;
+ case ir_unop_atan:
+ assert(!"not reached: should be handled by atan_to_arith");
+ break;
case ir_unop_saturate: {
glsl_to_tgsi_instruction *inst;
inst = emit(ir, TGSI_OPCODE_MOV, result_dst, op[0]);
-->8--
Note that this version also uses the new polynomial, rather than the
old approximation. But we get even better constant folding. Dunno if
it's worth it -- but it does get rid of that do_atan()-ugliness ;)
> I'm not worried about any potential performance impact because I've
> never seen any application use atan. Of course, that's also why that
> bug has sat open for long. :)
Yeah, atan is not very common, and when it's needed people seem to
tend to use 1d textures instead. So I think accuracy trumps
performance in this case.
> Reviewed-by: Ian Romanick <ian.d.romanick at intel.com>
>
>> + /* range-reduction fixup */
>> + body.emit(assign(tmp, add(tmp,
>> + csel(greater(abs(y_over_x),
>> + swizzle(imm(1.0f),
>> + SWIZZLE_XXXX,
>> + type->components())),
I guess this should work just as well (sorry for my lack of ir_builder-fu):
- swizzle(imm(1.0f),
- SWIZZLE_XXXX,
- type->components())),
+ imm(1.0f, type->components())),
>> + add(mul(tmp,
>> + imm(-2.0f)),
>> + imm(M_PI_2f)),
>> + imm(0.0f)))));
Hm. Don't I need to expand this last immediate to a vector of
type->components() size as well?
If so, this patch should go on top:
diff --git a/src/glsl/builtin_functions.cpp b/src/glsl/builtin_functions.cpp
index 1820dd5..bfa46eb 100644
--- a/src/glsl/builtin_functions.cpp
+++ b/src/glsl/builtin_functions.cpp
@@ -2749,13 +2749,11 @@ builtin_builder::do_atan(ir_factory &body,
const glsl_type *type, ir_variable *r
/* range-reduction fixup */
body.emit(assign(tmp, add(tmp,
csel(greater(abs(y_over_x),
- swizzle(imm(1.0f),
- SWIZZLE_XXXX,
- type->components())),
+ imm(1.0f, type->components())),
add(mul(tmp,
imm(-2.0f)),
imm(M_PI_2f)),
- imm(0.0f)))));
+ imm(0.0f, type->components())))));
/* sign fixup */
body.emit(assign(res, mul(tmp, sign(y_over_x))));
More information about the mesa-dev
mailing list