^(-n) is slow

[ntessore]

Related to #2741, raising to negative integer power does not seem to use the LLVM powi intrinsic and is much slower. This is somewhat unexpected (at least for me it was).

   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: http://docs.julialang.org
   _ _   _| |_  __ _   |  Type "help()" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.4.0-dev+1443 (2014-11-07 15:03 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit 0790c92* (0 days old master)
|__/                   |  x86_64-apple-darwin14.0.0

julia> f(x) = x^-2
f (generic function with 1 method)

julia> g(x) = 1/x^2
g (generic function with 1 method)

julia> h(x) = 1/x/x
h (generic function with 1 method)

julia> f(1.) == g(1.) == h(1.)
true

julia> code_llvm(f, (Float64,))

define double @"julia_f;40365"(double) {
top:
  %1 = call double @pow(double %0, double -2.000000e+00), !dbg !286
  ret double %1, !dbg !286
}

julia> code_llvm(g, (Float64,))

define double @"julia_g;40366"(double) {
top:
  %1 = fmul double %0, %0, !dbg !289
  %2 = fdiv double 1.000000e+00, %1, !dbg !289
  ret double %2, !dbg !289
}

julia> code_llvm(h, (Float64,))

define double @"julia_h;40367"(double) {
top:
  %1 = fdiv double 1.000000e+00, %0, !dbg !292
  %2 = fdiv double %1, %0, !dbg !292
  ret double %2, !dbg !292
}

julia> 

[ntessore]

At least I understood from this page that raising to negative integer power could be done without the floating-point pow call.

[vtjnash]

llvm lowers most powi calls to the pow function. it is unclear to me that 1/x/x and 1/x^2 and x^-2 would necessarily represent the same floating point operation, to 1 ulp, due to possible rounding differences

[ntessore]

Right, of course they don't, I should not have included these. I was merely expecting x^{-n} to be as "elementary" as x^n.

[ntessore]

Maybe I should replace #8878 with a version that says to be careful with negative powers, and use 1/x^n for better performance.

[stevengj]

The question is not whether pow(x, -n) and 1/x^n give exactly the same result (I don't think IEEE guarantees that exponentiation is accurate to 1 ulp), but rather (a) is one significantly more accurate than the other and (b) are there cases where the latter chooses the wrong branch cut (for complex x) or gives a similarly "wrong" answer?

If the answer to both questions is no, and in general if we don't violate any IEEE 754 guarantees about pow, then it should be a legitimate optimization for us to evaluate x^(-n) as 1/x^n.

[vtjnash]

not IEEE, but LLVM may be doing it anyways (configurable with http://llvm.org/docs/LangRef.html#fpmath-metadata)

[stevengj]

cc: @alanedelman

[simonster]

If we privilege accuracy over speed, we shouldn't use powi for small positive integer exponents either, per my comments in #2741. I'm not sure what IEEE guarantees, though.

[oscardssmith]

Fixed on master.