summaryrefslogtreecommitdiff
path: root/libm/powerpc/s_modf.c
blob: 403c54b79716d92ce8cec1c9443cfe39d5cf92ac (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
/*******************************************************************************
**      File:   rndint.c
**
**      Contains: C source code for implementations of floating-point
**                functions which round to integral value or format, as
**                defined in header <fp.h>.  In particular, this file
**                contains implementations of functions rint, nearbyint,
**                rinttol, round, roundtol, trunc, modf and modfl.  This file
**                targets PowerPC or Power platforms.
**
**      Written by: A. Sazegari, Apple AltiVec Group
**	    Created originally by Jon Okada, Apple Numerics Group
**
**      Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
**
**      Change History (most recent first):
**
**      13 Jul 01  ram  replaced --setflm calls with inline assembly
**      03 Mar 01  ali	first port to os x using gcc, added the crucial __setflm
**                      definition.
**				1. removed double_t, put in double for now.
**				2. removed iclass from nearbyint.
**				3. removed wrong comments intrunc.
**				4.
**      13 May 97  ali  made performance improvements in rint, rinttol, roundtol
**                      and trunc by folding some of the taligent ideas into this
**                      implementation.  nearbyint is faster than the one in taligent,
**                      rint is more elegant, but slower by %30 than the taligent one.
**      09 Apr 97  ali  deleted modfl and deferred to AuxiliaryDD.c
**      15 Sep 94  ali  Major overhaul and performance improvements of all functions.
**      20 Jul 94  PAF  New faster version
**      16 Jul 93  ali  Added the modfl function.
**      18 Feb 93  ali  Changed the return value of fenv functions
**                      feclearexcept and feraiseexcept to their new
**                      NCEG X3J11.1/93-001 definitions.
**      16 Dec 92  JPO  Removed __itrunc implementation to a
**                      separate file.
**      15 Dec 92  JPO  Added __itrunc implementation and modified
**                      rinttol to include conversion from double
**                      to long int format.  Modified roundtol to
**                      call __itrunc.
**      10 Dec 92  JPO  Added modf (double) implementation.
**      04 Dec 92  JPO  First created.
**
*******************************************************************************/

#include <limits.h>
#include <math.h>

#define      SET_INVALID      0x01000000UL

typedef union
      {
      struct {
#if defined(__BIG_ENDIAN__)
        unsigned long int hi;
        unsigned long int lo;
#else
        unsigned long int lo;
        unsigned long int hi;
#endif
      } words;
      double dbl;
      } DblInHex;

static const unsigned long int signMask = 0x80000000ul;
static const double twoTo52      = 4503599627370496.0;
static const double doubleToLong = 4503603922337792.0;	            // 2^52
static const DblInHex Huge       = {{ 0x7FF00000, 0x00000000 }};
static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};


/*******************************************************************************
*                                                                              *
*     The function nearbyint rounds its double argument to integral value      *
*     according to the current rounding direction and returns the result in    *
*     double format.  This function does not signal inexact.                   *
*                                                                              *
********************************************************************************
*                                                                              *
*     This function calls fabs and copysign.		                         *
*                                                                              *
*******************************************************************************/

double nearbyint ( double x )
      {
	double y;
	double OldEnvironment;

	y = twoTo52;

	asm ("mffs %0" : "=f" (OldEnvironment));	/* get the environement */

      if ( fabs ( x ) >= y )                          /* huge case is exact */
            return x;
      if ( x < 0 ) y = -y;                                   /* negative case */
      y = ( x + y ) - y;                                    /* force rounding */
      if ( y == 0.0 )                        /* zero results mirror sign of x */
            y = copysign ( y, x );
//	restore old flags
	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
      return ( y );
	}

/*******************************************************************************
*                                                                              *
*     The function rinttol converts its double argument to integral value      *
*     according to the current rounding direction and returns the result in    *
*     long int format.  This conversion signals invalid if the argument is a   *
*     NaN or the rounded intermediate result is out of range of the            *
*     destination long int format, and it delivers an unspecified result in    *
*     this case.  This function signals inexact if the rounded result is       *
*     within range of the long int format but unequal to the operand.          *
*                                                                              *
*******************************************************************************/

long int rinttol ( double x )
      {
      register double y;
      DblInHex argument, OldEnvironment;
      unsigned long int xHead;
      register long int target;

      argument.dbl = x;
      target = ( argument.words.hi < signMask );        // flag positive sign
      xHead = argument.words.hi & 0x7ffffffful;         // high 32 bits of x

      if ( target )
/*******************************************************************************
*    Sign of x is positive.                                                    *
*******************************************************************************/
            {
            if ( xHead < 0x41dffffful )
                  {                                    // x is safely in long range
                  y = ( x + twoTo52 ) - twoTo52;       // round at binary point
                  argument.dbl = y + doubleToLong;     // force result into argument.words.lo
                  return ( ( long ) argument.words.lo );
                  }

		asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment

            if ( xHead > 0x41dffffful )
                  {                                    // x is safely out of long range
                  OldEnvironment.words.lo |= SET_INVALID;
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
                  return ( LONG_MAX );
                  }

/*******************************************************************************
*     x > 0.0 and may or may not be out of range of long.                      *
*******************************************************************************/

            y = ( x + twoTo52 ) - twoTo52;             // do rounding
            if ( y > ( double ) LONG_MAX )
                  {                                    // out of range of long
                  OldEnvironment.words.lo |= SET_INVALID;
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
                  return ( LONG_MAX );
                  }
            argument.dbl = y + doubleToLong;           // in range
            return ( ( long ) argument.words.lo );      // return result & flags
            }

/*******************************************************************************
*    Sign of x is negative.                                                    *
*******************************************************************************/
      if ( xHead < 0x41e00000ul )
            {                                          // x is safely in long range
            y = ( x - twoTo52 ) + twoTo52;
            argument.dbl = y + doubleToLong;
            return ( ( long ) argument.words.lo );
            }

	asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment

      if ( xHead > 0x41e00000ul )
            {                                          // x is safely out of long range
            OldEnvironment.words.lo |= SET_INVALID;
		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
            return ( LONG_MIN );
            }

/*******************************************************************************
*    x < 0.0 and may or may not be out of range of long.                       *
*******************************************************************************/

      y = ( x - twoTo52 ) + twoTo52;                   // do rounding
      if ( y < ( double ) LONG_MIN )
            {                                          // out of range of long
            OldEnvironment.words.lo |= SET_INVALID;
		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
            return ( LONG_MIN );
            }
      argument.dbl = y + doubleToLong;                       // in range
      return ( ( long ) argument.words.lo );           // return result & flags
      }

/*******************************************************************************
*                                                                              *
*     The function round rounds its double argument to integral value          *
*     according to the "add half to the magnitude and truncate" rounding of    *
*     Pascal's Round function and FORTRAN's ANINT function and returns the     *
*     result in double format.  This function signals inexact if an ordered    *
*     return value is not equal to the operand.                                *
*                                                                              *
*******************************************************************************/

double round ( double x )
      {
      DblInHex argument, OldEnvironment;
      register double y, z;
      register unsigned long int xHead;
      register long int target;

      argument.dbl = x;
      xHead = argument.words.hi & 0x7fffffffUL;      // xHead <- high half of |x|
      target = ( argument.words.hi < signMask );     // flag positive sign

      if ( xHead < 0x43300000ul )
/*******************************************************************************
*     Is |x| < 2.0^52?                                                        *
*******************************************************************************/
            {
            if ( xHead < 0x3ff00000ul )
/*******************************************************************************
*     Is |x| < 1.0?                                                           *
*******************************************************************************/
                  {
			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
                  if ( xHead < 0x3fe00000ul )
/*******************************************************************************
*     Is |x| < 0.5?                                                           *
*******************************************************************************/
                        {
                        if ( ( xHead | argument.words.lo ) != 0ul )
                              OldEnvironment.words.lo |= 0x02000000ul;
				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
                        if ( target )
                              return ( 0.0 );
                        else
                              return ( -0.0 );
                        }
/*******************************************************************************
*     Is 0.5 ² |x| < 1.0?                                                      *
*******************************************************************************/
                  OldEnvironment.words.lo |= 0x02000000ul;
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
                  if ( target )
                        return ( 1.0 );
                  else
                        return ( -1.0 );
                  }
/*******************************************************************************
*     Is 1.0 < |x| < 2.0^52?                                                   *
*******************************************************************************/
            if ( target )
                  {                                     // positive x
                  y = ( x + twoTo52 ) - twoTo52;        // round at binary point
                  if ( y == x )                         // exact case
                        return ( x );
                  z = x + 0.5;                          // inexact case
                  y = ( z + twoTo52 ) - twoTo52;        // round at binary point
                  if ( y > z )
                        return ( y - 1.0 );
                  else
                        return ( y );
                  }

/*******************************************************************************
*     Is x < 0?                                                                *
*******************************************************************************/
            else
                  {
                  y = ( x - twoTo52 ) + twoTo52;        // round at binary point
                  if ( y == x )
                        return ( x );
                  z = x - 0.5;
                  y = ( z - twoTo52 ) + twoTo52;        // round at binary point
                  if ( y < z )
                        return ( y + 1.0 );
                  else
                  return ( y );
                  }
            }
/*******************************************************************************
*      |x| >= 2.0^52 or x is a NaN.                                            *
*******************************************************************************/
      return ( x );
      }

/*******************************************************************************
*                                                                              *
*     The function roundtol converts its double argument to integral format    *
*     according to the "add half to the magnitude and chop" rounding mode of   *
*     Pascal's Round function and FORTRAN's NINT function.  This conversion    *
*     signals invalid if the argument is a NaN or the rounded intermediate     *
*     result is out of range of the destination long int format, and it        *
*     delivers an unspecified result in this case.  This function signals      *
*     inexact if the rounded result is within range of the long int format but *
*     unequal to the operand.                                                  *
*                                                                              *
*******************************************************************************/

long int roundtol ( double x )
	{
	register double y, z;
	DblInHex argument, OldEnvironment;
	register unsigned long int xhi;
	register long int target;
	const DblInHex kTZ = {{ 0x0, 0x1 }};
	const DblInHex kUP = {{ 0x0, 0x2 }};

	argument.dbl = x;
	xhi = argument.words.hi & 0x7ffffffful;	        	// high 32 bits of x
	target = ( argument.words.hi < signMask );         	// flag positive sign

	if ( xhi > 0x41e00000ul )
/*******************************************************************************
*     Is x is out of long range or NaN?                                        *
*******************************************************************************/
		{
		asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// get environment
		OldEnvironment.words.lo |= SET_INVALID;
		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
		if ( target )			              	// pin result
			return ( LONG_MAX );
		else
			return ( LONG_MIN );
		}

	if ( target )
/*******************************************************************************
*     Is sign of x is "+"?                                                     *
*******************************************************************************/
		{
		if ( x < 2147483647.5 )
/*******************************************************************************
*     x is in the range of a long.                                             *
*******************************************************************************/
			{
			y = ( x + doubleToLong ) - doubleToLong; 	// round at binary point
			if ( y != x )
				{		                    	// inexact case
				asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// save environment
				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kTZ.dbl )); // truncate rounding
				z = x + 0.5;		        	// truncate x + 0.5
				argument.dbl = z + doubleToLong;
				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
				return ( ( long ) argument.words.lo );
				}

			argument.dbl = y + doubleToLong;		// force result into argument.words.lo
			return ( ( long ) argument.words.lo ); 	// return long result
			}
/*******************************************************************************
*     Rounded positive x is out of the range of a long.                        *
*******************************************************************************/
		asm ("mffs %0" : "=f" (OldEnvironment.dbl));
		OldEnvironment.words.lo |= SET_INVALID;
		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
		return ( LONG_MAX );		              	// return pinned result
		}
/*******************************************************************************
*     x < 0.0 and may or may not be out of the range of a long.                *
*******************************************************************************/
	if ( x > -2147483648.5 )
/*******************************************************************************
*     x is in the range of a long.                                             *
*******************************************************************************/
		{
		y = ( x + doubleToLong ) - doubleToLong;	  	// round at binary point
		if ( y != x )
			{			                    	// inexact case
			asm ("mffs %0" : "=f" (OldEnvironment.dbl));	// save environment
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( kUP.dbl )); // round up
			z = x - 0.5;		              	// truncate x - 0.5
			argument.dbl = z + doubleToLong;
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
			return ( ( long ) argument.words.lo );
			}

		argument.dbl = y + doubleToLong;
		return ( ( long ) argument.words.lo );	  	//  return long result
		}
/*******************************************************************************
*     Rounded negative x is out of the range of a long.                        *
*******************************************************************************/
	asm ("mffs %0" : "=f" (OldEnvironment.dbl));
	OldEnvironment.words.lo |= SET_INVALID;
	asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
	return ( LONG_MIN );			              	// return pinned result
	}

/*******************************************************************************
*                                                                              *
*     The function trunc truncates its double argument to integral value       *
*     and returns the result in double format.  This function signals          *
*     inexact if an ordered return value is not equal to the operand.          *
*                                                                              *
*******************************************************************************/

double trunc ( double x )
      {
	DblInHex argument,OldEnvironment;
	register double y;
	register unsigned long int xhi;
	register long int target;

	argument.dbl = x;
	xhi = argument.words.hi & 0x7fffffffUL;	      	// xhi <- high half of |x|
	target = ( argument.words.hi < signMask );	      	// flag positive sign

	if ( xhi < 0x43300000ul )
/*******************************************************************************
*     Is |x| < 2.0^53?                                                         *
*******************************************************************************/
		{
		if ( xhi < 0x3ff00000ul )
/*******************************************************************************
*     Is |x| < 1.0?                                                            *
*******************************************************************************/
			{
			if ( ( xhi | argument.words.lo ) != 0ul )
				{                             	// raise deserved INEXACT
				asm ("mffs %0" : "=f" (OldEnvironment.dbl));
				OldEnvironment.words.lo |= 0x02000000ul;
				asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
				}
			if ( target )	                  	// return properly signed zero
				return ( 0.0 );
			else
				return ( -0.0 );
			}
/*******************************************************************************
*     Is 1.0 < |x| < 2.0^52?                                                   *
*******************************************************************************/
		if ( target )
			{
			y = ( x + twoTo52 ) - twoTo52;      	// round at binary point
			if ( y > x )
				return ( y - 1.0 );
			else
				return ( y );
			}

		else
			{
			y = ( x - twoTo52 ) + twoTo52;      	// round at binary point.
			if ( y < x )
				return ( y + 1.0 );
			else
				return ( y );
			}
		}
/*******************************************************************************
*      Is |x| >= 2.0^52 or x is a NaN.                                         *
*******************************************************************************/
	return ( x );
	}

/*******************************************************************************
*     The modf family of functions separate a floating-point number into its   *
*     fractional and integral parts, returning the fractional part and writing *
*     the integral part in floating-point format to the object pointed to by a *
*     pointer argument.  If the input argument is integral or infinite in      *
*     value, the return value is a zero with the sign of the input argument.   *
*     The modf family of functions raises no floating-point exceptions. older  *
*     implemenation set the INVALID flag due to signaling NaN input.           *
*                                                                              *
*******************************************************************************/

/*******************************************************************************
*     modf is the double implementation.                                       *
*******************************************************************************/

double modf ( double x, double *iptr )
      {
      register double OldEnvironment, xtrunc;
      register unsigned long int xHead, signBit;
      DblInHex argument;

      argument.dbl = x;
      xHead = argument.words.hi & 0x7ffffffful;            // |x| high bit pattern
      signBit = ( argument.words.hi & 0x80000000ul );      // isolate sign bit
	  if (xHead == 0x7ff81fe0)
			signBit = signBit | 0;

      if ( xHead < 0x43300000ul )
/*******************************************************************************
*     Is |x| < 2.0^53?                                                         *
*******************************************************************************/
            {
            if ( xHead < 0x3ff00000ul )
/*******************************************************************************
*     Is |x| < 1.0?                                                            *
*******************************************************************************/
                  {
                  argument.words.hi = signBit;             // truncate to zero
                  argument.words.lo = 0ul;
                  *iptr = argument.dbl;
                  return ( x );
                  }
/*******************************************************************************
*     Is 1.0 < |x| < 2.0^52?                                                   *
*******************************************************************************/
			asm ("mffs %0" : "=f" (OldEnvironment));	// save environment
			// round toward zero
			asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( TOWARDZERO.dbl ));
            if ( signBit == 0ul )                         // truncate to integer
                  xtrunc = ( x + twoTo52 ) - twoTo52;
            else
                  xtrunc = ( x - twoTo52 ) + twoTo52;
		// restore caller's env
		asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment ));
            *iptr = xtrunc;                               // store integral part
            if ( x != xtrunc )                            // nonzero fraction
                  return ( x - xtrunc );
            else
                  {                                       // zero with x's sign
                  argument.words.hi = signBit;
                  argument.words.lo = 0ul;
                  return ( argument.dbl );
                  }
            }

      *iptr = x;                                          // x is integral or NaN
      if ( x != x )                                       // NaN is returned
            return x;
      else
            {                                             // zero with x's sign
            argument.words.hi = signBit;
            argument.words.lo = 0ul;
            return ( argument.dbl );
            }
      }