1/* 2 * Copyright (c) 1992, 1993 3 * The Regents of the University of California. All rights reserved. 4 * 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. 8 * 9 * %sccs.include.redist.c% 10 * 11 * from: $Header: modf.s,v 1.3 92/06/20 00:00:54 torek Exp $ 12 */ 13 14#if defined(LIBC_SCCS) && !defined(lint) 15 .asciz "@(#)modf.s 8.1 (Berkeley) 06/04/93" 16#endif /* LIBC_SCCS and not lint */ 17 18#include "DEFS.h" 19#include <machine/fsr.h> 20 21/* 22 * double modf(double val, double *iptr) 23 * 24 * Returns the fractional part of `val', storing the integer part of 25 * `val' in *iptr. Both *iptr and the return value have the same sign 26 * as `val'. 27 * 28 * Method: 29 * 30 * We use the fpu's normalization hardware to compute the integer portion 31 * of the double precision argument. Sun IEEE double precision numbers 32 * have 52 bits of mantissa, 11 bits of exponent, and one bit of sign, 33 * with the sign occupying bit 31 of word 0, and the exponent bits 30:20 34 * of word 0. Thus, values >= 2^52 are by definition integers. 35 * 36 * If we take a value that is in the range [+0..2^52) and add 2^52, all 37 * of the fractional bits fall out and all of the integer bits are summed 38 * with 2^52. If we then subtract 2^52, we get those integer bits back. 39 * This must be done with rounding set to `towards 0' or `towards -inf'. 40 * `Toward -inf' fails when the value is 0 (we get -0 back).... 41 * 42 * Note that this method will work anywhere, but is machine dependent in 43 * various aspects. 44 * 45 * Stack usage: 46 * 4@[%fp - 4] saved %fsr 47 * 4@[%fp - 8] new %fsr with rounding set to `towards 0' 48 * 8@[%fp - 16] space for moving between %i and %f registers 49 * Register usage: 50 * %i0%i1 double val; 51 * %l0 scratch 52 * %l1 sign bit (0x80000000) 53 * %i2 double *iptr; 54 * %f2:f3 `magic number' 2^52, in fpu registers 55 * %f4:f5 double v, in fpu registers 56 */ 57 58 .align 8 59Lmagic: 60 .word 0x43300000 ! sign = 0, exponent = 52 + 1023, mantissa = 0 61 .word 0 ! (i.e., .double 0r4503599627370496e+00) 62 63L0: 64 .word 0 ! 0.0 65 .word 0 66 67ENTRY(modf) 68 save %sp, -64-16, %sp 69 70 /* 71 * First, compute v = abs(val) by clearing sign bit, 72 * and then set up the fpu registers. This would be 73 * much easier if we could do alu operations on fpu registers! 74 */ 75 sethi 0x80000000, %l1 ! sign bit 76 andn %i0, %l1, %l0 77 st %l0, [%fp - 16] 78 sethi %hi(Lmagic), %l0 79 ldd [%l0 + %lo(Lmagic)], %f2 80 st %i1, [%fp - 12] 81 ldd [%fp - 16], %f4 ! %f4:f5 = v 82 83 /* 84 * Is %f4:f5 >= %f2:f3 ? If so, it is all integer bits. 85 * It is probably less, though. 86 */ 87 fcmped %f4, %f2 88 nop ! fpop2 delay 89 fbuge Lbig ! if >= (or unordered), go out 90 nop 91 92 /* 93 * v < 2^52, so add 2^52, then subtract 2^52, but do it all 94 * with rounding set towards zero. We leave any enabled 95 * traps enabled, but change the rounding mode. This might 96 * not be so good. Oh well.... 97 */ 98 st %fsr, [%fp - 4] ! %l5 = current FSR mode 99 set FSR_RD, %l3 ! %l3 = rounding direction mask 100 ld [%fp - 4], %l5 101 set FSR_RD_RZ << FSR_RD_SHIFT, %l4 102 andn %l5, %l3, %l6 103 or %l6, %l4, %l6 ! round towards zero, please 104 and %l5, %l3, %l5 ! save original rounding mode 105 st %l6, [%fp - 8] 106 ld [%fp - 8], %fsr 107 108 faddd %f4, %f2, %f4 ! %f4:f5 += 2^52 109 fsubd %f4, %f2, %f4 ! %f4:f5 -= 2^52 110 111 /* 112 * Restore %fsr, but leave exceptions accrued. 113 */ 114 st %fsr, [%fp - 4] 115 ld [%fp - 4], %l6 116 andn %l6, %l3, %l6 ! %l6 = %fsr & ~FSR_RD; 117 or %l5, %l6, %l5 ! %l5 |= %l6; 118 st %l5, [%fp - 4] 119 ld [%fp - 4], %fsr ! restore %fsr, leaving accrued stuff 120 121 /* 122 * Now insert the original sign in %f4:f5. 123 * This is a lot of work, so it is conditional here. 124 */ 125 btst %l1, %i0 126 be 1f 127 nop 128 st %f4, [%fp - 16] 129 ld [%fp - 16], %g1 130 or %l1, %g1, %g1 131 st %g1, [%fp - 16] 132 ld [%fp - 16], %f4 1331: 134 135 /* 136 * The value in %f4:f5 is now the integer portion of the original 137 * argument. We need to store this in *ival (%i2), subtract it 138 * from the original value argument (%i0:i1), and return the result. 139 */ 140 std %f4, [%i2] ! *ival = %f4:f5; 141 std %i0, [%fp - 16] 142 ldd [%fp - 16], %f0 ! %f0:f1 = val; 143 fsubd %f0, %f4, %f0 ! %f0:f1 -= %f4:f5; 144 ret 145 restore 146 147Lbig: 148 /* 149 * We get here if the original comparison of %f4:f5 (v) to 150 * %f2:f3 (2^52) came out `greater or unordered'. In this 151 * case the integer part is the original value, and the 152 * fractional part is 0. 153 */ 154 sethi %hi(L0), %l0 155 std %f0, [%i2] ! *ival = val; 156 ldd [%l0 + %lo(L0)], %f0 ! return 0.0; 157 ret 158 restore 159