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

   1 /*      $NetBSD: modf.S,v 1.4 2003/08/07 16:42:27 agc 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.4 2003/08/07 16:42:27 agc 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 + BIAS - 4]      saved %fsr
  75  *      4@[%fp + BIAS - 8]      new %fsr with rounding set to `towards 0'
  76  *      8@[%fp + BIAS - 16]     space for moving between %i and %f registers
  77  * Register usage:
  78  *      %f0:f1          double val;
  79  *      %l0             scratch
  80  *      %l1             sign bit (0x80000000)
  81  *      %i1             double *iptr;
  82  *      %f2:f3          `magic number' 2^52, in fpu registers
  83  *      %f4:f5          double v, in fpu registers
  84  *      %f6:f7          double temp.
  85  */
  86
  87         .align  8
  88 Lmagic:
  89         .word   0x43300000      ! sign = 0, exponent = 52 + 1023, mantissa = 0
  90         .word   0               ! (i.e., .double 0r4503599627370496e+00)
  91
  92 L0:
  93         .word   0               ! 0.0
  94         .word   0
  95
  96 ENTRY(modf)
  97         save    %sp, -CC64FSZ-16, %sp
  98
  99         /*
 100          * First, compute v = abs(val)
 101          */
 102         fabsd   %f0, %f4                ! %f4:f5 = v
 103         fcmped  %fcc1, %f0, %f4         ! %fcc1 = (val == abs(val))
 104 #ifdef PIC
 105         PICCY_SET(Lmagic, %l0, %o7)
 106         ldd     [%l0], %f2
 107 #else
 108         sethi   %hi(Lmagic), %l0
 109         ldd     [%l0 + %lo(Lmagic)], %f2
 110 #endif
 111
 112         /*
 113          * Is %f4:f5 >= %f2:f3 ?  If so, it is all integer bits.
 114          * It is probably less, though.
 115          */
 116         fcmped  %f4, %f2
 117         fbuge   Lbig                    ! if >= (or unordered), go out
 118         nop
 119
 120         /*
 121          * v < 2^52, so add 2^52, then subtract 2^52, but do it all
 122          * with rounding set towards zero.  We leave any enabled
 123          * traps enabled, but change the rounding mode.  This might
 124          * not be so good.  Oh well....
 125          */
 126         st      %fsr, [%fp + BIAS - 4]  ! %l5 = current FSR mode
 127         set     FSR_RD, %l3             ! %l3 = rounding direction mask
 128         ld      [%fp + BIAS - 4], %l5
 129         set     FSR_RD_RZ << FSR_RD_SHIFT, %l4
 130         andn    %l5, %l3, %l6
 131         or      %l6, %l4, %l6           ! round towards zero, please
 132         and     %l5, %l3, %l5           ! save original rounding mode
 133         st      %l6, [%fp + BIAS - 8]
 134         ld      [%fp + BIAS - 8], %fsr
 135
 136         faddd   %f4, %f2, %f4           ! %f4:f5 += 2^52
 137         fsubd   %f4, %f2, %f4           ! %f4:f5 -= 2^52
 138
 139         /*
 140          * Restore %fsr, but leave exceptions accrued.
 141          */
 142         st      %fsr, [%fp + BIAS - 4]
 143         ld      [%fp + BIAS - 4], %l6
 144         andn    %l6, %l3, %l6           ! %l6 = %fsr & ~FSR_RD;
 145         or      %l5, %l6, %l5           ! %l5 |= %l6;
 146         st      %l5, [%fp + BIAS - 4]
 147         ld      [%fp + BIAS - 4], %fsr  ! restore %fsr, leaving accrued stuff
 148
 149         /*
 150          * Now insert the original sign in %f4:f5.
 151          * %fcc1 should still have the reults of (val == abs(val))
 152          * from above, so we use a conditional move on %fcc1 to:
 153          *
 154          *      %f4 = (val == abs(val)) ? %f4 : -%f4
 155          *
 156          */
 157         fnegd   %f4, %f6
 158         fmovdnz %fcc1, %f6, %f4
 159 1:
 160
 161         /*
 162          * The value in %f4:f5 is now the integer portion of the original
 163          * argument.  We need to store this in *ival (%i1), subtract it
 164          * from the original value argument (%d0), and return the result.
 165          */
 166         std     %f4, [%i1]              ! *ival = %f4:f5;
 167         fsubd   %f0, %f4, %f0           ! %f0:f1 -= %f4:f5;
 168         ret
 169         restore
 170
 171 Lbig:
 172         /*
 173          * We get here if the original comparison of %f4:f5 (v) to
 174          * %f2:f3 (2^52) came out `greater or unordered'.  In this
 175          * case the integer part is the original value, and the
 176          * fractional part is 0.
 177          */
 178 #ifdef PIC
 179         PICCY_SET(L0, %l0, %o7)
 180         std     %f0, [%i1]              ! *ival = val;
 181         ldd     [%l0], %f0              ! return 0.0;
 182 #else
 183         sethi   %hi(L0), %l0
 184         std     %f0, [%i1]              ! *ival = val;
 185         ldd     [%l0 + %lo(L0)], %f0    ! return 0.0;
 186 #endif
 187         ret
 188         restore
 189