1/*- 2 * Copyright (c) 1991, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 5 * This code is derived from software contributed to Berkeley by 6 * Ralph Campbell. 7 * 8 * %sccs.include.redist.c% 9 */ 10 11#include <machine/machAsmDefs.h> 12 13#if defined(LIBC_SCCS) && !defined(lint) 14 ASMSTR("@(#)ldexp.s 8.1 (Berkeley) 06/04/93") 15#endif /* LIBC_SCCS and not lint */ 16 17#define DEXP_INF 0x7ff 18#define DEXP_BIAS 1023 19#define DEXP_MIN -1022 20#define DEXP_MAX 1023 21#define DFRAC_BITS 52 22#define DIMPL_ONE 0x00100000 23#define DLEAD_ZEROS 31 - 20 24#define STICKYBIT 1 25#define GUARDBIT 0x80000000 26#define DSIGNAL_NAN 0x00040000 27#define DQUIET_NAN0 0x0007ffff 28#define DQUIET_NAN1 0xffffffff 29 30/* 31 * double ldexp(x, N) 32 * double x; int N; 33 * 34 * Return x * (2**N), for integer values N. 35 */ 36LEAF(ldexp) 37 mfc1 v1, $f13 # get MSW of x 38 mfc1 t3, $f12 # get LSW of x 39 sll t1, v1, 1 # get x exponent 40 srl t1, t1, 32 - 11 41 beq t1, DEXP_INF, 9f # is it a NAN or infinity? 42 beq t1, zero, 1f # zero or denormalized number? 43 addu t1, t1, a2 # scale exponent 44 sll v0, a2, 20 # position N for addition 45 bge t1, DEXP_INF, 8f # overflow? 46 addu v0, v0, v1 # multiply by (2**N) 47 ble t1, zero, 4f # underflow? 48 mtc1 v0, $f1 # save MSW of result 49 mtc1 t3, $f0 # save LSW of result 50 j ra 511: 52 sll t2, v1, 32 - 20 # get x fraction 53 srl t2, t2, 32 - 20 54 srl t0, v1, 31 # get x sign 55 bne t2, zero, 1f 56 beq t3, zero, 9f # result is zero 571: 58/* 59 * Find out how many leading zero bits are in t2,t3 and put in t9. 60 */ 61 move v0, t2 62 move t9, zero 63 bne t2, zero, 1f 64 move v0, t3 65 addu t9, 32 661: 67 srl t4, v0, 16 68 bne t4, zero, 1f 69 addu t9, 16 70 sll v0, 16 711: 72 srl t4, v0, 24 73 bne t4, zero, 1f 74 addu t9, 8 75 sll v0, 8 761: 77 srl t4, v0, 28 78 bne t4, zero, 1f 79 addu t9, 4 80 sll v0, 4 811: 82 srl t4, v0, 30 83 bne t4, zero, 1f 84 addu t9, 2 85 sll v0, 2 861: 87 srl t4, v0, 31 88 bne t4, zero, 1f 89 addu t9, 1 90/* 91 * Now shift t2,t3 the correct number of bits. 92 */ 931: 94 subu t9, t9, DLEAD_ZEROS # dont count normal leading zeros 95 li t1, DEXP_MIN + DEXP_BIAS 96 subu t1, t1, t9 # adjust exponent 97 addu t1, t1, a2 # scale exponent 98 li v0, 32 99 blt t9, v0, 1f 100 subu t9, t9, v0 # shift fraction left >= 32 bits 101 sll t2, t3, t9 102 move t3, zero 103 b 2f 1041: 105 subu v0, v0, t9 # shift fraction left < 32 bits 106 sll t2, t2, t9 107 srl t4, t3, v0 108 or t2, t2, t4 109 sll t3, t3, t9 1102: 111 bge t1, DEXP_INF, 8f # overflow? 112 ble t1, zero, 4f # underflow? 113 sll t2, t2, 32 - 20 # clear implied one bit 114 srl t2, t2, 32 - 20 1153: 116 sll t1, t1, 31 - 11 # reposition exponent 117 sll t0, t0, 31 # reposition sign 118 or t0, t0, t1 # put result back together 119 or t0, t0, t2 120 mtc1 t0, $f1 # save MSW of result 121 mtc1 t3, $f0 # save LSW of result 122 j ra 1234: 124 li v0, 0x80000000 125 ble t1, -52, 7f # is result too small for denorm? 126 sll t2, v1, 31 - 20 # clear exponent, extract fraction 127 or t2, t2, v0 # set implied one bit 128 blt t1, -30, 2f # will all bits in t3 be shifted out? 129 srl t2, t2, 31 - 20 # shift fraction back to normal position 130 subu t1, t1, 1 131 sll t4, t2, t1 # shift right t2,t3 based on exponent 132 srl t8, t3, t1 # save bits shifted out 133 negu t1 134 srl t3, t3, t1 135 or t3, t3, t4 136 srl t2, t2, t1 137 bge t8, zero, 1f # does result need to be rounded? 138 addu t3, t3, 1 # round result 139 sltu t4, t3, 1 140 sll t8, t8, 1 141 addu t2, t2, t4 142 bne t8, zero, 1f # round result to nearest 143 and t3, t3, ~1 1441: 145 mtc1 t3, $f0 # save denormalized result (LSW) 146 mtc1 t2, $f1 # save denormalized result (MSW) 147 bge v1, zero, 1f # should result be negative? 148 neg.d $f0, $f0 # negate result 1491: 150 j ra 1512: 152 mtc1 zero, $f1 # exponent and upper fraction 153 addu t1, t1, 20 # compute amount to shift right by 154 sll t8, t2, t1 # save bits shifted out 155 negu t1 156 srl t3, t2, t1 157 bge t8, zero, 1f # does result need to be rounded? 158 addu t3, t3, 1 # round result 159 sltu t4, t3, 1 160 sll t8, t8, 1 161 mtc1 t4, $f1 # exponent and upper fraction 162 bne t8, zero, 1f # round result to nearest 163 and t3, t3, ~1 1641: 165 mtc1 t3, $f0 166 bge v1, zero, 1f # is result negative? 167 neg.d $f0, $f0 # negate result 1681: 169 j ra 1707: 171 mtc1 zero, $f0 # result is zero 172 mtc1 zero, $f1 173 beq t0, zero, 1f # is result positive? 174 neg.d $f0, $f0 # negate result 1751: 176 j ra 1778: 178 li t1, 0x7ff00000 # result is infinity (MSW) 179 mtc1 t1, $f1 180 mtc1 zero, $f0 # result is infinity (LSW) 181 bge v1, zero, 1f # should result be negative infinity? 182 neg.d $f0, $f0 # result is negative infinity 1831: 184 add.d $f0, $f0 # cause overflow faults if enabled 185 j ra 1869: 187 mov.d $f0, $f12 # yes, result is just x 188 j ra 189END(ldexp) 190