1/*	$NetBSD: n_cbrt.S,v 1.9 2024/05/07 15:15:31 riastradh Exp $	*/
2/*
3 * Copyright (c) 1985, 1993
4 *	The Regents of the University of California.  All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 * 1. Redistributions of source code must retain the above copyright
10 *    notice, this list of conditions and the following disclaimer.
11 * 2. Redistributions in binary form must reproduce the above copyright
12 *    notice, this list of conditions and the following disclaimer in the
13 *    documentation and/or other materials provided with the distribution.
14 * 3. Neither the name of the University nor the names of its contributors
15 *    may be used to endorse or promote products derived from this software
16 *    without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28 * SUCH DAMAGE.
29 *
30 *	@(#)cbrt.s	8.1 (Berkeley) 6/4/93
31 */
32
33#include <machine/asm.h>
34
35/*
36 * double cbrt(double arg)
37 * W. Kahan, 10/13/80. revised 1/13/84 for keeping sign symmetry
38 * error check by E LeBlanc, 8/18/82
39 * Revised and tested by K.C. Ng, 5/2/85
40 * Max error less than 0.667 ulps (unit in the last places)
41 */
42
43ENTRY(cbrtf, 0)
44	cvtfd	4(%ap),-(%sp)
45	calls	$2,_C_LABEL(cbrt)
46	cvtdf	%r0,%r0
47	ret
48END(cbrtf)
49
50#ifdef WEAK_ALIAS
51WEAK_ALIAS(cbrtl, cbrt)
52#endif
53
54STRONG_ALIAS(cbrt, d_cbrt)
55ENTRY(d_cbrt, 0x00c0)		# save %r6 & %r7
56	movq	4(%ap),%r0	# %r0 = argument x
57	jbr 	dcbrt2
58END(d_cbrt)
59
60ENTRY(dcbrt_, 0x00c0)		# save %r6 & %r7
61	movq	*4(%ap),%r0	# %r0 = argument x
62
63dcbrt2:	bicw3	$0x807f,%r0,%r2	# biased exponent of x
64	jeql	return		# dcbrt(0)=0  dcbrt(res)=res. operand
65	bicw3	$0x7fff,%r0,%ap	# ap has sign(x)
66	xorw2	%ap,%r0		# %r0 is abs(x)
67	movl	%r0,%r2		# %r2 has abs(x)
68	rotl	$16,%r2,%r2	# %r2 = |x| with bits unscrambled
69	divl2	$3,%r2		# rough dcbrt with bias/3
70	addl2	B,%r2		# restore bias, diminish fraction
71	rotl	$16,%r2,%r2	# %r2=|q|=|dcbrt| to 5 bits
72	mulf3	%r2,%r2,%r3	# %r3 =qq
73	divf2	%r0,%r3		# %r3 = qq/x
74	mulf2	%r2,%r3
75	addf2	C,%r3		# %r3 = s = C + qqq/x
76	divf3	%r3,D,%r4		# %r4 = D/s
77	addf2	E,%r4
78	addf2	%r4,%r3		# %r3 = s + E + D/s
79	divf3	%r3,F,%r3		# %r3 = F / (s + E + D/s)
80	addf2	G,%r3		# %r3 = G + F / (s + E + D/s)
81	mulf2	%r3,%r2		# %r2 = q%r3 = new q to 23 bits
82	clrl	%r3		# %r2:%r3 = q as double float
83	muld3	%r2,%r2,%r4	# %r4:%r5 = qq exactly
84	divd2	%r4,%r0		# %r0:%r1 = x/(q*q) rounded
85	subd3	%r2,%r0,%r6	# %r6:%r7 = x/(q*q) - q exactly
86	movq    %r2,%r4		# %r4:%r5 = q
87	addw2	$0x80,%r4	# %r4:%r5 = 2 * q
88	addd2	%r0,%r4		# %r4:%r5 = 2*q + x/(q*q)
89	divd2	%r4,%r6		# %r6:%r7 = (x/(q*q)-q)/(2*q+x/(q*q))
90	muld2	%r2,%r6		# %r6:%r7 = q*(x/(q*q)-q)/(2*q+x/(q*q))
91	addd3	%r6,%r2,%r0	# %r0:%r1 = q + %r6:%r7
92	bisw2	%ap,%r0		# restore the sign bit
93return:
94	ret			# error less than 0.667 ulps
95END(dcbrt_)
96
97	_ALIGN_TEXT
98B :	.long		 721142941		# (86-0.03306235651)*(2^23)
99C :	.float		0f0.5428571429		# 19/35
100D :	.float		0f-0.7053061224		# -864/1225
101E :	.float		0f1.414285714		# 99/70
102F :	.float		0f1.607142857		# 45/28
103G :	.float		0f0.3571428571		# 5/14
104