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