Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
4973 right-hear 1 
#include
2 
 
3 
/* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */
4 
	.data
5 
pinf:
6 
	.long	0xFF800000
7 
 
8 
NaN:
9 
	.long	0xFFC00000
10 
 
11 
temp:
12 
	.long	0, 0
13 
 
14 
onethird:
15 
	.long	1431655765
16 
 
17 
two54:
18 
	.long	0x5A800000
19 
 
20 
a0:
21 
	.float	+1.87277957900533E+00
22 
 
23 
a1:
24 
	.float	-1.87243905326548E+00
25 
 
26 
a2:
27 
	.float	+1.60286399719912E+00
28 
 
29 
a3:
30 
	.float	-7.46198924594210E-01
31 
 
32 
a4:
33 
	.float	+1.42994392730009E-01
34 
 
35 
b0:
36 
	.float	14.
37 
 
38 
b1:
39 
	.float	-7.
40 
 
41 
b2:
42 
	.float	+2.
43 
 
44 
one9th:
45 
	.tfloat +0.11111111111111111111
46 
 
47 
	.text
48 
MK_C_SYM(cbrt)
49 
 
50 
	movl	8(%esp), %eax
51 
	movl	%eax, %ecx		/* Save sign */
52 
 
53 
	andl	$0x7FFFFFFF, %eax	/* fabs */
54 
	movl	%eax, 8(%esp)
55 
 
56 
	cmpl	$0x7FF00000, %eax	/* Control flows straight through for */
57 
	jae	abarg			/* normal args: 0 < fabs(x) < +inf */
58 
	testl	$0x7FF00000, %eax
59 
	jz	verysmall
60 
 
61 
	mull	onethird
62 
	addl	$0x2A9F7893, %edx
63 
	movl	%edx, temp+4		/* First approximation good */
64 
					/* to 5.5 bits */
65 
 
66 
have55:
67 
	fldl	4(%esp)
68 
	fld1
69 
	fdivp				/* recip */
70 
 
71 
	fldl	temp			/* Load approximation */
72 
					/* 4rd-order minimax to 24 bits */
73 
	fld	%st(0)			/* x	x   recip	*/
74 
	fmul	%st(1)			/* x^2	x   recip	*/
75 
	fmul	%st(1)			/* x^3	x   recip	*/
76 
	fmul	%st(2)			/* y	x   recip	*/
77 
	fld	%st(0)			/* y	y   x	  recip */
78 
	fmuls	a4			/* P1'	y   x	  recip */
79 
	fadds	a3			/* P1	y   x	  recip */
80 
	fmul	%st(1)			/* P2'	y   x	  recip */
81 
	fadds	a2			/* P2	y   x	  recip */
82 
	fmul	%st(1)			/* P3'	y   x	  recip */
83 
	fadds	a1			/* P3	y   x	  recip */
84 
	fmulp				/* P4'	x   recip	*/
85 
	fadds	a0			/* P4	x   recip	*/
86 
	fmulp				/* x'	recip		*/
87 
					/* 2nd-order Taylor to 64 bits */
88 
	fld	%st(0)			/* x	x   recip */
89 
	fmul	%st(1)			/* x^2	x   recip */
90 
	fmul	%st(1)			/* x^3	x   recip */
91 
	fmul	%st(2)			/* y	x   recip */
92 
	ffree	%st(2)			/* y	x	  */
93 
	fld	%st(0)			/* y	y   x	  */
94 
	fmuls	b2
95 
	fadds	b1
96 
	fmulp				/* ccc	x */
97 
	fadds	b0			/* P(y) x */
98 
	fmulp				/* x''	  */
99 
	fldt	one9th
100 
	fmulp
101 
 
102 
cleanup:				/* Restore sign */
103 
	testl	%ecx, %ecx
104 
	jns	end
105 
	fchs
106 
 
107 
end:
108 
	ret
109 
 
110 
verysmall:				/* Exponent is 0 */
111 
	movl	8(%esp), %eax
112 
	testl	%eax, %eax
113 
	jnz	denormal
114 
	movl	4(%esp), %eax
115 
	testl	%eax, %eax
116 
	jz	special			/* x = 0 */
117 
 
118 
denormal:
119 
	fldl	4(%esp)
120 
	fmuls	two54			/* Multiply by 2^54 to normalize */
121 
	fstpl	temp
122 
 
123 
	movl	temp+4, %eax
124 
	mull	onethird
125 
	addl	$0x297F7893, %edx	/* Undo 2^54 multiplier */
126 
	movl	%edx, temp+4		/* First approximation to 5.5 bits */
127 
	movl	$0, temp
128 
 
129 
	jmp	have55
130 
 
131 
abarg:					/* x = inf, or NaN */
132 
	testl	$0x000FFFFF, %eax
133 
	jnz	badarg
134 
	movl	4(%esp), %eax
135 
	testl	%eax, %eax
136 
	jz	special
137 
 
138 
badarg:					/* arg is negative or NaN */
139 
	movl	$1, C_SYM(errno)
140 
	flds	NaN
141 
	ret
142 
 
143 
special:
144 
	fldl	4(%esp)			/* x = 0 or inf: just load x */
145 
	jmp	cleanup