Home | History | Annotate | Line # | Download | only in gen
      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