Added run-time and ISO C library support for C99 complex numbers.

Version 5.52. Tagged as 'RISC_OSLib-5_52'

Makefile

@@ -166,6 +166,7 @@ A_ABSSYM  = RISC_OSLib:o.a_abssym
 #
 HEADERS =\
  CLIB:h.assert \
+ CLIB:h.complex \
  CLIB:h.ctype \
  CLIB:h.errno \
  CLIB:h.float \
@@ -285,7 +286,7 @@ ANSI_OBJS =\
  o.alloc o.armprof o.armsys o.clib o.ctype o.error o.fpprintf \
  o.kernel o.locale o.math o.memcpset o.overmgr o.printf o.scanf \
  o.signal o.sort o.stdio o.stdlib o.string o.swiv o.time o.fenv \
- o.longlong o.mathl o.mathasm
+ o.longlong o.mathl o.mathasm o.complex o.cxsupport
 
 #
 # Objects for Clib:o.ansilibm
@@ -294,7 +295,7 @@ ANSI_MOD_OBJS =\
  m_o.alloc m_o.armprof m_o.armsys m_o.clib m_o.ctype m_o.error m_o.fpprintf \
  m_o.kernel m_o.locale m_o.math m_o.memcpset m_o.overmgr m_o.printf m_o.scanf \
  m_o.signal m_o.sort m_o.stdio m_o.stdlib m_o.string m_o.swiv m_o.time m_o.fenv \
- m_o.longlong m_o.mathl m_o.mathasm m_o.cl_stub_m
+ m_o.longlong m_o.mathl m_o.mathasm m_o.complex m_o.cxsupport m_o.cl_stub_m
 
 #
 # Objects for softload shared library
@@ -304,7 +305,9 @@ RM_OBJS =\
  rm_o.cl_modbody \
  rm_o.alloc \
  rm_o.armsys \
+ rm_o.complex \
  rm_o.ctype \
+ rm_o.cxsupport \
  rm_o.error \
  rm_o.fpprintf \
  rm_o.locale \
@@ -662,6 +665,7 @@ CLIB:h.swis: derived.swis
 
 # Headers:
 CLIB:h.assert: clib.h.assert  ;   ${CP} clib.h.assert $@ ${CPFLAGS}
+CLIB:h.complex: clib.h.complex;   ${CP} clib.h.complex $@ ${CPFLAGS}
 CLIB:h.ctype: clib.h.ctype;       ${CP} clib.h.ctype $@ ${CPFLAGS}
 CLIB:h.errno: clib.h.errno;       ${CP} clib.h.errno $@ ${CPFLAGS}
 CLIB:h.fenv: clib.h.fenv;         ${CP} clib.h.fenv $@ ${CPFLAGS}

VersionASM

@@ -11,13 +11,13 @@
                         GBLS    Module_HelpVersion
                         GBLS    Module_ComponentName
                         GBLS    Module_ComponentPath
-Module_MajorVersion     SETS    "5.51"
-Module_Version          SETA    551
+Module_MajorVersion     SETS    "5.52"
+Module_Version          SETA    552
 Module_MinorVersion     SETS    ""
-Module_Date             SETS    "27 Oct 2004"
-Module_ApplicationDate  SETS    "27-Oct-04"
+Module_Date             SETS    "07 Mar 2005"
+Module_ApplicationDate  SETS    "07-Mar-05"
 Module_ComponentName    SETS    "RISC_OSLib"
 Module_ComponentPath    SETS    "RiscOS/Sources/Lib/RISC_OSLib"
-Module_FullVersion      SETS    "5.51"
-Module_HelpVersion      SETS    "5.51 (27 Oct 2004)"
+Module_FullVersion      SETS    "5.52"
+Module_HelpVersion      SETS    "5.52 (07 Mar 2005)"
                         END

VersionNum

-/* (5.51)
+/* (5.52)
  *
  * This file is automatically maintained by srccommit, do not edit manually.
  * Last processed by srccommit version: 1.2.
  *
  */
-#define Module_MajorVersion_CMHG        5.51
+#define Module_MajorVersion_CMHG        5.52
 #define Module_MinorVersion_CMHG        
-#define Module_Date_CMHG                27 Oct 2004
+#define Module_Date_CMHG                07 Mar 2005
 
-#define Module_MajorVersion             "5.51"
-#define Module_Version                  551
+#define Module_MajorVersion             "5.52"
+#define Module_Version                  552
 #define Module_MinorVersion             ""
-#define Module_Date                     "27 Oct 2004"
+#define Module_Date                     "07 Mar 2005"
 
-#define Module_ApplicationDate          "27-Oct-04"
+#define Module_ApplicationDate          "07-Mar-05"
 
 #define Module_ComponentName            "RISC_OSLib"
 #define Module_ComponentPath            "RiscOS/Sources/Lib/RISC_OSLib"
 
-#define Module_FullVersion              "5.51"
-#define Module_HelpVersion              "5.51 (27 Oct 2004)"
-#define Module_LibraryVersionInfo       "5:51"
+#define Module_FullVersion              "5.52"
+#define Module_HelpVersion              "5.52 (07 Mar 2005)"
+#define Module_LibraryVersionInfo       "5:52"

c/complex

+/* Copyright 2005 Castle Technology Ltd
+ *
+ * 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.
+ */
+
+/* complex.c: ISO 'C' (9899:1999) library code, sections 7.3 and G.6 */
+/* Copyright (C) Acorn Computers Ltd. 2005 */
+
+/* version 1.00 */
+
+#include <float.h>
+#include <fenv.h>
+#include <math.h>
+#include <complex.h>
+
+/* All this code uses Annex F/G-style infinity/nan/exception handling. */
+/* To achieve this we use the built-in (ie FPA) versions of sin/cos etc, */
+/* because the <math.h> ones still do the old HUGE_VAL/errno nonsense. */
+/* When this is changed, we'll tidy this up. */
+
+#pragma STDC FENV_ACCESS ON
+
+#define _pi_    0x1.921FB54442D18p+1      // ..4698998C
+#define _pi_2   0x1.921FB54442D18p0
+#define _pi_4   0x1.921FB54442D18p-1
+#define _3pi_4  0x1.2D97C7F3321D2p+1
+
+extern double __log_cabs(double complex);
+extern double _new_cosh(double x);
+extern double _new_sinh(double x);
+extern double _new_tanh(double x);
+
+static complex double cfetidyexcept(const fenv_t *envp, complex double x)
+{
+    int exc = fetestexcept(FE_ALL_EXCEPT);
+    if (exc & (FE_DIVBYZERO|FE_INVALID))
+    {   if (exc & (FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT))
+            feclearexcept(FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT);
+    }
+    else if (exc & FE_UNDERFLOW)
+    {   if (!(isless(fabs(creal(x)), DBL_MIN) ||
+              isless(fabs(cimag(x)), DBL_MIN)))
+            feclearexcept(FE_UNDERFLOW);
+    }
+    feupdateenv(envp);
+    return x;
+}
+
+static complex float cfetidyexceptf(const fenv_t *envp, complex float x)
+{
+    int exc = fetestexcept(FE_ALL_EXCEPT);
+    if (exc & (FE_DIVBYZERO|FE_INVALID))
+    {   if (exc & (FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT))
+            feclearexcept(FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT);
+    }
+    else if (exc & FE_UNDERFLOW)
+    {   if (!(isless(fabsf(crealf(x)), FLT_MIN) ||
+              isless(fabsf(cimagf(x)), FLT_MIN)))
+            feclearexcept(FE_UNDERFLOW);
+    }
+    feupdateenv(envp);
+    return x;
+}
+
+double complex cacos(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    int cx, cy;
+
+    if (signbit(y)) return conj(cacos(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cx == FP_ZERO)
+    {   if (cy == FP_ZERO) return _pi_2 - I*0;
+        if (cy == FP_NAN) return _pi_2 - I*y;
+    }
+    else if (cx == FP_INFINITE)
+    {   if (cy == FP_NAN) return y - I * INFINITY;
+        if (cy == FP_INFINITE) return (x < 0 ? _3pi_4 : _pi_4) - I * INFINITY;
+        return (x < 0 ? _pi_ : +0.0) - I * INFINITY;
+    }
+    else if (cy == FP_INFINITE)
+    {   if (cx == FP_NAN) return x - I * INFINITY;
+        return _pi_2 - I * INFINITY;
+    }
+
+    return (2/I)*clog(csqrt(0.5*(1+z))+I*csqrt(0.5*(1-z)));
+}
+
+static float complex narrow(double complex func(double complex), float complex z)
+{
+    fenv_t env;
+    feholdexcept(&env);
+    return cfetidyexceptf(&env, (float complex) func(z));
+}
+
+float complex cacosf(float complex z)
+{
+    return narrow(cacos, z);
+}
+
+double complex casin(double complex z)
+{
+    return -I*casinh(I*z);
+}
+
+float complex casinf(float complex z)
+{
+    return -I*casinhf(I*z);
+}
+
+double complex catan(double complex z)
+{
+    return -I*catanh(I*z);
+}
+
+float complex catanf(float complex z)
+{
+    return -I*catanhf(I*z);
+}
+
+double complex ccos(double complex z)
+{
+    return ccosh(I*z);
+}
+
+float complex ccosf(float complex z)
+{
+    return ccoshf(I*z);
+}
+
+double complex csin(double complex z)
+{
+    return -I*csinh(I*z);
+}
+
+float complex csinf(float complex z)
+{
+    return -I*csinhf(I*z);
+}
+
+double complex ctan(double complex z)
+{
+    return -I*ctanh(I*z);
+}
+
+float complex ctanf(float complex z)
+{
+    return -I*ctanhf(I*z);
+}
+
+double complex cacosh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    int cx, cy;
+
+    if (signbit(y)) return conj(cacosh(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cy == FP_INFINITE)
+    {
+        if (cx == FP_INFINITE) return +INFINITY + (x > 0 ? I*_pi_4 : I*_3pi_4);
+        if (cx == FP_NAN) return +INFINITY + I*y;
+    }
+
+    return 2*clog(csqrt(0.5*(z+1))+csqrt(0.5*(z-1)));
+}
+
+float complex cacoshf(float complex z)
+{
+    return narrow(cacosh, z);
+}
+
+double complex casinh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    int cx, cy;
+
+    if (signbit(x)) return -casinh(-z);
+    if (signbit(y)) return conj(casinh(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cy == FP_INFINITE)
+    {   if (cx == FP_NAN) return y + I*x;
+        if (cx == FP_INFINITE) return +INFINITY + I*_pi_4;
+        return +INFINITY + I*_pi_2;
+    }
+    else if (cy == FP_NAN)
+    {   if (cx == FP_NAN || cx == FP_INFINITE) return z;
+        return y + I*y;
+    }
+    else if (cx == FP_INFINITE)
+    {   return x;
+    }
+    else if (cx == FP_NAN)
+    {   if (cy == FP_ZERO) return z;
+        return x + I*x;
+    }
+
+    return clog(z+csqrt(1+z*z));
+}
+
+float complex casinhf(float complex z)
+{
+    return narrow(casinh, z);
+}
+
+double complex catanh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    int cx, cy;
+
+    if (signbit(x)) return -casinh(-z);
+    if (signbit(y)) return conj(catanh(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cy == FP_NAN)
+    {   if (cx == FP_ZERO || cx == FP_INFINITE) return 0 + I*y;
+        if (cx == FP_NAN) return z;
+        return y + I*y;
+    }
+    else if (cx == FP_INFINITE || cy == FP_INFINITE)
+    {   return 0 + I*_pi_2;
+    }
+    else if (cx == FP_NAN)
+    {   return x + I*x;
+    }
+    else if (cy == FP_ZERO)
+    {   if (x == 1) return x / y; /* +INFINITY; DIVBYZERO */
+    }
+
+    return 0.5*(clog(1+z)-clog(1-z));
+}
+
+float complex catanhf(float complex z)
+{
+    return narrow(catanh, z);
+}
+
+double complex ccosh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    double rx, ry;
+    int cx, cy;
+    fenv_t fe;
+
+    if (signbit(x) != signbit(y)) return conj(ccosh(conj(z)));
+
+    x = fabs(x); y = fabs(y);
+    cx = fpclassify(x), cy = fpclassify(y);
+
+    if (cx == FP_NAN)
+    {   if (cy == FP_ZERO || cy == FP_NAN) return x + I*y;
+        return x + I*x;
+    }
+    else if (cy == FP_NAN)
+    {   if (cx == FP_ZERO) return y + I*0;
+        if (cx == FP_INFINITE) return x + I*y;
+        return y + I*y;
+    }
+    else if (cx == FP_INFINITE)
+    {   if (cy == FP_ZERO) return x + I*y;
+        if (cy == FP_INFINITE) return x + I*__d_sin(y);
+    }
+    else if (cy == FP_INFINITE)
+    {   if (cx == FP_ZERO) return __d_cos(y) + I*0;
+    }
+
+    /* Remaining cases:                                                      */
+    /* finite          + I*finite            just need to check for overflow */
+    /* finite non-zero + I*infinity          will generate NaN+i*NaN +IVO    */
+    /* infinity        + I*finite non-zero   inf*cis(y), clear INX after     */
+
+    feholdexcept(&fe);
+
+    rx = _new_cosh(x) * __d_cos(y);
+    ry = _new_sinh(x) * __d_sin(y);
+
+    if (isnan(ry) && y == 0)
+    {   ry = 0;
+        feclearexcept(FE_INVALID);
+    }
+    if (cx == FP_INFINITE)
+        feclearexcept(FE_INEXACT);
+
+    return cfetidyexcept(&fe, rx + I * ry);
+}
+
+float complex ccoshf(float complex z)
+{
+    return narrow(ccosh, z);
+}
+
+double complex csinh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    double rx, ry;
+    int cx, cy;
+    fenv_t fe;
+
+    if (signbit(x)) return -csinh(-z);
+    if (signbit(y)) return conj(csinh(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cx == FP_NAN)
+    {   if (cy == FP_ZERO || cy == FP_NAN) return x + I*y;
+        return x + I*x;
+    }
+    else if (cy == FP_NAN)
+    {   if (cx == FP_ZERO || cx == FP_INFINITE) return x + I*y;
+        return y + I*y;
+    }
+    else if (cy == FP_INFINITE)
+    {   if (cx == FP_ZERO || cx == FP_INFINITE) return x + I*__d_sin(y);
+    }
+    else if (cx == FP_INFINITE)
+    {   if (cy == FP_ZERO) return x + I*y;
+    }
+
+    /* Remaining cases:                                                      */
+    /* finite          + I*finite            just need to check for overflow */
+    /* finite non-zero + I*infinity          will generate NaN+i*NaN +IVO    */
+    /* infinity        + I*finite non-zero   inf*cis(y), clear INX after     */
+
+    feholdexcept(&fe);
+
+    rx = _new_sinh(x) * __d_cos(y);
+    ry = _new_cosh(x) * __d_sin(y);
+
+    if (isnan(ry) && y == 0)
+    {   ry = 0;
+        feclearexcept(FE_INVALID);
+    }
+    if (cx == FP_INFINITE)
+        feclearexcept(FE_INEXACT);
+
+    return cfetidyexcept(&fe, rx + I * ry);
+}
+
+float complex csinhf(float complex z)
+{
+    return narrow(csinh, z);
+}
+
+double complex ctanh(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    double complex r;
+    int cx, cy;
+    fenv_t fe;
+
+    if (signbit(x)) return -ctanh(-z);
+    if (signbit(y)) return conj(ctanh(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+    if (cx == FP_NAN)
+    {   if (cy == FP_ZERO || cy == FP_NAN) return x + I*y;
+        return x + I*x;
+    }
+    else if (cy == FP_NAN)
+    {   if (cx == FP_INFINITE) return 1 + I*0;
+        return y + I*y;
+    }
+    else if (cy == FP_INFINITE)
+    {   if (cx == FP_INFINITE) return 1 + I*0;
+        return __d_tan(y) * (1+I);
+    }
+
+    /* Remaining cases:                                                      */
+    /* finite          + I*finite            just need to check for overflow */
+    /* infinity        + I*finite            1 + i0 sin(2y), clear INX after */
+
+    feholdexcept(&fe);
+
+    #define _asinh_DBL_MAX__4 0x1.633CE8FB9F87Ep+7 /* asinh(DBL_MAX)/4 */
+    if (x > _asinh_DBL_MAX__4)
+    {   r = I*0*__d_sin(2*y);
+        feclearexcept(FE_INEXACT);
+        if (fetestexcept(FE_OVERFLOW|FE_INVALID))
+        {   /* Don't really want to return a NAN for large y - act as for */
+            /* infinite y */
+            feclearexcept(FE_ALL_EXCEPT);
+            r = I*0;
+        }
+        r += _new_tanh(x);
+    }
+    else
+    {   double t = __d_tan(y);
+        /* Note that t will be finite if in range - overflow is impossible. */
+        /* A range error leading to a NaN is possible, of course. */
+        double beta = 1 + t*t;
+        double s = _new_sinh(x);
+        double c = __d_sqrt(1 + s*s);
+        r = (beta*c*s + I*t) / (1 + beta*s*s);
+    }
+
+    return cfetidyexcept(&fe, r);
+}
+
+float complex ctanhf(float complex z)
+{
+    return narrow(ctanh, z);
+}
+
+double complex cexp(double complex z)
+{
+    double x = creal(z), y = cimag(z);
+    double ex, rx, ry;
+    int cx, cy;
+    fenv_t fe;
+    if (signbit(y)) return conj(cexp(conj(z)));
+
+    cx = fpclassify(x), cy = fpclassify(y);
+
+    if (cx == FP_NAN)
+    {   if (cy == FP_ZERO || cy == FP_NAN) return z;
+        return x + I*x;
+    }
+    else if (cy == FP_NAN)
+    {   if (cx == FP_INFINITE) return signbit(x) ? 0 : z;
+        return y + I*y;
+    }
+    else if (cy == FP_INFINITE)
+    {   if (cx == FP_INFINITE) return signbit(x) ? 0 : INFINITY + I*__d_sin(y);
+        return __d_cos(y) + I*__d_sin(y); /* Just to get right sort of NaN */
+    }
+
+    feholdexcept(&fe);
+
+    ex = __d_exp(x);
+    if (isinf(ex) && y == 0) /* Either infinite x, or overflow */
+    {   rx = ex;
+        ry = 0;
+    }
+    else
+    {   rx = ex * __d_cos(y);
+        ry = ex * __d_sin(y);
+        /* If x was infinite, result is exact - clear inexact from cos/sin */
+        if (isinf(x)) feclearexcept(FE_INEXACT);
+    }
+    feupdateenv(&fe);
+    return rx + I*ry;
+}
+
+float complex cexpf(float complex z)
+{
+    return narrow(cexp, z);
+}
+
+double complex clog(double complex z)
+{
+    /* For once, the standard formula gets all exceptional cases right */
+    /* We use a special log_cabs function to avoid overflow though */
+    return __log_cabs(z) + I*carg(z);
+}
+
+float complex clogf(float complex z)
+{
+    return narrow(clog, z);
+}
+
+double complex cpow(double complex x, double complex y)
+{
+    /* Simplistic implementation - the standard does allow it */
+    return cexp(y * clog(x));
+}
+
+float complex cpowf(float complex x, float complex y)
+{
+    fenv_t env;
+    feholdexcept(&env);
+    return cfetidyexceptf(&env, (float complex) cpow(x,y));
+}
+
+/* end of complex.c */

c/math

@@ -79,6 +79,10 @@
 /* so that it still works with the Shared C Library module.       */
 /* const double __huge_val = 1.79769313486231571e+308; */
 
+/* New versions with Annex F behaviour; exported for complex.c */
+double _new_cosh(double x);
+double _new_sinh(double x);
+double _new_tanh(double x);
 
 double frexp(double d, int *lvn)
 {
@@ -384,6 +388,23 @@ float logbf(float s)
 #define _ln2hi_    0x1.62E42FEEp-1
 #define _ln2lo_    0x1.A39EF35793C76p-33
 
+static double fetidyexcept(const fenv_t *envp, double x)
+{
+    int exc = fetestexcept(FE_ALL_EXCEPT);
+    if (exc & (FE_DIVBYZERO|FE_INVALID))
+    {   if (exc & (FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT))
+            feclearexcept(FE_UNDERFLOW|FE_OVERFLOW|FE_INEXACT);
+    }
+    else if (exc & FE_UNDERFLOW)
+    {   /* Got to check for spurious underflow. We are allowed to raise underflow */
+        /* whenever result is tiny (even if exact), so this is sufficient. */
+        if (!isless(fabs(x), DBL_MIN))
+            feclearexcept(FE_UNDERFLOW);
+    }
+    feupdateenv(envp);
+    return x;
+}
+
 static float fetidyexceptf(const fenv_t *envp, float x)
 {
     int exc = fetestexcept(FE_ALL_EXCEPT);
@@ -482,6 +503,37 @@ double sinh(double x)
     }
 }
 
+double _new_sinh(double x)
+{
+    int sign;
+    double y;
+    if (x==0) return x;
+    if (isless(x,0.0)) y = -x, sign = 1; else y = x, sign = 0;
+    if (isgreater(y,1.0))
+    {   fenv_t env;
+        feholdexcept(&env);
+        double w = y - _sinh_lnv, z, r;
+        z = __d_exp(w);  /* may overflow */
+        if (z < 1.0e10) z = z - _sinh_vm2/z;
+        r = z + _sinh_v2m1 * z;
+        return fetidyexcept(&env, sign ? -r : +r);
+    }
+    else if (islessequal(y,1.0e-10))
+    {   feraiseexcept(FE_INEXACT);
+        return x;
+    }
+    else /* NaN case comes here */
+    {
+        double g = x*x;
+        double r;
+        /* I assume this will always raise inexact (except for NaN). */
+        /* Use a (minimax) rational approximation. See Cody & Waite.     */
+        r = ((((_sinh_p3*g + _sinh_p2)*g + _sinh_p1)*g + _sinh_p0)*g) /
+             (((g + _sinh_q2)*g + _sinh_q1)*g + _sinh_q0);
+        return x + x*r;
+    }
+}
+
 double cosh(double x)
 {
     x = fabs(x);
@@ -504,6 +556,25 @@ double cosh(double x)
     return 0.5*(x + 1/x);
 }
 
+double _new_cosh(double x)
+{
+    fenv_t env;
+    feholdexcept(&env);
+    x = fabs(x);
+    if (isgreater(x,1.0))
+    {   x = __d_exp(x - _sinh_lnv);
+        if (x < 1.0e10) x = x + _sinh_vm2/x;
+        x = x + _sinh_v2m1 * x;
+    }
+    else
+    {
+    /* This second part is satisfactory, even though it is simple!           */
+        x = __d_exp(x);
+        x = 0.5*(x + 1/x);
+    }
+    return fetidyexcept(&env, x);
+}
+
 double tanh(double x)
 {
 /* The first two exits avoid premature overflow as well as needless use  */
@@ -537,6 +608,42 @@ double tanh(double x)
     else return x;
 }
 
+double _new_tanh(double x)
+{
+/* The first two exits avoid premature overflow as well as needless use  */
+/* of the exp() function.                                                */
+    int sign;
+    if (isgreater(x,27.0))
+    {   if (!isinf(x)) feraiseexcept(FE_INEXACT);
+        return 1.0;  /* here exp(2x) dominates 1.0    */
+    }
+    else if (isless(x,-27.0))
+    {   if (!isinf(x)) feraiseexcept(FE_INEXACT);
+        return -1.0;
+    }
+    if (x == 0.0) return x;
+    if (isless(x,0.0)) x = -x, sign = 1;
+    else sign = 0;
+    if (isgreater(x,0.549306144334054846)) /* ln(3)/2 is crossover point */
+    {   x = exp(2.0*x);
+        x = 1.0 - 2.0/(x + 1.0);
+        if (sign) return -x;
+        else return x;
+    }
+    if (!islessequal(x,1.0e-10))
+    {   double y = x*x; /* NaN case comes here */
+        /* I assume this will raise inexact (except for NaN) */
+        /* minimax rational approximation                                */
+        y = (((_tanh_p2*y + _tanh_p1)*y + _tanh_p0)*y) /
+             (((y + _tanh_q2)*y + _tanh_q1)*y + _tanh_q0);
+        x = x + x*y;
+    }
+    else
+        feraiseexcept(FE_INEXACT);
+    if (sign) return -x;
+    else return x;
+}
+
 /* Simple forms for now */
 float coshf(float x)
 {

clib/h/complex

+/* Copyright 2005 Castle Technology Ltd
+ *
+ * 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.
+ */
+#pragma force_top_level
+#pragma include_only_once
+
+/* complex.h: ISO 'C' (9899:1999) library header, sections 7.3 and G.6 */
+/* Copyright (C) Acorn Computers Ltd. 2005 */
+/* version 1.00 */
+
+/*
+ * The header <complex.h> defines macros and declares functions that support
+ * complex arithmetic. Each synopsis specifies a family of functions
+ * consistinc of a principal function with one or more double complex
+ * parameters and a double complex or double return value; and other
+ * functions with the same name but with f and l suffixes which are
+ * corresponding functions with float and long double parameters and return
+ * value.
+ *
+ * This implementation supports Annex G of ISO/IEC 9899:1999, which adds
+ * imaginary types, and defines "I" to be _Imaginary_I rather than
+ * _Complex_I.
+ */
+
+#ifndef __complex_h
+#define __complex_h
+
+#ifdef __cplusplus
+#  include <complex++.h>
+#else
+
+#undef imaginary
+#undef complex
+#undef I
+#define imaginary    _Imaginary
+#define complex      _Complex
+   /* a program may undefine and perhaps then redefine the macros complex,
+   /* imaginary, and I. */
+
+#define _Imaginary_I ___i
+   /* a constant expression of type const float _Imaginary, with the value */
+   /* of the imaginary unit. */
+#define _Complex_I   (0.0F+___i)
+   /* a constant expression of type const float _Complex, with the value */
+   /* of the imaginary unit. */
+
+#define I            _Imaginary_I
+   /* I is imaginary. On an implementation without imaginary types, I */
+   /* would expand to _Complex_I. */
+
+/* #pragma STDC CX_LIMITED_RANGE ON|OFF|DEFAULT */
+   /* The usual mathmatical formulas for complex multiply, divide and */
+   /* absolute value are problematic because of their treatment of */
+   /* infinities and because of undue overflow and underflow. The */
+   /* CX_LIMITED_RANGE pragma can be used to inform the implementation */
+   /* that (where the state is "on") the usual mathmatical formulas are */
+   /* acceptable. The pragma can occur either outside external */
+   /* declarations or preceding all explicit declarations and statements */
+   /* inside a compound statement. When outside external declarations, the */
+   /* pragma takes effect from its occurrence until another */
+   /* CX_LIMITED_RANGE pragma is encountered, or until the end of the */
+   /* translation unit. When inside a compound statement, the pragma takes */
+   /* effect from its occurrence until another CX_LIMITED_RANGE pragma is */
+   /* encountered (including within a nested compound statement), or until */
+   /* the end of the compound statement; at the end of a compound */
+   /*  statement the state for the pragma is restored to its condition */
+   /* just before the compound statement. If this pragma is used in any */
+   /* other context, the behaviour is undefined. The default state for */
+   /* this pragma is "off". */
+
+double complex cacos(double complex /*z*/);
+float complex cacosf(float complex /*z*/);
+long double complex cacosl(long double complex /*z*/);
+   /* computes the complex arc cosine of z, with branch cuts outside the */
+   /* interval [-1,+1] along the real axis. */
+   /* Returns: the complex arc cosine value, in the range of a strip */
+   /*          mathematically unbounded along the imaginary axis and in */
+   /*          the interval [0,Pi] along the real axis. */
+double complex casin(double complex /*z*/);
+float complex casinf(float complex /*z*/);
+long double complex casinl(long double complex /*z*/);
+   /* computes the complex arc sine of z, with branch cuts outside the */
+   /* interval [-1,+1] along the real axis. */
+   /* Returns: the complex arc sine value, in the range of a strip */
+   /*          mathematically unbounded along the imaginary axis and in */
+   /*          the interval [-Pi/2,+Pi/2] along the real axis. */
+double complex catan(double complex /*z*/);
+float complex catanf(float complex /*z*/);
+long double complex catanl(long double complex /*z*/);
+   /* computes the complex arc tangent of z, with branch cuts outside the */
+   /* interval [-i,+i] along the imaginary axis. */
+   /* Returns: the complex arc tangent value, in the range of a strip */
+   /*          mathematically unbounded along the imaginary axis and in */
+   /*          the interval [-Pi/2,+Pi/2] along the real axis. */
+double complex ccos(double complex /*z*/);
+float complex ccosf(float complex /*z*/);
+long double complex ccosl(long double complex /*z*/);
+   /* computes the complex cosine of z. */
+   /* Returns: the complex cosine value. */
+double complex csin(double complex /*z*/);
+float complex csinf(float complex /*z*/);
+long double complex csinl(long double complex /*z*/);
+   /* computes the complex sine of z. */
+   /* Returns: the complex sine value. */
+double complex ctan(double complex /*z*/);
+float complex ctanf(float complex /*z*/);
+long double complex ctanl(long double complex /*z*/);
+   /* computes the complex tangent of z. */
+   /* Returns: the complex tangent value. */
+
+double complex cacosh(double complex /*z*/);
+float complex cacoshf(float complex /*z*/);
+long double complex cacoshl(long double complex /*z*/);
+   /* computes the complex arc hyperbolic cosine of z, with a branch cut */
+   /* at values less than 1 along the real axis. */
+   /* Returns: the complex arc hyperbolic cosine value, in the range of a */
+   /*          half-strip of non-negative values along the real axis and */
+   /*          in the interval [-i*Pi,+i*Pi] along the imaginary axis. */
+double complex casinh(double complex /*z*/);
+float complex casinhf(float complex /*z*/);
+long double complex casinhl(long double complex /*z*/);
+   /* computes the complex arc hyperbolic sine of z, with branch cuts */
+   /* outside the interval [-i,+i] along the imaginary axis. */
+   /* Returns: the complex arc hyperbolic sine value, in the range of a */
+   /*          strip mathematically unbounded along the real axis and in */
+   /*          the interval [-i*Pi/2,+i*Pi/2] along the imaginary axis. */
+double complex catanh(double complex /*z*/);
+float complex catanhf(float complex /*z*/);
+long double complex catanhl(long double complex /*z*/);
+   /* computes the complex arc hyperbolic tangent of z, with branch cuts */
+   /* outside the interval [-1,+1] along the real axis. */
+   /* Returns: the complex arc hyperbolic tangent value, in the range of a */
+   /*          strip mathematically unbounded along the real axis and in */
+   /*          the interval [-i*Pi/2,+i*Pi/2] along the imaginary axis. */
+double complex ccosh(double complex /*z*/);
+float complex ccoshf(float complex /*z*/);
+long double complex ccoshl(long double complex /*z*/);
+   /* computes the complex hyperbolic cosine of z. */
+   /* Returns: the complex hyperbolic cosine value. */
+double complex csinh(double complex /*z*/);
+float complex csinhf(float complex /*z*/);
+long double complex csinhl(long double complex /*z*/);
+   /* computes the complex hyperbolic sine of z. */
+   /* Returns: the complex hyperbolic sine value. */
+double complex ctanh(double complex /*z*/);
+float complex ctanhf(float complex /*z*/);
+long double complex ctanhl(long double complex /*z*/);
+   /* computes the complex hyperbolic tangent of z. */
+   /* Returns: the complex hyperbolic tangent value. */
+
+double complex cexp(double complex /*z*/);
+float complex cexpf(float complex /*z*/);
+long double complex cexpl(long double complex /*z*/);
+   /* computes the complex base-e exponential of z. */
+   /* Returns: the complex base-e exponential value */
+double complex clog(double complex /*z*/);
+float complex clogf(float complex /*z*/);
+long double complex clogl(long double complex /*z*/);
+   /* computes the complex natural (base-e) logarithm of z, with a branch */
+   /* cut along the negative real axis. */
+   /* Returns: the complex natural logarithm value, in the range of a */
+   /*          strip mathematically unbounded along the real axis and in */
+   /*          the interval [-i*Pi,+i*Pi] along the imaginary axis. */
+
+double cabs(double complex /*z*/);
+float cabsf(float complex /*z*/);
+long double cabsl(long double complex /*z*/);
+   /* computes the complex absolute value (also called norm, modulus, or */
+   /* magnitude) of z. */
+   /* Returns: the complex absolute value */
+double complex cpow(double complex /*x*/, double complex /*y*/);
+float complex cpowf(float complex /*x*/, float complex /*y*/);
+long double complex cpowl(long double complex /*x*/, long double complex /*y*/);
+   /* computes the complex power function x^y, with a branch cut for the */
+   /* first parameter along the negative real axis. */
+   /* Returns: the complex power function value */
+double complex csqrt(double complex /*z*/);
+float complex csqrtf(float complex /*z*/);
+long double complex csqrtl(long double complex /*z*/);
+   /* computes the complex square root of z, with a branch cut along the */
+   /* negative real axis. */
+   /* Returns: the complex square root value, in the range of the right */
+   /*          half-plane (including the imaginary axis). */
+
+double carg(double complex /*z*/);
+float cargf(float complex /*z*/);
+long double cargl(long double complex /*z*/);
+   /* computes the argument (also called phase angle) of z, with a branch */
+   /* cut along the negative real axis. */
+   /* Returns: the value of the argument in the interval [-Pi,+Pi] */
+double cimag(double complex /*z*/);
+float cimagf(float complex /*z*/);
+long double cimagl(long double complex /*z*/);
+#define cimag(z) ((void) sizeof cimag(z), (double _Imaginary)(z)/_Imaginary_I)
+#define cimagf(z) ((void) sizeof cimagf(z), (float _Imaginary)(z)/_Imaginary_I)
+#define cimagl(z) ((void) sizeof cimagl(z), (long double _Imaginary)(z)/_Imaginary_I)
+   /* computes the imaginary part of z. */
+   /* Returns: the imaginary part value (as a real) */
+double complex conj(double complex /*z*/);
+float complex conjf(float complex /*z*/);
+long double complex conjl(long double complex /*z*/);
+#define conj(z) ((void) sizeof conj(z), __conj (double _Complex)(z))
+#define conjf(z) ((void) sizeof conjf(z), __conj (float _Complex)(z))
+#define conjl(z) ((void) sizeof conjl(z), __conj (long double _Complex)(z))
+   /* computes the complex conjugate of z, by reversing the sign of its */
+   /* imaginary part. */
+   /* Returns: the complex conjugate value */
+double complex cproj(double complex /*z*/);
+float complex cprojf(float complex /*z*/);
+long double complex cprojl(long double complex /*z*/);
+   /* computes a projection of z onto the Riemann sphere: z projects to z */
+   /* except that all complex infinities (even those with one infinite */
+   /* part and one NaN part) project to positive infinity on the real */
+   /* axis. If z has an infinite part, then cproj(z) is equivalent to */
+   /*          INFINITY + I * copysign(0.0, cimag(z))           */
+   /* Returns: the value of the projection onto the Riemann sphere */
+double creal(double complex /*z*/);
+float crealf(float complex /*z*/);
+long double creall(long double complex /*z*/);
+#define creal(z) ((void) sizeof creal(z), (double)(z))
+#define crealf(z) ((void) sizeof crealf(z), (float)(z))
+#define creall(z) ((void) sizeof creall(z), (long double)(z))
+   /* computes the real part of z. */
+   /* Returns: the real part value */
+
+#endif /* __cplusplus */
+
+#endif
+
+/* end of complex.h */

clib/h/fenv

@@ -17,13 +17,35 @@
 
 /* fenv.h: ISO 'C' (9899:1999) library header, section 7.6 */
 /* Copyright (C) Acorn Computers Ltd. 2002 */
-/* version 1.00 */
+/* version 1.01 */
 
 #ifndef __fenv_h
 #define __fenv_h
 
-typedef unsigned int fexcept_t;
-  /* represents the floating-point status flags collectively */
+/*
+ * The header <fenv.h> declares two types and several macros and functions
+ * to provide access to the floating-point environment. The floating-point
+ * environment refers to any floating-point status flags and contral modes
+ * supported by the implementation. A floating-point status flag is a
+ * system variable whose value is set (but never cleared) when a floating-
+ * point exception is raised, which occurs as a side effect of exceptional
+ * floating-point arithmetic to provide auxiliary information. A floating-
+ * point control mode is a system variable whose value may be set by the
+ * user to affect the subsequent behaviour of floating-point arithmetic.
+ *
+ * Certain programming conventions support the intended model of use for
+ * the programming environment:
+ *
+ *   - a function call does not alter its caller's floating-point control
+ *     modes, clear its caller's floating-point status flags, nor depend
+ *     on the state of its caller's floating-point status flags unless the
+ *     function is so documented;
+ *   - a function call is assumed to require default floating-point control
+ *     modes, unless its documentation promises otherwise;
+ *   - a function call is assumed to have the potential for raising
+ *     floating-point exceptions, unless its documentation promises
+ *     otherwise.
+ */
 
 typedef struct __fenv_t_struct
 { unsigned __status;
@@ -31,6 +53,9 @@ typedef struct __fenv_t_struct
 } fenv_t;
   /* represents the entire floating-point environment */
 
+typedef unsigned int fexcept_t;
+  /* represents the floating-point status flags collectively */
+
 #define FE_INVALID    0x01
 #define FE_DIVBYZERO  0x02
 #define FE_OVERFLOW   0x04
@@ -47,6 +72,32 @@ typedef struct __fenv_t_struct
   /* at program startup. It can be used as an argument to <fenv.h> */
   /* functions that manage the floating point environment. */
 
+/* #pragma STDC FENV_ACCESS ON|OFF|DEFAULT */
+  /* The FENV_ACCESS pragma provides a means to inform the implementation */
+  /* when a program might access the floating-point environment to test */
+  /* floating-point status flags or run under non-default floating-point */
+  /* control modes. The pragma shall occur either outside external */
+  /* declarations or preceding all explicit declarations and statements */
+  /* inside a compound statement. When outside external declarations, the */
+  /* pragma takes effect from its occurrence until another FENV_ACCESS */
+  /* pragma is encountered, or until the end of the translation unit. When */
+  /* inside a compound statement, the pragma takes effect from its */
+  /* occurrence until another FENV_ACCESS pragma is encountered (including */
+  /* within a nested compound statement), or until the end of the compound */
+  /* statement; at the end of a compound statement the state for the */
+  /* pragma is restored to its condition just before the compound */
+  /* statement. If this pragma is used in any other context, the behaviour */
+  /* is undefined. If part of a program tests floating-point status flags, */
+  /* sets floating-point control modes, or runs under non-default mode */
+  /* settings, but was translated with the state for the FENV_ACCESS */
+  /* pragma "off", the behaviour is undefined. The default state ("on" or */
+  /* "off") for the pragma is implementation-defined ("off" in  */
+  /* Norcroft C). (When execution passes from a part of the program */
+  /* translated with FENV_ACCESS "off" to a part translated with */
+  /* FENV_ACCESS "on", the state of the floating-point status flags is */
+  /* unspecified and the floatng-point control modes have their default */
+  /* settings.) */
+
 #ifdef __cplusplus
 extern "C" {
 #endif

clib/h/math

@@ -18,7 +18,7 @@
 /* math.h: ISO 'C' (9899:1999) library header, section 7.12 */
 /* Copyright (C) Codemist Ltd. */
 /* Copyright (C) Acorn Computers Ltd. 1991 */
-/* version 0.06 */
+/* version 0.07 */
 
 #ifndef __math_h
 #define __math_h
@@ -70,39 +70,57 @@ extern const double HUGE_VAL;
 #define FP_NORMAL    2
 #define FP_INFINITE  3
 #define FP_NAN       4
-   /* the mutually exclusive kinds of floating-point values for fpclassify() */
+   /* the mutually exclusive kinds of floating-point values for fpclassify */
 #endif
 
 #define FP_FAST_FMAF 1
    /* indicates that the fmaf function generally executes about as fast as */
-   /* a multiply and an add of float operands. */
+   /* a multiply and an add of float operands. fma and fmal are slow. */
 
 #define FP_ILOGB0    (-0x7FFFFFFF) /* -INT_MAX */
 #define FP_ILOGBNAN  (~0x7FFFFFFF) /* INT_MIN */
-   /* integer constant expressions whose values are returned by ilogb(x) if */
-   /* x is zero or NaN, respectively. */
+   /* integer constant expressions whose values are returned by ilogb(x) */
+   /* if x is zero or NaN, respectively. */
 
 #define MATH_ERRNO 1
 #define MATH_ERREXCEPT 2
 
 #undef math_errhandling
-   /* math_errhandling is not currently provided, as an indication that this */
-   /* is in flux. At present, "invalid", "divide-by-zero" and "overflow" */
-   /* traps are enabled by default, so to raise these would cause SIGFPE. */
-   /* C90 does not allow <math.h> functions to do this, so the functions */
-   /* that were present in C90 do not raise these exceptions. They instead */
-   /* set errno to EDOM or ERANGE, as required by C90. */
+   /* math_errhandling is not currently provided, as an indication that */
+   /* this is in flux. At present, "invalid", "divide-by-zero" and */
+   /* "overflow" traps are enabled by default, so to raise these would */
+   /* cause SIGFPE. C90 does not allow <math.h> functions to do this, so */
+   /* the functions present in C90 do not raise these exceptions. */
+   /* They instead set errno to EDOM or ERANGE, as required by C90. */
    /* The new C99 functions do raise exceptions, and don't set errno. But */
    /* because this will (by default) lead to SIGFPE being raised, the */
-   /* behaviour does not currently conform to the C99 standard. To avoid the */
-   /* traps you can use feholdexcept(). A future C99 version of the library */
-   /* will have all traps disabled by default, making these functions */
-   /* conforming, and will have new versions of the C90 functions which do */
-   /* raise exceptions. Note that some C99 functions may still be using */
-   /* errno and DBL_MAX, but this will be changed (even for existing */
-   /* clients) by newer Shared C libraries. The original C90 functions will */
-   /* not change, to preserve compatibility for older clients. */
-   /* See also HUGE_VAL above. */
+   /* behaviour does not currently conform to the C99 standard. To avoid */
+   /* the traps you can use feholdexcept(). A future C99 version of the */
+   /* library will have all traps disabled by default, making these */
+   /* functions conforming, and will have new versions of the C90 */
+   /* functions which do raise exceptions. Note that some C99 functions */
+   /* may still be using errno and DBL_MAX, but this will be changed (even */
+   /* for existing clients) by newer Shared C libraries. The original C90 */
+   /* functions will not change, to preserve compatibility for older */
+   /* clients. See also HUGE_VAL above. */
+
+/* #pragma STDC FP_CONTRACT ON|OFF|DEFAULT */
+   /* The FP_CONTRACT pragma can be used to allow (if the state is "on") */
+   /* or disallow (if the state is "off") the implementation to contract */
+   /* expressions (refer to section 6.5 of the ISO standard). Each pragma */
+   /* can occur either outside external declarations or preceding all */
+   /* explicit declarations and statements inside a compound statement. */
+   /* When outside external declarations, the pragma takes effect from its */
+   /* occurrence until another FP_CONTRACT pragma is encountered, or until */
+   /* the end of the translation unit. When inside a compound statement, */
+   /* the pragma takes effect from its occurrence until another */
+   /* FP_CONTRACT pragma is encountered (including within a nested */
+   /* compound statement), or until the end of the compound statement; at */
+   /* the end of a compound statement the state for the pragma is restored */
+   /* to its condition just before the compound statement. If this pragma */
+   /* is used in any other context, the behaviour is undefined. The */
+   /* default state for this pragma is implementation-defined ("off" in */
+   /* Norcroft C.). */
 
 #ifdef __cplusplus
 extern "C" {
@@ -173,7 +191,8 @@ double atan2(double /*y*/, double /*x*/);
 float atan2f(float /*y*/, float /*x*/);
 long double atan2l(long double /*y*/, long double /*x*/);
    /* computes the principal value of the arc tangent of y/x, using the */
-   /* signs of both arguments to determine the quadrant of the return value */
+   /* signs of both arguments to determine the quadrant of the return */
+   /* value */
    /* Returns: the arc tangent of y/x, in the range -Pi to Pi. */
 double cos(double /*x*/);
 float cosf(float /*x*/);
@@ -262,9 +281,9 @@ int ilogbf(float /*x*/);
 int ilogbl(long double /*x*/);
    /* extracts the exponent of x as a signed int value. If x is zero it */
    /* computes the value FP_ILOGB0; if x is infinite it computes the value */
-   /* INT_MAX; if x is a NaN it computes the value FP_ILOGBNAN; otherwise it */
-   /* is equivalent to calling the corresponding logb function and casting */
-   /* the returned value to type int. */
+   /* INT_MAX; if x is a NaN it computes the value FP_ILOGBNAN; otherwise */
+   /* it is equivalent to calling the corresponding logb function and */
+   /* casting the returned value to type int. */
    /* Returns: the exponent of x as a signed int value. */
 double ldexp(double /*x*/, int /*exp*/);
 float ldexpf(float /*x*/, int /*exp*/);
@@ -275,14 +294,16 @@ long double ldexpl(long double /*x*/, int /*exp*/);
 double log(double /*x*/);
 float logf(float /*x*/);
 long double logl(long double /*x*/);
-   /* computes the base-e (natural) logarithm of x. A domain error occurs if */
-   /* the argument is negative. A range error occurs if the argument is zero. */
+   /* computes the base-e (natural) logarithm of x. A domain error occurs */
+   /* if the argument is negative. A range error occurs if the argument is */
+   /* zero. */
    /* Returns: the natural logarithm. */
 double log10(double /*x*/);
 float log10f(float /*x*/);
 long double log10l(long double /*x*/);
-   /* computes the base-ten (common) logarithm of x. A domain error occurs if */
-   /* the argument is negative. A range error occurs if the argument is zero. */
+   /* computes the base-ten (common) logarithm of x. A domain error occurs */
+   /* if the argument is negative. A range error occurs if the argument is */
+   /* zero. */
    /* Returns: the base-ten logarithm. */
 double log1p(double /*x*/);
 float log1pf(float /*x*/);
@@ -337,8 +358,8 @@ long double fabsl(long double /*x*/);
 double hypot(double /*x*/, double /*y*/);
 float hypotf(float /*x*/, float /*y*/);
 long double hypotl(long double /*x*/, long double /*y*/);
-   /* computes the square root of the sum of the squares of x and y, without */
-   /* undue overflow or underflow. A range error may occur. */
+   /* computes the square root of the sum of the squares of x and y, */
+   /* without undue overflow or underflow. A range error may occur. */
    /* Returns: sqrt(x^2 + y^2) */
 double pow(double /*x*/, double /*y*/);
 float powf(float /*x*/, float /*y*/);
@@ -386,12 +407,14 @@ double ceil(double /*x*/);
 float ceilf(float /*x*/);
 long double ceill(long double /*x*/);
    /* computes the smallest integer not less than x. */
-   /* Returns: the smallest integer not less than x, expressed as a double. */
+   /* Returns: the smallest integer not less than x, expressed as a */
+   /* floating-point number. */
 double floor(double /*x*/);
 float floorf(float /*x*/);
 long double floorl(long double /*x*/);
    /* computes the largest integer not greater than x. */
-   /* Returns: the largest integer not greater than x, expressed as a double */
+   /* Returns: the largest integer not greater than x, expressed as a */
+   /* floating-point number. */
 #pragma no_side_effects
 double nearbyint(double /*x*/);
 float nearbyintf(float /*x*/);
@@ -409,9 +432,13 @@ long double rintl(long double /*x*/);
 long int lrint(double /*x*/);
 long int lrintf(float /*x*/);
 long int lrintl(long double /*x*/);
+#ifdef __STDC_VERSION__
+#if __STDC_VERSION__ >= 199901
 long long int llrint(double /*x*/);
 long long int llrintf(float /*x*/);
 long long int llrintl(long double /*x*/);
+#endif
+#endif
    /* rounds its argument to an integer value, using the current rounding */
    /* direction. Raises "inexact" if the result differs from the argument. */
    /* Returns: the rounded integer value. */
@@ -426,17 +453,21 @@ long double roundl(long double /*x*/);
 long int lround(double /*x*/);
 long int lroundf(float /*x*/);
 long int lroundl(long double /*x*/);
+#ifdef __STDC_VERSION__
+#if __STDC_VERSION__ >= 199901
 long long int llround(double /*x*/);
 long long int llroundf(float /*x*/);
 long long int llroundl(long double /*x*/);
+#endif
+#endif
    /* rounds its argument to the nearest integer value, rounding halfway */
    /* cases away from zero. */
    /* Returns: the rounded integer value. */
 double trunc(double /*x*/);
 float truncf(float /*x*/);
 long double truncl(long double /*x*/);
-   /* rounds its argument to the integer value, nearest to but no larger in */
-   /* magnitude than the argument. */
+   /* rounds its argument to the integer value, nearest to but no larger */
+   /* in magnitude than the argument. */
    /* Returns: the truncated integer value. */
 double fmod(double /*x*/, double /*y*/);
 float fmodf(float /*x*/, float /*y*/);
@@ -455,10 +486,10 @@ double remquo(double /*x*/, double /*y*/, int * /*quo*/);
 float remquof(float /*x*/, float /*y*/, int * /*quo*/);
 long double remquol(long double /*x*/, long double /*y*/, int * /*quo*/);
    /* compute the same remainder as the remainder functions. In the object */
-   /* pointed to by quo they store a value whose sign is the sign of x/y and */
-   /* whose magnitude is congruent modulo 2^n to the magnitude of the */
-   /* integral quotient of x/y, where n is an implementation-defined integer */
-   /* greater than or equal to 3 (n=31 under RISC OS) */
+   /* pointed to by quo they store a value whose sign is the sign of x/y */
+   /* and whose magnitude is congruent modulo 2^n to the magnitude of the */
+   /* integral quotient of x/y, where n is an implementation-defined */
+   /* integer greater than or equal to 3 (n=31 under RISC OS) */
    /* Returns: x REM y */
 
 #pragma no_side_effects
@@ -481,8 +512,8 @@ long double nextafterl(long double /*x*/, long double /*y*/);
 double nexttoward(double /*x*/, long double /*y*/);
 float nexttowardf(float /*x*/, long double /*y*/);
 long double nexttowardl(long double /*x*/, long double /*y*/);
-   /* equivalent to the nextafter functions except that the second parameter */
-   /* has type long double. */
+   /* equivalent to the nextafter functions except that the second */
+   /* parameter has type long double. */
 double fdim(double /*x*/, double /*y*/);
 float fdimf(float /*x*/, float /*y*/);
 long double fdiml(long double /*x*/, long double /*y*/);
@@ -506,9 +537,9 @@ double fma(double /*x*/, double /*y*/, double /*z*/);
 float fmaf(float /*x*/, float /*y*/, float /*z*/);
 long double fmal(long double /*x*/, long double /*y*/, long double /*z*/);
    /* computes (x*y)+z, rounded as one ternary operation: it computes */
-   /* the value (as if) to infinite precision and rounds once to the result */
-   /* format, according to the rounding mode characterised by the value of */
-   /* FLT_ROUNDS. */
+   /* the value (as if) to infinite precision and rounds once to the */
+   /* result format, according to the rounding mode characterised by the */
+   /* value of FLT_ROUNDS. */
    /* Returns: (x*y)+z, rounded as one ternary operation. */
 
 #ifndef __cplusplus
@@ -525,10 +556,11 @@ long double fmal(long double /*x*/, long double /*y*/, long double /*z*/);
 #endif
 
 /* Some functions can be safely inlined - appropriate macros defined here. */
-/* float ones can't be inlined in C++ because of the odd calling convention */
-/* Some can't be inlined because they won't set errno - in future they */
-/* will be inlined in C99 mode, once the real functions are consistent. */
-/* Define __MATH_FORCE_INLINING if you don't care about errno vs exceptions. */
+/* float ones can't be inlined in C++ because of the odd calling */
+/* convention. Some can't be inlined because they won't set errno - in */
+/* future they will be inlined in C99 mode, once the real functions are */
+/* consistent. Define __MATH_FORCE_INLINING if you don't care about errno */
+/* vs exceptions. */
 double __d_acos(double);
 double __d_asin(double);
 double __d_atan(double);
@@ -538,6 +570,7 @@ double __d_tan(double);
 double __d_log(double);
 double __d_lg10(double);
 double __d_exp(double);
+double __d_abs(double);
 double __d_sqrt(double);
 double __d_floor(double);
 double __d_ceil(double);
@@ -576,6 +609,7 @@ __caller_narrow float __r_tan(float);
 __caller_narrow float __r_log(float);
 __caller_narrow float __r_lg10(float);
 __caller_narrow float __r_exp(float);
+__caller_narrow float __r_abs(float);
 __caller_narrow float __r_sqrt(float);
 __caller_narrow float __r_floor(float);
 __caller_narrow float __r_ceil(float);

clib/h/stdbool

@@ -17,15 +17,11 @@
 
 /* stdbool.h: ISO 'C' (9899:1999) library header, section 7.16 */
 /* Copyright (C) Acorn Computers Ltd. 2002 */
-/* version 1.01 */
+/* version 1.02 */
 
 #ifndef __stdbool_h
 #define __stdbool_h
 
-#if __STDC_VERSION__ < 199901
-#  error <stdbool.h> can only be used in C99
-#endif
-
 #define false 0
 #define true  1
 

clib/h/tgmath

@@ -17,16 +17,34 @@
 
 /* tgmath.h: ISO 'C' (9899:1999) library header, sections 7.22 and G.7 */
 /* Copyright (C) Acorn Computers Ltd. 2003 */
-/* version 1.01 */
+/* version 1.02 */
 
 /*
- * tgmath.h includes the headers <math.h> and <complex.h> and defines
- * several type-generic macros.
+ * tgmath.h includes the headers <math.h> and <complex.h> and defines several
+ * type-generic macros.
+ *
+ * Of the <math.h> and <complex.h> functions without an f or l suffix, several
+ * have one or more parameters whose corresponding real type is double. For each
+ * such function, except modf, there is a corresponding type-generic macro. The
+ * parameters whose corresponding real type is double in the function synopsis
+ * are generic parameters. Use of the macro invokes a function whose
+ * corresponding real type and domain are determined by the arguments for the
+ * generic parameters.
+ *
+ * Use of the macro invokes a function whose generic parameters have the
+ * corresponding real type as follows:
+ *
+ *   - First, if any argument for generic parameters has type long double,
+ *     the type determined is long double.
+ *   - Otherwise, if any argument for generic parameters has type double or
+ *     is of integer type, the type determined is double.
+ *   - Otherwise, the type determined is float.
  */
 #ifndef __tgmath_h
 #define __tgmath_h
 
 #include <math.h>
+#include <complex.h>
 
 /* Internal macros to simplify public macros */
 #define _TG(o, ...)      (___select(0, o##f, o, o##l, __VA_ARGS__))(__VA_ARGS__)
@@ -35,11 +53,54 @@
 
 #define _TGx(f,d,l, ...) (___select(0, f, d, l, __VA_ARGS__))(__VA_ARGS__)
 
-#define _TGC(o, ...)     _TG(o, __VA_ARGS__)
+#define _TGZ(r,i,c, ...) (___select(1, r##f, r, r##l, \
+                                       i##f, i, i##l, \
+                                       c##f, c, c##l, __VA_ARGS__))(__VA_ARGS__)
+
+#define _TGC(o, ...)     _TGZ(o,   c##o, c##o, __VA_ARGS__)
+#define _TGI(o, ...)     _TGZ(o, __i##o, c##o, __VA_ARGS__)
+
+
+/* Private inline functions to provide the special imaginary forms as
+ * per G.7. This could be done using an expression inside ___select,
+ * but that leads to explosive macro expansion with nested macros.
+ * Inlining will occur for these as per <math.h>.
+ */
+#define _DEF_IFN1(fn, val, t, im) \
+  static inline __caller_narrow t im __i##fn(t imaginary z) { return val(z/I); }
+
+#define _DEF_IFN(im, fn, val) \
+    _DEF_IFN1(fn,    val,    double, im) \
+    _DEF_IFN1(fn##f, val##f, float, im) \
+    _DEF_IFN1(fn##l, val##l, long double, im)
 
-#define _TGI(o,i, z)     _TG(o, z)
+_DEF_IFN(imaginary, asin,  I*asinh)
+_DEF_IFN(imaginary, atan,  I*atanh)
+_DEF_IFN(imaginary, asinh, I*asin)
+_DEF_IFN(imaginary, atanh, I*atan)
+_DEF_IFN(         , cos,   cosh)
+_DEF_IFN(imaginary, sin,   I*sinh)
+_DEF_IFN(imaginary, tan,   I*tanh)
+_DEF_IFN(         , cosh,  cos)
+_DEF_IFN(imaginary, sinh,  I*sin)
+_DEF_IFN(imaginary, tanh,  I*tan)
+_DEF_IFN(         , abs,   fabs)
 
-/* Type-generic macros for functions with both real and complex forms */
+#undef _DEF_IFN
+#undef _DEF_IFN1
+
+/*
+ * For each unsuffixed function in <math.h> for which there is a function in
+ * <complex.h> with the same name except for a c prefix, the corresponding
+ * type-generic macro (for both functions) has the same name as the function in
+ * <math.h>. The corresponding type-generic macro for fabs and cabs is fabs.
+ *
+ * If at least one argument for a generic parameter is complex, then use of
+ * the macro invokes a complex function; otherwise, if an argument is imaginary,
+ * the macro expands to an expression whose type is real, imaginary or complex,
+ * as appropriate for the particular function; otherwise, use of the macro
+ * invokes a real function.
+ */
 #undef acos
 #undef asin
 #undef atan
@@ -59,24 +120,32 @@
 #undef fabs
 
 #define acos(z)  _TGC(acos, z)
-#define asin(z)  _TGI(asin, I*asinh, z)
-#define atan(z)  _TGI(atan, I*atanh, z)
+#define asin(z)  _TGI(asin, z)
+#define atan(z)  _TGI(atan, z)
 #define acosh(z) _TGC(acosh, z)
-#define asinh(z) _TGI(asinh, I*asin, z)
-#define atanh(z) _TGI(atanh, I*atan, z)
-#define cos(z)   _TGI(cos, cosh, z)
-#define sin(z)   _TGI(sin, I*sinh, z)
-#define tan(z)   _TGI(tan, I*tanh, z)
-#define cosh(z)  _TGI(cosh, cos, z)
-#define sinh(z)  _TGI(sinh, I*sin, z)
-#define tanh(z)  _TGI(tanh, I*tan, z)
+#define asinh(z) _TGI(asinh, z)
+#define atanh(z) _TGI(atanh, z)
+#define cos(z)   _TGI(cos, z)
+#define sin(z)   _TGI(sin, z)
+#define tan(z)   _TGI(tan, z)
+#define cosh(z)  _TGI(cosh, z)
+#define sinh(z)  _TGI(sinh, z)
+#define tanh(z)  _TGI(tanh, z)
 #define exp(z)   _TGC(exp, z)
 #define log(z)   _TGC(log, z)
 #define pow(x,y) _TGC(pow, x,y)
 #define sqrt(z)  _TGC(sqrt, z)
-#define fabs(z)  _TG(fabs, z)
+#define fabs(z)  _TGZ(fabs,__iabs,cabs, z)
 
-/* Type-generic macros for real-only functions */
+/*
+ * For each unsuffixed function in <math.h> without a c-prefixed counterpart in
+ * <complex.h>, the corresponding type-generic macro has the same name as the
+ * function.
+ *
+ * If all arguments for generic parameters are real, then use of the macro
+ * invokes a real function; otherwise, use of the macro results in undefined
+ * behaviour.
+ */
 #undef atan2
 #undef cbrt
 #undef ceil
@@ -153,6 +222,26 @@
 #define tgamma(x)       _TG(tgamma, x)
 #define trunc(x)        _TG(trunc, x)
 
+/*
+ * For each unsuffixed function in <complex.h> that is not a c-prefixed
+ * counterpart to a function in <math.h>, the corresponding type-generic macro
+ * has the same name as the function.
+ *
+ * Use of the macro with any real, imaginary or complex argument invokes a
+ * complex function.
+ */
+#undef carg
+#undef cimag
+#undef conj
+#undef cproj
+#undef creal
+
+#define carg(z)         _TG(carg, z)
+#define cimag(z)        _TG(cimag, z)
+#define conj(z)         _TG(conj, z)
+#define cproj(z)        _TG(cproj, z)
+#define creal(z)        _TG(creal, z)
+
 /* Ensure inlining still happens, unless requested otherwise */
 #ifndef __TGMATH_NO_INLINING
 #undef atan
@@ -163,16 +252,75 @@
 #undef rint
 #undef lrint
 #undef remainder
-#undef fma
-#define atan(x)  _TGx(__r_atan, __d_atan, atanl, x)
-#define fabs(x)  ___select(0, __abs (x), __abs (double) (x), __abs (x), x)
+#undef cimag
+#undef conj
+#undef creal
+#define atan(z)  (___select(1, __r_atan, __d_atan, atanl, \
+                               __iatanf, __iatan, __iatanl,\
+                               catanf, catan, catanl, z))(z)
+#define fabs(z)  (___select(1, __r_abs, __d_abs, fabsl, \
+                               __iabsf, __iabs, __iabsl,\
+                               cabsf, cabs, cabsl, z))(z)
 #define floor(x) _TGx(__r_floor, __d_floor, floorl, x)
 #define ceil(x)  _TGx(__r_ceil, __d_ceil, ceill, x)
 #define trunc(x) _TGx(__r_trunc, __d_trunc, truncl, x)
 #define rint(x)  _TGx(__r_rint, __d_rint, rintl, x)
 #define lrint(x) _TGx(__r_lrint, __d_lrint, lrintl, x)
 #define remainder(x,y) _TGx(__r_rem, __d_rem, remainderl, x,y)
-#define fma(x,y,z) _TGx(__r_fma, fma, fmal, x,y,z)
+#define cimag(z) ___select(0, (float _Imaginary)(z)/_Imaginary_I, \
+                              (double _Imaginary)(z)/_Imaginary_I, \
+                              (long double _Imaginary)(z)/_Imaginary_I, z)
+#define conj(z)  ___select(0, __conj (float _Complex) (z), \
+                              __conj (double _Complex) (z), \
+                              __conj (long double _Complex) (z), z)
+#define creal(z) ___select(0, (float)(z), \
+                              (double)(z), \
+                              (long double)(z), z)
+#ifdef __MATH_FORCE_INLINING
+#undef acos
+#undef asin
+#undef cos
+#undef sin
+#undef tan
+#undef sqrt
+#undef log
+#undef exp
+#undef log10
+#define acos(z)  (___select(1, __r_acos, __d_acos, acosl, \
+                               __iacosf, __iacos, __iacosl,\
+                               cacosf, cacos, cacosl, z))(z)
+#define asin(z)  (___select(1, __r_asin, __d_asin, asinl, \
+                               __iasinf, __iasin, __iasinl,\
+                               casinf, casin, casinl, z))(z)
+#define cos(z)   (___select(1, __r_cos, __d_cos, cosl, \
+                               __icosf, __icos, __icosl,\
+                               ccosf, ccos, ccosl, z))(z)
+#define sin(z)   (___select(1, __r_sin, __d_sin, sinl, \
+                               __isinf, __isin, __isinl,\
+                               csinf, csin, csinl, z))(z)
+#define tan(z)   (___select(1, __r_tan, __d_tan, tanl, \
+                               __itanf, __itan, __itanl,\
+                               ctanf, ctan, ctanl, z))(z)
+#define sqrt(z)  (___select(1, __r_sqrt, __d_sqrt, sqrtl, \
+                               csqrtf, csqrt, csqrtl,\
+                               csqrtf, csqrt, csqrtl, z))(z)
+#define log(z)   (___select(1, __r_log, __d_log, logl, \
+                                clogf, clog, clogl,\
+                                clogf, clog, clogl, z))(z)
+#define exp(z)   (___select(1, __r_exp, __d_exp, expl, \
+                                cexpf, cexp, cexpl,\
+                                cexpf, cexp, cexpl, z))(z)
+#define log10(z)  _TGx(__r_lg10, __d_lg10, log10l, x)
+#endif
+#define cimag(z) ___select(0, (float _Imaginary)(z)/_Imaginary_I, \
+                              (double _Imaginary)(z)/_Imaginary_I, \
+                              (long double _Imaginary)(z)/_Imaginary_I, z)
+#define conj(z)  ___select(0, __conj (float _Complex) (z), \
+                              __conj (double _Complex) (z), \
+                              __conj (long double _Complex) (z), z)
+#define creal(z) ___select(0, (float)(z), \
+                              (double)(z), \
+                              (long double)(z), z)
 #endif
 
 #endif

clib/s/cl_entry2

@@ -153,4 +153,56 @@
         Entry   llround, imported, , unveneered
         Entry   llroundf, imported, , unveneered
 
+        Entry   _cxd_mul, imported, , unveneered
+        Entry   _cxf_mul, imported, , unveneered
+        Entry   _cxd_div, imported, , unveneered
+        Entry   _cxf_div, imported, , unveneered
+        Entry   _cxd_rdv, imported, , unveneered
+        Entry   _cxf_rdv, imported, , unveneered
+
+        Entry   cacos, imported, , unveneered
+        Entry   cacosf, imported, , unveneered
+        Entry   casin, imported, , unveneered
+        Entry   casinf, imported, , unveneered
+        Entry   catan, imported, , unveneered
+        Entry   catanf, imported, , unveneered
+        Entry   ccos, imported, , unveneered
+        Entry   ccosf, imported, , unveneered
+        Entry   csin, imported, , unveneered
+        Entry   csinf, imported, , unveneered
+        Entry   ctan, imported, , unveneered
+        Entry   ctanf, imported, , unveneered
+        Entry   cacosh, imported, , unveneered
+        Entry   cacoshf, imported, , unveneered
+        Entry   casinh, imported, , unveneered
+        Entry   casinhf, imported, , unveneered
+        Entry   catanh, imported, , unveneered
+        Entry   catanhf, imported, , unveneered
+        Entry   ccosh, imported, , unveneered
+        Entry   ccoshf, imported, , unveneered
+        Entry   csinh, imported, , unveneered
+        Entry   csinhf, imported, , unveneered
+        Entry   ctanh, imported, , unveneered
+        Entry   ctanhf, imported, , unveneered
+        Entry   cexp, imported, , unveneered
+        Entry   cexpf, imported, , unveneered
+        Entry   clog, imported, , unveneered
+        Entry   clogf, imported, , unveneered
+        Entry   cabs, imported, , unveneered
+        Entry   cabsf, imported, , unveneered
+        Entry   cpow, imported, , unveneered
+        Entry   cpowf, imported, , unveneered
+        Entry   csqrt, imported, , unveneered
+        Entry   csqrtf, imported, , unveneered
+        Entry   carg, imported, , unveneered
+        Entry   cargf, imported, , unveneered
+        Entry   cimag, imported, , unveneered
+        Entry   cimagf, imported, , unveneered
+        Entry   conj, imported, , unveneered
+        Entry   conjf, imported, , unveneered
+        Entry   cproj, imported, , unveneered
+        Entry   cprojf, imported, , unveneered
+        Entry   creal, imported, , unveneered
+        Entry   crealf, imported, , unveneered
+
         END

s/cxsupport

+; Copyright 2005 Castle Technology Ltd
+;
+; 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.
+;
+        GET objmacs.s
+
+        CodeArea
+
+        EXPORT  _cxd_mul
+        EXPORT  _cxf_mul
+        EXPORT  _cxd_div
+        EXPORT  _cxf_div
+        EXPORT  _cxd_rdv
+        EXPORT  _cxf_rdv
+
+        EXPORT  __log_cabs
+
+        EXPORT  cabs
+        EXPORT  cabsf
+        EXPORT  csqrt
+        EXPORT  csqrtf
+
+        EXPORT  carg
+        EXPORT  cargf
+        EXPORT  cimag
+        EXPORT  cimagf
+        EXPORT  conj
+        EXPORT  conjf
+        EXPORT  cproj
+        EXPORT  cprojf
+        EXPORT  creal
+        EXPORT  crealf
+
+        IMPORT  hypot
+        IMPORT  hypotf
+        IMPORT  _new_atan2
+        IMPORT  atan2f
+
+FP_ZERO * 0
+FP_SUBNORMAL * 1
+FP_NORMAL * 2
+FP_INFINITE * 3
+FP_NAN * 4
+
+|__fpclassifyl|
+        STFE    f0,[sp,#-12]!
+        LDMFD   sp!,{a1-a3}
+        MOVS    a4,a1,LSL #17
+        BEQ     %FT02
+        ADDS    a4,a4,#1:SHL:17
+        MOVCC   a1,#FP_NORMAL
+        MOVCC   pc,lr
+        ORRS    a4,a3,a2,LSL #1
+        MOVEQ   a1,#FP_INFINITE
+        MOVNE   a1,#FP_NAN
+        MOV     pc,lr
+
+02      ORRS    a4,a2,a3
+        MOVEQ   a1,#FP_ZERO
+        MOVEQ   pc,lr
+        TEQ     a2,#0
+        MOVMI   a1,#FP_NORMAL
+        MOVPL   a1,#FP_SUBNORMAL
+        MOV     pc,lr
+
+copysignl
+        STFE    f1,[sp,#-12]!
+        STFE    f0,[sp,#-12]!
+        LDMFD   sp,{a1-a4}
+        BIC     a1,a1,#&80000000
+        AND     a4,a4,#&80000000
+        ORR     a1,a1,a4
+        STR     a1,[sp,#0]
+        LDFE    f0,[sp],#24
+        MOV     pc,lr
+
+; (f4+i*f5) * (f6+i*f7) -> (f0+i*f1) in full extended precision
+_cx_mul
+        MUFE    f0,f4,f6
+        MUFE    f1,f5,f7
+        SUFE    f0,f0,f1
+        MUFE    f1,f5,f6
+        MUFE    f2,f4,f7
+        ADFE    f1,f1,f2
+        CMF     f0,f0           ; if not (NaN+i*NaN) return now
+        CMFVS   f1,f1
+        MOVVC   pc,lr
+
+        ; Recovery mechanism from Annex G
+
+        STMFD   sp!,{v1-v5,lr}
+        MOV     v5,#0           ; v5 = recalc flag
+        SFMFD   f0,2,[sp]!      ; remember original attempt
+
+        MVFE    f0,f4           ; v1-v4 = classification of f4-f7
+        BL      __fpclassifyl
+        MOV     v1,a1
+        MVFE    f0,f5
+        BL      __fpclassifyl
+        MOV     v2,a1
+        MVFE    f0,f6
+        BL      __fpclassifyl
+        MOV     v3,a1
+        MVFE    f0,f7
+        BL      __fpclassifyl
+        MOV     v4,a1
+
+        TEQ     v1,#FP_INFINITE
+        TEQNE   v2,#FP_INFINITE
+        BNE     %FT20
+
+        ; (f4+i*f5) is infinite - box the infinity and remove NaNs from other operand
+        TEQ     v1,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f4
+        BL      copysignl
+        MVFE    f4,f0
+
+        TEQ     v2,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f5
+        BL      copysignl
+        MVFE    f5,f0
+
+        TEQ     v3,#FP_NAN
+        BNE     %FT01
+        MVFE    f0,#0
+        MVFE    f1,f6
+        BL      copysignl
+        MVFE    f6,f0
+01
+        TEQ     v4,#FP_NAN
+        BNE     %FT01
+        MVFE    f0,#0
+        MVFE    f1,f7
+        BL      copysignl
+        MVFE    f7,f0
+01
+        MOV     v5,#1
+
+20      TEQ     v3,#FP_INFINITE
+        TEQNE   v4,#FP_INFINITE
+        BNE     %FT40
+
+        ; (f6+i*f7) is infinite - box the infinity and remove NaNs from other operand
+        TEQ     v3,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f6
+        BL      copysignl
+        MVFE    f6,f0
+
+        TEQ     v4,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f7
+        BL      copysignl
+        MVFE    f7,f0
+
+        TEQ     v1,#FP_NAN
+        BNE     %FT01
+        MVFE    f0,#0
+        MVFE    f1,f4
+        BL      copysignl
+        MVFE    f4,f0
+01
+        TEQ     v2,#FP_NAN
+        BNE     %FT01
+        MVFE    f0,#0
+        MVFE    f1,f5
+        BL      copysignl
+        MVFE    f5,f0
+01
+        MOV     v5,#1
+
+40      ; Annex G has guff to clear up NaNs generated by overflow
+        ; (eg big * big - big * big -> inf - inf -> NaN), but that
+        ; can't happen here because we know arguments were no larger
+        ; than double.
+
+        TEQ     v5,#1
+        LFMNEFD f0,2,[sp]!
+        LDMNEFD sp!,{v1-v5,pc}
+
+        LDFE    f3,ExtInf
+        MUFE    f0,f4,f6
+        MUFE    f1,f5,f7
+        SUFE    f0,f0,f1
+        MUFE    f1,f5,f6
+        MUFE    f2,f4,f7
+        ADFE    f1,f1,f2
+        MUFE    f0,f0,f3
+        MUFE    f1,f1,f3
+        ADD     sp,sp,#12*2
+        LDMFD   sp!,{v1-v5,pc}
+
+; (f4+i*f5) / (f6+i*f7) -> (f0+i*f1) in full extended precision
+_cx_div
+        MUFE    f0,f6,f6
+        MUFE    f1,f7,f7
+        ADFE    f3,f0,f1        ; f3 = denom
+        MUFE    f0,f4,f6
+        MUFE    f1,f5,f7
+        ADFE    f0,f0,f1
+        DVFE    f0,f0,f3
+        MUFE    f1,f5,f6
+        MUFE    f2,f4,f7
+        SUFE    f1,f1,f2
+        DVFE    f1,f1,f3
+        CMF     f0,f0
+        CMFVS   f1,f1
+        MOVVC   pc,lr
+
+        ; Based on Annex G implementation again. Recover infinities and zeros
+        ; from NaN+i*NaN results. Three cases...
+
+        STMFD   sp!,{v1-v4,lr}
+        MVFE    f0,f4           ; v1-v4 = classification of f4-f7
+        BL      __fpclassifyl
+        MOV     v1,a1
+        MVFE    f0,f5
+        BL      __fpclassifyl
+        MOV     v2,a1
+        MVFE    f0,f6
+        BL      __fpclassifyl
+        MOV     v3,a1
+        MVFE    f0,f7
+        BL      __fpclassifyl
+        MOV     v4,a1
+
+        CMF     f3,#0
+        BNE     %FT20
+        TEQ     v1,#FP_NAN
+        TEQEQ   v2,#FP_NAN
+        BEQ     %FT20
+
+        ; Case 1: non-zero / zero
+
+        LDFE    f0,ExtInf
+        MVFE    f1,f4
+        BL      copysignl
+        MUFE    f1,f0,f5
+        MUFE    f0,f0,f4
+        LDMFD   sp!,{v1-v4,pc}
+
+20      TEQ     v1,#FP_INFINITE
+        TEQNE   v2,#FP_INFINITE
+        BNE     %FT40
+        CMP     v3,#FP_NORMAL
+        CMPLS   v4,#FP_NORMAL
+        BHI     %FT40
+
+        ; Case 2: infinity / finite
+
+        TEQ     v1,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f4
+        BL      copysignl
+        MVFE    f4,f0
+
+        TEQ     v2,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f5
+        BL      copysignl
+        MVFE    f5,f0
+
+        LDFE    f3,ExtInf
+        MUFE    f0,f4,f6
+        MUFE    f1,f5,f7
+        ADFE    f0,f0,f1
+        MUFE    f0,f0,f3
+        MUFE    f1,f5,f6
+        MUFE    f2,f4,f7
+        SUFE    f1,f1,f2
+        MUFE    f1,f1,f3
+        LDMFD   sp!,{v1-v4,pc}
+
+40      TEQ     v3,#FP_INFINITE
+        TEQNE   v4,#FP_INFINITE
+        BNE     %FT60
+        CMP     v1,#FP_NORMAL
+        CMPLS   v2,#FP_NORMAL
+        BHI     %FT60
+
+        ; Case 3: finite / infinity
+
+        TEQ     v3,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f6
+        BL      copysignl
+        MVFE    f6,f0
+
+        TEQ     v4,#FP_INFINITE
+        MVFEQE  f0,#1
+        MVFNEE  f0,#0
+        MVFE    f1,f7
+        BL      copysignl
+        MVFE    f7,f0
+
+        MUFE    f0,f4,f6
+        MUFE    f1,f5,f7
+        ADFE    f0,f0,f1
+        MUFE    f0,f0,#0
+        MUFE    f1,f5,f6
+        MUFE    f2,f4,f7
+        SUFE    f1,f1,f2
+        MUFE    f1,f1,#0
+60      LDMFD   sp!,{v1-v4,pc}
+
+
+ExtInf  DCD     &00007FFF, &00000000, &00000000
+
+
+_cxd_mul
+        SFMFD   f4,4,[sp]!
+        MVFD    f4,f0
+        MVFD    f5,f1
+        MVFD    f6,f2
+        MVFD    f7,f3
+        MOV     ip,lr
+        BL      _cx_mul
+        MVFD    f0,f0
+        MVFD    f1,f1
+        LFMFD   f4,4,[sp]!
+        MOV     pc,ip
+
+_cxf_mul
+        SFMFD   f4,4,[sp]!
+        MVFS    f4,f0
+        MVFS    f5,f1
+        MVFS    f6,f2
+        MVFS    f7,f3
+        MOV     ip,lr
+        BL      _cx_mul
+        MVFS    f0,f0
+        MVFS    f1,f1
+        LFMFD   f4,4,[sp]!
+        MOV     pc,ip
+
+_cxd_div
+        SFMFD   f4,4,[sp]!
+        MVFD    f4,f0
+        MVFD    f5,f1
+        MVFD    f6,f2
+        MVFD    f7,f3
+_cxd_div1
+        MOV     ip,lr
+        BL      _cx_div
+        MVFD    f0,f0
+        MVFD    f1,f1
+        LFMFD   f4,4,[sp]!
+        MOV     pc,ip
+
+_cxd_rdv
+        SFMFD   f4,4,[sp]!
+        MVFD    f4,f2
+        MVFD    f5,f3
+        MVFD    f6,f0
+        MVFD    f7,f1
+        B       _cxd_div1
+
+_cxf_div
+        SFMFD   f4,4,[sp]!
+        MVFS    f4,f0
+        MVFS    f5,f1
+        MVFS    f6,f2
+        MVFS    f7,f3
+_cxf_div1
+        MOV     ip,lr
+        BL      _cx_div
+        MVFS    f0,f0
+        MVFS    f1,f1
+        LFMFD   f4,4,[sp]!
+        MOV     pc,ip
+
+_cxf_rdv
+        SFMFD   f4,4,[sp]!
+        MVFS    f4,f2
+        MVFS    f5,f3
+        MVFS    f6,f0
+        MVFS    f7,f1
+        B       _cxf_div1
+
+; calculate log(cabs(z)). Use extended intermediates to avoid over/underflow
+
+__log_cabs
+        BIC     a1,a1,#&80000000        ; Fast way to do fabs
+        BIC     a3,a3,#&80000000
+        STMFD   sp!,{a1-a4}
+        LDFD    f2,[sp],#8
+        LDFD    f3,[sp],#8
+        CMF     f2,f3
+        MUFVCE  f0,f2,f2
+        MUFVCE  f1,f3,f3
+        ADFVCE  f0,f0,f1
+        LGNVCD  f0,f0
+        MUFVCD  f0,f0,#0.5
+        MOVVC   pc,lr
+; cabs(nan,inf)==cabs(inf,nan)==inf; otherwise cabs(nan,x)=nan
+        LDFD    f0,DblInf
+        CMF     f2,f0
+        CMFNE   f3,f0
+        ADFNED  f0,f2,f3
+        MOV     pc,lr
+
+cabs
+        B       hypot
+
+cabsf
+        STMFD   sp!,{a1-a2}
+        LDFS    f0,[sp],#4
+        LDFS    f1,[sp],#4
+        STFD    f0,[sp,#-16]!
+        STFD    f1,[sp,#8]
+        LDMFD   sp!,{a1-a4}
+        B       hypotf
+
+csqrtx
+        CMF     f1,#0
+        BMI     csqrt_negi
+        BEQ     csqrt_real
+        BVS     csqrt_nani
+csqrt_posi
+        LDFE    f3,ExtInf
+        CMF     f1,f3
+        BEQ     csqrt_infi
+        MUFE    f2,f0,f0
+        MUFE    f3,f1,f1
+        ADFE    f2,f2,f3
+        SQTE    f2,f2           ; f2 = |z|
+        ABSE    f3,f0
+        ADFE    f2,f3,f2
+        MUFE    f2,f2,#0.5
+        SQTE    f2,f2           ; f2 = w = sqrt((|x|+|z|)/2)
+        CMF     f0,#0
+        BMI     csqrt_negr
+        MVFE    f0,f2
+        DVFE    f1,f1,f2
+        MUFE    f1,f1,#0.5
+        MOV     pc,lr
+csqrt_negr
+        DVFE    f0,f1,f2
+        MUFE    f0,f0,#0.5
+        MVFE    f1,f2
+        CMF     f0,f1
+        MOV     pc,lr
+
+csqrt_infi
+        MVFE    f0,f1
+        MOV     pc,lr
+
+csqrt_negi
+        MOV     a4,lr
+        MNFE    f1,f1
+        BL      csqrt_posi
+        MNFE    f1,f1
+        MOV     pc,a4
+
+csqrt_real
+        CMF     f0,#0
+        BMI     csqrt_negreal
+        BEQ     csqrt_zero
+        MVFVSE  f1,f0           ; NaN propagation
+        SQTE    f0,f0
+        MOV     pc,lr
+
+csqrt_negreal
+        STFE    f1,[sp,#-12]!
+        LDR     a3,[sp],#12
+        MNFE    f0,f0
+        SQTE    f1,f0
+        MVFE    f0,#0
+        TEQ     a3,#0
+        MNFMIE  f1,f1
+        MOV     pc,lr
+
+csqrt_zero
+        MVFE    f0,#0
+        MOV     pc,lr
+
+csqrt_nani
+        LDFE    f3,ExtInf
+        CMF     f0,f3
+        MOVEQ   pc,lr
+        CNF     f0,f3
+        ABSNEE  f0,f1
+        MOVNE   pc,lr
+        ; csqrt(-inf + I*NaN) - muck around for conj() identity
+        STFE    f1,[sp,#-12]!
+        LDR     a3,[sp],#12
+        TEQ     a3,#0
+        ABSE    f0,f1
+        MVFPLE  f1,f3
+        MNFMIE  f1,f3
+        MOV     pc,lr
+
+
+; (pos+I*0) -> (sqrt(pos)+I*0)
+; (+0+I*0)  -> (+0+I*0)
+; (-0+I*0)  -> (+0+I*0)
+; (NAN+I*0) -> (NaN+I*NaN)
+; (neg+I*0) -> (0+I*sqrt(-neg))
+
+csqrt
+        STMFD   sp!,{a1-a4}
+        LDFD    f0,[sp],#8
+        LDFD    f1,[sp],#8
+        MOV     ip,lr
+        BL      csqrtx
+        MVFD    f0,f0
+        MVFD    f1,f1
+        MOV     pc,ip
+
+csqrtf
+        STMFD   sp!,{a1-a2}
+        LDFS    f0,[sp],#4
+        LDFS    f1,[sp],#4
+        MOV     ip,lr
+        BL      csqrtx
+        MVFS    f0,f0
+        MVFS    f1,f1
+        MOV     pc,ip
+
+carg
+        MOV     ip,a1
+        MOV     a1,a3
+        MOV     a3,ip
+        MOV     ip,a2
+        MOV     a2,a4
+        MOV     a4,ip
+        B       _new_atan2
+
+cargf
+        STMFD   sp!,{a1-a2}
+        LDFS    f0,[sp],#4
+        LDFS    f1,[sp],#4
+        STFD    f0,[sp,#-8]!
+        STFD    f1,[sp,#-8]!
+        LDMIA   sp!,{a1-a4}
+        B       atan2f
+
+cimag
+        STMFD   sp!,{a3-a4}
+        LDFD    f0,[sp],#8
+        MOV     pc,lr
+
+cimagf
+        STR     a2,[sp,#-4]!
+        LDFS    f0,[sp],#4
+        MOV     pc,lr
+
+conj
+        EOR     a3,a3,#&80000000
+        STMFD   sp!,{a1-a4}
+        LDFD    f0,[sp],#8
+        LDFD    f1,[sp],#8
+        MOV     pc,lr
+
+conjf
+        EOR     a2,a2,#&80000000
+        STMFD   sp!,{a1-a2}
+        LDFS    f0,[sp],#4
+        LDFS    f1,[sp],#4
+        MOV     pc,lr
+
+DblInf
+        DCD     &7FF00000,&00000000
+cproj
+        LDFD    f2,DblInf
+        ABSD    f3,f0
+        CMF     f3,f2
+        ABSNED  f3,f1
+        CMFNE   f3,f2
+        MOVNE   pc,lr
+        MVFD    f0,f2
+        STFD    f1,[sp,#-8]!
+        LDMIA   sp!,{a1,a2}
+        TST     a1,#0
+        MVFPLD  f1,#0
+        MNFMID  f1,#0
+        MOV     pc,lr
+
+SglInf
+        DCD     &7F800000
+cprojf
+        LDFS    f2,SglInf
+        ABSS    f3,f0
+        CMF     f3,f2
+        ABSNES  f3,f1
+        CMFNE   f3,f2
+        MOVNE   pc,lr
+        MVFS    f0,f2
+        STFS    f1,[sp,#-4]!
+        LDR     a1,[sp],#4
+        TST     a1,#0
+        MVFPLS  f1,#0
+        MNFMIS  f1,#0
+        MOV     pc,lr
+
+creal
+        STMFD   sp!,{a1-a2}
+        LDFD    f0,[sp],#8
+        MOV     pc,lr
+
+crealf
+        STR     a1,[sp,#-4]!
+        LDFS    f0,[sp],#4
+        MOV     pc,lr
+
+        END

s/mathasm

@@ -22,7 +22,7 @@ FloatingPointArgsInRegs SETL {FALSE}
 
         CodeArea
 
-;        EXPORT  atan2
+        EXPORT  _new_atan2
         EXPORT  atan2f
         EXPORT  fma
         EXPORT  fmaf
@@ -499,7 +499,7 @@ $name   Atan2 $p
         Return  ,,LinkNotStacked
         MEND
 
-;atan2   Atan2   D
+_new_atan2   Atan2   D
 atan2f  Atan2   S
 
         END

s/mathl

@@ -12,9 +12,10 @@
 ; See the License for the specific language governing permissions and
 ; limitations under the License.
 ;
-; As long double is the same as double, all the long double <math.h> functions
-; can just bounce to the double ones. These bounces are included in ANSILib
-; and the stubs as is - the Shared C Library doesn't contain or export them.
+; As long double is the same as double, all the long double <math.h> and
+; <complex.h> functions can just bounce to the double ones. These bounces are
+; included in ANSILib and the stubs as is - the Shared C Library doesn't
+; contain or export them.
 
         AREA    |C$$code|, CODE, READONLY
 
@@ -92,6 +93,33 @@ $func.l B       $func
 
         DoL     fma
 
+        DoL     cacos
+        DoL     casin
+        DoL     catan
+        DoL     ccos
+        DoL     csin
+        DoL     ctan
+
+        DoL     cacosh
+        DoL     casinh
+        DoL     catanh
+        DoL     ccosh
+        DoL     csinh
+        DoL     ctanh
+
+        DoL     cexp
+        DoL     clog
+
+        DoL     cabs
+        DoL     cpow
+        DoL     csqrt
+
+        DoL     carg
+        DoL     cimag
+        DoL     conj
+        DoL     cproj
+        DoL     creal
+
         EXPORT  strtold
         IMPORT  strtod
 strtold B       strtod