Mephistopheles Faust
Mephistopheles Faust

Reputation: 75

Fastest way to find the absolute value of a double

Working on a Finite element solver for the fun of it, I encountered a math problem where the area of a triangle has to be solved multiple times for different triangular elements on a mesh.

Note: This is all being done in Cython:

For the 2D case, the determinant formula allows us to do this by the least number of arithmetic operations.

A = 0.5|[x_1(y_2-y_3) + x_2(y_3-y_1) + x_3(y_1-y_2)]|

However this will be called once for each element and, in the case of adaptive meshing, for every iteration too. So finding the absolute value becomes a problem. Going with a function that does comparisons and if-else branching seems wrong. So i thought that bit manipulation could be a solution here maybe and went with this approach.

from libc.stdint cimport uint64_t

cdef union DoubleIntUnion:
    double double_nr
    uint64_t int_nr

cdef double fastABS(double number) nogil:
    cdef DoubleIntUnion u_number
    u_number.double_nr = number
    u_number.int_nr &= <uint64_t>0x7FFFFFFFFFFFFFFF
    return u_number.double_nr

this works. But is there any way to do it better in this case or any suggestions ?

Upvotes: 0

Views: 225

Answers (1)

Maxim Egorushkin
Maxim Egorushkin

Reputation: 136495

If you just call std::abs, which inlines floating point number binary and to mask the sign bit, it would be challenging to beat it.

from libcpp.cmath cimport abs as std_abs

cdef double fastABS(double number) noexcept:
    return std_abs(number)

Your solution with union cast involves loading a SIMD register into a general purpose register, masking the sign bit, and loading the general purpose register back into SIMD register. Such loads are not cheap and best avoided.

Benchmark your variants to get full understanding.

Upvotes: 3

Related Questions