Blame view

build2/epsilon-master/liba/src/external/openbsd/s_rint.c 1.97 KB
6663b6c9   adorian   projet complet av...
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
  /* @(#)s_rint.c 5.1 93/09/24 */
  /*
   * ====================================================
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice 
   * is preserved.
   * ====================================================
   */
  
  /*
   * rint(x)
   * Return x rounded to integral value according to the prevailing
   * rounding mode.
   * Method:
   *	Using floating addition.
   * Exception:
   *	Inexact flag raised if x not equal to rint(x).
   */
  
  #include <sys/cdefs.h>
  #include <float.h>
  #include <math.h>
  
  #include "math_private.h"
  
  static const double
  TWO52[2]={
    4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
   -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  };
  
  double
  rint(double x)
  {
  	int32_t i0,jj0,sx;
  	u_int32_t i,i1;
  	double w,t;
  	EXTRACT_WORDS(i0,i1,x);
  	sx = (i0>>31)&1;
  	jj0 = ((i0>>20)&0x7ff)-0x3ff;
  	if(jj0<20) {
  	    if(jj0<0) { 	
  		if(((i0&0x7fffffff)|i1)==0) return x;
  		i1 |= (i0&0x0fffff);
  		i0 &= 0xfffe0000;
  		i0 |= ((i1|-i1)>>12)&0x80000;
  		SET_HIGH_WORD(x,i0);
  	        w = TWO52[sx]+x;
  	        t =  w-TWO52[sx];
  		GET_HIGH_WORD(i0,t);
  		SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
  	        return t;
  	    } else {
  		i = (0x000fffff)>>jj0;
  		if(((i0&i)|i1)==0) return x; /* x is integral */
  		i>>=1;
  		if(((i0&i)|i1)!=0) {
  		    if(jj0==19) i1 = 0x40000000; else
  		    i0 = (i0&(~i))|((0x20000)>>jj0);
  		}
  	    }
  	} else if (jj0>51) {
  	    if(jj0==0x400) return x+x;	/* inf or NaN */
  	    else return x;		/* x is integral */
  	} else {
  	    i = ((u_int32_t)(0xffffffff))>>(jj0-20);
  	    if((i1&i)==0) return x;	/* x is integral */
  	    i>>=1;
  	    if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(jj0-20));
  	}
  	INSERT_WORDS(x,i0,i1);
  	w = TWO52[sx]+x;
  	return w-TWO52[sx];
  }
  
  #if LDBL_MANT_DIG == 53
  #ifdef __weak_alias
  __weak_alias(rintl, rint);
  #endif /* __weak_alias */
  #endif /* LDBL_MANT_DIG == 53 */