# SIMD and Inverse Tangent

I'm calculating spherical UV coordinates in a loop, trying to get GCC to vectorize the code.

Here's a compile-able example of what it looks like.

#include <cmath>
#include <cstddef>

/// Calculates UV coordinates from positions on a sphere.
/// The sphere has a radius of one.
//
/// @param pos_array The 3D sphere positions to get the UV coordinates of.
/// The positions are in the format: x, y, z, w (where w is unused)
//
/// @param uv_array This is where the UV coordinates go.
/// A series of 2D vectors (no padding)
void calc_uv_coordinates(const float* pos_array, float* uv_array) noexcept {

constexpr std::size_t chunk_size = 256;

for (std::size_t i = 0; i < chunk_size; i++) {

const float* pos = pos_array + (i * 4);

float* uv = uv_array + (i * 2);

uv = 0.5 + (std::atan2(pos, pos) / (2 * M_PI));
uv = 0.5 - (std::asin(pos) / M_PI);
}
}


But the inverse trig functions are keeping it from being vectorized, as per GCC's report:

gcc -c -O3 -fopt-info-vec-missed test.cc


Prints:

test.cc:16:29: missed: couldn't vectorize loop
/usr/include/c++/9/cmath:145:28: missed: statement clobbers memory: _25 = __builtin_atan2f (_4, _3);
/usr/include/c++/9/cmath:145:28: missed: statement clobbers memory: _25 = __builtin_atan2f (_4, _3);
/usr/include/c++/9/cmath:107:27: missed: statement clobbers memory: _24 = __builtin_asinf (_9);


Which is understandable, I don't think there's any hardware implementations of those functions.

Edit: There is for Intel.

In game development, how is this kind of issue dealt with? I could think of two approaches:

• Write a vectorized atan2 and asin that compute many inverse trig values at once
• Write a SIMD-friendly approximation for each function

A quick google search led me here, where I found an atan approximation.

inline constexpr float fast_atan_approx2(float x) noexcept {
return x / (1 + (0.28 * x * x));
}

inline constexpr float fast_atan_approx2(float x) noexcept {
return (0.9724 * x) - (0.1919 * x * x * x);
}


And from the comments, another approximation.

But looking at the graphs, I don't think they're going to be suitable for texture mapping. Here's a picture, where:

• blue is true arctan
• orange is the first approximation
• red is the second approximation The approximations are only decent in the domain [-1, 1]. So just to restate my question, how is this dealt with? Is there a way to vectorize these trig functions? Do people generally just accept the bottleneck?

• How do you compute your position array? There may be opportunities there to share some knowledge and avoid re-computing it later. – DMGregory Nov 24 '19 at 20:35
• The position values are all in the range of [-1, 1] and they're computed using the quadratic equation (ray-sphere intersections). – tay10r Nov 24 '19 at 20:37
• Those are atan() approximations, you show, not atan2() approximations. For atan2(y,x) approximations, see here: gist.github.com/volkansalma/2972237 – Bram Nov 25 '19 at 2:42
• My first question is whether you really need to have the angles. Instead of spherical UV you can use a cubemap so texture the sphere. – ratchet freak Nov 25 '19 at 17:17
• @ratchetfreak no I don't technically need the angle, just the UV coordinates. I don't know much about cube mapping. I guess I would have to look into to cube mapping with SIMD and see what the image differences are. – tay10r Nov 25 '19 at 19:02

First, you comment that your position values are in [-1,1] range, so that approximation would be just fine, no? Also, that approximation functions you show are for atan(x) and not for atan2(y,x). See here for atan2(y,x) approximations.

Second, you could nudge the compiler even more towards SIMD by supplying your input as structure-of-arrays.

void calc_uv_coordinates
(
const float* __restrict__ posx,
const float* __restrict__ posy,
const float* __restrict__ posz
float* __restrict__ u,
float* __restrict__ v
) noexcept
{
for (int i = 0; i < 256; i++)
u[i] = 0.5f + fast_atan2( posz[i], posx[i] ) / ( 2 * M_PI );

for (int i = 0; i < 256; ++i )
v[i] = 0.5f - fast_asin2( posy[i] ) / M_PI;
}


Note the use of __restrict__ to tell the compiler there is no overlap between the arrays.

And if the compiler still fails to vectorize, there are always intrinsics if you don't mind hardcoding for the CPU architecture.

Oh, and it also helps to make sure your arrays are aligned on 256 bits, if you can. I use __attribute__ ((aligned(32))) for this, or if the memory is dynamic: posix_memalign()

• The approximation doesn't work for the domain [-1, 1] because (for example) atan2(1, 0.1) is evaluated as atan(4) and that's out of the domain of the approximations. I appreciate the other points you're making but a) I'm already using SOA and b) the restrict qualifier isn't technically legal in C++. Maybe since you brought it up, I'll add it as a conditional macro, but the original problem is still unsolved. – tay10r Nov 25 '19 at 1:32
• @tay10r To go fully SoA, you would have the vertices stored in memory as x,x,x,...x, y,y,y,...y, z,z,z,...,z as in my example, and not as x,y,z,w, x,y,z,w, .... x,y,z,w as you have in your example. – Bram Nov 25 '19 at 2:47
• Ah, I didn't know that. What's the difference in the code generation? Edit never mind I'll just mess around with it and see. Thanks for the tip! – tay10r Nov 25 '19 at 3:45