1 files changed, 106 insertions, 0 deletions
diff --git a/math/test/trigpi_references.h b/math/test/trigpi_references.h
new file mode 100644
index 000000000000..3dc5a3173436
--- /dev/null
+++ b/math/test/trigpi_references.h
@@ -0,0 +1,106 @@
+/*
+ * Extended precision scalar reference functions for trigpi.
+ *
+ * Copyright (c) 2023-2024, Arm Limited.
+ * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
+ */
+
+#include "math_config.h"
+
+#ifndef M_PIl
+#  define M_PIl 3.141592653589793238462643383279502884l
+#endif
+
+long double
+arm_math_sinpil (long double x)
+{
+  /* sin(inf) should return nan, as defined by C23.  */
+  if (isinf (x))
+    return __math_invalid (x);
+
+  long double ax = fabsl (x);
+
+  /* Return 0 for all values above 2^64 to prevent
+     overflow when casting to uint64_t.  */
+  if (ax >= 0x1p64)
+    return x < 0 ? -0.0l : 0.0l;
+
+  /* All integer cases should return 0, with unchanged sign for zero.  */
+  if (x == 0.0l)
+    return x;
+  if (ax == (uint64_t) ax)
+    return x < 0 ? -0.0l : 0.0l;
+
+  return sinl (x * M_PIl);
+}
+
+long double
+arm_math_cospil (long double x)
+{
+  /* cos(inf) should return nan, as defined by C23.  */
+  if (isinf (x))
+    return __math_invalid (x);
+
+  long double ax = fabsl (x);
+
+  if (ax >= 0x1p64)
+    return 1;
+
+  uint64_t m = (uint64_t) ax;
+
+  /* Integer values of cospi(x) should return +/-1.
+    The sign depends on if x is odd or even.  */
+  if (m == ax)
+    return (m & 1) ? -1 : 1;
+
+  /* Values of Integer + 0.5 should always return 0.  */
+  if (ax - 0.5 == m || ax + 0.5 == m)
+    return 0;
+
+  return cosl (ax * M_PIl);
+}
+
+long double
+arm_math_tanpil (long double x)
+{
+  /* inf and x = n + 0.5 for any integral n should return nan.  */
+  if (fabsl (x) >= 0x1p54l)
+    {
+      if (isinf (x))
+	return __math_invalid (x);
+      return x < 0 ? -0.0l : 0.0l;
+    }
+
+  long double i = roundl (x);
+  long double f = x - i;
+  int64_t m = (int64_t) i;
+
+  if (x == 0)
+    {
+      return x;
+    }
+  else if (x == i)
+    {
+      if (x < 0)
+	{
+	  return m & 1 ? 0.0l : -0.0l;
+	}
+      else
+	{
+	  return m & 1 ? -0.0l : 0.0l;
+	}
+    }
+  else if (fabsl (f) == 0.5l)
+    {
+      if (x < 0)
+	{
+	  return m & 1 ? -1.0l / 0.0l : 1.0l / 0.0l;
+	}
+      else
+	{
+	  return m & 1 ? 1.0l / 0.0l : -1.0l / 0.0l;
+	}
+    }
+
+  return tanl (f * M_PIl);
+}