Adjust complex sqrt in compiler to use the C functions

frontend/lib/immediates/CMakeLists.txt

@@ -18,10 +18,14 @@
 target_sources(ChplFrontend-obj
                PRIVATE
 
-               ifa_vars.cpp
+               # C sources
+               complex-support.h
+               complex-support.c
+
+               # C++ sources
                hash_multipliers.cpp
+               ifa_vars.cpp
                num.cpp
                num.h
                prim_data.h
-
               )

frontend/lib/immediates/complex-support.c

@@ -0,0 +1,79 @@
+/*
+ * Copyright 2024 Hewlett Packard Enterprise Development LP
+ * Other additional copyright holders may be indicated within.
+ *
+ * The entirety of this work is licensed under the Apache License,
+ * Version 2.0 (the "License"); you may not use this file except
+ * in compliance with the License.
+ *
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "complex-support.h"
+
+#include <complex.h>
+
+static double complex makeDoubleComplex(double re, double im) {
+  // Some test environments don't have a working CMPLX
+  // so use another way to set if needed.
+#if defined(CMPLX)
+  return CMPLX(re, im);
+#else
+#ifndef CHPL_DONT_USE_CMPLX_PTR_ALIASING
+#define cmplx_re64(c) (((double *)&(c))[0])
+#define cmplx_im64(c) (((double *)&(c))[1])
+  double complex val;
+  cmplx_re64(val) = re;
+  cmplx_im64(val) = im;
+  return val;
+#else
+  // This can generate bad values in the face of inf/nan values
+  return re + im*_Complex_I;
+#endif
+#endif
+}
+
+static float complex makeFloatComplex(float re, float im) {
+  // Some test environments don't have a working CMPLXF
+  // so use another way to set if needed.
+#if defined(CMPLXF)
+  return CMPLXF(re, im);
+#else
+#ifndef CHPL_DONT_USE_CMPLX_PTR_ALIASING
+#define cmplx_re32(c) (((float *)&(c))[0])
+#define cmplx_im32(c) (((float *)&(c))[1])
+  float complex val;
+  cmplx_re32(val) = re;
+  cmplx_im32(val) = im;
+  return val;
+#else
+  // This can generate bad values in the face of inf/nan values
+  return re + im*_Complex_I;
+#endif
+#endif
+}
+
+struct complex64 complexSqrt64(struct complex64 x) {
+  float complex c = makeFloatComplex(x.r, x.i);
+  float complex n = csqrtf(c);
+  struct complex64 ret;
+  ret.r = crealf(n);
+  ret.i = cimagf(n);
+  return ret;
+}
+struct complex128 complexSqrt128(struct complex128 x) {
+  double complex c = makeDoubleComplex(x.r, x.i);
+  double complex n = csqrt(c);
+  struct complex128 ret;
+  ret.r = creal(n);
+  ret.i = cimag(n);
+  return ret;
+}

frontend/lib/immediates/complex-support.h

@@ -0,0 +1,42 @@
+/*
+ * Copyright 2024 Hewlett Packard Enterprise Development LP
+ * Other additional copyright holders may be indicated within.
+ *
+ * The entirety of this work is licensed under the Apache License,
+ * Version 2.0 (the "License"); you may not use this file except
+ * in compliance with the License.
+ *
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef IMMEDIATES_COMPLEX_SQRT_H
+#define IMMEDIATES_COMPLEX_SQRT_H
+
+/*
+
+ This header exists to work around a problem with complex square root being
+ innacurate on libc++ by using C functions to do the complex square root.
+
+ */
+
+struct complex64 {
+  float r;
+  float i;
+};
+struct complex128 {
+  double r;
+  double i;
+};
+
+struct complex64 complexSqrt64(struct complex64 x);
+struct complex128 complexSqrt128(struct complex128 x);
+
+#endif

frontend/lib/immediates/num.cpp

@@ -31,7 +31,6 @@
 #include <cmath>
 #include <cstdio>
 #include <cstring>
-#include <complex>
 
 #include "num.h"
 #include "prim_data.h"
@@ -489,23 +488,6 @@ coerce_immediate(chpl::Context* context, Immediate *from, Immediate *to) {
           break; \
       }
 
-static complex64 complexSqrt64(complex64 x) {
-  auto c = std::complex<float>(x.r, x.i);
-  auto n = std::sqrt(c);
-  complex64 ret;
-  ret.r = std::real(n);
-  ret.i = std::imag(n);
-  return ret;
-}
-static complex128 complexSqrt128(complex128 x) {
-  auto c = std::complex<double>(x.r, x.i);
-  auto n = std::sqrt(c);
-  complex128 ret;
-  ret.r = std::real(n);
-  ret.i = std::imag(n);
-  return ret;
-}
-
 static void doFoldSqrt(chpl::Context* context,
                        Immediate &im1, /* input */
                        Immediate *imm  /* output */) {

frontend/lib/immediates/num.h

@@ -24,6 +24,10 @@
 #include "chpl/framework/Context.h"
 #include "chpl/framework/UniqueString.h"
 
+extern "C" {
+  #include "complex-support.h"
+}
+
 #include <cstdlib>
 #include <cstring>
 #include <sstream>
@@ -38,15 +42,6 @@
 
 extern unsigned int open_hash_multipliers[256];
 
-struct complex64 {
-  float r;
-  float i;
-};
-struct complex128 {
-  double r;
-  double i;
-};
-
 using ImmString = chpl::detail::PODUniqueString;
 
 //