xref: /src/lib/libm/arch/vax/n_support.S

/*	$NetBSD: n_support.S,v 1.13 2024/07/17 12:00:48 riastradh Exp $	*/
/*
 * Copyright (c) 1985, 1993
 *	The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 *	@(#)support.s	8.1 (Berkeley) 6/4/93
 */
#include <machine/asm.h>

WEAK_ALIAS(logbl,logb)
WEAK_ALIAS(copysignl, _copysignl)
WEAK_ALIAS(_copysignl, copysign)

	.text
_sccsid:
	.asciz "@(#)support.s\t1.3 (Berkeley) 8/21/85; 8.1 (ucb.elefunt) 6/4/93"

/*
 * copysign(x,y),
 * logb(x),
 * scalb(x,N),
 * finite(x),
 * drem(x,y),
 * Coded in vax assembly language by K.C. Ng,  3/14/85.
 * Revised by K.C. Ng on 4/9/85.
 */

/*
 * double copysign(double x,double y)
 */

ENTRY(copysign, 0)
	movq	4(%ap),%r0		# load x into %r0
	bicw3	$0x807f,%r0,%r2		# mask off the exponent of x
	beql	Lz			# if zero or reserved op then return x
	bicw3	$0x7fff,12(%ap),%r2	# copy the sign bit of y into %r2
	bicw2	$0x8000,%r0		# replace x by |x|
	bisw2	%r2,%r0			# copy the sign bit of y to x
Lz:	ret
END(copysign)

ENTRY(copysignf, 0)
	movl	4(%ap),%r0		# load x into %r0
	bicw3	$0x807f,%r0,%r2		# mask off the exponent of x
	beql	1f			# if zero or reserved op then return x
	bicw3	$0x7fff,8(%ap),%r2	# copy the sign bit of y into %r2
	bicw2	$0x8000,%r0		# replace x by |x|
	bisw2	%r2,%r0			# copy the sign bit of y to x
1:	ret
END(copysignf)

/*
 * float logbf(float x);
 */
ENTRY(logbf, 0)
	cvtfd   4(%ap),-(%sp)
	calls   $2,_C_LABEL(logb)
	cvtdf   %r0,%r0
	ret
END(logbf)

/*
 * double logb(double x);
 */
ENTRY(logb, 0)
	bicl3	$0xffff807f,4(%ap),%r0	# mask off the exponent of x
	beql    Ln
	ashl	$-7,%r0,%r0		# get the bias exponent
	subl2	$129,%r0			# get the unbias exponent
	cvtld	%r0,%r0			# return the answer in double
	ret
Ln:	movq	4(%ap),%r0		# %r0:1 = x (zero or reserved op)
	bneq	1f			# simply return if reserved op
	movq 	$0x0000fe00ffffcfff,%r0  # -2147483647.0
1:	ret
END(logb)

/*
 * long finite(double x);
 */
#ifndef __GFLOAT__
WEAK_ALIAS(finitef, _finitef)
STRONG_ALIAS(_finitef, _finite)
#endif
WEAK_ALIAS(finite, _finite)
ENTRY(_finite, 0)
	bicw3	$0x7f,4(%ap),%r0	# mask off the significand
	cmpw	%r0,$0x8000		# to see if x is the reserved op
	beql	1f			# if so, return FALSE (0)
	movl	$1,%r0			# else return TRUE (1)
	ret
1:	clrl	%r0
	ret
END(_finite)

/* int isnan(double x);
 */
#if 0
ENTRY(isnan, 0)
	clrl	%r0
	ret
#endif

/* int isnanf(float x);
 */
ENTRY(isnanf, 0)
	clrl	%r0
	ret
END(isnanf)

/*
 * double scalb(x,N)
 * double x; double N;
 */
	.set	ERANGE,34

ENTRY(scalb, 0)
	movq	4(%ap),%r0
	bicl3	$0xffff807f,%r0,%r3
	beql	ret1			# 0 or reserved operand
	movq	12(%ap),%r4
	cvtdl	%r4, %r2
	cmpl	%r2,$0x12c
	bgeq	ovfl
	cmpl	%r2,$-0x12c
	bleq	unfl
	ashl	$7,%r2,%r2
	addl2	%r2,%r3
	bleq	unfl
	cmpl	%r3,$0x8000
	bgeq	ovfl
	addl2	%r2,%r0
	ret
ovfl:	pushl	$ERANGE
	calls	$1,_C_LABEL(infnan)	# if it returns
	bicw3	$0x7fff,4(%ap),%r2	# get the sign of input arg
	bisw2	%r2,%r0			# re-attach the sign to %r0/1
	ret
unfl:	movq	$0,%r0
ret1:	ret
END(scalb)

/*
 * DREM(X,Y)
 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
 * DOUBLE PRECISION (VAX D format 56 bits)
 * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85.
 */
	.set	EDOM,33

ENTRY(drem, 0x0fc0)
	subl2	$12,%sp
	movq	4(%ap),%r0		#%r0=x
	movq	12(%ap),%r2		#%r2=y
	jeql	Rop			#if y=0 then generate reserved op fault
	bicw3	$0x007f,%r0,%r4		#check if x is Rop
	cmpw	%r4,$0x8000
	jeql	Ret			#if x is Rop then return Rop
	bicl3	$0x007f,%r2,%r4		#check if y is Rop
	cmpw	%r4,$0x8000
	jeql	Ret			#if y is Rop then return Rop
	bicw2	$0x8000,%r2		#y  := |y|
	movw	$0,-4(%fp)		#-4(%fp) = nx := 0
	cmpw	%r2,$0x1c80		#yexp ? 57
	bgtr	C1			#if yexp > 57 goto C1
	addw2	$0x1c80,%r2		#scale up y by 2**57
	movw	$0x1c80,-4(%fp)		#nx := 57 (exponent field)
C1:
	movw	-4(%fp),-8(%fp)		#-8(%fp) = nf := nx
	bicw3	$0x7fff,%r0,-12(%fp)	#-12(%fp) = sign of x
	bicw2	$0x8000,%r0		#x  := |x|
	movq	%r2,%r10			#y1 := y
	bicl2	$0xffff07ff,%r11		#clear the last 27 bits of y1
loop:
	cmpd	%r0,%r2			#x ? y
	bleq	E1			#if x <= y goto E1
 /* begin argument reduction */
	movq	%r2,%r4			#t =y
	movq	%r10,%r6			#t1=y1
	bicw3	$0x807f,%r0,%r8		#xexp= exponent of x
	bicw3	$0x807f,%r2,%r9		#yexp= exponent fo y
	subw2	%r9,%r8			#xexp-yexp
	subw2	$0x0c80,%r8		#k=xexp-yexp-25(exponent bit field)
	blss	C2			#if k<0 goto C2
	addw2	%r8,%r4			#t +=k
	addw2	%r8,%r6			#t1+=k, scale up t and t1
C2:
	divd3	%r4,%r0,%r8		#x/t
	cvtdl	%r8,%r8			#n=[x/t] truncated
	cvtld	%r8,%r8			#float(n)
	subd2	%r6,%r4			#t:=t-t1
	muld2	%r8,%r4			#n*(t-t1)
	muld2	%r8,%r6			#n*t1
	subd2	%r6,%r0			#x-n*t1
	subd2	%r4,%r0			#(x-n*t1)-n*(t-t1)
	jbr	loop
E1:
	movw	-4(%fp),%r6		#%r6=nx
	beql	C3			#if nx=0 goto C3
	addw2	%r6,%r0			#x:=x*2**57 scale up x by nx
	movw	$0,-4(%fp)		#clear nx
	jbr	loop
C3:
	movq	%r2,%r4			#%r4 = y
	subw2	$0x80,%r4		#%r4 = y/2
	cmpd	%r0,%r4			#x:y/2
	blss	E2			#if x < y/2 goto E2
	bgtr	C4			#if x > y/2 goto C4
	cvtdl	%r8,%r8			#ifix(float(n))
	blbc	%r8,E2			#if the last bit is zero, goto E2
C4:
	subd2	%r2,%r0			#x-y
E2:
	xorw2	-12(%fp),%r0		#x^sign (exclusive or)
	movw	-8(%fp),%r6		#%r6=nf
	bicw3	$0x807f,%r0,%r8		#%r8=exponent of x
	bicw2	$0x7f80,%r0		#clear the exponent of x
	subw2	%r6,%r8			#%r8=xexp-nf
	bgtr	C5			#if xexp-nf is positive goto C5
	movw	$0,%r8			#clear %r8
	movq	$0,%r0			#x underflow to zero
C5:
	bisw2	%r8,%r0			/* put %r8 into x's exponent field */
	ret
Rop:					#Reserved operand
	pushl	$EDOM
	calls	$1,_C_LABEL(infnan)	#generate reserved op fault
	ret
Ret:
	movq	$0x8000,%r0		#propagate reserved op
	ret
END(drem)