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