1/* $OpenBSD: ldexp.S,v 1.9 2019/01/05 12:16:59 visa Exp $ */ 2/*- 3 * Copyright (c) 1991, 1993 4 * The Regents of the University of California. All rights reserved. 5 * 6 * This code is derived from software contributed to Berkeley by 7 * Ralph Campbell. 8 * 9 * Redistribution and use in source and binary forms, with or without 10 * modification, are permitted provided that the following conditions 11 * are met: 12 * 1. Redistributions of source code must retain the above copyright 13 * notice, this list of conditions and the following disclaimer. 14 * 2. Redistributions in binary form must reproduce the above copyright 15 * notice, this list of conditions and the following disclaimer in the 16 * documentation and/or other materials provided with the distribution. 17 * 3. Neither the name of the University nor the names of its contributors 18 * may be used to endorse or promote products derived from this software 19 * without specific prior written permission. 20 * 21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 31 * SUCH DAMAGE. 32 */ 33 34#include "SYS.h" 35 36#define DEXP_INF 0x7ff 37#define DEXP_BIAS 1023 38#define DEXP_MIN -1022 39#define DEXP_MAX 1023 40#define DFRAC_BITS 52 41#define DIMPL_ONE 0x00100000 42#define DLEAD_ZEROS 63 - 52 43#define STICKYBIT 1 44#define GUARDBIT 0x80000000 45#define DSIGNAL_NAN 0x00040000 46#define DQUIET_NAN0 0x0007ffff 47#define DQUIET_NAN1 0xffffffff 48 49/* 50 * double ldexp(x, N) 51 * double x; int N; 52 * 53 * Return x * (2**N), for integer values N. 54 */ 55LEAF(ldexp, 0) 56 .set reorder 57 dmfc1 t3, $f12 # get x 58 dsll t1, t3, 1 # get x exponent 59 dsrl t1, t1, 64 - 11 60 beq t1, DEXP_INF, 9f # is it a NAN or infinity? 61 beq t1, zero, 1f # zero or denormalized number? 62 daddu t1, a1 # scale exponent 63 dsll v0, a1, 52 # position N for addition 64 bge t1, DEXP_INF, 8f # overflow? 65 daddu v0, t3, v0 # multiply by (2**N) 66 ble t1, zero, 4f # underflow? 67 dmtc1 v0, $f0 # save result 68 j ra 691: 70 dsll t2, t3, 64 - 52 # get x fraction 71 dsrl t2, t2, 64 - 52 72 dsrl t0, t3, 63 # get x sign 73 beq t2, zero, 9f # result is zero 74/* 75 * Find out how many leading zero bits are in t2 and put in t9. 76 */ 77 move v0, t2 78 move t9, zero 79 dsrl ta0, v0, 32 80 bne ta0, zero, 1f 81 daddu t9, 32 82 dsll v0, 32 831: 84 dsrl ta0, v0, 16 85 bne ta0, zero, 1f 86 daddu t9, 16 87 dsll v0, 16 881: 89 dsrl ta0, v0, 24 90 bne ta0, zero, 1f 91 daddu t9, 8 92 dsll v0, 8 931: 94 dsrl ta0, v0, 28 95 bne ta0, zero, 1f 96 daddu t9, 4 97 dsll v0, 4 981: 99 dsrl ta0, v0, 30 100 bne ta0, zero, 1f 101 daddu t9, 2 102 dsll v0, 2 1031: 104 dsrl ta0, v0, 31 105 bne ta0, zero, 1f 106 daddu t9, 1 107/* 108 * Now shift t2 the correct number of bits. 109 */ 1101: 111 dsubu t9, t9, DLEAD_ZEROS # dont count normal leading zeros 112 li t1, DEXP_MIN + DEXP_BIAS 113 subu t1, t1, t9 # adjust exponent 114 addu t1, t1, a2 # scale exponent 115 dsll t2, t9 116 117 bge t1, DEXP_INF, 8f # overflow? 118 ble t1, zero, 4f # underflow? 119 dsll t2, t2, 64 - 52 # clear implied one bit 120 dsrl t2, t2, 64 - 52 121 122 dsll t1, t1, 63 - 11 # reposition exponent 123 dsll t0, t0, 63 # reposition sign 124 or t0, t0, t1 # put result back together 125 or t0, t0, t2 126 dmtc1 t0, $f0 # save result 127 j ra 1284: 129 dli v0, 0x8000000000000000 130 ble t1, -52, 7f # is result too small for denorm? 131 dsll t2, t3, 63 - 52 # clear exponent, extract fraction 132 or t2, t2, v0 # set implied one bit 133 dsrl t2, t2, 63 - 52 # shift fraction back to normal position 134 subu t1, t1, 1 135 dsrl t8, t2, t1 # save bits shifted out 136 negu t1 137 dsrl t2, t2, t1 138 bge t8, zero, 1f # does result need to be rounded? 139 daddu t2, t2, 1 # round result 140 dsll t8, t8, 1 141 bne t8, zero, 1f # round result to nearest 142 ori t2, 1 143 xori t2, 1 1441: 145 dmtc1 t2, $f0 # save denormalized result (LSW) 146 bge t3, zero, 1f # should result be negative? 147 neg.d $f0, $f0 # negate result 1481: 149 j ra 1507: 151 dmtc1 zero, $f0 # result is zero 152 bge t3, zero, 1f # is result positive? 153 neg.d $f0, $f0 # negate result 1541: 155 j ra 1568: 157 dli t1, 0x7ff0000000000000 # result is infinity (MSW) 158 dmtc1 t1, $f0 159 bge t3, zero, 1f # should result be negative infinity? 160 neg.d $f0, $f0 # result is negative infinity 1611: 162 add.d $f0, $f0, $f0 # cause overflow faults if enabled 163 j ra 1649: 165 mov.d $f0, $f12 # yes, result is just x 166 j ra 167END_STRONG(ldexp) 168