1/*
2 * Copyright © 2015 Intel Corporation
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 * and/or sell copies of the Software, and to permit persons to whom the
9 * Software is furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice (including the next
12 * paragraph) shall be included in all copies or substantial portions of the
13 * Software.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
18 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
20 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
21 * IN THE SOFTWARE.
22 */
23
24#ifndef _ROUNDING_H
25#define _ROUNDING_H
26
27#include "c99_math.h"
28
29#include <limits.h>
30#include <stdint.h>
31
32#if defined(__SSE__) || defined(_MSC_VER)
33/* MSVC always has SSE nowadays */
34#include <xmmintrin.h>
35#include <emmintrin.h>
36#endif
37
38#ifdef __SSE4_1__
39#include <smmintrin.h>
40#endif
41
42/* The C standard library has functions round()/rint()/nearbyint() that round
43 * their arguments according to the rounding mode set in the floating-point
44 * control register. While there are trunc()/ceil()/floor() functions that do
45 * a specific operation without modifying the rounding mode, there is no
46 * roundeven() in any version of C.
47 *
48 * Technical Specification 18661 (ISO/IEC TS 18661-1:2014) adds roundeven(),
49 * but it's unfortunately not implemented by glibc.
50 *
51 * This implementation differs in that it does not raise the inexact exception.
52 *
53 * We use rint() to implement these functions, with the assumption that the
54 * floating-point rounding mode has not been changed from the default Round
55 * to Nearest.
56 */
57
58/**
59 * \brief Rounds \c x to the nearest integer, with ties to the even integer.
60 */
61static inline float
62_mesa_roundevenf(float x)
63{
64#ifdef __SSE4_1__
65   float ret;
66   __m128 m = _mm_load_ss(&x);
67   m = _mm_round_ss(m, m, _MM_FROUND_CUR_DIRECTION | _MM_FROUND_NO_EXC);
68   _mm_store_ss(&ret, m);
69   return ret;
70#else
71   return rintf(x);
72#endif
73}
74
75/**
76 * \brief Rounds \c x to the nearest integer, with ties to the even integer.
77 */
78static inline double
79_mesa_roundeven(double x)
80{
81#ifdef __SSE4_1__
82   double ret;
83   __m128d m = _mm_load_sd(&x);
84   m = _mm_round_sd(m, m, _MM_FROUND_CUR_DIRECTION | _MM_FROUND_NO_EXC);
85   _mm_store_sd(&ret, m);
86   return ret;
87#else
88   return rint(x);
89#endif
90}
91
92/**
93 * \brief Rounds \c x to the nearest integer, with ties to the even integer,
94 * and returns the value as a long int.
95 */
96static inline long
97_mesa_lroundevenf(float x)
98{
99#if defined(__SSE__) || defined(_MSC_VER)
100#if LONG_MAX == INT64_MAX
101   return _mm_cvtss_si64(_mm_load_ss(&x));
102#elif LONG_MAX == INT32_MAX
103   return _mm_cvtss_si32(_mm_load_ss(&x));
104#else
105#error "Unsupported long size"
106#endif
107#else
108   return lrintf(x);
109#endif
110}
111
112/**
113 * \brief Rounds \c x to the nearest integer, with ties to the even integer,
114 * and returns the value as a long int.
115 */
116static inline long
117_mesa_lroundeven(double x)
118{
119#if defined(__SSE2__) || defined(_MSC_VER)
120#if LONG_MAX == INT64_MAX
121   return _mm_cvtsd_si64(_mm_load_sd(&x));
122#elif LONG_MAX == INT32_MAX
123   return _mm_cvtsd_si32(_mm_load_sd(&x));
124#else
125#error "Unsupported long size"
126#endif
127#else
128   return lrint(x);
129#endif
130}
131
132/**
133 * \brief Rounds \c x to the nearest integer, with ties to the even integer,
134 * and returns the value as an int64_t.
135 */
136static inline int64_t
137_mesa_i64roundevenf(float x)
138{
139#if LONG_MAX == INT64_MAX
140   return _mesa_lroundevenf(x);
141#elif LONG_MAX == INT32_MAX
142   return llrintf(x);
143#else
144#error "Unsupported long size"
145#endif
146}
147
148#endif
149