summaryrefslogtreecommitdiff
path: root/sysdeps/ieee754/ldbl-128/s_lrintl.c
blob: 23f828f8623bd9ba38efb07de7e0df216feedb61 (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
/* Round argument to nearest integral value according to current rounding
   direction.
   Copyright (C) 1997-2016 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and
		  Jakub Jelinek <jj@ultra.linux.cz>, 1999.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, see
   <http://www.gnu.org/licenses/>.  */

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

#include <math_private.h>
#include <fix-fp-int-convert-overflow.h>

static const long double two112[2] =
{
  5.19229685853482762853049632922009600E+33L, /* 0x406F000000000000, 0 */
 -5.19229685853482762853049632922009600E+33L  /* 0xC06F000000000000, 0 */
};

long int
__lrintl (long double x)
{
  int32_t j0;
  u_int64_t i0,i1;
  long double w;
  long double t;
  long int result;
  int sx;

  GET_LDOUBLE_WORDS64 (i0, i1, x);
  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
  sx = i0 >> 63;
  i0 &= 0x0000ffffffffffffLL;
  i0 |= 0x0001000000000000LL;

  if (j0 < (int32_t) (8 * sizeof (long int)) - 1)
    {
      if (j0 < 48)
	{
#if defined FE_INVALID || defined FE_INEXACT
	  /* X < LONG_MAX + 1 implied by J0 < 31.  */
	  if (sizeof (long int) == 4
	      && x > (long double) LONG_MAX)
	    {
	      /* In the event of overflow we must raise the "invalid"
		 exception, but not "inexact".  */
	      t = __nearbyintl (x);
	      feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID);
	    }
	  else
#endif
	    {
	      w = two112[sx] + x;
	      t = w - two112[sx];
	    }
	  GET_LDOUBLE_WORDS64 (i0, i1, t);
	  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
	  i0 &= 0x0000ffffffffffffLL;
	  i0 |= 0x0001000000000000LL;

	  result = (j0 < 0 ? 0 : i0 >> (48 - j0));
	}
      else if (j0 >= 112)
	result = ((long int) i0 << (j0 - 48)) | (i1 << (j0 - 112));
      else
	{
#if defined FE_INVALID || defined FE_INEXACT
	  /* X < LONG_MAX + 1 implied by J0 < 63.  */
	  if (sizeof (long int) == 8
	      && x > (long double) LONG_MAX)
	    {
	      /* In the event of overflow we must raise the "invalid"
		 exception, but not "inexact".  */
	      t = __nearbyintl (x);
	      feraiseexcept (t == LONG_MAX ? FE_INEXACT : FE_INVALID);
	    }
	  else
#endif
	    {
	      w = two112[sx] + x;
	      t = w - two112[sx];
	    }
	  GET_LDOUBLE_WORDS64 (i0, i1, t);
	  j0 = ((i0 >> 48) & 0x7fff) - 0x3fff;
	  i0 &= 0x0000ffffffffffffLL;
	  i0 |= 0x0001000000000000LL;

	  if (j0 == 48)
	    result = (long int) i0;
	  else
	    result = ((long int) i0 << (j0 - 48)) | (i1 >> (112 - j0));
	}
    }
  else
    {
      /* The number is too large.  Unless it rounds to LONG_MIN,
	 FE_INVALID must be raised and the return value is
	 unspecified.  */
#if defined FE_INVALID || defined FE_INEXACT
      if (x < (long double) LONG_MIN
	  && x > (long double) LONG_MIN - 1.0L)
	{
	  /* If truncation produces LONG_MIN, the cast will not raise
	     the exception, but may raise "inexact".  */
	  t = __nearbyintl (x);
	  feraiseexcept (t == LONG_MIN ? FE_INEXACT : FE_INVALID);
	  return LONG_MIN;
	}
      else if (FIX_LDBL_LONG_CONVERT_OVERFLOW && x != (long double) LONG_MIN)
	{
	  feraiseexcept (FE_INVALID);
	  return sx == 0 ? LONG_MAX : LONG_MIN;
	}

#endif
      return (long int) x;
    }

  return sx ? -result : result;
}

weak_alias (__lrintl, lrintl)