/* Copyright (C) 2011 DJ Delorie, see COPYING.DJ for details */
/* Copyright (C) 1999 DJ Delorie, see COPYING.DJ for details */
#include <libc/asmdefs.h>

	.data
pinf:
	.long	0xFF800000

NaN:
	.long	0xFFC00000

temp:
	.long	0, 0

onethird:
	.long	1431655765

two54:
	.long	0x5A800000

a0:
	.float	+1.87277957900533E+00

a1:
	.float	-1.87243905326548E+00

a2:
	.float	+1.60286399719912E+00

a3:
	.float	-7.46198924594210E-01

a4:
	.float	+1.42994392730009E-01

b0:
	.float	14.

b1:
	.float	-7.

b2:
	.float	+2.

one9th:
	.tfloat +0.11111111111111111111

	.text
	.globl	_cbrt
_cbrt:

	movl	8(%esp), %eax
	movl	%eax, %ecx		/* Save sign */

	andl	$0x7FFFFFFF, %eax	/* fabs */
	movl	%eax, 8(%esp)

	cmpl	$0x7FF00000, %eax	/* Control flows straight through for */
	jae	abarg			/* normal args: 0 < fabs(x) < +inf */
	testl	$0x7FF00000, %eax
	jz	verysmall

	mull	onethird
	addl	$0x2A9F7893, %edx
	movl	%edx, temp+4		/* First approximation good */
					/* to 5.5 bits */

have55:
	fldl	4(%esp)
	fld1
	fdivp				/* recip */

	fldl	temp			/* Load approximation */
					/* 4rd-order minimax to 24 bits */
	fld	%st(0)			/* x	x   recip	*/
	fmul	%st(1)			/* x^2	x   recip	*/
	fmul	%st(1)			/* x^3	x   recip	*/
	fmul	%st(2)			/* y	x   recip	*/
	fld	%st(0)			/* y	y   x	  recip */
	fmuls	a4			/* P1'	y   x	  recip */
	fadds	a3			/* P1	y   x	  recip */
	fmul	%st(1)			/* P2'	y   x	  recip */
	fadds	a2			/* P2	y   x	  recip */
	fmul	%st(1)			/* P3'	y   x	  recip */
	fadds	a1			/* P3	y   x	  recip */
	fmulp				/* P4'	x   recip	*/
	fadds	a0			/* P4	x   recip	*/
	fmulp				/* x'	recip		*/
					/* 2nd-order Taylor to 64 bits */
	fld	%st(0)			/* x	x   recip */
	fmul	%st(1)			/* x^2	x   recip */
	fmul	%st(1)			/* x^3	x   recip */
	fmul	%st(2)			/* y	x   recip */
	ffree	%st(2)			/* y	x	  */
	fld	%st(0)			/* y	y   x	  */
	fmuls	b2
	fadds	b1
	fmulp				/* ccc	x */
	fadds	b0			/* P(y) x */
	fmulp				/* x''	  */
	fldt	one9th
	fmulp

cleanup:				/* Restore sign */
	testl	%ecx, %ecx
	jns	end
	fchs

end:
	ret

verysmall:				/* Exponent is 0 */
	movl	8(%esp), %eax
	testl	%eax, %eax
	jnz	denormal
	movl	4(%esp), %eax
	testl	%eax, %eax
	jz	special			/* x = 0 */
	
denormal:
	fldl	4(%esp)
	fmuls	two54			/* Multiply by 2^54 to normalize */
	fstpl	temp

	movl	temp+4, %eax
	mull	onethird
	addl	$0x297F7893, %edx	/* Undo 2^54 multiplier */
	movl	%edx, temp+4		/* First approximation to 5.5 bits */
	movl	$0, temp

	jmp	have55

abarg:					/* x = inf, or NaN */
	testl	$0x000FFFFF, %eax
	jnz	badarg
	movl	4(%esp), %eax
	testl	%eax, %eax
	jz	special

badarg:					/* arg is negative or NaN */
	movl	$1, _errno
	flds	NaN
	ret

special:
	fldl	4(%esp)			/* x = 0 or inf: just load x */
	jmp	cleanup