modf.S revision 1.7 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
198