|
 |
8ae002 |
From 1ddf0e88cfea375eee1c7da2127a8520c02862b3 Mon Sep 17 00:00:00 2001
|
|
|
8ae002 |
From: Anton Blanchard <anton@samba.org>
|
|
|
8ae002 |
Date: Tue, 28 Jun 2016 21:59:40 +1000
|
|
|
8ae002 |
Subject: [PATCH] powerpc: Add a POWER8-optimized version of sinf()
|
|
|
8ae002 |
|
|
|
8ae002 |
This uses the implementation of sinf() in sysdeps/x86_64/fpu/s_sinf.S
|
|
|
8ae002 |
as inspiration.
|
|
|
8ae002 |
|
|
|
8ae002 |
(cherry picked from commit aa95fc13f5b02044eadc3af3d9e1c025f2e1edda)
|
|
|
8ae002 |
---
|
|
|
8ae002 |
ChangeLog | 10 +
|
|
|
8ae002 |
sysdeps/powerpc/powerpc64/fpu/multiarch/Makefile | 3 +-
|
|
|
8ae002 |
.../powerpc64/fpu/multiarch/s_sinf-power8.S | 26 ++
|
|
|
8ae002 |
.../powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c | 26 ++
|
|
|
8ae002 |
sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf.c | 31 ++
|
|
|
8ae002 |
sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S | 519 +++++++++++++++++++++
|
|
|
8ae002 |
6 files changed, 614 insertions(+), 1 deletion(-)
|
|
|
8ae002 |
create mode 100644 sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-power8.S
|
|
|
8ae002 |
create mode 100644 sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c
|
|
|
8ae002 |
create mode 100644 sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf.c
|
|
|
8ae002 |
create mode 100644 sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S
|
|
|
8ae002 |
|
|
|
8ae002 |
diff --git a/ChangeLog b/ChangeLog
|
|
|
8ae002 |
index 6cb2578..6d6aab3 100644
|
|
|
8ae002 |
diff --git a/sysdeps/powerpc/powerpc64/fpu/multiarch/Makefile b/sysdeps/powerpc/powerpc64/fpu/multiarch/Makefile
|
|
|
8ae002 |
index 331763e..add1fb8 100644
|
|
|
8ae002 |
--- a/sysdeps/powerpc/powerpc64/fpu/multiarch/Makefile
|
|
|
8ae002 |
+++ b/sysdeps/powerpc/powerpc64/fpu/multiarch/Makefile
|
|
|
8ae002 |
@@ -25,7 +25,8 @@ libm-sysdep_routines += s_isnan-power7 s_isnan-power6x s_isnan-power6 \
|
|
|
8ae002 |
e_hypot-power7 e_hypotf-ppc64 e_hypotf-power7 \
|
|
|
8ae002 |
s_isnan-power8 s_isinf-power8 s_finite-power8 \
|
|
|
8ae002 |
s_llrint-power8 s_llround-power8 \
|
|
|
8ae002 |
- e_expf-power8 e_expf-ppc64
|
|
|
8ae002 |
+ e_expf-power8 e_expf-ppc64 \
|
|
|
8ae002 |
+ s_sinf-ppc64 s_sinf-power8
|
|
|
8ae002 |
|
|
|
8ae002 |
CFLAGS-s_logbf-power7.c = -mcpu=power7
|
|
|
8ae002 |
CFLAGS-s_logbl-power7.c = -mcpu=power7
|
|
|
8ae002 |
diff --git a/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-power8.S b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-power8.S
|
|
|
8ae002 |
new file mode 100644
|
|
|
8ae002 |
index 0000000..579019c
|
|
|
8ae002 |
--- /dev/null
|
|
|
8ae002 |
+++ b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-power8.S
|
|
|
8ae002 |
@@ -0,0 +1,26 @@
|
|
|
8ae002 |
+/* sinf(). PowerPC64/POWER8 version.
|
|
|
8ae002 |
+ Copyright (C) 2016 Free Software Foundation, Inc.
|
|
|
8ae002 |
+ This file is part of the GNU C Library.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is free software; you can redistribute it and/or
|
|
|
8ae002 |
+ modify it under the terms of the GNU Lesser General Public
|
|
|
8ae002 |
+ License as published by the Free Software Foundation; either
|
|
|
8ae002 |
+ version 2.1 of the License, or (at your option) any later version.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is distributed in the hope that it will be useful,
|
|
|
8ae002 |
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
8ae002 |
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
8ae002 |
+ Lesser General Public License for more details.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ You should have received a copy of the GNU Lesser General Public
|
|
|
8ae002 |
+ License along with the GNU C Library; if not, see
|
|
|
8ae002 |
+ <http://www.gnu.org/licenses/>. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <sysdep.h>
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#undef weak_alias
|
|
|
8ae002 |
+#define weak_alias(a, b)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#define __sinf __sinf_power8
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S>
|
|
|
8ae002 |
diff --git a/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c
|
|
|
8ae002 |
new file mode 100644
|
|
|
8ae002 |
index 0000000..eaf83fa
|
|
|
8ae002 |
--- /dev/null
|
|
|
8ae002 |
+++ b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf-ppc64.c
|
|
|
8ae002 |
@@ -0,0 +1,26 @@
|
|
|
8ae002 |
+/* sinf(). PowerPC64 default version.
|
|
|
8ae002 |
+ Copyright (C) 2016 Free Software Foundation, Inc.
|
|
|
8ae002 |
+ This file is part of the GNU C Library.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is free software; you can redistribute it and/or
|
|
|
8ae002 |
+ modify it under the terms of the GNU Lesser General Public
|
|
|
8ae002 |
+ License as published by the Free Software Foundation; either
|
|
|
8ae002 |
+ version 2.1 of the License, or (at your option) any later version.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is distributed in the hope that it will be useful,
|
|
|
8ae002 |
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
8ae002 |
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
8ae002 |
+ Lesser General Public License for more details.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ You should have received a copy of the GNU Lesser General Public
|
|
|
8ae002 |
+ License along with the GNU C Library; if not, see
|
|
|
8ae002 |
+ <http://www.gnu.org/licenses/>. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <sysdep.h>
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#undef weak_alias
|
|
|
8ae002 |
+#define weak_alias(a, b)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#define __sinf __sinf_ppc64
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <sysdeps/ieee754/flt-32/s_sinf.c>
|
|
|
8ae002 |
diff --git a/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf.c b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf.c
|
|
|
8ae002 |
new file mode 100644
|
|
|
8ae002 |
index 0000000..4269d58
|
|
|
8ae002 |
--- /dev/null
|
|
|
8ae002 |
+++ b/sysdeps/powerpc/powerpc64/fpu/multiarch/s_sinf.c
|
|
|
8ae002 |
@@ -0,0 +1,31 @@
|
|
|
8ae002 |
+/* Multiple versions of sinf.
|
|
|
8ae002 |
+ Copyright (C) 2016 Free Software Foundation, Inc.
|
|
|
8ae002 |
+ This file is part of the GNU C Library.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is free software; you can redistribute it and/or
|
|
|
8ae002 |
+ modify it under the terms of the GNU Lesser General Public
|
|
|
8ae002 |
+ License as published by the Free Software Foundation; either
|
|
|
8ae002 |
+ version 2.1 of the License, or (at your option) any later version.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is distributed in the hope that it will be useful,
|
|
|
8ae002 |
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
8ae002 |
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
8ae002 |
+ Lesser General Public License for more details.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ You should have received a copy of the GNU Lesser General Public
|
|
|
8ae002 |
+ License along with the GNU C Library; if not, see
|
|
|
8ae002 |
+ <http://www.gnu.org/licenses/>. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <math.h>
|
|
|
8ae002 |
+#include <shlib-compat.h>
|
|
|
8ae002 |
+#include "init-arch.h"
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+extern __typeof (__sinf) __sinf_ppc64 attribute_hidden;
|
|
|
8ae002 |
+extern __typeof (__sinf) __sinf_power8 attribute_hidden;
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+libc_ifunc (__sinf,
|
|
|
8ae002 |
+ (hwcap2 & PPC_FEATURE2_ARCH_2_07)
|
|
|
8ae002 |
+ ? __sinf_power8
|
|
|
8ae002 |
+ : __sinf_ppc64);
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+weak_alias (__sinf, sinf)
|
|
|
8ae002 |
diff --git a/sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S b/sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S
|
|
|
8ae002 |
new file mode 100644
|
|
|
8ae002 |
index 0000000..3b8f5af
|
|
|
8ae002 |
--- /dev/null
|
|
|
8ae002 |
+++ b/sysdeps/powerpc/powerpc64/power8/fpu/s_sinf.S
|
|
|
8ae002 |
@@ -0,0 +1,519 @@
|
|
|
8ae002 |
+/* Optimized sinf(). PowerPC64/POWER8 version.
|
|
|
8ae002 |
+ Copyright (C) 2016 Free Software Foundation, Inc.
|
|
|
8ae002 |
+ This file is part of the GNU C Library.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is free software; you can redistribute it and/or
|
|
|
8ae002 |
+ modify it under the terms of the GNU Lesser General Public
|
|
|
8ae002 |
+ License as published by the Free Software Foundation; either
|
|
|
8ae002 |
+ version 2.1 of the License, or (at your option) any later version.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ The GNU C Library is distributed in the hope that it will be useful,
|
|
|
8ae002 |
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
8ae002 |
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
8ae002 |
+ Lesser General Public License for more details.
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ You should have received a copy of the GNU Lesser General Public
|
|
|
8ae002 |
+ License along with the GNU C Library; if not, see
|
|
|
8ae002 |
+ <http://www.gnu.org/licenses/>. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#include <sysdep.h>
|
|
|
8ae002 |
+#define _ERRNO_H 1
|
|
|
8ae002 |
+#include <bits/errno.h>
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#define FRAMESIZE (FRAME_MIN_SIZE+16)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#define FLOAT_EXPONENT_SHIFT 23
|
|
|
8ae002 |
+#define FLOAT_EXPONENT_BIAS 127
|
|
|
8ae002 |
+#define INTEGER_BITS 3
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+#define PI_4 0x3f490fdb /* PI/4 */
|
|
|
8ae002 |
+#define NINEPI_4 0x40e231d6 /* 9 * PI/4 */
|
|
|
8ae002 |
+#define TWO_PN5 0x3d000000 /* 2^-5 */
|
|
|
8ae002 |
+#define TWO_PN27 0x32000000 /* 2^-27 */
|
|
|
8ae002 |
+#define INFINITY 0x7f800000
|
|
|
8ae002 |
+#define TWO_P23 0x4b000000 /* 2^27 */
|
|
|
8ae002 |
+#define FX_FRACTION_1_28 0x9249250 /* 0x100000000 / 28 + 1 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Implements the function
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ float [fp1] sinf (float [fp1] x) */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .machine power8
|
|
|
8ae002 |
+EALIGN(__sinf, 4, 0)
|
|
|
8ae002 |
+ addis r9,r2,L(anchor)@toc@ha
|
|
|
8ae002 |
+ addi r9,r9,L(anchor)@toc@l
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lis r4,PI_4@h
|
|
|
8ae002 |
+ ori r4,r4,PI_4@l
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ xscvdpspn v0,v1
|
|
|
8ae002 |
+ mfvsrd r8,v0
|
|
|
8ae002 |
+ rldicl r3,r8,32,33 /* Remove sign bit. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ bge L(greater_or_equal_pio4)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lis r4,TWO_PN5@h
|
|
|
8ae002 |
+ ori r4,r4,TWO_PN5@l
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ blt L(less_2pn5)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Chebyshev polynomial of the form:
|
|
|
8ae002 |
+ * x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp9,(L(S0)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp10,(L(S1)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp11,(L(S2)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp12,(L(S3)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp13,(L(S4)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp2,fp1,fp1 /* x^2 */
|
|
|
8ae002 |
+ fmul fp3,fp2,fp1 /* x^3 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp13,fp12 /* S3+x^2*S4 */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp11 /* S2+x^2*(S3+x^2*S4) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp10 /* S1+x^2*(S2+x^2*(S3+x^2*S4)) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp9 /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))) */
|
|
|
8ae002 |
+ fmadd fp1,fp3,fp4,fp1 /* x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))) */
|
|
|
8ae002 |
+ frsp fp1,fp1 /* Round to single precision. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(greater_or_equal_pio4):
|
|
|
8ae002 |
+ lis r4,NINEPI_4@h
|
|
|
8ae002 |
+ ori r4,r4,NINEPI_4@l
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ bge L(greater_or_equal_9pio4)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Calculate quotient of |x|/(PI/4). */
|
|
|
8ae002 |
+ lfd fp2,(L(invpio4)-L(anchor))(r9)
|
|
|
8ae002 |
+ fabs fp1,fp1 /* |x| */
|
|
|
8ae002 |
+ fmul fp2,fp1,fp2 /* |x|/(PI/4) */
|
|
|
8ae002 |
+ fctiduz fp2,fp2
|
|
|
8ae002 |
+ mfvsrd r3,v2 /* n = |x| mod PI/4 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Now use that quotient to find |x| mod (PI/2). */
|
|
|
8ae002 |
+ addi r7,r3,1
|
|
|
8ae002 |
+ rldicr r5,r7,2,60 /* ((n+1) >> 1) << 3 */
|
|
|
8ae002 |
+ addi r6,r9,(L(pio2_table)-L(anchor))
|
|
|
8ae002 |
+ lfdx fp4,r5,r6
|
|
|
8ae002 |
+ fsub fp1,fp1,fp4
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(reduced):
|
|
|
8ae002 |
+ /* Now we are in the range -PI/4 to PI/4. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Work out if we are in a positive or negative primary interval. */
|
|
|
8ae002 |
+ rldicl r4,r7,62,63 /* ((n+1) >> 2) & 1 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* We are operating on |x|, so we need to add back the original
|
|
|
8ae002 |
+ sign. */
|
|
|
8ae002 |
+ rldicl r8,r8,33,63 /* (x >> 31) & 1, ie the sign bit. */
|
|
|
8ae002 |
+ xor r4,r4,r8 /* 0 if result should be positive,
|
|
|
8ae002 |
+ 1 if negative. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Load a 1.0 or -1.0. */
|
|
|
8ae002 |
+ addi r5,r9,(L(ones)-L(anchor))
|
|
|
8ae002 |
+ sldi r4,r4,3
|
|
|
8ae002 |
+ lfdx fp0,r4,r5
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Are we in the primary interval of sin or cos? */
|
|
|
8ae002 |
+ andi. r4,r7,0x2
|
|
|
8ae002 |
+ bne L(cos)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Chebyshev polynomial of the form:
|
|
|
8ae002 |
+ x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp9,(L(S0)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp10,(L(S1)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp11,(L(S2)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp12,(L(S3)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp13,(L(S4)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp2,fp1,fp1 /* x^2 */
|
|
|
8ae002 |
+ fmul fp3,fp2,fp1 /* x^3 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp13,fp12 /* S3+x^2*S4 */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp11 /* S2+x^2*(S3+x^2*S4) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp10 /* S1+x^2*(S2+x^2*(S3+x^2*S4)) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp9 /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))) */
|
|
|
8ae002 |
+ fmadd fp4,fp3,fp4,fp1 /* x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))) */
|
|
|
8ae002 |
+ fmul fp4,fp4,fp0 /* Add in the sign. */
|
|
|
8ae002 |
+ frsp fp1,fp4 /* Round to single precision. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(cos):
|
|
|
8ae002 |
+ /* Chebyshev polynomial of the form:
|
|
|
8ae002 |
+ 1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))). */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp9,(L(C0)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp10,(L(C1)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp11,(L(C2)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp12,(L(C3)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp13,(L(C4)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp2,fp1,fp1 /* x^2 */
|
|
|
8ae002 |
+ lfd fp3,(L(DPone)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp13,fp12 /* C3+x^2*C4 */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp11 /* C2+x^2*(C3+x^2*C4) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp10 /* C1+x^2*(C2+x^2*(C3+x^2*C4)) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp9 /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))) */
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp4,fp3 /* 1.0 + x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))) */
|
|
|
8ae002 |
+ fmul fp4,fp4,fp0 /* Add in the sign. */
|
|
|
8ae002 |
+ frsp fp1,fp4 /* Round to single precision. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(greater_or_equal_9pio4):
|
|
|
8ae002 |
+ lis r4,INFINITY@h
|
|
|
8ae002 |
+ ori r4,r4,INFINITY@l
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ bge L(inf_or_nan)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lis r4,TWO_P23@h
|
|
|
8ae002 |
+ ori r4,r4,TWO_P23@l
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ bge L(greater_or_equal_2p23)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fabs fp1,fp1 /* |x| */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Calculate quotient of |x|/(PI/4). */
|
|
|
8ae002 |
+ lfd fp2,(L(invpio4)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp3,(L(DPone)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp4,(L(DPhalf)-L(anchor))(r9)
|
|
|
8ae002 |
+ fmul fp2,fp1,fp2 /* |x|/(PI/4) */
|
|
|
8ae002 |
+ friz fp2,fp2 /* n = floor(|x|/(PI/4)) */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Calculate (n + 1) / 2. */
|
|
|
8ae002 |
+ fadd fp2,fp2,fp3 /* n + 1 */
|
|
|
8ae002 |
+ fmul fp3,fp2,fp4 /* (n + 1) / 2 */
|
|
|
8ae002 |
+ friz fp3,fp3
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp4,(L(pio2hi)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp5,(L(pio2lo)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp6,fp4,fp3
|
|
|
8ae002 |
+ fadd fp6,fp6,fp1
|
|
|
8ae002 |
+ fmadd fp1,fp5,fp3,fp6
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fctiduz fp2,fp2
|
|
|
8ae002 |
+ mfvsrd r7,v2 /* n + 1 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ b L(reduced)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(inf_or_nan):
|
|
|
8ae002 |
+ bne L(skip_errno_setting) /* Is a NAN? */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* We delayed the creation of the stack frame, as well as the saving of
|
|
|
8ae002 |
+ the link register, because only at this point, we are sure that
|
|
|
8ae002 |
+ doing so is actually needed. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ stfd fp1,-8(r1)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Save the link register. */
|
|
|
8ae002 |
+ mflr r0
|
|
|
8ae002 |
+ std r0,16(r1)
|
|
|
8ae002 |
+ cfi_offset(lr, 16)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Create the stack frame. */
|
|
|
8ae002 |
+ stdu r1,-FRAMESIZE(r1)
|
|
|
8ae002 |
+ cfi_adjust_cfa_offset(FRAMESIZE)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ bl JUMPTARGET(__errno_location)
|
|
|
8ae002 |
+ nop
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Restore the stack frame. */
|
|
|
8ae002 |
+ addi r1,r1,FRAMESIZE
|
|
|
8ae002 |
+ cfi_adjust_cfa_offset(-FRAMESIZE)
|
|
|
8ae002 |
+ /* Restore the link register. */
|
|
|
8ae002 |
+ ld r0,16(r1)
|
|
|
8ae002 |
+ mtlr r0
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp1,-8(r1)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* errno = EDOM */
|
|
|
8ae002 |
+ li r4,EDOM
|
|
|
8ae002 |
+ stw r4,0(r3)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(skip_errno_setting):
|
|
|
8ae002 |
+ fsub fp1,fp1,fp1 /* x - x */
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(greater_or_equal_2p23):
|
|
|
8ae002 |
+ fabs fp1,fp1
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ srwi r4,r3,FLOAT_EXPONENT_SHIFT
|
|
|
8ae002 |
+ subi r4,r4,FLOAT_EXPONENT_BIAS
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* We reduce the input modulo pi/4, so we need 3 bits of integer
|
|
|
8ae002 |
+ to determine where in 2*pi we are. Index into our array
|
|
|
8ae002 |
+ accordingly. */
|
|
|
8ae002 |
+ addi r4,r4,INTEGER_BITS
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* To avoid an expensive divide, for the range we care about (0 - 127)
|
|
|
8ae002 |
+ we can transform x/28 into:
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ x/28 = (x * ((0x100000000 / 28) + 1)) >> 32
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ mulhwu returns the top 32 bits of the 64 bit result, doing the
|
|
|
8ae002 |
+ shift for us in the same instruction. The top 32 bits are undefined,
|
|
|
8ae002 |
+ so we have to mask them. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lis r6,FX_FRACTION_1_28@h
|
|
|
8ae002 |
+ ori r6,r6,FX_FRACTION_1_28@l
|
|
|
8ae002 |
+ mulhwu r5,r4,r6
|
|
|
8ae002 |
+ clrldi r5,r5,32
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Get our pointer into the invpio4_table array. */
|
|
|
8ae002 |
+ sldi r4,r5,3
|
|
|
8ae002 |
+ addi r6,r9,(L(invpio4_table)-L(anchor))
|
|
|
8ae002 |
+ add r4,r4,r6
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp2,0(r4)
|
|
|
8ae002 |
+ lfd fp3,8(r4)
|
|
|
8ae002 |
+ lfd fp4,16(r4)
|
|
|
8ae002 |
+ lfd fp5,24(r4)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp6,fp2,fp1
|
|
|
8ae002 |
+ fmul fp7,fp3,fp1
|
|
|
8ae002 |
+ fmul fp8,fp4,fp1
|
|
|
8ae002 |
+ fmul fp9,fp5,fp1
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Mask off larger integer bits in highest double word that we don't
|
|
|
8ae002 |
+ care about to avoid losing precision when combining with smaller
|
|
|
8ae002 |
+ values. */
|
|
|
8ae002 |
+ fctiduz fp10,fp6
|
|
|
8ae002 |
+ mfvsrd r7,v10
|
|
|
8ae002 |
+ rldicr r7,r7,0,(63-INTEGER_BITS)
|
|
|
8ae002 |
+ mtvsrd v10,r7
|
|
|
8ae002 |
+ fcfidu fp10,fp10 /* Integer bits. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fsub fp6,fp6,fp10 /* highest -= integer bits */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Work out the integer component, rounded down. Use the top two
|
|
|
8ae002 |
+ limbs for this. */
|
|
|
8ae002 |
+ fadd fp10,fp6,fp7 /* highest + higher */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fctiduz fp10,fp10
|
|
|
8ae002 |
+ mfvsrd r7,v10
|
|
|
8ae002 |
+ andi. r0,r7,1
|
|
|
8ae002 |
+ fcfidu fp10,fp10
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Subtract integer component from highest limb. */
|
|
|
8ae002 |
+ fsub fp12,fp6,fp10
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ beq L(even_integer)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Our integer component is odd, so we are in the -PI/4 to 0 primary
|
|
|
8ae002 |
+ region. We need to shift our result down by PI/4, and to do this
|
|
|
8ae002 |
+ in the mod (4/PI) space we simply subtract 1. */
|
|
|
8ae002 |
+ lfd fp11,(L(DPone)-L(anchor))(r9)
|
|
|
8ae002 |
+ fsub fp12,fp12,fp11
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Now add up all the limbs in order. */
|
|
|
8ae002 |
+ fadd fp12,fp12,fp7
|
|
|
8ae002 |
+ fadd fp12,fp12,fp8
|
|
|
8ae002 |
+ fadd fp12,fp12,fp9
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* And finally multiply by pi/4. */
|
|
|
8ae002 |
+ lfd fp13,(L(pio4)-L(anchor))(r9)
|
|
|
8ae002 |
+ fmul fp1,fp12,fp13
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ addi r7,r7,1
|
|
|
8ae002 |
+ b L(reduced)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(even_integer):
|
|
|
8ae002 |
+ lfd fp11,(L(DPone)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Now add up all the limbs in order. */
|
|
|
8ae002 |
+ fadd fp12,fp12,fp7
|
|
|
8ae002 |
+ fadd fp12,r12,fp8
|
|
|
8ae002 |
+ fadd fp12,r12,fp9
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* We need to check if the addition of all the limbs resulted in us
|
|
|
8ae002 |
+ overflowing 1.0. */
|
|
|
8ae002 |
+ fcmpu 0,fp12,fp11
|
|
|
8ae002 |
+ bgt L(greater_than_one)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* And finally multiply by pi/4. */
|
|
|
8ae002 |
+ lfd fp13,(L(pio4)-L(anchor))(r9)
|
|
|
8ae002 |
+ fmul fp1,fp12,fp13
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ addi r7,r7,1
|
|
|
8ae002 |
+ b L(reduced)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(greater_than_one):
|
|
|
8ae002 |
+ /* We did overflow 1.0 when adding up all the limbs. Add 1.0 to our
|
|
|
8ae002 |
+ integer, and subtract 1.0 from our result. Since that makes the
|
|
|
8ae002 |
+ integer component odd, we need to subtract another 1.0 as
|
|
|
8ae002 |
+ explained above. */
|
|
|
8ae002 |
+ addi r7,r7,1
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp11,(L(DPtwo)-L(anchor))(r9)
|
|
|
8ae002 |
+ fsub fp12,fp12,fp11
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* And finally multiply by pi/4. */
|
|
|
8ae002 |
+ lfd fp13,(L(pio4)-L(anchor))(r9)
|
|
|
8ae002 |
+ fmul fp1,fp12,fp13
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ addi r7,r7,1
|
|
|
8ae002 |
+ b L(reduced)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(less_2pn5):
|
|
|
8ae002 |
+ lis r4,TWO_PN27@h
|
|
|
8ae002 |
+ ori r4,r4,TWO_PN27@l
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ cmpw r3,r4
|
|
|
8ae002 |
+ blt L(less_2pn27)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* A simpler Chebyshev approximation is close enough for this range:
|
|
|
8ae002 |
+ x+x^3*(SS0+x^2*SS1). */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp10,(L(SS0)-L(anchor))(r9)
|
|
|
8ae002 |
+ lfd fp11,(L(SS1)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp2,fp1,fp1 /* x^2 */
|
|
|
8ae002 |
+ fmul fp3,fp2,fp1 /* x^3 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmadd fp4,fp2,fp11,fp10 /* SS0+x^2*SS1 */
|
|
|
8ae002 |
+ fmadd fp1,fp3,fp4,fp1 /* x+x^3*(SS0+x^2*SS1) */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ frsp fp1,fp1 /* Round to single precision. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(less_2pn27):
|
|
|
8ae002 |
+ cmpwi r3,0
|
|
|
8ae002 |
+ beq L(zero)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Handle some special cases:
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ sinf(subnormal) raises inexact/underflow
|
|
|
8ae002 |
+ sinf(min_normalized) raises inexact/underflow
|
|
|
8ae002 |
+ sinf(normalized) raises inexact. */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ lfd fp2,(L(small)-L(anchor))(r9)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ fmul fp2,fp1,fp2 /* x * small */
|
|
|
8ae002 |
+ fsub fp1,fp1,fp2 /* x - x * small */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ frsp fp1,fp1
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 16
|
|
|
8ae002 |
+L(zero):
|
|
|
8ae002 |
+ blr
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+END (__sinf)
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .section .rodata, "a"
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ .balign 8
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(anchor):
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Chebyshev constants for sin, range -PI/4 - PI/4. */
|
|
|
8ae002 |
+L(S0): .8byte 0xbfc5555555551cd9
|
|
|
8ae002 |
+L(S1): .8byte 0x3f81111110c2688b
|
|
|
8ae002 |
+L(S2): .8byte 0xbf2a019f8b4bd1f9
|
|
|
8ae002 |
+L(S3): .8byte 0x3ec71d7264e6b5b4
|
|
|
8ae002 |
+L(S4): .8byte 0xbe5a947e1674b58a
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Chebyshev constants for sin, range 2^-27 - 2^-5. */
|
|
|
8ae002 |
+L(SS0): .8byte 0xbfc555555543d49d
|
|
|
8ae002 |
+L(SS1): .8byte 0x3f8110f475cec8c5
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+ /* Chebyshev constants for cos, range -PI/4 - PI/4. */
|
|
|
8ae002 |
+L(C0): .8byte 0xbfdffffffffe98ae
|
|
|
8ae002 |
+L(C1): .8byte 0x3fa55555545c50c7
|
|
|
8ae002 |
+L(C2): .8byte 0xbf56c16b348b6874
|
|
|
8ae002 |
+L(C3): .8byte 0x3efa00eb9ac43cc0
|
|
|
8ae002 |
+L(C4): .8byte 0xbe923c97dd8844d7
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(invpio2):
|
|
|
8ae002 |
+ .8byte 0x3fe45f306dc9c883 /* 2/PI */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(invpio4):
|
|
|
8ae002 |
+ .8byte 0x3ff45f306dc9c883 /* 4/PI */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(invpio4_table):
|
|
|
8ae002 |
+ .8byte 0x0000000000000000
|
|
|
8ae002 |
+ .8byte 0x3ff45f306c000000
|
|
|
8ae002 |
+ .8byte 0x3e3c9c882a000000
|
|
|
8ae002 |
+ .8byte 0x3c54fe13a8000000
|
|
|
8ae002 |
+ .8byte 0x3aaf47d4d0000000
|
|
|
8ae002 |
+ .8byte 0x38fbb81b6c000000
|
|
|
8ae002 |
+ .8byte 0x3714acc9e0000000
|
|
|
8ae002 |
+ .8byte 0x3560e4107c000000
|
|
|
8ae002 |
+ .8byte 0x33bca2c756000000
|
|
|
8ae002 |
+ .8byte 0x31fbd778ac000000
|
|
|
8ae002 |
+ .8byte 0x300b7246e0000000
|
|
|
8ae002 |
+ .8byte 0x2e5d2126e8000000
|
|
|
8ae002 |
+ .8byte 0x2c97003248000000
|
|
|
8ae002 |
+ .8byte 0x2ad77504e8000000
|
|
|
8ae002 |
+ .8byte 0x290921cfe0000000
|
|
|
8ae002 |
+ .8byte 0x274deb1cb0000000
|
|
|
8ae002 |
+ .8byte 0x25829a73e0000000
|
|
|
8ae002 |
+ .8byte 0x23fd1046be000000
|
|
|
8ae002 |
+ .8byte 0x2224baed10000000
|
|
|
8ae002 |
+ .8byte 0x20709d338e000000
|
|
|
8ae002 |
+ .8byte 0x1e535a2f80000000
|
|
|
8ae002 |
+ .8byte 0x1cef904e64000000
|
|
|
8ae002 |
+ .8byte 0x1b0d639830000000
|
|
|
8ae002 |
+ .8byte 0x1964ce7d24000000
|
|
|
8ae002 |
+ .8byte 0x17b908bf16000000
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(pio4):
|
|
|
8ae002 |
+ .8byte 0x3fe921fb54442d18 /* PI/4 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+/* PI/2 as a sum of two doubles. We only use 32 bits of the upper limb
|
|
|
8ae002 |
+ to avoid losing significant bits when multiplying with up to
|
|
|
8ae002 |
+ (2^22)/(pi/2). */
|
|
|
8ae002 |
+L(pio2hi):
|
|
|
8ae002 |
+ .8byte 0xbff921fb54400000
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(pio2lo):
|
|
|
8ae002 |
+ .8byte 0xbdd0b4611a626332
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(pio2_table):
|
|
|
8ae002 |
+ .8byte 0
|
|
|
8ae002 |
+ .8byte 0x3ff921fb54442d18 /* 1 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x400921fb54442d18 /* 2 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x4012d97c7f3321d2 /* 3 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x401921fb54442d18 /* 4 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x401f6a7a2955385e /* 5 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x4022d97c7f3321d2 /* 6 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x4025fdbbe9bba775 /* 7 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x402921fb54442d18 /* 8 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x402c463abeccb2bb /* 9 * PI/2 */
|
|
|
8ae002 |
+ .8byte 0x402f6a7a2955385e /* 10 * PI/2 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(small):
|
|
|
8ae002 |
+ .8byte 0x3cd0000000000000 /* 2^-50 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(ones):
|
|
|
8ae002 |
+ .8byte 0x3ff0000000000000 /* +1.0 */
|
|
|
8ae002 |
+ .8byte 0xbff0000000000000 /* -1.0 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(DPhalf):
|
|
|
8ae002 |
+ .8byte 0x3fe0000000000000 /* 0.5 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(DPone):
|
|
|
8ae002 |
+ .8byte 0x3ff0000000000000 /* 1.0 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+L(DPtwo):
|
|
|
8ae002 |
+ .8byte 0x4000000000000000 /* 2.0 */
|
|
|
8ae002 |
+
|
|
|
8ae002 |
+weak_alias(__sinf, sinf)
|
|
|
8ae002 |
--
|
|
|
8ae002 |
2.1.0
|
|
|
8ae002 |
|