numpy：np.abs到底是如何工作的？

Question

I'm trying to implement my own absolute function for gonum dense vectors in Go. I'm wandering if there's a better way of getting the absolute value of an array than squaring and then square rooting?

My main issue is that I've had to implement my own element wise Newtonian square-root function on these vectors and there's a balance between implementation speed and accuracy. If I could avoid using this square-root function I'd be happy.

Answer 1

NumPy source code can be tricky to navigate, because it has so many functions for so many data types. You can find the C-level source code for the absolute value function in the file scalarmath.c.src. This file is actually a template with function definitions that are later replicated by the build system for several data types. Note each function is the "kernel" that is run for each element of the array (looping through the array is done somewhere else). The functions are always called <name of the type>_ctype_absolute, where <name of the type> is the data type it applies to and is generally templated. Let's go through them.

/**begin repeat
 * #name = ubyte, ushort, uint, ulong, ulonglong#
 */

#define @name@_ctype_absolute @name@_ctype_positive

/**end repeat**/

This one is for unsigned types. In this case, the absolute value is the same as np.positive, which just copies the value without doing anything (it is what you get if you have an array a and you do +a).

/**begin repeat
 * #name = byte, short, int, long, longlong#
 * #type = npy_byte, npy_short, npy_int, npy_long, npy_longlong#
 */
static void
@name@_ctype_absolute(@type@ a, @type@ *out)
{
    *out = (a < 0 ? -a : a);
}
/**end repeat**/

This one is for signed integers. Pretty straightforward.

/**begin repeat
 * #name = float, double, longdouble#
 * #type = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */
static void
@name@_ctype_absolute(@type@ a, @type@ *out)
{
    *out = npy_fabs@c@(a);
}
/**end repeat**/

This is for floating-point values. Here npy_fabsf, npy_fabs and npy_fabsl functions are used. These are declared in npy_math.h, but defined through templated C code in npy_math_internal.h.src, essentially calling the C/C99 counterparts (unless C99 is not available, in which case fabsf and fabsl are emulated with fabs). You might think that the previous code should work as well for floating-point types, but actually these are more complicated, since they have things like NaN, infinity or signed zeros, so it is better to use the standard C functions that deal with everything reliably.

static void
half_ctype_absolute(npy_half a, npy_half *out)
{
    *out = a&0x7fffu;
}

This is actually not templated, it is the absolute value function for half-precision floating-point values. Turns out you can change sign by just doing that bitwise operation (set the first bit to 0), since half-precision is simpler (if more limited) than other floating-point types (it's usually the same for those, but with special cases).

/**begin repeat
 * #name = cfloat, cdouble, clongdouble#
 * #type = npy_cfloat, npy_cdouble, npy_clongdouble#
 * #rtype = npy_float, npy_double, npy_longdouble#
 * #c = f,,l#
 */
static void
@name@_ctype_absolute(@type@ a, @rtype@ *out)
{
    *out = npy_cabs@c@(a);
}
/**end repeat**/

This last one is for complex types. These use npy_cabsf, npycabs and npy_cabsl functions, again declared in npy_math.h but in this case template-implemented in npy_math_complex.c.src using C99 functions (unless that is not available, in which case it is emulated with np.hypot).