1 /* Compute cubic root of float value.
2 Copyright (C) 1997-2021 Free Software Foundation, Inc.
3 This file is part of the GNU C Library.
5 The GNU C Library is free software; you can redistribute it and/or
6 modify it under the terms of the GNU Lesser General Public
7 License as published by the Free Software Foundation; either
8 version 2.1 of the License, or (at your option) any later version.
10 The GNU C Library is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 Lesser General Public License for more details.
15 You should have received a copy of the GNU Lesser General Public
16 License along with the GNU C Library; if not, see
17 <https://www.gnu.org/licenses/>. */
19 #include <machine/asm.h>
20 #include <libm-alias-float.h>
26 f3: .double 0.191502161678719066
27 ASM_SIZE_DIRECTIVE(f3)
29 f2: .double 0.697570460207922770
30 ASM_SIZE_DIRECTIVE(f2)
32 f1: .double 0.492659620528969547
33 ASM_SIZE_DIRECTIVE(f1)
35 #define CBRT2 1.2599210498948731648
36 #define ONE_CBRT2 0.793700525984099737355196796584
37 #define SQR_CBRT2 1.5874010519681994748
38 #define ONE_SQR_CBRT2 0.629960524947436582364439673883
42 factor: .double ONE_SQR_CBRT2
47 ASM_SIZE_DIRECTIVE(factor)
50 two25: .byte 0, 0, 0, 0x4c
51 ASM_SIZE_DIRECTIVE(two25)
54 #define MO(op) op##@GOTOFF(%ebx)
55 #define MOX(op,x) op##@GOTOFF(%ebx,x,1)
58 #define MOX(op,x) op(x)
66 andl $0x7fffffff, %eax
68 cmpl $0x7f800000, %eax
73 cfi_adjust_cfa_offset (4)
74 cfi_rel_offset (ebx, 0)
78 cmpl $0x00800000, %eax
96 andl $0x7fffffff, %eax
99 andl $0x807fffff, %edx
101 orl $0x3f000000, %edx
106 flds 8(%esp) /* xm */
110 flds 4(%esp) /* xm */
114 /* The following code has two tracks:
115 a) compute the normalized cbrt value
116 b) compute xe/3 and xe%3
117 The right track computes the value for b) and this is done
118 in an optimized way by avoiding division.
120 But why two tracks at all? Very easy: efficiency. Some FP
121 instruction can overlap with a certain amount of integer (and
122 FP) instructions. So we get (except for the imull) all
123 instructions for free. */
125 fld %st(0) /* xm : xm */
126 fmull MO(f3) /* f3*xm : xm */
127 movl $1431655766, %eax
128 fsubrl MO(f2) /* f2-f3*xm : xm */
130 fmul %st(1) /* (f2-f3*xm)*xm : xm */
132 faddl MO(f1) /* u:=f1+(f2-f3*xm)*xm : xm */
134 fld %st /* u : u : xm */
136 fmul %st(1) /* u*u : u : xm */
137 fld %st(2) /* xm : u*u : u : xm */
138 fadd %st /* 2*xm : u*u : u : xm */
139 fxch %st(1) /* u*u : 2*xm : u : xm */
140 fmul %st(2) /* t2:=u*u*u : 2*xm : u : xm */
142 fadd %st, %st(1) /* t2 : t2+2*xm : u : xm */
143 leal (%edx,%edx,2),%edx
144 fadd %st(0) /* 2*t2 : t2+2*xm : u : xm */
146 faddp %st, %st(3) /* t2+2*xm : u : 2*t2+xm */
148 fmulp /* u*(t2+2*xm) : 2*t2+xm */
149 fdivp %st, %st(1) /* u*(t2+2*xm)/(2*t2+xm) */
150 fmull MOX(16+factor,%ecx) /* u*(t2+2*xm)/(2*t2+xm)*FACT */
152 cfi_adjust_cfa_offset (4)
153 fildl (%esp) /* xe/3 : u*(t2+2*xm)/(2*t2+xm)*FACT */
154 fxch /* u*(t2+2*xm)/(2*t2+xm)*FACT : xe/3 */
155 fscale /* u*(t2+2*xm)/(2*t2+xm)*FACT*2^xe/3 */
157 cfi_adjust_cfa_offset (-4)
161 cfi_adjust_cfa_offset (-4)
172 /* Return the argument. */
176 libm_alias_float (__cbrt, cbrt)