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