1*05a0b428SJohn Marino /* $OpenBSD: s_lrint.c,v 1.10 2013/07/03 04:46:36 espie Exp $ */
2*05a0b428SJohn Marino /* $NetBSD: lrint.c,v 1.3 2004/10/13 15:18:32 drochner Exp $ */
3*05a0b428SJohn Marino
4*05a0b428SJohn Marino /*-
5*05a0b428SJohn Marino * Copyright (c) 2004
6*05a0b428SJohn Marino * Matthias Drochner. All rights reserved.
7*05a0b428SJohn Marino *
8*05a0b428SJohn Marino * Redistribution and use in source and binary forms, with or without
9*05a0b428SJohn Marino * modification, are permitted provided that the following conditions
10*05a0b428SJohn Marino * are met:
11*05a0b428SJohn Marino * 1. Redistributions of source code must retain the above copyright
12*05a0b428SJohn Marino * notice, this list of conditions and the following disclaimer.
13*05a0b428SJohn Marino * 2. Redistributions in binary form must reproduce the above copyright
14*05a0b428SJohn Marino * notice, this list of conditions and the following disclaimer in the
15*05a0b428SJohn Marino * documentation and/or other materials provided with the distribution.
16*05a0b428SJohn Marino *
17*05a0b428SJohn Marino * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
18*05a0b428SJohn Marino * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19*05a0b428SJohn Marino * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20*05a0b428SJohn Marino * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
21*05a0b428SJohn Marino * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22*05a0b428SJohn Marino * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23*05a0b428SJohn Marino * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24*05a0b428SJohn Marino * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25*05a0b428SJohn Marino * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26*05a0b428SJohn Marino * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27*05a0b428SJohn Marino * SUCH DAMAGE.
28*05a0b428SJohn Marino */
29*05a0b428SJohn Marino
30*05a0b428SJohn Marino #include <sys/types.h>
31*05a0b428SJohn Marino #include <sys/limits.h>
32*05a0b428SJohn Marino #include <float.h>
33*05a0b428SJohn Marino #include <math.h>
34*05a0b428SJohn Marino #include <ieeefp.h>
35*05a0b428SJohn Marino #include <machine/ieee.h>
36*05a0b428SJohn Marino
37*05a0b428SJohn Marino #include "math_private.h"
38*05a0b428SJohn Marino
39*05a0b428SJohn Marino #ifndef LRINTNAME
40*05a0b428SJohn Marino #define LRINTNAME lrint
41*05a0b428SJohn Marino #define RESTYPE long int
42*05a0b428SJohn Marino #define RESTYPE_MIN LONG_MIN
43*05a0b428SJohn Marino #define RESTYPE_MAX LONG_MAX
44*05a0b428SJohn Marino #endif
45*05a0b428SJohn Marino
46*05a0b428SJohn Marino #define RESTYPE_BITS (sizeof(RESTYPE) * 8)
47*05a0b428SJohn Marino
48*05a0b428SJohn Marino static const double
49*05a0b428SJohn Marino TWO52[2]={
50*05a0b428SJohn Marino 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
51*05a0b428SJohn Marino -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
52*05a0b428SJohn Marino };
53*05a0b428SJohn Marino
54*05a0b428SJohn Marino RESTYPE
LRINTNAME(double x)55*05a0b428SJohn Marino LRINTNAME(double x)
56*05a0b428SJohn Marino {
57*05a0b428SJohn Marino u_int32_t i0, i1;
58*05a0b428SJohn Marino int e, s, shift;
59*05a0b428SJohn Marino RESTYPE res;
60*05a0b428SJohn Marino
61*05a0b428SJohn Marino GET_HIGH_WORD(i0, x);
62*05a0b428SJohn Marino e = i0 >> DBL_FRACHBITS;
63*05a0b428SJohn Marino s = e >> DBL_EXPBITS;
64*05a0b428SJohn Marino e = (e & 0x7ff) - DBL_EXP_BIAS;
65*05a0b428SJohn Marino
66*05a0b428SJohn Marino /* 1.0 x 2^31 (or 2^63) is already too large */
67*05a0b428SJohn Marino if (e >= (int)RESTYPE_BITS - 1)
68*05a0b428SJohn Marino return (s ? RESTYPE_MIN : RESTYPE_MAX); /* ??? unspecified */
69*05a0b428SJohn Marino
70*05a0b428SJohn Marino /* >= 2^52 is already an exact integer */
71*05a0b428SJohn Marino if (e < DBL_FRACBITS) {
72*05a0b428SJohn Marino volatile double t = x; /* clip extra precision */
73*05a0b428SJohn Marino /* round, using current direction */
74*05a0b428SJohn Marino t += TWO52[s];
75*05a0b428SJohn Marino t -= TWO52[s];
76*05a0b428SJohn Marino x = t;
77*05a0b428SJohn Marino }
78*05a0b428SJohn Marino
79*05a0b428SJohn Marino EXTRACT_WORDS(i0, i1, x);
80*05a0b428SJohn Marino e = ((i0 >> DBL_FRACHBITS) & 0x7ff) - DBL_EXP_BIAS;
81*05a0b428SJohn Marino i0 &= 0xfffff;
82*05a0b428SJohn Marino i0 |= (1 << DBL_FRACHBITS);
83*05a0b428SJohn Marino
84*05a0b428SJohn Marino if (e < 0)
85*05a0b428SJohn Marino return (0);
86*05a0b428SJohn Marino
87*05a0b428SJohn Marino shift = e - DBL_FRACBITS;
88*05a0b428SJohn Marino if (shift >=0)
89*05a0b428SJohn Marino res = (shift < RESTYPE_BITS ? (RESTYPE)i1 << shift : 0);
90*05a0b428SJohn Marino else
91*05a0b428SJohn Marino res = (shift > -RESTYPE_BITS ? (RESTYPE)i1 >> -shift : 0);
92*05a0b428SJohn Marino shift += 32;
93*05a0b428SJohn Marino if (shift >=0)
94*05a0b428SJohn Marino res |= (shift < RESTYPE_BITS ? (RESTYPE)i0 << shift : 0);
95*05a0b428SJohn Marino else
96*05a0b428SJohn Marino res |= (shift > -RESTYPE_BITS ? (RESTYPE)i0 >> -shift : 0);
97*05a0b428SJohn Marino
98*05a0b428SJohn Marino return (s ? -res : res);
99*05a0b428SJohn Marino }
100*05a0b428SJohn Marino
101*05a0b428SJohn Marino #if LDBL_MANT_DIG == DBL_MANT_DIG
102*05a0b428SJohn Marino __strong_alias(lrintl, lrint);
103*05a0b428SJohn Marino #endif /* LDBL_MANT_DIG == DBL_MANT_DIG */
104