Merge branch 'master' into pz_cbrtl

README.md

@@ -74,6 +74,7 @@ to get latency instead of reciprocal throughput.
 
 When you run ./perf.sh acosf, it does the following:
 
+   $ export OPENMP=-fopenmp
    $ cd src/binary32/acos
    $ make clean
    $ make CFLAGS="-O3 -march=native"

ci.sh

@@ -11,7 +11,7 @@ fi
 
 # use the same order as on https://core-math.gitlabpages.inria.fr/
 FUNCTIONS_EXHAUSTIVE=(acosf acoshf acospif asinf asinhf asinpif atanf atanhf atanpif cbrtf cosf coshf cospif erff erfcf expf exp10f exp10m1f exp2f exp2m1f expm1f lgammaf logf log10f log10p1f log1pf log2f log2p1f rsqrtf sinf sinhf sinpif tanf tanhf tanpif tgammaf)
-FUNCTIONS_WORST=(acos acosh acospi asin asinh asinpi atan atan2 atan2f atan2pi atan2pif atanh atanpi cbrt cos cosh cospi erf erfc exp exp10 exp10m1 exp2 exp2m1 hypot hypotf log log10 log10p1 log1p log2 log2p1 pow powf rsqrt sin sinh sinpi tan tanh tanpi)
+FUNCTIONS_WORST=(acos acosh acospi asin asinh asinpi atan atan2 atan2f atan2pi atan2pif atanh atanpi cbrt cos cosh cospi erf erfc exp exp10 exp10m1 exp2 exp2l exp2m1 hypot hypotf log log10 log10p1 log1p log2 log2p1 pow powf rsqrt sin sinh sinpi tan tanh tanpi)
 FUNCTIONS_SPECIAL=(atan2pif hypotf)
 
 echo "Reference commit is $LAST_COMMIT"

src/binary32/atan2/Makefile

@@ -1,3 +1,14 @@
 FUNCTION_UNDER_TEST := atan2f
 
 include ../support/Makefile.bivariate
+
+all:: check_special
+
+check_special.o: check_special.c
+	$(CC) $(CFLAGS) $(CORE_MATH_DEFINES) -c -fopenmp -o $@ $<
+
+check_special: check_special.o $(CORE_MATH_OBJS) $(FUNCTION_UNDER_TEST)_mpfr.o
+	$(CC) $(LDFLAGS) -fopenmp -o $@ $^ -lmpfr -lm
+
+clean::
+	rm -f check_special

src/binary32/atan2/check_special.c

@@ -0,0 +1,160 @@
+/* Generate special cases for atan2f testing.
+
+Copyright (c) 2022-2024 Stéphane Glondu and Paul Zimmermann, Inria.
+
+This file is part of the CORE-MATH project
+(https://core-math.gitlabpages.inria.fr/).
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdint.h>
+#include <unistd.h>
+#include <fenv.h>
+#include <math.h>
+#include <omp.h>
+#include <mpfr.h>
+
+float cr_atan2f (float, float);
+void ref_init (void);
+
+int rnd1[] = { FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD };
+int rnd2[] = { MPFR_RNDN, MPFR_RNDZ, MPFR_RNDU, MPFR_RNDD };
+
+int rnd = 0;
+
+int verbose = 0;
+
+/* reference code using MPFR */
+static float
+ref_atan2 (float y, float x, int rnd)
+{
+  mpfr_t xi, yi;
+  mpfr_inits2 (24, xi, yi, NULL);
+  mpfr_set_flt (xi, x, MPFR_RNDN);
+  mpfr_set_flt (yi, y, MPFR_RNDN);
+  int inex = mpfr_atan2 (xi, yi, xi, rnd2[rnd]);
+  mpfr_subnormalize (xi, inex, rnd2[rnd]);
+  float ret = mpfr_get_flt (xi, MPFR_RNDN);
+  mpfr_clears (xi, yi, NULL);
+  return ret;
+}
+
+typedef union { uint32_t n; float x; } union_t;
+
+static float
+asfloat (uint32_t n)
+{
+  union_t u;
+  u.n = n;
+  return u.x;
+}
+
+static void
+check (float y, float x)
+{
+  float z, t;
+  t = ref_atan2 (y, x, rnd);
+  z = cr_atan2f (y, x);
+  if ((isnan (t) && !isnan(z)) || (!isnan (t) && isnan(z)) ||
+      (!isnan (t) && !isnan(z) && z != t))
+  {
+    printf ("FAIL y=%a x=%a ref=%a z=%a\n", y, x, t, z);
+    exit (1);
+  }
+}
+
+#define N 1000000000
+
+static void
+check_random (int i)
+{
+  long l;
+  float x, y;
+  struct drand48_data buffer[1];
+  ref_init ();
+  fesetround (rnd1[rnd]);
+  srand48_r (i, buffer);
+  for (int n = 0; n < N; n++)
+  {
+    lrand48_r (buffer, &l);
+    y = asfloat (l);
+    lrand48_r (buffer, &l);
+    x = asfloat (l);
+    check (y, x);
+    check (y, -x);
+    check (-y, x);
+    check (-y, -x);
+  }
+}
+
+int
+main (int argc, char *argv[])
+{
+  while (argc >= 2)
+    {
+      if (strcmp (argv[1], "--rndn") == 0)
+        {
+          rnd = 0;
+          argc --;
+          argv ++;
+        }
+      else if (strcmp (argv[1], "--rndz") == 0)
+        {
+          rnd = 1;
+          argc --;
+          argv ++;
+        }
+      else if (strcmp (argv[1], "--rndu") == 0)
+        {
+          rnd = 2;
+          argc --;
+          argv ++;
+        }
+      else if (strcmp (argv[1], "--rndd") == 0)
+        {
+          rnd = 3;
+          argc --;
+          argv ++;
+        }
+      else if (strcmp (argv[1], "--verbose") == 0)
+        {
+          verbose = 1;
+          argc --;
+          argv ++;
+        }
+      else
+        {
+          fprintf (stderr, "Error, unknown option %s\n", argv[1]);
+          exit (1);
+        }
+    }
+
+  int nthreads;
+#pragma omp parallel
+  nthreads = omp_get_num_threads ();
+  /* check random values */
+#pragma omp parallel for
+  for (int i = 0; i < nthreads; i++)
+    check_random (getpid () + i);
+  return 0;
+}