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