Fix T78324: Different Sky Texture results between CPU and GPU

The problem here was numerical precision: The code calculates the angle between
sun and view direction, and the usual acos(dot(a, b)) approach for that has
poor numerical performance for almost parallel angles.

As a result, the generally tiny difference between floating point computation
between CPU and GPU was enough to make the sun vanish at different radii,
causing different results.

The new version fixes the difference by making the computation much more robust
on both platforms.

intern/cycles/kernel/shaders/node_sky_texture.osl

@@ -122,6 +122,11 @@ vector geographical_to_direction(float lat, float lon)
   return vector(cos(lat) * cos(lon), cos(lat) * sin(lon), sin(lat));
 }
 
+float precise_angle(vector a, vector b)
+{
+  return 2.0 * atan2(length(a - b), length(a + b));
+}
+
 color sky_radiance_nishita(vector dir, float nishita_data[9], string filename)
 {
   /* definitions */
@@ -138,7 +143,7 @@ color sky_radiance_nishita(vector dir, float nishita_data[9], string filename)
   if (dir[2] >= 0.0) {
     /* definitions */
     vector sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2);
-    float sun_dir_angle = acos(dot(dir, sun_dir));
+    float sun_dir_angle = precise_angle(dir, sun_dir);
     float half_angular = angular_diameter / 2.0;
     float dir_elevation = M_PI_2 - direction[0];
 

intern/cycles/kernel/svm/svm_sky.h

@@ -145,7 +145,7 @@ ccl_device float3 sky_radiance_nishita(KernelGlobals *kg,
   if (dir.z >= 0.0f) {
     /* definitions */
     float3 sun_dir = geographical_to_direction(sun_elevation, sun_rotation + M_PI_2_F);
-    float sun_dir_angle = acos(dot(dir, sun_dir));
+    float sun_dir_angle = precise_angle(dir, sun_dir);
     float half_angular = angular_diameter / 2.0f;
     float dir_elevation = M_PI_2_F - direction.x;
 

intern/cycles/util/util_math.h

@@ -787,6 +787,16 @@ ccl_device_inline float compare_floats(float a, float b, float abs_diff, int ulp
   return (abs(__float_as_int(a) - __float_as_int(b)) < ulp_diff);
 }
 
+/* Calculate the angle between the two vectors a and b.
+ * The usual approach acos(dot(a, b)) has severe precision issues for small angles,
+ * which are avoided by this method.
+ * Based on "Mangled Angles" from https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf
+ */
+ccl_device_inline float precise_angle(float3 a, float3 b)
+{
+  return 2.0f * atan2f(len(a - b), len(a + b));
+}
+
 CCL_NAMESPACE_END
 
 #endif /* __UTIL_MATH_H__ */