2 * Copyright (c) 1992, 1993
3 * The Regents of the University of California. All rights reserved.
5 * This software was developed by the Computer Systems Engineering group
6 * at Lawrence Berkeley Laboratory under DARPA contract BG 91-66 and
7 * contributed to Berkeley.
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
12 * 1. Redistributions of source code must retain the above copyright
13 * notice, this list of conditions and the following disclaimer.
14 * 2. Redistributions in binary form must reproduce the above copyright
15 * notice, this list of conditions and the following disclaimer in the
16 * documentation and/or other materials provided with the distribution.
17 * 4. Neither the name of the University nor the names of its contributors
18 * may be used to endorse or promote products derived from this software
19 * without specific prior written permission.
21 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
22 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24 * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
25 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
27 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
30 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
33 * from: Header: modf.s,v 1.3 92/06/20 00:00:54 torek Exp
36 #if defined(LIBC_SCCS) && !defined(lint)
37 .asciz "@(#)modf.s 8.1 (Berkeley) 6/4/93"
39 .asciz "$NetBSD: modf.S,v 1.2 2000/07/23 07:12:22 eeh Exp $"
41 #endif /* LIBC_SCCS and not lint */
42 #include <machine/asm.h>
43 __FBSDID("$FreeBSD$");
48 * double modf(double val, double *iptr)
50 * Returns the fractional part of `val', storing the integer part of
51 * `val' in *iptr. Both *iptr and the return value have the same sign
56 * We use the fpu's normalization hardware to compute the integer portion
57 * of the double precision argument. Sun IEEE double precision numbers
58 * have 52 bits of mantissa, 11 bits of exponent, and one bit of sign,
59 * with the sign occupying bit 31 of word 0, and the exponent bits 30:20
60 * of word 0. Thus, values >= 2^52 are by definition integers.
62 * If we take a value that is in the range [+0..2^52) and add 2^52, all
63 * of the fractional bits fall out and all of the integer bits are summed
64 * with 2^52. If we then subtract 2^52, we get those integer bits back.
65 * This must be done with rounding set to `towards 0' or `towards -inf'.
66 * `Toward -inf' fails when the value is 0 (we get -0 back)....
68 * Note that this method will work anywhere, but is machine dependent in
72 * 4@[%fp + SPOFF - 4] saved %fsr
73 * 4@[%fp + SPOFF - 8] new %fsr with rounding set to `towards 0'
74 * 8@[%fp + SPOFF - 16] space for moving between %i and %f registers
78 * %l1 sign bit (0x80000000)
80 * %f2:f3 `magic number' 2^52, in fpu registers
81 * %f4:f5 double v, in fpu registers
87 .word 0x43300000 ! sign = 0, exponent = 52 + 1023, mantissa = 0
88 .word 0 ! (i.e., .double 0r4503599627370496e+00)
95 save %sp, -CCFSZ - 16, %sp
96 PIC_PROLOGUE(%l6, %l7)
99 * First, compute v = abs(val)
101 fabsd %f0, %f4 ! %f4:f5 = v
102 fcmped %fcc1, %f0, %f4 ! %fcc1 = (val == abs(val))
103 SET(.Lmagic, %l7, %l0)
107 * Is %f4:f5 >= %f2:f3 ? If so, it is all integer bits.
108 * It is probably less, though.
111 fbuge .Lbig ! if >= (or unordered), go out
115 * v < 2^52, so add 2^52, then subtract 2^52, but do it all
116 * with rounding set towards zero. We leave any enabled
117 * traps enabled, but change the rounding mode. This might
118 * not be so good. Oh well....
120 st %fsr, [%fp + SPOFF - 4] ! %l5 = current FSR mode
121 set FSR_RD_MASK, %l3 ! %l3 = rounding direction mask
122 ld [%fp + SPOFF - 4], %l5
125 or %l6, %l4, %l6 ! round towards zero, please
126 and %l5, %l3, %l5 ! save original rounding mode
127 st %l6, [%fp + SPOFF - 8]
128 ld [%fp + SPOFF - 8], %fsr
130 faddd %f4, %f2, %f4 ! %f4:f5 += 2^52
131 fsubd %f4, %f2, %f4 ! %f4:f5 -= 2^52
134 * Restore %fsr, but leave exceptions accrued.
136 st %fsr, [%fp + SPOFF - 4]
137 ld [%fp + SPOFF - 4], %l6
138 andn %l6, %l3, %l6 ! %l6 = %fsr & ~FSR_RD_MASK;
139 or %l5, %l6, %l5 ! %l5 |= %l6;
140 st %l5, [%fp + SPOFF - 4]
141 ld [%fp + SPOFF - 4], %fsr ! restore %fsr, leaving accrued stuff
144 * Now insert the original sign in %f4:f5.
145 * %fcc1 should still have the results of (val == abs(val))
146 * from above, so we use a conditional move on %fcc1 to:
148 * %f4 = (val == abs(val)) ? %f4 : -%f4
152 fmovdnz %fcc1, %f6, %f4
156 * The value in %f4:f5 is now the integer portion of the original
157 * argument. We need to store this in *ival (%i1), subtract it
158 * from the original value argument (%d0), and return the result.
160 std %f4, [%i1] ! *ival = %f4:f5;
161 fsubd %f0, %f4, %f0 ! %f0:f1 -= %f4:f5;
167 * We get here if the original comparison of %f4:f5 (v) to
168 * %f2:f3 (2^52) came out `greater or unordered'. In this
169 * case the integer part is the original value, and the
170 * fractional part is 0.
173 std %f0, [%i1] ! *ival = val;
174 ldd [%l0], %f0 ! return 0.0;