/*
   libffm	- Free, pretty fast replacement for some math (libm) routines 
			on Linux/AXP, optimized for the 21164

   Copyright (C) 1998  Joachim Wesner <joachim.wesner@frankfurt.netsurf.de>
                  and  Kazushige Goto <goto@statabo.rim.or.jp>

   This library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Library General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU Library General Public
   License along with this library (see file COPYING.LIB); if not, write 
   to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, 
   MA 02139, USA.
*/
/*
   Fast sin/cos approximation including range reduction by Joachim Wesner,
   <joachim.wesner@frankfurt.netsurf.de>, see also the work published in 
   mc 8/1991 p. 78-93. Done in June 1998.

   No special handling of illegal arguments or NANs, yet !!!
*/

/*
  improved and modified by Kazushige Goto <goto@statabo.rim.or.jp>

  note:

   This program's target is Alpha 21164.
   It takes only 88 clocks at most in any parameters.
   But this reduction routine can not accept NaN, Inf or extra
   Huge values.  To get faster code(to fit 21164), I modified
   last routine.

  originally,

  x = (((((((C[0]*yy + C[1])*yy + C[2])*yy + C[3])*yy +
          C[4])*yy + C[5])*yy + C[6])*yy + C[7])*x;

  I modified to

  r1 = K[8] * yy;  r1 = r1 + K[12]; z1 = yy * yy;
  r2 = K[9] * yy;  r2 = r2 + K[13]; z2 = z1 * z1;
  r3 = K[10] * yy; r3 = r3 + K[14]; z3 = z2 * z1;
  r4 = K[11] * yy; r4 = r4 + K[15];

  r1 = z3 * r1;  r2 = z2 * r2;  r3 = z1 * r3;
  r1 = r1 + r2;  r1 = r1 + r3;  r1 = r1 + r4;
  x = r1 * x;

  This routine has slightly error.  My aim is to make routine as fast
  as I can.  If you worry about such error, please use glibc's generic
  sin/cos routine.  But its error is only 1-2 bit of LSB. Most people
  need not worry about.  This technique is provided by
  Naohiko Shimizu <nshimizu@et.u-tokai.ac.jp>).

  Also, I programmed to speed up for first access.  Different from
  hardware function routine(i.e. intel x86, Mips),  we must take
  into instruction and data movement from main memory.  At a first
  access,  alpha waits very-very long time(about 1000 to 10000 clocks,
  it depends on the coding).  This is not owing to me, the users must
  take consideration. You must try to call continuously.

*/

/*
   further modified to use a higher order approximation by Joachim Wesner
   (joachim.wesner@frankfurt.netsurf.de).
*/

	.set noreorder
	.set noat

#ifdef __ELF__
	.section .rodata
#else
	.rdata
#endif

	.align 5

        /* Derived from Hart & Cheney, "Computer Approximations" */
        /* Approximation #3345, rescaled,  18.6 digits precision */
        /*   old # 3344 only had Precision 15.9 digits */

K:
	.t_floating  3.18309886183790691216e-1
	.t_floating  1.00000000000000000000e0
	.t_floating  5.00000000000000000000e-1
	.t_floating  3.14159202575683593750e0
	.t_floating  6.27832832833519205451e-7
	.t_floating  1.24467443437932267676e-13
	.t_floating  2.72047909578888629421e-15
	.t_floating  1.60589364903715969799e-10
	.t_floating  2.75573192101527645153e-6
	.t_floating  8.33333333333316668423e-3
	.t_floating -7.64291780689105089624e-13
	.t_floating -2.50521067982745964283e-8
	.t_floating -1.98412698412018482239e-4
	.t_floating -1.66666666666666657415e-1
	.t_floating  1.00000000000000000000e0

.text
	.align 5

	.globl sin
	.ent sin
sin:
	lda	$30, -16($30)
	ldgp	$29,.-sin($27)
	lda	$1,K			# address of K
	.frame	$30,16,$26,0

#ifdef PROF
	lda	$28, _mcount
	jsr	$28, ($28), _mcount
	unop
	unop
#endif
	.prologue 1

 	ldt	$f10,   0($1)		# K[0]
	ldt	$f24,  64($1)		# K[8]
	ldt	$f14,  32($1)		# K[4]
	ldt	$f28,  96($1)		# K[12]

	mult	$f16, $f10, $f0		# z = x * K[0]
	ldt	$f11,   8($1)		# K[1]
	clr	$0
	br	$31, $common
	.end sin

	.align 4

	.globl cos
	.ent cos
cos:
	lda	$30, -16($30)
	ldgp	$29,.-cos($27)
	lda	$1,K			# address of K
	.frame	$30,16,$26,0

#ifdef PROF
	lda	$28, _mcount
	jsr	$28, ($28), _mcount
	unop
	unop
#endif
	.prologue 1

 	ldt	$f10,   0($1)		# K[0]
	ldt	$f24,  64($1)		# K[8]
	ldt	$f14,  32($1)		# K[4]
	ldt	$f28,  96($1)		# K[12]

	mult	$f16, $f10, $f0		# z = x * K[0]
	ldt	$f11,   8($1)		# K[1]
	mov	1, $0
	unop
	.align 4

$common:
	mult	$f16,$f16,$f17		# yy = x * x
	ldt	$f12,  16($1)		# K[2]
	ldt	$f13,  24($1)		# K[3]
	ldt	$f15,  40($1)		# K[5]

	addt	$f0,  $f11, $f19	# r1 = z + K[1]
 	ldt	$f23,  56($1)		# K[7]
	addt	$f0,  $f12, $f20 	# r2 = z + K[2]
	ldt	$f25,  72($1)		# K[9]

	ldt	$f22,  48($1)		# K[6]
	subt	$f0,  $f12, $f21	# r3 = z - K[2]
	fblt	$f16, $negative

	cvttqc	$f19, $f18

	cvttqc	$f20, $f20
	br	$31, $skip

$negative:
	cvttqc	$f0,  $f18
	cvttqc	$f21, $f20

$skip:
	stt	$f18, 8($30)
	cvtqt	$f18, $f18
	stt	$f20, 0($30)
	cvtqt	$f20, $f0

	beq	$0, $sine
	subt	$f18, $f12, $f0
$sine:
	ldt	$f26,  80($1)		# K[10]
	ldt	$f27,  88($1)		# K[11]

	ldt	$f29, 104($1)		# K[13]
	mult	$f0,  $f13, $f13	# z * K[3]
	ldt	$f30, 112($1)		# K[14]
	fbeq	$f0, $skip_reduction

	mult	$f0,  $f14, $f14	# z * K[4]
	mult	$f0,  $f15, $f15	# z * K[5]
	subt	$f16, $f13, $f13	# x - z*K[3]
	subt	$f13, $f14, $f13	# (x-z*K[3]) - z*K[4]

	subt	$f13, $f15, $f16	# ((x-z*K[3]) - z*K[4]) - z*K[5]
	mult	$f16,$f16,$f17		# yy = x * x
	unop
	unop
	.align 4

$skip_reduction:
	mult	$f17, $f17, $f18	# z1 = yy * yy
	mult	$f24, $f17, $f24	# r3 = K[8] * yy
	mult	$f23, $f17, $f23	# r2 = K[7] * yy
	mult	$f25, $f17, $f25	# r4 = K[9] * yy

	mult	$f18, $f17, $f19	# z2 = z1 * yy
	mult	$f22, $f17, $f22	# r1 = K[6] * yy
	addt	$f24, $f28, $f24	# r3 = r3 + K[12]
	addt	$f23, $f27, $f23	# r2 = r2 + K[11]

	addq	$30, 8, $2
	addt	$f25, $f29, $f25	# r4 = r4 + K[13]
	mult	$f19, $f18, $f20	# z3 = z2 * z1
	addt	$f22, $f26, $f22	# r1 = r1 + K[10]

	cmoveq	$0, $30, $2
	mult	$f19, $f24, $f24	# r3 = z2 * r3
	mult	$f20, $f18, $f21	# z4 = z3 * z1
	mult	$f20, $f23, $f23	# r2 = z3 * r2

	ldq	$4,  0($2)
	mult	$f17, $f25, $f25	# r4 = yy * r4
	mult	$f21, $f22, $f22	# r1 = z4 * r1
	addt	$f22, $f23, $f22	# r1 = r1 + r2
	addt	$f25, $f30, $f25	# r4 = r4 + K[14]

	addt	$f22, $f24, $f22	# r1 = r1 + r3
	blbc	$4,$skip_negative
	fneg	$f16, $f16		# x = -x
$skip_negative:
	addt	$f22, $f25, $f22	# r1 = r1 + r4

	mult	$f22, $f16, $f0		# r1 * x
	lda	$30, 16($30)
	ret	$31, ($26), 1
	.end cos