29 files changed, 322 insertions, 318 deletions
diff --git a/Makefile b/Makefile
index a9c1128..7336a0e 100644
--- a/Makefile
+++ b/Makefile
@@ -21,9 +21,7 @@ bindir = $(prefix)/bin
 libdir = $(prefix)/lib
 includedir = $(prefix)/include
 
-HACK = $(srcdir)/math/s_sincosf.c
-
-MATH_SRCS = $(filter-out $(HACK),$(wildcard $(srcdir)/math/*.[cS]))
+MATH_SRCS = $(wildcard $(srcdir)/math/*.[cS])
 MATH_BASE = $(basename $(MATH_SRCS))
 MATH_OBJS = $(MATH_BASE:$(srcdir)/%=build/%.o)
 RTEST_SRCS = $(wildcard $(srcdir)/test/rtest/*.[cS])
diff --git a/math/cosf.c b/math/cosf.c
new file mode 100644
index 0000000..6b5284c
--- /dev/null
+++ b/math/cosf.c
@@ -0,0 +1 @@
+#include "single/s_cosf.c"
diff --git a/math/e_exp2f.c b/math/exp2f.c
index fca66fe..fca66fe 100644
--- a/math/e_exp2f.c
+++ b/math/exp2f.c
diff --git a/math/e_exp2f_data.c b/math/exp2f_data.c
index a74844b..a74844b 100644
--- a/math/e_exp2f_data.c
+++ b/math/exp2f_data.c
diff --git a/math/e_expf.c b/math/expf.c
index 6b81310..6b81310 100644
--- a/math/e_expf.c
+++ b/math/expf.c
diff --git a/math/funder.c b/math/funder.c
index 5e5f6b2..6208e12 100644
--- a/math/funder.c
+++ b/math/funder.c
@@ -1,63 +1 @@
-/*
- * funder.c - manually provoke SP exceptions for mathlib
- *
- * Copyright (c) 2009-2015, Arm Limited.
- * SPDX-License-Identifier: Apache-2.0
- *
- * 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 "math_private.h"
-#include <fenv.h>
-
-__inline float __mathlib_flt_infnan2(float x, float y)
-{
-  return x+y;
-}
-
-__inline float __mathlib_flt_infnan(float x)
-{
-  return x+x;
-}
-
-float __mathlib_flt_underflow(void)
-{
-#ifdef CLANG_EXCEPTIONS
-  feraiseexcept(FE_UNDERFLOW);
-#endif
-  return 0x1p-95F * 0x1p-95F;
-}
-
-float __mathlib_flt_overflow(void)
-{
-#ifdef CLANG_EXCEPTIONS
-  feraiseexcept(FE_OVERFLOW);
-#endif
-  return 0x1p+97F * 0x1p+97F;
-}
-
-float __mathlib_flt_invalid(void)
-{
-#ifdef CLANG_EXCEPTIONS
-  feraiseexcept(FE_INVALID);
-#endif
-  return 0.0f / 0.0f;
-}
-
-float __mathlib_flt_divzero(void)
-{
-#ifdef CLANG_EXCEPTIONS
-  feraiseexcept(FE_DIVBYZERO);
-#endif
-  return 1.0f / 0.0f;
-}
+#include "single/funder.c"
diff --git a/math/e_log2f.c b/math/log2f.c
index dab7005..dab7005 100644
--- a/math/e_log2f.c
+++ b/math/log2f.c
diff --git a/math/e_log2f_data.c b/math/log2f_data.c
index c68c670..c68c670 100644
--- a/math/e_log2f_data.c
+++ b/math/log2f_data.c
diff --git a/math/e_logf.c b/math/logf.c
index eb06a63..eb06a63 100644
--- a/math/e_logf.c
+++ b/math/logf.c
diff --git a/math/e_logf_data.c b/math/logf_data.c
index 521d413..521d413 100644
--- a/math/e_logf_data.c
+++ b/math/logf_data.c
diff --git a/math/e_powf.c b/math/powf.c
index 8ffbe84..8ffbe84 100644
--- a/math/e_powf.c
+++ b/math/powf.c
diff --git a/math/e_powf_log2_data.c b/math/powf_log2_data.c
index c21406b..c21406b 100644
--- a/math/e_powf_log2_data.c
+++ b/math/powf_log2_data.c
diff --git a/math/rem_pio2.c b/math/rem_pio2.c
new file mode 100644
index 0000000..16edd74
--- /dev/null
+++ b/math/rem_pio2.c
@@ -0,0 +1 @@
+#include "single/e_rem_pio2.c"
diff --git a/math/rredf.c b/math/rredf.c
index b463480..c96fee4 100644
--- a/math/rredf.c
+++ b/math/rredf.c
@@ -1,251 +1 @@
-/*
- * rredf.c - trigonometric range reduction function
- *
- * Copyright (c) 2009-2015, Arm Limited.
- * SPDX-License-Identifier: Apache-2.0
- *
- * 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.
- */
-
-/*
- * This code is intended to be used as the second half of a range
- * reducer whose first half is an inline function defined in
- * rredf.h. Each trig function performs range reduction by invoking
- * that, which handles the quickest and most common cases inline
- * before handing off to this function for everything else. Thus a
- * reasonable compromise is struck between speed and space. (I
- * hope.) In particular, this approach avoids a function call
- * overhead in the common case.
- */
-
-#include "math_private.h"
-
-#ifdef __cplusplus
-extern "C" {
-#endif /* __cplusplus */
-
-/*
- * Input values to this function:
- *  - x is the original user input value, unchanged by the
- *    first-tier reducer in the case where it hands over to us.
- *  - q is still the place where the caller expects us to leave the
- *    quadrant code.
- *  - k is the IEEE bit pattern of x (which it would seem a shame to
- *    recompute given that the first-tier reducer already went to
- *    the effort of extracting it from the VFP). FIXME: in softfp,
- *    on the other hand, it's unconscionably wasteful to replicate
- *    this value into a second register and we should change the
- *    prototype!
- */
-float ARM__mathlib_rredf2(float x, int *q, unsigned k)
-{
-    /*
-     * First, weed out infinities and NaNs, and deal with them by
-     * returning a negative q.
-     */
-    if ((k << 1) >= 0xFF000000) {
-        *q = -1;
-        return x;
-    }
-    /*
-     * We do general range reduction by multiplying by 2/pi, and
-     * retaining the bottom two bits of the integer part and an
-     * initial chunk of the fraction below that. The integer bits
-     * are directly output as *q; the fraction is then multiplied
-     * back up by pi/2 before returning it.
-     *
-     * To get this right, we don't have to multiply by the _whole_
-     * of 2/pi right from the most significant bit downwards:
-     * instead we can discard any bit of 2/pi with a place value
-     * high enough that multiplying it by the LSB of x will yield a
-     * place value higher than 2. Thus we can bound the required
-     * work by a reasonably small constant regardless of the size of
-     * x (unlike, for instance, the IEEE remainder operation).
-     *
-     * At the other end, however, we must take more care: it isn't
-     * adequate just to acquire two integer bits and 24 fraction
-     * bits of (2/pi)x, because if a lot of those fraction bits are
-     * zero then we will suffer significance loss. So we must keep
-     * computing fraction bits as far down as 23 bits below the
-     * _highest set fraction bit_.
-     *
-     * The immediate question, therefore, is what the bound on this
-     * end of the job will be. In other words: what is the smallest
-     * difference between an integer multiple of pi/2 and a
-     * representable IEEE single precision number larger than the
-     * maximum size handled by rredf.h?
-     *
-     * The most difficult cases for each exponent can readily be
-     * found by Tim Peters's modular minimisation algorithm, and are
-     * tabulated in mathlib/tests/directed/rredf.tst. The single
-     * worst case is the IEEE single-precision number 0x6F79BE45,
-     * whose numerical value is in the region of 7.7*10^28; when
-     * reduced mod pi/2, it attains the value 0x30DDEEA9, or about
-     * 0.00000000161. The highest set bit of this value is the one
-     * with place value 2^-30; so its lowest is 2^-53. Hence, to be
-     * sure of having enough fraction bits to output at full single
-     * precision, we must be prepared to collect up to 53 bits of
-     * fraction in addition to our two bits of integer part.
-     *
-     * To begin with, this means we must store the value of 2/pi to
-     * a precision of 128+53 = 181 bits. That's six 32-bit words.
-     * (Hardly a chore, unlike the equivalent problem in double
-     * precision!)
-     */
-    {
-        static const unsigned twooverpi[] = {
-            /* We start with a zero word, because that takes up less
-             * space than the array bounds checking and special-case
-             * handling that would have to occur in its absence. */
-            0,
-            /* 2/pi in hex is 0.a2f9836e... */
-            0xa2f9836e, 0x4e441529, 0xfc2757d1,
-            0xf534ddc0, 0xdb629599, 0x3c439041,
-            /* Again, to avoid array bounds overrun, we store a spare
-             * word at the end. And it would be a shame to fill it
-             * with zeroes when we could use more bits of 2/pi... */
-            0xfe5163ab
-        };
-
-        /*
-         * Multiprecision multiplication of this nature is more
-         * readily done in integers than in VFP, since we can use
-         * UMULL (on CPUs that support it) to multiply 32 by 32 bits
-         * at a time whereas the VFP would only be able to do 12x12
-         * without losing accuracy.
-         *
-         * So extract the mantissa of the input number as a 32-bit
-         * integer.
-         */
-        unsigned mantissa = 0x80000000 | (k << 8);
-
-        /*
-         * Now work out which part of our stored value of 2/pi we're
-         * supposed to be multiplying by.
-         *
-         * Let the IEEE exponent field of x be e. With its bias
-         * removed, (e-127) is the index of the set bit at the top
-         * of 'mantissa' (i.e. that set bit has real place value
-         * 2^(e-127)). So the lowest set bit in 'mantissa', 23 bits
-         * further down, must have place value 2^(e-150).
-         *
-         * We begin taking an interest in the value of 2/pi at the
-         * bit which multiplies by _that_ to give something with
-         * place value at most 2. In other words, the highest bit of
-         * 2/pi we're interested in is the one with place value
-         * 2/(2^(e-150)) = 2^(151-e).
-         *
-         * The bit at the top of the first (zero) word of the above
-         * array has place value 2^31. Hence, the bit we want to put
-         * at the top of the first word we extract from that array
-         * is the one at bit index n, where 31-n = 151-e and hence
-         * n=e-120.
-         */
-        int topbitindex = ((k >> 23) & 0xFF) - 120;
-        int wordindex = topbitindex >> 5;
-        int shiftup = topbitindex & 31;
-        int shiftdown = 32 - shiftup;
-        unsigned word1, word2, word3;
-        if (shiftup) {
-            word1 = (twooverpi[wordindex] << shiftup) | (twooverpi[wordindex+1] >> shiftdown);
-            word2 = (twooverpi[wordindex+1] << shiftup) | (twooverpi[wordindex+2] >> shiftdown);
-            word3 = (twooverpi[wordindex+2] << shiftup) | (twooverpi[wordindex+3] >> shiftdown);
-        } else {
-            word1 = twooverpi[wordindex];
-            word2 = twooverpi[wordindex+1];
-            word3 = twooverpi[wordindex+2];
-        }
-
-        /*
-         * Do the multiplications, and add them together.
-         */
-        unsigned long long mult1 = (unsigned long long)word1 * mantissa;
-        unsigned long long mult2 = (unsigned long long)word2 * mantissa;
-        unsigned long long mult3 = (unsigned long long)word3 * mantissa;
-
-        unsigned /* bottom3 = (unsigned)mult3, */ top3 = (unsigned)(mult3 >> 32);
-        unsigned bottom2 = (unsigned)mult2, top2 = (unsigned)(mult2 >> 32);
-        unsigned bottom1 = (unsigned)mult1, top1 = (unsigned)(mult1 >> 32);
-
-        unsigned out3, out2, out1, carry;
-
-        out3 = top3 + bottom2; carry = (out3 < top3);
-        out2 = top2 + bottom1 + carry; carry = carry ? (out2 <= top2) : (out2 < top2);
-        out1 = top1 + carry;
-
-        /*
-         * The two words we multiplied to get mult1 had their top
-         * bits at (respectively) place values 2^(151-e) and
-         * 2^(e-127). The value of those two bits multiplied
-         * together will have ended up in bit 62 (the
-         * topmost-but-one bit) of mult1, i.e. bit 30 of out1.
-         * Hence, that bit has place value 2^(151-e+e-127) = 2^24.
-         * So the integer value that we want to output as q,
-         * consisting of the bits with place values 2^1 and 2^0,
-         * must be 23 and 24 bits below that, i.e. in bits 7 and 6
-         * of out1.
-         *
-         * Or, at least, it will be once we add 1/2, to round to the
-         * _nearest_ multiple of pi/2 rather than the next one down.
-         */
-        *q = (out1 + (1<<5)) >> 6;
-
-        /*
-         * Now we construct the output fraction, which is most
-         * simply done in the VFP. We just extract three consecutive
-         * bit strings from our chunk of binary data, convert them
-         * to integers, equip each with an appropriate FP exponent,
-         * add them together, and (don't forget) multiply back up by
-         * pi/2. That way we don't have to work out ourselves where
-         * the highest fraction bit ended up.
-         *
-         * Since our displacement from the nearest multiple of pi/2
-         * can be positive or negative, the topmost of these three
-         * values must be arranged with its 2^-1 bit at the very top
-         * of the word, and then treated as a _signed_ integer.
-         */
-        {
-            int i1 = (out1 << 26) | ((out2 >> 19) << 13);
-            unsigned i2 = out2 << 13;
-            unsigned i3 = out3;
-            float f1 = i1, f2 = i2 * (1.0f/524288.0f), f3 = i3 * (1.0f/524288.0f/524288.0f);
-
-            /*
-             * Now f1+f2+f3 is a representation, potentially to
-             * twice double precision, of 2^32 times ((2/pi)*x minus
-             * some integer). So our remaining job is to multiply
-             * back down by (pi/2)*2^-32, and convert back to one
-             * single-precision output number.
-             */
-
-            /* Normalise to a prec-and-a-half representation... */
-            float ftop = CLEARBOTTOMHALF(f1+f2+f3), fbot = f3-((ftop-f1)-f2);
-
-            /* ... and multiply by a prec-and-a-half value of (pi/2)*2^-32. */
-            float ret = (ftop * 0x1.92p-32F) + (ftop * 0x1.fb5444p-44F + fbot * 0x1.921fb6p-32F);
-
-            /* Just before we return, take the input sign into account. */
-            if (k & 0x80000000) {
-                *q = 0x10000000 - *q;
-                ret = -ret;
-            }
-            return ret;
-        }
-    }
-}
-
-#ifdef __cplusplus
-} /* end of extern "C" */
-#endif /* __cplusplus */
-
-/* end of rredf.c */
+#include "single/rredf.c"
diff --git a/math/sinf.c b/math/sinf.c
new file mode 100644
index 0000000..859e9ba
--- /dev/null
+++ b/math/sinf.c
@@ -0,0 +1 @@
+#include "single/s_sinf.c"
diff --git a/math/dunder.c b/math/single/dunder.c
index 07af6d2..07af6d2 100644
--- a/math/dunder.c
+++ b/math/single/dunder.c
diff --git a/math/e_rem_pio2.c b/math/single/e_rem_pio2.c
index ec18f35..ec18f35 100644
--- a/math/e_rem_pio2.c
+++ b/math/single/e_rem_pio2.c
diff --git a/math/single/funder.c b/math/single/funder.c
new file mode 100644
index 0000000..5e5f6b2
--- /dev/null
+++ b/math/single/funder.c
@@ -0,0 +1,63 @@
+/*
+ * funder.c - manually provoke SP exceptions for mathlib
+ *
+ * Copyright (c) 2009-2015, Arm Limited.
+ * SPDX-License-Identifier: Apache-2.0
+ *
+ * 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 "math_private.h"
+#include <fenv.h>
+
+__inline float __mathlib_flt_infnan2(float x, float y)
+{
+  return x+y;
+}
+
+__inline float __mathlib_flt_infnan(float x)
+{
+  return x+x;
+}
+
+float __mathlib_flt_underflow(void)
+{
+#ifdef CLANG_EXCEPTIONS
+  feraiseexcept(FE_UNDERFLOW);
+#endif
+  return 0x1p-95F * 0x1p-95F;
+}
+
+float __mathlib_flt_overflow(void)
+{
+#ifdef CLANG_EXCEPTIONS
+  feraiseexcept(FE_OVERFLOW);
+#endif
+  return 0x1p+97F * 0x1p+97F;
+}
+
+float __mathlib_flt_invalid(void)
+{
+#ifdef CLANG_EXCEPTIONS
+  feraiseexcept(FE_INVALID);
+#endif
+  return 0.0f / 0.0f;
+}
+
+float __mathlib_flt_divzero(void)
+{
+#ifdef CLANG_EXCEPTIONS
+  feraiseexcept(FE_DIVBYZERO);
+#endif
+  return 1.0f / 0.0f;
+}
diff --git a/math/ieee_status.c b/math/single/ieee_status.c
index e9f4d16..e9f4d16 100644
--- a/math/ieee_status.c
+++ b/math/single/ieee_status.c
diff --git a/math/math_private.h b/math/single/math_private.h
index 0b57072..0b57072 100644
--- a/math/math_private.h
+++ b/math/single/math_private.h
diff --git a/math/poly.c b/math/single/poly.c
index 6f25bf5..6f25bf5 100644
--- a/math/poly.c
+++ b/math/single/poly.c
diff --git a/math/single/rredf.c b/math/single/rredf.c
new file mode 100644
index 0000000..b463480
--- /dev/null
+++ b/math/single/rredf.c
@@ -0,0 +1,251 @@
+/*
+ * rredf.c - trigonometric range reduction function
+ *
+ * Copyright (c) 2009-2015, Arm Limited.
+ * SPDX-License-Identifier: Apache-2.0
+ *
+ * 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.
+ */
+
+/*
+ * This code is intended to be used as the second half of a range
+ * reducer whose first half is an inline function defined in
+ * rredf.h. Each trig function performs range reduction by invoking
+ * that, which handles the quickest and most common cases inline
+ * before handing off to this function for everything else. Thus a
+ * reasonable compromise is struck between speed and space. (I
+ * hope.) In particular, this approach avoids a function call
+ * overhead in the common case.
+ */
+
+#include "math_private.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif /* __cplusplus */
+
+/*
+ * Input values to this function:
+ *  - x is the original user input value, unchanged by the
+ *    first-tier reducer in the case where it hands over to us.
+ *  - q is still the place where the caller expects us to leave the
+ *    quadrant code.
+ *  - k is the IEEE bit pattern of x (which it would seem a shame to
+ *    recompute given that the first-tier reducer already went to
+ *    the effort of extracting it from the VFP). FIXME: in softfp,
+ *    on the other hand, it's unconscionably wasteful to replicate
+ *    this value into a second register and we should change the
+ *    prototype!
+ */
+float ARM__mathlib_rredf2(float x, int *q, unsigned k)
+{
+    /*
+     * First, weed out infinities and NaNs, and deal with them by
+     * returning a negative q.
+     */
+    if ((k << 1) >= 0xFF000000) {
+        *q = -1;
+        return x;
+    }
+    /*
+     * We do general range reduction by multiplying by 2/pi, and
+     * retaining the bottom two bits of the integer part and an
+     * initial chunk of the fraction below that. The integer bits
+     * are directly output as *q; the fraction is then multiplied
+     * back up by pi/2 before returning it.
+     *
+     * To get this right, we don't have to multiply by the _whole_
+     * of 2/pi right from the most significant bit downwards:
+     * instead we can discard any bit of 2/pi with a place value
+     * high enough that multiplying it by the LSB of x will yield a
+     * place value higher than 2. Thus we can bound the required
+     * work by a reasonably small constant regardless of the size of
+     * x (unlike, for instance, the IEEE remainder operation).
+     *
+     * At the other end, however, we must take more care: it isn't
+     * adequate just to acquire two integer bits and 24 fraction
+     * bits of (2/pi)x, because if a lot of those fraction bits are
+     * zero then we will suffer significance loss. So we must keep
+     * computing fraction bits as far down as 23 bits below the
+     * _highest set fraction bit_.
+     *
+     * The immediate question, therefore, is what the bound on this
+     * end of the job will be. In other words: what is the smallest
+     * difference between an integer multiple of pi/2 and a
+     * representable IEEE single precision number larger than the
+     * maximum size handled by rredf.h?
+     *
+     * The most difficult cases for each exponent can readily be
+     * found by Tim Peters's modular minimisation algorithm, and are
+     * tabulated in mathlib/tests/directed/rredf.tst. The single
+     * worst case is the IEEE single-precision number 0x6F79BE45,
+     * whose numerical value is in the region of 7.7*10^28; when
+     * reduced mod pi/2, it attains the value 0x30DDEEA9, or about
+     * 0.00000000161. The highest set bit of this value is the one
+     * with place value 2^-30; so its lowest is 2^-53. Hence, to be
+     * sure of having enough fraction bits to output at full single
+     * precision, we must be prepared to collect up to 53 bits of
+     * fraction in addition to our two bits of integer part.
+     *
+     * To begin with, this means we must store the value of 2/pi to
+     * a precision of 128+53 = 181 bits. That's six 32-bit words.
+     * (Hardly a chore, unlike the equivalent problem in double
+     * precision!)
+     */
+    {
+        static const unsigned twooverpi[] = {
+            /* We start with a zero word, because that takes up less
+             * space than the array bounds checking and special-case
+             * handling that would have to occur in its absence. */
+            0,
+            /* 2/pi in hex is 0.a2f9836e... */
+            0xa2f9836e, 0x4e441529, 0xfc2757d1,
+            0xf534ddc0, 0xdb629599, 0x3c439041,
+            /* Again, to avoid array bounds overrun, we store a spare
+             * word at the end. And it would be a shame to fill it
+             * with zeroes when we could use more bits of 2/pi... */
+            0xfe5163ab
+        };
+
+        /*
+         * Multiprecision multiplication of this nature is more
+         * readily done in integers than in VFP, since we can use
+         * UMULL (on CPUs that support it) to multiply 32 by 32 bits
+         * at a time whereas the VFP would only be able to do 12x12
+         * without losing accuracy.
+         *
+         * So extract the mantissa of the input number as a 32-bit
+         * integer.
+         */
+        unsigned mantissa = 0x80000000 | (k << 8);
+
+        /*
+         * Now work out which part of our stored value of 2/pi we're
+         * supposed to be multiplying by.
+         *
+         * Let the IEEE exponent field of x be e. With its bias
+         * removed, (e-127) is the index of the set bit at the top
+         * of 'mantissa' (i.e. that set bit has real place value
+         * 2^(e-127)). So the lowest set bit in 'mantissa', 23 bits
+         * further down, must have place value 2^(e-150).
+         *
+         * We begin taking an interest in the value of 2/pi at the
+         * bit which multiplies by _that_ to give something with
+         * place value at most 2. In other words, the highest bit of
+         * 2/pi we're interested in is the one with place value
+         * 2/(2^(e-150)) = 2^(151-e).
+         *
+         * The bit at the top of the first (zero) word of the above
+         * array has place value 2^31. Hence, the bit we want to put
+         * at the top of the first word we extract from that array
+         * is the one at bit index n, where 31-n = 151-e and hence
+         * n=e-120.
+         */
+        int topbitindex = ((k >> 23) & 0xFF) - 120;
+        int wordindex = topbitindex >> 5;
+        int shiftup = topbitindex & 31;
+        int shiftdown = 32 - shiftup;
+        unsigned word1, word2, word3;
+        if (shiftup) {
+            word1 = (twooverpi[wordindex] << shiftup) | (twooverpi[wordindex+1] >> shiftdown);
+            word2 = (twooverpi[wordindex+1] << shiftup) | (twooverpi[wordindex+2] >> shiftdown);
+            word3 = (twooverpi[wordindex+2] << shiftup) | (twooverpi[wordindex+3] >> shiftdown);
+        } else {
+            word1 = twooverpi[wordindex];
+            word2 = twooverpi[wordindex+1];
+            word3 = twooverpi[wordindex+2];
+        }
+
+        /*
+         * Do the multiplications, and add them together.
+         */
+        unsigned long long mult1 = (unsigned long long)word1 * mantissa;
+        unsigned long long mult2 = (unsigned long long)word2 * mantissa;
+        unsigned long long mult3 = (unsigned long long)word3 * mantissa;
+
+        unsigned /* bottom3 = (unsigned)mult3, */ top3 = (unsigned)(mult3 >> 32);
+        unsigned bottom2 = (unsigned)mult2, top2 = (unsigned)(mult2 >> 32);
+        unsigned bottom1 = (unsigned)mult1, top1 = (unsigned)(mult1 >> 32);
+
+        unsigned out3, out2, out1, carry;
+
+        out3 = top3 + bottom2; carry = (out3 < top3);
+        out2 = top2 + bottom1 + carry; carry = carry ? (out2 <= top2) : (out2 < top2);
+        out1 = top1 + carry;
+
+        /*
+         * The two words we multiplied to get mult1 had their top
+         * bits at (respectively) place values 2^(151-e) and
+         * 2^(e-127). The value of those two bits multiplied
+         * together will have ended up in bit 62 (the
+         * topmost-but-one bit) of mult1, i.e. bit 30 of out1.
+         * Hence, that bit has place value 2^(151-e+e-127) = 2^24.
+         * So the integer value that we want to output as q,
+         * consisting of the bits with place values 2^1 and 2^0,
+         * must be 23 and 24 bits below that, i.e. in bits 7 and 6
+         * of out1.
+         *
+         * Or, at least, it will be once we add 1/2, to round to the
+         * _nearest_ multiple of pi/2 rather than the next one down.
+         */
+        *q = (out1 + (1<<5)) >> 6;
+
+        /*
+         * Now we construct the output fraction, which is most
+         * simply done in the VFP. We just extract three consecutive
+         * bit strings from our chunk of binary data, convert them
+         * to integers, equip each with an appropriate FP exponent,
+         * add them together, and (don't forget) multiply back up by
+         * pi/2. That way we don't have to work out ourselves where
+         * the highest fraction bit ended up.
+         *
+         * Since our displacement from the nearest multiple of pi/2
+         * can be positive or negative, the topmost of these three
+         * values must be arranged with its 2^-1 bit at the very top
+         * of the word, and then treated as a _signed_ integer.
+         */
+        {
+            int i1 = (out1 << 26) | ((out2 >> 19) << 13);
+            unsigned i2 = out2 << 13;
+            unsigned i3 = out3;
+            float f1 = i1, f2 = i2 * (1.0f/524288.0f), f3 = i3 * (1.0f/524288.0f/524288.0f);
+
+            /*
+             * Now f1+f2+f3 is a representation, potentially to
+             * twice double precision, of 2^32 times ((2/pi)*x minus
+             * some integer). So our remaining job is to multiply
+             * back down by (pi/2)*2^-32, and convert back to one
+             * single-precision output number.
+             */
+
+            /* Normalise to a prec-and-a-half representation... */
+            float ftop = CLEARBOTTOMHALF(f1+f2+f3), fbot = f3-((ftop-f1)-f2);
+
+            /* ... and multiply by a prec-and-a-half value of (pi/2)*2^-32. */
+            float ret = (ftop * 0x1.92p-32F) + (ftop * 0x1.fb5444p-44F + fbot * 0x1.921fb6p-32F);
+
+            /* Just before we return, take the input sign into account. */
+            if (k & 0x80000000) {
+                *q = 0x10000000 - *q;
+                ret = -ret;
+            }
+            return ret;
+        }
+    }
+}
+
+#ifdef __cplusplus
+} /* end of extern "C" */
+#endif /* __cplusplus */
+
+/* end of rredf.c */
diff --git a/math/rredf.h b/math/single/rredf.h
index d888487..d888487 100644
--- a/math/rredf.h
+++ b/math/single/rredf.h
diff --git a/math/s_cosf.c b/math/single/s_cosf.c
index f8adad6..f8adad6 100644
--- a/math/s_cosf.c
+++ b/math/single/s_cosf.c
diff --git a/math/s_sincosf.c b/math/single/s_sincosf.c
index 73712eb..73712eb 100644
--- a/math/s_sincosf.c
+++ b/math/single/s_sincosf.c
diff --git a/math/s_sinf.c b/math/single/s_sinf.c
index 1640329..1640329 100644
--- a/math/s_sinf.c
+++ b/math/single/s_sinf.c
diff --git a/math/s_tanf.c b/math/single/s_tanf.c
index f628165..f628165 100644
--- a/math/s_tanf.c
+++ b/math/single/s_tanf.c
diff --git a/math/tanf.c b/math/tanf.c
new file mode 100644
index 0000000..36ecb4f
--- /dev/null
+++ b/math/tanf.c
@@ -0,0 +1 @@
+#include "single/s_tanf.c"
diff --git a/test/mathtest.c b/test/mathtest.c
index 7d77260..cb4e5d1 100644
--- a/test/mathtest.c
+++ b/test/mathtest.c
@@ -56,7 +56,7 @@ _Pragma(STR(import IMPORT_SYMBOL))
 #endif
 
 EXTERN_C int ARM__ieee754_rem_pio2(double, double *);
-#include "../math/rredf.h"
+#include "../math/single/rredf.h"
 
 int sp_rem_pio2(float x, float *y) {
   int q;