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