1/* 2 * Copyright (c) 1985 Regents of the University of California. 3 * All rights reserved. 4 * 5 * Redistribution and use in source and binary forms are permitted 6 * provided that this notice is preserved and that due credit is given 7 * to the University of California at Berkeley. The name of the University 8 * may not be used to endorse or promote products derived from this 9 * software without specific prior written permission. This software 10 * is provided ``as is'' without express or implied warranty. 11 * 12 * All recipients should regard themselves as participants in an ongoing 13 * research project and hence should feel obligated to report their 14 * experiences (good or bad) with these elementary function codes, using 15 * the sendbug(8) program, to the authors. 16 * 17 * @(#)support.s 5.2 (Berkeley) 04/29/88 18 */ 19 .data 20 .align 2 21_sccsid: 22.asciz "@(#)support.s 1.3 (Berkeley) 8/21/85; 5.2 (ucb.elefunt) 04/29/88" 23 24/* 25 * copysign(x,y), 26 * logb(x), 27 * scalb(x,N), 28 * finite(x), 29 * drem(x,y), 30 * Coded in vax assembly language by K.C. Ng, 3/14/85. 31 * Revised by K.C. Ng on 4/9/85. 32 */ 33 34/* 35 * double copysign(x,y) 36 * double x,y; 37 */ 38 .globl _copysign 39 .text 40 .align 1 41_copysign: 42 .word 0x4 43 movq 4(ap),r0 # load x into r0 44 bicw3 $0x807f,r0,r2 # mask off the exponent of x 45 beql Lz # if zero or reserved op then return x 46 bicw3 $0x7fff,12(ap),r2 # copy the sign bit of y into r2 47 bicw2 $0x8000,r0 # replace x by |x| 48 bisw2 r2,r0 # copy the sign bit of y to x 49Lz: ret 50 51/* 52 * double logb(x) 53 * double x; 54 */ 55 .globl _logb 56 .text 57 .align 1 58_logb: 59 .word 0x0 60 bicl3 $0xffff807f,4(ap),r0 # mask off the exponent of x 61 beql Ln 62 ashl $-7,r0,r0 # get the bias exponent 63 subl2 $129,r0 # get the unbias exponent 64 cvtld r0,r0 # return the answer in double 65 ret 66Ln: movq 4(ap),r0 # r0:1 = x (zero or reserved op) 67 bneq 1f # simply return if reserved op 68 movq $0x0000fe00ffffcfff,r0 # -2147483647.0 691: ret 70 71/* 72 * long finite(x) 73 * double x; 74 */ 75 .globl _finite 76 .text 77 .align 1 78_finite: 79 .word 0x0000 80 bicw3 $0x7f,4(ap),r0 # mask off the mantissa 81 cmpw r0,$0x8000 # to see if x is the reserved op 82 beql 1f # if so, return FALSE (0) 83 movl $1,r0 # else return TRUE (1) 84 ret 851: clrl r0 86 ret 87 88/* 89 * double scalb(x,N) 90 * double x; int N; 91 */ 92 .globl _scalb 93 .set ERANGE,34 94 .text 95 .align 1 96_scalb: 97 .word 0xc 98 movq 4(ap),r0 99 bicl3 $0xffff807f,r0,r3 100 beql ret1 # 0 or reserved operand 101 movl 12(ap),r2 102 cmpl r2,$0x12c 103 bgeq ovfl 104 cmpl r2,$-0x12c 105 bleq unfl 106 ashl $7,r2,r2 107 addl2 r2,r3 108 bleq unfl 109 cmpl r3,$0x8000 110 bgeq ovfl 111 addl2 r2,r0 112 ret 113ovfl: pushl $ERANGE 114 calls $1,_infnan # if it returns 115 bicw3 $0x7fff,4(ap),r2 # get the sign of input arg 116 bisw2 r2,r0 # re-attach the sign to r0/1 117 ret 118unfl: movq $0,r0 119ret1: ret 120 121/* 122 * DREM(X,Y) 123 * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE) 124 * DOUBLE PRECISION (VAX D format 56 bits) 125 * CODED IN VAX ASSEMBLY LANGUAGE BY K.C. NG, 4/8/85. 126 */ 127 .globl _drem 128 .set EDOM,33 129 .text 130 .align 1 131_drem: 132 .word 0xffc 133 subl2 $12,sp 134 movq 4(ap),r0 #r0=x 135 movq 12(ap),r2 #r2=y 136 jeql Rop #if y=0 then generate reserved op fault 137 bicw3 $0x007f,r0,r4 #check if x is Rop 138 cmpw r4,$0x8000 139 jeql Ret #if x is Rop then return Rop 140 bicl3 $0x007f,r2,r4 #check if y is Rop 141 cmpw r4,$0x8000 142 jeql Ret #if y is Rop then return Rop 143 bicw2 $0x8000,r2 #y := |y| 144 movw $0,-4(fp) #-4(fp) = nx := 0 145 cmpw r2,$0x1c80 #yexp ? 57 146 bgtr C1 #if yexp > 57 goto C1 147 addw2 $0x1c80,r2 #scale up y by 2**57 148 movw $0x1c80,-4(fp) #nx := 57 (exponent field) 149C1: 150 movw -4(fp),-8(fp) #-8(fp) = nf := nx 151 bicw3 $0x7fff,r0,-12(fp) #-12(fp) = sign of x 152 bicw2 $0x8000,r0 #x := |x| 153 movq r2,r10 #y1 := y 154 bicl2 $0xffff07ff,r11 #clear the last 27 bits of y1 155loop: 156 cmpd r0,r2 #x ? y 157 bleq E1 #if x <= y goto E1 158 /* begin argument reduction */ 159 movq r2,r4 #t =y 160 movq r10,r6 #t1=y1 161 bicw3 $0x807f,r0,r8 #xexp= exponent of x 162 bicw3 $0x807f,r2,r9 #yexp= exponent fo y 163 subw2 r9,r8 #xexp-yexp 164 subw2 $0x0c80,r8 #k=xexp-yexp-25(exponent bit field) 165 blss C2 #if k<0 goto C2 166 addw2 r8,r4 #t +=k 167 addw2 r8,r6 #t1+=k, scale up t and t1 168C2: 169 divd3 r4,r0,r8 #x/t 170 cvtdl r8,r8 #n=[x/t] truncated 171 cvtld r8,r8 #float(n) 172 subd2 r6,r4 #t:=t-t1 173 muld2 r8,r4 #n*(t-t1) 174 muld2 r8,r6 #n*t1 175 subd2 r6,r0 #x-n*t1 176 subd2 r4,r0 #(x-n*t1)-n*(t-t1) 177 brb loop 178E1: 179 movw -4(fp),r6 #r6=nx 180 beql C3 #if nx=0 goto C3 181 addw2 r6,r0 #x:=x*2**57 scale up x by nx 182 movw $0,-4(fp) #clear nx 183 brb loop 184C3: 185 movq r2,r4 #r4 = y 186 subw2 $0x80,r4 #r4 = y/2 187 cmpd r0,r4 #x:y/2 188 blss E2 #if x < y/2 goto E2 189 bgtr C4 #if x > y/2 goto C4 190 cvtdl r8,r8 #ifix(float(n)) 191 blbc r8,E2 #if the last bit is zero, goto E2 192C4: 193 subd2 r2,r0 #x-y 194E2: 195 xorw2 -12(fp),r0 #x^sign (exclusive or) 196 movw -8(fp),r6 #r6=nf 197 bicw3 $0x807f,r0,r8 #r8=exponent of x 198 bicw2 $0x7f80,r0 #clear the exponent of x 199 subw2 r6,r8 #r8=xexp-nf 200 bgtr C5 #if xexp-nf is positive goto C5 201 movw $0,r8 #clear r8 202 movq $0,r0 #x underflow to zero 203C5: 204 bisw2 r8,r0 #put r8 into x's exponent field 205 ret 206Rop: #Reserved operand 207 pushl $EDOM 208 calls $1,_infnan #generate reserved op fault 209 ret 210Ret: 211 movq $0x8000,r0 #propagate reserved op 212 ret 213