lib/libc/arch/sparc/gen/modf.S

   1 /*      $NetBSD: modf.S,v 1.7 2013/09/12 15:36:16 joerg Exp $   */
   2
   3 /*
   4  * Copyright (c) 1992, 1993
   5  *      The Regents of the University of California.  All rights reserved.
   6  *
   7  * This software was developed by the Computer Systems Engineering group
   8  * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
   9  * contributed to Berkeley.
  10  *
  11  * Redistribution and use in source and binary forms, with or without
  12  * modification, are permitted provided that the following conditions
  13  * are met:
  14  * 1. Redistributions of source code must retain the above copyright
  15  *    notice, this list of conditions and the following disclaimer.
  16  * 2. Redistributions in binary form must reproduce the above copyright
  17  *    notice, this list of conditions and the following disclaimer in the
  18  *    documentation and/or other materials provided with the distribution.
  19  * 3. Neither the name of the University nor the names of its contributors
  20  *    may be used to endorse or promote products derived from this software
  21  *    without specific prior written permission.
  22  *
  23  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  24  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  25  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  26  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  27  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  28  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  29  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  30  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  31  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  32  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  33  * SUCH DAMAGE.
  34  *
  35  * from: Header: modf.s,v 1.3 92/06/20 00:00:54 torek Exp
  36  */
  37
  38 #include <machine/asm.h>
  39 #if defined(LIBC_SCCS) && !defined(lint)
  40 #if 0
  41         .asciz "@(#)modf.s      8.1 (Berkeley) 6/4/93"
  42 #else
  43         RCSID("$NetBSD: modf.S,v 1.7 2013/09/12 15:36:16 joerg Exp $")
  44 #endif
  45 #endif /* LIBC_SCCS and not lint */
  46
  47 #include <machine/fsr.h>
  48
  49 /*
  50  * double modf(double val, double *iptr)
  51  *
  52  * Returns the fractional part of `val', storing the integer part of
  53  * `val' in *iptr.  Both *iptr and the return value have the same sign
  54  * as `val'.
  55  *
  56  * Method:
  57  *
  58  * We use the fpu's normalization hardware to compute the integer portion
  59  * of the double precision argument.  Sun IEEE double precision numbers
  60  * have 52 bits of mantissa, 11 bits of exponent, and one bit of sign,
  61  * with the sign occupying bit 31 of word 0, and the exponent bits 30:20
  62  * of word 0.  Thus, values >= 2^52 are by definition integers.
  63  *
  64  * If we take a value that is in the range [+0..2^52) and add 2^52, all
  65  * of the fractional bits fall out and all of the integer bits are summed
  66  * with 2^52.  If we then subtract 2^52, we get those integer bits back.
  67  * This must be done with rounding set to `towards 0' or `towards -inf'.
  68  * `Toward -inf' fails when the value is 0 (we get -0 back)....
  69  *
  70  * Note that this method will work anywhere, but is machine dependent in
  71  * various aspects.
  72  *
  73  * Stack usage:
  74  *      4@[%fp - 4]     saved %fsr
  75  *      4@[%fp - 8]     new %fsr with rounding set to `towards 0'
  76  *      8@[%fp - 16]    space for moving between %i and %f registers
  77  * Register usage:
  78  *      %i0%i1          double val;
  79  *      %l0             scratch
  80  *      %l1             sign bit (0x80000000)
  81  *      %i2             double *iptr;
  82  *      %f2:f3          `magic number' 2^52, in fpu registers
  83  *      %f4:f5          double v, in fpu registers
  84  */
  85
  86         .align  8
  87 Lmagic:
  88         .word   0x43300000      ! sign = 0, exponent = 52 + 1023, mantissa = 0
  89         .word   0               ! (i.e., .double 0r4503599627370496e+00)
  90
  91 L0:
  92         .word   0               ! 0.0
  93         .word   0
  94
  95 ENTRY(modf)
  96         save    %sp, -64-16, %sp
  97
  98         /*
  99          * First, compute v = abs(val) by clearing sign bit,
 100          * and then set up the fpu registers.  This would be
 101          * much easier if we could do alu operations on fpu registers!
 102          */
 103         sethi   %hi(0x80000000), %l1    ! sign bit
 104         andn    %i0, %l1, %l0
 105         st      %l0, [%fp - 16]
 106 #ifdef __PIC__
 107         PICCY_SET(Lmagic, %l0, %o7)
 108         ldd     [%l0], %f2
 109 #else
 110         sethi   %hi(Lmagic), %l0
 111         ldd     [%l0 + %lo(Lmagic)], %f2
 112 #endif
 113         st      %i1, [%fp - 12]
 114         ldd     [%fp - 16], %f4         ! %f4:f5 = v
 115
 116         /*
 117          * Is %f4:f5 >= %f2:f3 ?  If so, it is all integer bits.
 118          * It is probably less, though.
 119          */
 120         fcmped  %f4, %f2
 121         nop                             ! fpop2 delay
 122         fbuge   Lbig                    ! if >= (or unordered), go out
 123         nop
 124
 125         /*
 126          * v < 2^52, so add 2^52, then subtract 2^52, but do it all
 127          * with rounding set towards zero.  We leave any enabled
 128          * traps enabled, but change the rounding mode.  This might
 129          * not be so good.  Oh well....
 130          */
 131         st      %fsr, [%fp - 4]         ! %l5 = current FSR mode
 132         set     FSR_RD, %l3             ! %l3 = rounding direction mask
 133         ld      [%fp - 4], %l5
 134         set     FSR_RD_RZ << FSR_RD_SHIFT, %l4
 135         andn    %l5, %l3, %l6
 136         or      %l6, %l4, %l6           ! round towards zero, please
 137         and     %l5, %l3, %l5           ! save original rounding mode
 138         st      %l6, [%fp - 8]
 139         ld      [%fp - 8], %fsr
 140
 141         faddd   %f4, %f2, %f4           ! %f4:f5 += 2^52
 142         fsubd   %f4, %f2, %f4           ! %f4:f5 -= 2^52
 143
 144         /*
 145          * Restore %fsr, but leave exceptions accrued.
 146          */
 147         st      %fsr, [%fp - 4]
 148         ld      [%fp - 4], %l6
 149         andn    %l6, %l3, %l6           ! %l6 = %fsr & ~FSR_RD;
 150         or      %l5, %l6, %l5           ! %l5 |= %l6;
 151         st      %l5, [%fp - 4]
 152         ld      [%fp - 4], %fsr         ! restore %fsr, leaving accrued stuff
 153
 154         /*
 155          * Now insert the original sign in %f4:f5.
 156          * This is a lot of work, so it is conditional here.
 157          */
 158         btst    %l1, %i0
 159         be      1f
 160         nop
 161         st      %f4, [%fp - 16]
 162         ld      [%fp - 16], %g1
 163         or      %l1, %g1, %g1
 164         st      %g1, [%fp - 16]
 165         ld      [%fp - 16], %f4
 166 1:
 167
 168         /*
 169          * The value in %f4:f5 is now the integer portion of the original
 170          * argument.  We need to store this in *ival (%i2), subtract it
 171          * from the original value argument (%i0:i1), and return the result.
 172          */
 173         std     %f4, [%i2]              ! *ival = %f4:f5;
 174         std     %i0, [%fp - 16]
 175         ldd     [%fp - 16], %f0         ! %f0:f1 = val;
 176         fsubd   %f0, %f4, %f0           ! %f0:f1 -= %f4:f5;
 177         ret
 178         restore
 179
 180 Lbig:
 181         /*
 182          * We get here if the original comparison of %f4:f5 (v) to
 183          * %f2:f3 (2^52) came out `greater or unordered'.  In this
 184          * case the integer part is the original value, and the
 185          * fractional part is 0.
 186          */
 187 #ifdef __PIC__
 188         PICCY_SET(L0, %l0, %o7)
 189         std     %f0, [%i2]              ! *ival = val;
 190         ldd     [%l0], %f0              ! return 0.0;
 191 #else
 192         sethi   %hi(L0), %l0
 193         std     %f0, [%i2]              ! *ival = val;
 194         ldd     [%l0 + %lo(L0)], %f0    ! return 0.0;
 195 #endif
 196         ret
 197         restore