/*---------------------------------------------------------------------------+
 |  polynomial.S                                                             |
 |                                                                           |
 | Fixed point arithmetic polynomial evaluation.                             |
 |                                                                           |
 | Copyright (C) 1992,1993                                                   |
 |                       W. Metzenthen, 22 Parker St, Ormond, Vic 3163,      |
 |                       Australia.  E-mail   billm@vaxc.cc.monash.edu.au    |
 |                                                                           |
 | Call from C as:                                                           |
 |   void polynomial(unsigned accum[], unsigned x[], unsigned terms[][2],    |
 |                   int n)                                                  |
 |                                                                           |
 | Computes:                                                                 |
 | terms[0] + (terms[1] + (terms[2] + ... + (terms[n-1]*x)*x)*x)*x) ... )*x  |
 | The result is returned in accum.                                          |
 |                                                                           |
 +---------------------------------------------------------------------------*/

	.file	"fpolynom.s"

#include "fpu_asm.h"


/*	#define	EXTRA_PRECISE	// Do not use: not complete */

#define	TERM_SIZE	$8
#define	SUM_MS		-20(%ebp)	/* sum ms long */
#define SUM_MIDDLE	-24(%ebp)	/* sum middle long */
#define	SUM_LS		-28(%ebp)	/* sum ls long */
#define	SUM_LS_HI	-25(%ebp)	/* high byte of sum ls */
#define	ACCUM_MS	-4(%ebp)	/* accum ms long */
#define	ACCUM_MIDDLE	-8(%ebp)	/* accum middle long */
#define	ACCUM_LS	-12(%ebp)	/* accum ls long */
#define ACCUM_LS_HI	-9(%ebp)	/* high byte of accum ls */

.text
	.align 2,144
.globl _polynomial
_polynomial:
	pushl	%ebp
	movl	%esp,%ebp
	subl	$32,%esp
	pushl	%esi
	pushl	%edi
	pushl	%ebx

	movl	PARAM2,%esi		/* x */
	movl	PARAM3,%edi		/* terms */

	movl	TERM_SIZE,%eax
	mull	PARAM4			/* n */
	addl	%eax,%edi

	movl	4(%edi),%edx		/* terms[n] */
	movl	%edx,SUM_MS
	movl	(%edi),%edx		/* terms[n] */
	movl	%edx,SUM_MIDDLE
	xor	%eax,%eax
	movl	%eax,SUM_LS

	subl	TERM_SIZE,%edi
	decl	PARAM4
	js	L_accum_done

L_accum_loop:
	xor	%eax,%eax
	movl	%eax,ACCUM_MS
	movl	%eax,ACCUM_MIDDLE

	movl	SUM_MIDDLE,%eax
	mull	(%esi)			/* x ls long */
/*	movl	%eax,-16(%ebp)		// Not needed */
	movl	%edx,ACCUM_LS

	movl	SUM_MIDDLE,%eax
	mull	4(%esi)			/* x ms long */
	addl	%eax,ACCUM_LS
	adcl	%edx,ACCUM_MIDDLE
	adcl	$0,ACCUM_MS

	movl	SUM_MS,%eax
	mull	(%esi)			/* x ls long */
	addl	%eax,ACCUM_LS
	adcl	%edx,ACCUM_MIDDLE
	adcl	$0,ACCUM_MS

	movl	SUM_MS,%eax
	mull	4(%esi)			/* x ms long */
	addl	%eax,ACCUM_MIDDLE
	adcl	%edx,ACCUM_MS

/*
 * Now put the sum of next term and the accumulator
 * into the sum register
 */
	movl	ACCUM_MIDDLE,%eax
	addl	(%edi),%eax		/* term ls long */
	movl	%eax,SUM_MIDDLE
	movl	ACCUM_MS,%eax
	adcl	4(%edi),%eax		/* term ms long */
	movl	%eax,SUM_MS

#ifdef EXTRA_PRECISE
	movl	ACCUM_LS,%eax
	movl	%eax,SUM_LS
#else
	testb	$0x80,ACCUM_LS_HI	/* ms bit of ACCUM_LS */
	je	L_no_poly_round

	addl	$1,SUM_MIDDLE
	adcl	$0,SUM_MS
L_no_poly_round:
#endif EXTRA_PRECISE

	subl	TERM_SIZE,%edi
	decl	PARAM4
	jns	L_accum_loop

L_accum_done:
#ifdef EXTRA_PRECISE
/* Round the result */
	testb	$128,SUM_LS_HI
	je	L_poly_done

	addl	$1,SUM_MIDDLE
	adcl	$0,SUM_MS
#endif EXTRA_PRECISE

L_poly_done:
	movl	PARAM1,%edi		/* accum */
	movl	SUM_MIDDLE,%eax
	movl	%eax,(%edi)
	movl	SUM_MS,%eax
	movl	%eax,4(%edi)

	popl	%ebx
	popl	%edi
	popl	%esi
	leave
	ret