emufncalib.s

Vxworks OS source code

/* Copyright 1991-2001 Wind River Systems, Inc. */

/*
modification history
--------------------
01e,28aug01,hdn  added FUNC/FUNC_LABEL,	replaced .align with .balign
01d,16apr98,hdn  fixed comments which cause assembler error.
01c,03feb98,hdn  fixed sin/cos problem (SPR7764) w USSW v.311
01b,24oct96,yp   stuck in a # in USSC mailing addr so cpp will work
01a,24jan95,hdn  original US Software version.
*/

#
#       Filename:   EMUFNC.ASM
#
#       Copyright (C) 1990,1991 By
#       # United States Software Corporation
#       # 14215 N.W. Science Park Drive
#       # Portland, Oregon 97229
#
#       This software is furnished under a license and may be used
#       and copied only in accordance with the terms of such license
#       and with the inclusion of the above copyright notice.
#       This software or any other copies thereof may not be provided
#       or otherwise made available to any other person.  No title to
#       and ownership of the software is hereby transferred.
#
#       The information in this software is subject to change without
#       notice and should not be construed as a commitment by United
#       States Software Corporation.
#
#       Version:        V3.09  01/19/95 GOFAST  WRS GNU version of GF-PROT
#       Released:       1 March 1991
#
#

#define _ASMLANGUAGE
#include "vxWorks.h"
#include "asm.h"
	
	.data
	.globl	FUNC(copyright_wind_river)
	.long	FUNC(copyright_wind_river)

	.include	"emuIncALib.s"
	.text

#        extrn   wiep:near
#        extrn   tnormal:near,shiftr:near,exproc:near,exprop:near
#        extrn   epmul:near
#        extrn   epadd1:near,epmul1:near,epdiv1:near
#        extrn   retnan:near,getqn:near,chkarg:near,exret:near

#DGROUP  GROUP   hwseg
#CGROUP  GROUP   emuseg

#hwseg   segment RW para public 'DATA'
#        assume  ds:DGROUP, es:DGROUP, ss:DGROUP
#        extrn   h_ctrl:word,h_stat:word
#hwseg   ends

#emuseg  segment public USE32 'CODE'
#        assume  cs:CGROUP

__lib:

#        subttl  (2 to power x) - 1
#        page

# raise 2 to power x and subtract 1 (x between -1 and 1)
# input:        EP in cx:edx:eax
# output:       result in cx:edx:eax

	.globl	epf2xm1

	.balign	16
expcon:	.word	0,0xf717,0x0001,0x0000,0x0000,0x0000 # ln2^18/18!
	.word	0,0x0889,0x0033,0x0000,0x0000,0x0000 # ln2^17/17!
	.word	0,0xa26b,0x04e3,0x0000,0x0000,0x0000 # ln2^16/16!
	.word	0,0xa10e,0x70db,0x0000,0x0000,0x0000 # ln2^15/15!
	.word	0,0x26ac,0x8a4b,0x0009,0x0000,0x0000 # ln2^14/14!
	.word	0,0x8b36,0xb0c9,0x00c0,0x0000,0x0000 # ln2^13/13!
	.word	0,0x7e14,0xeb28,0x0e1d,0x0000,0x0000 # ln2^12/12!
	.word	0,0x8dd9,0x639a,0xf465,0x0000,0x0000 # ln2^11/11!
	.word	0,0xc764,0x8ac5,0x267a,0x000f,0x0000 # ln2^10/10!
	.word	0,0x3e1e,0x9caf,0x929e,0x00da,0x0000 # ln2^9/9!
	.word	0,0x11fe,0xd2e4,0x0111,0x0b16,0x0000 # ln2^8/8!
	.word	0,0x1a1b,0x22c3,0xff16,0x7ff2,0x0000 # ln2^7/7!
	.word	0,0xdbd2,0xb1e1,0x4be1,0x0c24,0x0005 # ln2^6/6!
	.word	0,0x20e2,0xce62,0xcf14,0xb0ff,0x002b # ln2^5/5!
	.word	0,0x9ccb,0xe772,0xfba4,0x2ab6,0x013b # ln2^4/4!
	.word	0,0xcce9,0x2fe2,0xc128,0xc235,0x071a # ln2^3/3!
	.word	0,0x6f16,0x8ea8,0x82c5,0xbdff,0x1ebf # ln2^2/2!
	.word	0,0xe4f1,0xbcd5,0xe8e7,0x0bfb,0x58b9 # ln2^1/1!
	.word	-1

exp87:	orl	%ecx,%ecx	#argument + or - 1
	jns	exp89		#   +1 returns +1
	decl	%ecx		#   -1 returns -1/2
exp89:	ret			

expw:	btl	$16,%ecx	#zero returns as such
	jc	exp89		

	call	chkarg		#weed out illegals
	jc	exp89		

	cmpw	$0x7fff,%cx	#-INF returns -1, +INF as such
	jnz	exp2		
	orl	%ecx,%ecx	
	jns	exp89		
	movl	$0x80003fff,%ecx	
	ret			

epf2xm1:			
	testl	$0x00030000,%ecx	#jump if special value
	jnz	expw		

exp2:	#PREEXS  precis          #set precision exception
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL1		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL1:				

	cmpw	$0x3fff,%cx	#handle limit cases +1 and -1
	jz	exp87		

	movl	$expcon,%ebx	#call the polynomial routine
	call	sigmax		
	addw	$0x3ffe,%cx	

	call	epmul		#then multiply with x
	ret			

# sin: trigonometric sine of x
# cos: trigonometric cosine of x
# tan: trigonometric tangent of x
# input:        x in ecx:edx:eax and ST
# return:       answer in ecx:edx:eax

	.globl	epsin,epcos,eptan

	.balign	16
sincon:				#coefficients for sin(x)
	.word	0,0x1dc0,0x654b,0x0000,0x0000,0x0000
	.word	0,0x6332,0x6030,0xff94,0xffff,0xffff
	.word	0,0xa1b4,0x184e,0x5849,0x0000,0x0000
	.word	0,0x763a,0x3015,0x3375,0xffca,0xffff
	.word	0,0x338e,0x56c7,0xe3a5,0x171d,0x0000
	.word	0,0x7f98,0x97f9,0xf97f,0x7f97,0xfff9
	.word	0,0x1111,0x1111,0x1111,0x1111,0x0111
	.word	0,0xaaab,0xaaaa,0xaaaa,0xaaaa,0xeaaa
	.word	0,0x0000,0x0000,0x0000,0x0000,0x8000
	.word	-1

	.balign	16
coscon:				#coefficients for cos(x)
	.word	0,0x61b8,0xfa5f,0xffff,0xffff,0xffff
	.word	0,0xf9cc,0xb9fc,0x0006,0x0000,0x0000
	.word	0,0xcfe2,0xa2d5,0xf9b1,0xffff,0xffff
	.word	0,0x3624,0x3bfe,0x7bb6,0x0004,0x0000
	.word	0,0x1472,0x10ec,0x3609,0xfdb0,0xffff
	.word	0,0xd00d,0x0d00,0x00d0,0xd00d,0x0000
	.word	0,0x7d28,0x27d2,0xd27d,0x7d27,0xffd2
	.word	0,0x5555,0x5555,0x5555,0x5555,0x0555
	.word	0,0x0000,0x0000,0x0000,0x0000,0xc000
	.word	0,0x0000,0x0000,0x0000,0x0000,0x8000
	.word	-1

	.balign	16
tancon:				#coefficients for tan(x)
	.word	0,0x49fa,0xa843,0xffff,0xffff,0xffff
	.word	0,0xd55c,0x9e11,0xfffc,0xffff,0xffff
	.word	0,0x5074,0x9d9a,0xffde,0xffff,0xffff
	.word	0,0xcc5e,0x827f,0xfeb6,0xffff,0xffff
	.word	0,0xb1e6,0x0f07,0xf34c,0xffff,0xffff
	.word	0,0xd54e,0x9a5c,0x82a0,0xffff,0xffff
	.word	0,0xb45a,0x38a6,0x2a9e,0xfffb,0xffff
	.word	0,0xd548,0xc10c,0x4b63,0xffd0,0xffff
	.word	0,0xe9e2,0xb1f6,0x24d3,0xfe29,0xffff
	.word	0,0x810a,0xbbda,0xfa29,0xedd7,0xffff
	.word	0,0x12ec,0x478a,0x86a0,0x4cab,0xffff
	.word	0,0x54e8,0xc43d,0x1b32,0x10a2,0xfff9
	.word	0,0x5100,0xaa65,0x0ffb,0xa655,0xffba
	.word	0,0xd27d,0x7d27,0x27d2,0xd27d,0xfd27
	.word	0,0x5555,0x5555,0x5555,0x5555,0xd555
	.word	0,0x0000,0x0000,0x0000,0x0000,0x8000
	.word	-1

epsin:				
	#PREEXS  precis          #not precise
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL2		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL2:				

	cmpb	$4,%bl		#sin(x+pi) = -sinx
	jc	sin2		#   so if octant > 3 we flip the sign
	xorb	$0x80,11(%edi)	#   and subtract 4 from octant number
	andb	$-5,%bl		

sin2:	jpo	cos3		#if octant is 1 or 2 we take cosx

sin3:	movl	%eax,(%edi)	#store reduced value, leave sign
	movl	%edx,4(%edi)	
	movw	%cx,8(%edi)	

	call	square		#we use x^2 in the polynomial

	lea	sincon,%ebx	#do the polynomial
	call	sigmax		
	addw	$0x3ffe,%cx	

	call	epmul1		#multiply with x
	ret			

epcos:				
	#PREEXS  precis          #not precise
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL3		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL3:				

	andl	$0xffff,%ecx	#store for possible use, sign cleared
	movl	%ecx,8(%edi)	

	testb	$2,%bl		#we flip the sign in octants 2 and 3
	jz	cos2		
	xorb	$0x80,11(%edi)	

cos2:	cmpb	$4,%bl		#cos(x+pi) = -cosx
	jc	cos2b		#   so if octant > 3 we flip the sign
	xorb	$0x80,11(%edi)	#   and subtract 4 from octant number
	andb	$-5,%bl		

cos2b:	jpo	sin3		#if octant is 1 or 2 we take sinx

cos3:	call	square		#we use x^2 in the polynomial

	lea	coscon,%ebx	#do the polynomial
	call	sigmax		
	addw	$0x3ffe,%cx	

	movw	$0,8(%edi)	#put in the sign
	orl	8(%edi),%ecx	
	ret			

eptan:				
	#PREEXS  precis          #not precise
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL4		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL4:				
	movl	%eax,(%edi)	
	movl	%edx,4(%edi)	
	movw	%cx,8(%edi)	

	andb	$-5,%bl		#tan(x+pi) = tan(x) 

	testb	$2,%bl		#if octant is 2 or 3 flip sign
	jz	tan3		
	xorb	$0x80,11(%edi)	

tan3:	call	square		#we use x^2 in the polynomial

	pushl	%ebx		#save the octant

	lea	tancon,%ebx	#do the polynomial
	call	sigmax		
	addw	$0x3ffe,%cx	

	popl	%ebx		#tan(x) = x/SUM
	orb	%bl,%bl		#   if octant is 1 or 2 we take 1/tan(x)
	jpo	tan7		
	xchgw	8(%edi),%cx	
	xchgl	4(%edi),%edx	
	xchgl	0(%edi),%eax	
tan7:	call	epdiv1		

	ret			#return

#        subttl   base 2 logarithm of x
#        page

# take base 2 logarithm of x
# input:        x in cx:edx:eax and ST
# output:       log2(x) in cx:edx:eax

	.globl	eplog,eplog1

	.balign	16
logcon:	.word	0,0x26b2,0xff99,0xd343,0x0f07,0x04dc #log2(e)/19
	.word	0,0xd0e5,0xc350,0xdd0f,0x6b26,0x056e #log2(e)/17
	.word	0,0xb98d,0xff7d,0xa533,0xcec5,0x0627 #log2(e)/15
	.word	0,0x388f,0xd807,0x34c5,0x3d5a,0x071a #log2(e)/13
	.word	0,0x2b91,0x1694,0xca01,0xd424,0x0864 #log2(e)/11
	.word	0,0x3540,0x54c7,0xbe01,0x589e,0x0a42 #log2(e)/9
	.word	0,0xb22e,0x6c9f,0x3d6f,0xbb15,0x0d30 #log2(e)/7
	.word	0,0x2ca7,0xfe79,0xef9b,0x6c50,0x1277 #log2(e)/5
	.word	0,0x9fc1,0xfd74,0x3a03,0x09dc,0x1ec7 #log2(e)/3
	.word	0,0xdf44,0xf85d,0xae0b,0x1d94,0x5c55 #log2(e)
	.word	-1

	.balign	16
flnlog:	.word	0x0000,0x0000,0x0000,0x8000,0x3fff,0x8000 #log2(0.5) 
	.word	0xf0c1,0x0852,0xcb8c,0xd47f,0x3ffe,0x8000 #log2(0.5625)
	.word	0xd407,0xcb91,0x1ed0,0xad96,0x3ffe,0x8000 #log2(0.625)
	.word	0xc407,0x3457,0xb07f,0x8a62,0x3ffe,0x8000 #log2(0.6875)
	.word	0xf0c1,0x0852,0xcb8c,0xd47f,0x3ffd,0x8000 #log2(0.75)
	.word	0x432d,0x8773,0xf71b,0x995f,0x3ffd,0x8000 #log2(0.8125)
	.word	0x9333,0xfde9,0xc055,0xc544,0x3ffc,0x8000 #log2(0.875)
	.word	0x633a,0x7dda,0x24b6,0xbeb0,0x3ffb,0x8000 #log2(0.9375)
	.word	0x0000,0x0000,0x0000,0x0000,0x0000,1 #log2(1)
	.word	0x3cfd,0xdeb4,0xd1cf,0xae00,0x3ffc,0 #log2(1.125)
	.word	0x57f2,0x68dc,0xc25e,0xa4d3,0x3ffd,0 #log2(1.25)
	.word	0x77f2,0x9750,0x9f01,0xeb3a,0x3ffd,0 #log2(1.375)
	.word	0x87a0,0xfbd6,0x1a39,0x95c0,0x3ffe,0 #log2(1.5)
	.word	0x5e69,0x3c46,0x0472,0xb350,0x3ffe,0 #log2(1.625)
	.word	0x9b33,0x8085,0xcfea,0xceae,0x3ffe,0 #log2(1.75)
	.word	0xb399,0x3044,0xfb69,0xe829,0x3ffe,0 #log2(1.875)

logzer:	orb	$zerod,h_stat	#return -INF
	call	exproc		
	btsl	$31,%edx	
	jmp	log88		

lognan:	orb	$invop,h_stat	#return NaN
	call	exproc		
	movl	$0xc0000000,%edx	
log88:	movl	$0x80027fff,%ecx	
	xorl	%eax,%eax	
log89:	ret			

logw:	btl	$16,%ecx	#zero returns -INF
	jc	logzer		

	call	chkarg		#weed out illegals
	jc	log89		

	orl	%ecx,%ecx	#negative illegal
	js	lognan		

	cmpw	$0x7fff,%cx	#+INF returns as such
	jz	log89		

	call	tnormal		#handle denormal
	#PUTEP   
	movl	%eax,(%edi)	
	movl	%edx,4(%edi)	
	movl	%ecx,8(%edi)	

eplog:				
	testl	$0x80030000,%ecx	#jump if special value
	jnz	logw		

log2:	movl	%edx,%ebx	#precision exception unless power of 2
	addl	%ebx,%ebx	
	orl	%eax,%ebx	
	jz	log3		
	#PREEXS  precis
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL5		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL5:				

log3:	movl	%ecx,%ebx	#we make the exponent 0 or -1
	cmpw	$0x3fff,%cx	#   depending on if x > 1 or < 1
	movw	$0x3fff,%cx	#   push the needed adjustment
	jge	log31		
	decl	%ecx		
log31:	movw	%cx,8(%edi)	
	subl	%ecx,%ebx	
	movswl	%bx,%ebx	
	pushl	%ebx		

	xorl	%ebx,%ebx	#calculate area number (0-15) 
	shldl	$4,%edx,%ebx	
	testb	$1,%cl		
	jnz	log5		
	andb	$7,%bl		

log5:	cmpb	$7,%bl		#handle 0.9375 < x < specially
	jnz	log21		
	incl	%ebx		
	pushl	%ebx		
	stc			
	rcrl	$1,%edx		
	rcrl	$1,%eax		
	incw	%cx		
	xchgl	(%edi),%eax	
	xchgl	4(%edi),%edx	
	xchgw	8(%edi),%cx	
	notl	%edx		
	negl	%eax		
	sbbl	$-1,%edx	
	xorl	$0x80000000,%ecx	
	jmp	log51		

log21:	pushl	%ebx		#interpolation area number

	movl	%edx,%ebx	#registers = x+C
	andl	$0xf0000000,%ebx	
	addl	%ebx,%edx	
	rcrl	$1,%edx		
	rcrl	$1,%eax		
	incw	%cx		

	xchgl	(%edi),%eax	#swap arguments
	xchgl	4(%edi),%edx	
	xchgw	8(%edi),%cx	

	andl	$0x0fffffff,%edx	#registers = x-C
log54:	jnz	log51		
	orl	%eax,%eax	
	jnz	log51		
	xorl	%ecx,%ecx	
	jmp	log57		
log51:	decw	%cx		
	addl	%eax,%eax	
	adcl	%edx,%edx	
	jns	log51		
				
log52:	call	epdiv1		#calculate d = (M-C)/(M+C)
	addb	%bl,%bl		
	adcl	$0,%eax		
	adcl	$0,%edx		
	#PUTEP   
	movl	%eax,(%edi)	
	movl	%edx,4(%edi)	
	movl	%ecx,8(%edi)	

	call	square		#square that

	lea	logcon,%ebx	#calculate polynomial for d^2
	call	sigmax		
	addb	$0x40,%ch	

	call	epmul1		#then multiply by d
	addb	%bl,%bl		
	adcl	$0,%eax		
	adcl	$0,%edx		

log57:	popl	%ebx		#add base log2(C)
	shll	$2,%ebx		
	lea	flnlog(%ebx,%ebx,2),%ebx	
	movl	%cs:(%ebx),%esi	
	movl	%esi,(%edi)	
	movl	%cs:4(%ebx),%esi	
	movl	%esi,4(%edi)	
	movl	%cs:8(%ebx),%esi	
	movl	%esi,8(%edi)	
	call	epadd1		
	addb	%bl,%bl		
	adcl	$0,%eax		
	adcl	$0,%edx		
	#PUTEP   
	movl	%eax,(%edi)	
	movl	%edx,4(%edi)	
	movl	%ecx,8(%edi)	

	popl	%edx		#get exponent adjustment
	call	wiep		#   convert to EP, add to log2(d)
	call	epadd1		
	addb	%bl,%bl		
	adcl	$0,%eax		
	adcl	$0,%edx		
log9:	ret			

lgpw61:	orl	%ecx,%ecx	#-INF returns NaN
	js	lognan		
	ret			

lgpw:	btl	$16,%ecx	#zero returns 0
	jc	log9		

	call	chkarg		#weed out illegals
	jc	log9		

	cmpw	$0x7fff,%cx	#INF jumps
	jz	lgpw61		

lgpw8:	#PREEXS  precis          #very small multiplied by log2(e)
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL6		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL6:				
	movl	$0x5C17F0BC,0(%edi)	
	movl	$0xB8AA3B29,4(%edi)	
	movl	$0x3fff,8(%edi)	
	call	epmul		
	ret			

eplog1:				
	testl	$0x00030000,%ecx	#jump if special value
	jnz	lgpw		

	cmpb	$0x10,%ch	#very small skip some code
	jc	lgpw8		

	#PREEXS  precis
	orb	$precis,h_stat	
	testb	$precis,h_ctrl	
	jnz	LL7		
	orw	$0x8080,h_stat	
	movl	$exret,-4(%ebp)	
LL7:				

	xorl	%ebx,%ebx	#push zero exponent adjustment
	pushl	%ebx		

	orl	%ecx,%ecx	#go handle x < 0 separately
	js	lgp30		

	subw	$tBIAS-4,%cx	#calculate area number (8-15) 
	cmpw	$9,%cx		
	jnc	lgp3		
	shldl	%cl,%edx,%ebx	
lgp3:	addb	$8,%bl		
	pushl	%ebx		

	movw	$0x4000,%cx	#calculate x+1+C
	subw	8(%edi),%cx	
	call	shiftr		
	btsl	$30,%edx	
	movl	%edx,%ebx	
	andl	$0x78000000,%ebx	
	addl	%ebx,%edx	
	movl	$0x4000,%ecx	

	xchgl	(%edi),%eax	#swap arguments
	xchgl	4(%edi),%edx	
	xchgl	8(%edi),%ecx	

	movl	%ecx,%ebx	#registers = x-C
	addw	$4-tBIAS,%bx	
	jle	log52		
	xchgl	%ecx,%ebx	
	shll	%cl,%edx	
	shrl	%cl,%edx	
	xchgl	%ecx,%ebx	
	jmp	log54		
				
lgp30:	subw	$tBIAS-5,%cx	#calculate area number (0-7) 
	cmpw	$9,%cx		
	jnc	lgp33		
	shldl	%cl,%edx,%ebx	
	negl	%ebx		
lgp33:	addl	$8,%ebx		
	pushl	%ebx		

	movw	$0x3fff,%cx	#calculate x+1+C
	subw	8(%edi),%cx	
	call	shiftr