lib/libc/mips/gen/ldexp.S

   1 /*      $NetBSD: ldexp.S,v 1.8 2003/08/07 16:42:15 agc Exp $    */
   2
   3 /*-
   4  * Copyright (c) 1991, 1993
   5  *      The Regents of the University of California.  All rights reserved.
   6  *
   7  * This code is derived from software contributed to Berkeley by
   8  * Ralph Campbell.
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  * 3. Neither the name of the University nor the names of its contributors
  19  *    may be used to endorse or promote products derived from this software
  20  *    without specific prior written permission.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  23  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  24  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  25  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  26  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  27  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  28  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  29  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  30  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  31  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  32  * SUCH DAMAGE.
  33  */
  34
  35 #include <machine/asm.h>
  36 __FBSDID("$FreeBSD$");
  37
  38 #if defined(LIBC_SCCS) && !defined(lint)
  39         ASMSTR("from: @(#)ldexp.s       8.1 (Berkeley) 6/4/93")
  40         ASMSTR("$NetBSD: ldexp.S,v 1.8 2003/08/07 16:42:15 agc Exp $")
  41 #endif /* LIBC_SCCS and not lint */
  42
  43 #ifdef __ABICALLS__
  44         .abicalls
  45 #endif
  46
  47 #define DEXP_INF        0x7ff
  48 #define DEXP_BIAS       1023
  49 #define DEXP_MIN        -1022
  50 #define DEXP_MAX        1023
  51 #define DFRAC_BITS      52
  52 #define DIMPL_ONE       0x00100000
  53 #define DLEAD_ZEROS     31 - 20
  54 #define STICKYBIT       1
  55 #define GUARDBIT        0x80000000
  56 #define DSIGNAL_NAN     0x00040000
  57 #define DQUIET_NAN0     0x0007ffff
  58 #define DQUIET_NAN1     0xffffffff
  59
  60 /*
  61  * double ldexp(x, N)
  62  *      double x; int N;
  63  *
  64  * Return x * (2**N), for integer values N.
  65  */
  66 LEAF(ldexp)
  67         mfc1    v1, $f13                # get MSW of x
  68         mfc1    t3, $f12                # get LSW of x
  69         sll     t1, v1, 1               # get x exponent
  70         srl     t1, t1, 32 - 11
  71         beq     t1, DEXP_INF, 9f        # is it a NAN or infinity?
  72         beq     t1, zero, 1f            # zero or denormalized number?
  73         addu    t1, t1, a2              # scale exponent
  74         sll     v0, a2, 20              # position N for addition
  75         bge     t1, DEXP_INF, 8f        # overflow?
  76         addu    v0, v0, v1              # multiply by (2**N)
  77         ble     t1, zero, 4f            # underflow?
  78         mtc1    v0, $f1                 # save MSW of result
  79         mtc1    t3, $f0                 # save LSW of result
  80         j       ra
  81 1:
  82         sll     t2, v1, 32 - 20         # get x fraction
  83         srl     t2, t2, 32 - 20
  84         srl     t0, v1, 31              # get x sign
  85         bne     t2, zero, 1f
  86         beq     t3, zero, 9f            # result is zero
  87 1:
  88 /*
  89  * Find out how many leading zero bits are in t2,t3 and put in t9.
  90  */
  91         move    v0, t2
  92         move    t9, zero
  93         bne     t2, zero, 1f
  94         move    v0, t3
  95         addu    t9, 32
  96 1:
  97         srl     ta0, v0, 16
  98         bne     ta0, zero, 1f
  99         addu    t9, 16
 100         sll     v0, 16
 101 1:
 102         srl     ta0, v0, 24
 103         bne     ta0, zero, 1f
 104         addu    t9, 8
 105         sll     v0, 8
 106 1:
 107         srl     ta0, v0, 28
 108         bne     ta0, zero, 1f
 109         addu    t9, 4
 110         sll     v0, 4
 111 1:
 112         srl     ta0, v0, 30
 113         bne     ta0, zero, 1f
 114         addu    t9, 2
 115         sll     v0, 2
 116 1:
 117         srl     ta0, v0, 31
 118         bne     ta0, zero, 1f
 119         addu    t9, 1
 120 /*
 121  * Now shift t2,t3 the correct number of bits.
 122  */
 123 1:
 124         subu    t9, t9, DLEAD_ZEROS     # dont count normal leading zeros
 125         li      t1, DEXP_MIN + DEXP_BIAS
 126         subu    t1, t1, t9              # adjust exponent
 127         addu    t1, t1, a2              # scale exponent
 128         li      v0, 32
 129         blt     t9, v0, 1f
 130         subu    t9, t9, v0              # shift fraction left >= 32 bits
 131         sll     t2, t3, t9
 132         move    t3, zero
 133         b       2f
 134 1:
 135         subu    v0, v0, t9              # shift fraction left < 32 bits
 136         sll     t2, t2, t9
 137         srl     ta0, t3, v0
 138         or      t2, t2, ta0
 139         sll     t3, t3, t9
 140 2:
 141         bge     t1, DEXP_INF, 8f        # overflow?
 142         ble     t1, zero, 4f            # underflow?
 143         sll     t2, t2, 32 - 20         # clear implied one bit
 144         srl     t2, t2, 32 - 20
 145 3:
 146         sll     t1, t1, 31 - 11         # reposition exponent
 147         sll     t0, t0, 31              # reposition sign
 148         or      t0, t0, t1              # put result back together
 149         or      t0, t0, t2
 150         mtc1    t0, $f1                 # save MSW of result
 151         mtc1    t3, $f0                 # save LSW of result
 152         j       ra
 153 4:
 154         li      v0, 0x80000000
 155         ble     t1, -52, 7f             # is result too small for denorm?
 156         sll     t2, v1, 31 - 20         # clear exponent, extract fraction
 157         or      t2, t2, v0              # set implied one bit
 158         blt     t1, -30, 2f             # will all bits in t3 be shifted out?
 159         srl     t2, t2, 31 - 20         # shift fraction back to normal position
 160         subu    t1, t1, 1
 161         sll     ta0, t2, t1             # shift right t2,t3 based on exponent
 162         srl     t8, t3, t1              # save bits shifted out
 163         negu    t1
 164         srl     t3, t3, t1
 165         or      t3, t3, ta0
 166         srl     t2, t2, t1
 167         bge     t8, zero, 1f            # does result need to be rounded?
 168         addu    t3, t3, 1               # round result
 169         sltu    ta0, t3, 1
 170         sll     t8, t8, 1
 171         addu    t2, t2, ta0
 172         bne     t8, zero, 1f            # round result to nearest
 173         and     t3, t3, ~1
 174 1:
 175         mtc1    t3, $f0                 # save denormalized result (LSW)
 176         mtc1    t2, $f1                 # save denormalized result (MSW)
 177         bge     v1, zero, 1f            # should result be negative?
 178         neg.d   $f0, $f0                # negate result
 179 1:
 180         j       ra
 181 2:
 182         mtc1    zero, $f1               # exponent and upper fraction
 183         addu    t1, t1, 20              # compute amount to shift right by
 184         sll     t8, t2, t1              # save bits shifted out
 185         negu    t1
 186         srl     t3, t2, t1
 187         bge     t8, zero, 1f            # does result need to be rounded?
 188         addu    t3, t3, 1               # round result
 189         sltu    ta0, t3, 1
 190         sll     t8, t8, 1
 191         mtc1    ta0, $f1                        # exponent and upper fraction
 192         bne     t8, zero, 1f            # round result to nearest
 193         and     t3, t3, ~1
 194 1:
 195         mtc1    t3, $f0
 196         bge     v1, zero, 1f            # is result negative?
 197         neg.d   $f0, $f0                # negate result
 198 1:
 199         j       ra
 200 7:
 201         mtc1    zero, $f0               # result is zero
 202         mtc1    zero, $f1
 203         beq     t0, zero, 1f            # is result positive?
 204         neg.d   $f0, $f0                # negate result
 205 1:
 206         j       ra
 207 8:
 208         li      t1, 0x7ff00000          # result is infinity (MSW)
 209         mtc1    t1, $f1
 210         mtc1    zero, $f0               # result is infinity (LSW)
 211         bge     v1, zero, 1f            # should result be negative infinity?
 212         neg.d   $f0, $f0                # result is negative infinity
 213 1:
 214         add.d   $f0, $f0                # cause overflow faults if enabled
 215         j       ra
 216 9:
 217         mov.d   $f0, $f12               # yes, result is just x
 218         j       ra
 219 END(ldexp)