lib/libcr/mips/gen/ldexp.S

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