s_rintl.c (f839bac29c5d8167fb7a25138583006684bf9751) s_rintl.c (d3f9671a7d1d27db78ea3d24b9edaf22e8751b70)
1/*-
2 * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright

--- 18 unchanged lines hidden (view full) ---

27#include <sys/cdefs.h>
28__FBSDID("$FreeBSD$");
29
30#include <float.h>
31#include <math.h>
32
33#include "fpmath.h"
34
1/*-
2 * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright

--- 18 unchanged lines hidden (view full) ---

27#include <sys/cdefs.h>
28__FBSDID("$FreeBSD$");
29
30#include <float.h>
31#include <math.h>
32
33#include "fpmath.h"
34
35#if LDBL_MAX_EXP != 0x4000
36/* We also require the usual bias, min exp and expsign packing. */
37#error "Unsupported long double format"
38#endif
39
40#define BIAS (LDBL_MAX_EXP - 1)
41
42static const float
43shift[2] = {
35static const long double
36shift[2]={
44#if LDBL_MANT_DIG == 64
45 0x1.0p63, -0x1.0p63
46#elif LDBL_MANT_DIG == 113
47 0x1.0p112, -0x1.0p112
48#else
49#error "Unsupported long double format"
50#endif
51};
37#if LDBL_MANT_DIG == 64
38 0x1.0p63, -0x1.0p63
39#elif LDBL_MANT_DIG == 113
40 0x1.0p112, -0x1.0p112
41#else
42#error "Unsupported long double format"
43#endif
44};
52static const float zero[2] = { 0.0, -0.0 };
53
54long double
55rintl(long double x)
56{
57 union IEEEl2bits u;
45
46long double
47rintl(long double x)
48{
49 union IEEEl2bits u;
58 uint32_t expsign;
59 int ex, sign;
50 short sign;
60
61 u.e = x;
51
52 u.e = x;
62 expsign = u.xbits.expsign;
63 ex = expsign & 0x7fff;
64
53
65 if (ex >= BIAS + LDBL_MANT_DIG - 1) {
66 if (ex == BIAS + LDBL_MAX_EXP)
67 return (x + x); /* Inf, NaN, or unsupported format */
68 return (x); /* finite and already an integer */
54 if (u.bits.exp >= LDBL_MANT_DIG + LDBL_MAX_EXP - 2) {
55 /*
56 * The biased exponent is greater than the number of digits
57 * in the mantissa, so x is inf, NaN, or an integer.
58 */
59 if (u.bits.exp == 2 * LDBL_MAX_EXP - 1)
60 return (x + x); /* inf or NaN */
61 else
62 return (x);
69 }
63 }
70 sign = expsign >> 15;
71
72 /*
73 * The following code assumes that intermediate results are
74 * evaluated in long double precision. If they are evaluated in
64
65 /*
66 * The following code assumes that intermediate results are
67 * evaluated in long double precision. If they are evaluated in
75 * greater precision, double rounding may occur, and if they are
68 * greater precision, double rounding will occur, and if they are
76 * evaluated in less precision (as on i386), results will be
77 * wildly incorrect.
78 */
69 * evaluated in less precision (as on i386), results will be
70 * wildly incorrect.
71 */
79 x += shift[sign];
80 x -= shift[sign];
81
82 /*
83 * If the result is +-0, then it must have the same sign as x, but
84 * the above calculation doesn't always give this. Fix up the sign.
85 */
86 if (ex < BIAS && x == 0.0L)
87 return (zero[sign]);
88
89 return (x);
72 sign = u.bits.sign;
73 u.e = shift[sign] + x;
74 u.e -= shift[sign];
75 u.bits.sign = sign;
76 return (u.e);
90}
77}