1/*
2 * Copyright (c) 2014 Advanced Micro Devices, Inc.
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a copy
5 * of this software and associated documentation files (the "Software"), to deal
6 * in the Software without restriction, including without limitation the rights
7 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8 * copies of the Software, and to permit persons to whom the Software is
9 * furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
17 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
18 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
19 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
20 * THE SOFTWARE.
21 */
22
23#include <clc/clc.h>
24
25#include <math/clc_remainder.h>
26#include "../clcmacro.h"
27#include "config.h"
28#include "math.h"
29
30_CLC_DEF _CLC_OVERLOAD float __clc_remquo(float x, float y, __private int *quo)
31{
32    x = __clc_flush_denormal_if_not_supported(x);
33    y = __clc_flush_denormal_if_not_supported(y);
34    int ux = as_int(x);
35    int ax = ux & EXSIGNBIT_SP32;
36    float xa = as_float(ax);
37    int sx = ux ^ ax;
38    int ex = ax >> EXPSHIFTBITS_SP32;
39
40    int uy = as_int(y);
41    int ay = uy & EXSIGNBIT_SP32;
42    float ya = as_float(ay);
43    int sy = uy ^ ay;
44    int ey = ay >> EXPSHIFTBITS_SP32;
45
46    float xr = as_float(0x3f800000 | (ax & 0x007fffff));
47    float yr = as_float(0x3f800000 | (ay & 0x007fffff));
48    int c;
49    int k = ex - ey;
50
51    uint q = 0;
52
53    while (k > 0) {
54        c = xr >= yr;
55        q = (q << 1) | c;
56        xr -= c ? yr : 0.0f;
57        xr += xr;
58	--k;
59    }
60
61    c = xr > yr;
62    q = (q << 1) | c;
63    xr -= c ? yr : 0.0f;
64
65    int lt = ex < ey;
66
67    q = lt ? 0 : q;
68    xr = lt ? xa : xr;
69    yr = lt ? ya : yr;
70
71    c = (yr < 2.0f * xr) | ((yr == 2.0f * xr) & ((q & 0x1) == 0x1));
72    xr -= c ? yr : 0.0f;
73    q += c;
74
75    float s = as_float(ey << EXPSHIFTBITS_SP32);
76    xr *= lt ? 1.0f : s;
77
78    int qsgn = sx == sy ? 1 : -1;
79    int quot = (q & 0x7f) * qsgn;
80
81    c = ax == ay;
82    quot = c ? qsgn : quot;
83    xr = c ? 0.0f : xr;
84
85    xr = as_float(sx ^ as_int(xr));
86
87    c = ax > PINFBITPATT_SP32 | ay > PINFBITPATT_SP32 | ax == PINFBITPATT_SP32 | ay == 0;
88    quot = c ? 0 : quot;
89    xr = c ? as_float(QNANBITPATT_SP32) : xr;
90
91    *quo = quot;
92
93    return xr;
94}
95// remquo singature is special, we don't have macro for this
96#define __VEC_REMQUO(TYPE, VEC_SIZE, HALF_VEC_SIZE) \
97_CLC_DEF _CLC_OVERLOAD TYPE##VEC_SIZE __clc_remquo(TYPE##VEC_SIZE x, TYPE##VEC_SIZE y, __private int##VEC_SIZE *quo) \
98{ \
99	int##HALF_VEC_SIZE lo, hi; \
100	TYPE##VEC_SIZE ret; \
101	ret.lo = __clc_remquo(x.lo, y.lo, &lo); \
102	ret.hi = __clc_remquo(x.hi, y.hi, &hi); \
103	(*quo).lo = lo; \
104	(*quo).hi = hi; \
105	return ret; \
106}
107__VEC_REMQUO(float, 2,)
108__VEC_REMQUO(float, 3, 2)
109__VEC_REMQUO(float, 4, 2)
110__VEC_REMQUO(float, 8, 4)
111__VEC_REMQUO(float, 16, 8)
112
113#ifdef cl_khr_fp64
114_CLC_DEF _CLC_OVERLOAD double __clc_remquo(double x, double y, __private int *pquo)
115{
116    ulong ux = as_ulong(x);
117    ulong ax = ux & ~SIGNBIT_DP64;
118    ulong xsgn = ux ^ ax;
119    double dx = as_double(ax);
120    int xexp = convert_int(ax >> EXPSHIFTBITS_DP64);
121    int xexp1 = 11 - (int) clz(ax & MANTBITS_DP64);
122    xexp1 = xexp < 1 ? xexp1 : xexp;
123
124    ulong uy = as_ulong(y);
125    ulong ay = uy & ~SIGNBIT_DP64;
126    double dy = as_double(ay);
127    int yexp = convert_int(ay >> EXPSHIFTBITS_DP64);
128    int yexp1 = 11 - (int) clz(ay & MANTBITS_DP64);
129    yexp1 = yexp < 1 ? yexp1 : yexp;
130
131    int qsgn = ((ux ^ uy) & SIGNBIT_DP64) == 0UL ? 1 : -1;
132
133    // First assume |x| > |y|
134
135    // Set ntimes to the number of times we need to do a
136    // partial remainder. If the exponent of x is an exact multiple
137    // of 53 larger than the exponent of y, and the mantissa of x is
138    // less than the mantissa of y, ntimes will be one too large
139    // but it doesn't matter - it just means that we'll go round
140    // the loop below one extra time.
141    int ntimes = max(0, (xexp1 - yexp1) / 53);
142    double w =  ldexp(dy, ntimes * 53);
143    w = ntimes == 0 ? dy : w;
144    double scale = ntimes == 0 ? 1.0 : 0x1.0p-53;
145
146    // Each time round the loop we compute a partial remainder.
147    // This is done by subtracting a large multiple of w
148    // from x each time, where w is a scaled up version of y.
149    // The subtraction must be performed exactly in quad
150    // precision, though the result at each stage can
151    // fit exactly in a double precision number.
152    int i;
153    double t, v, p, pp;
154
155    for (i = 0; i < ntimes; i++) {
156        // Compute integral multiplier
157        t = trunc(dx / w);
158
159        // Compute w * t in quad precision
160        p = w * t;
161        pp = fma(w, t, -p);
162
163        // Subtract w * t from dx
164        v = dx - p;
165        dx = v + (((dx - v) - p) - pp);
166
167        // If t was one too large, dx will be negative. Add back one w.
168        dx += dx < 0.0 ? w : 0.0;
169
170        // Scale w down by 2^(-53) for the next iteration
171        w *= scale;
172    }
173
174    // One more time
175    // Variable todd says whether the integer t is odd or not
176    t = floor(dx / w);
177    long lt = (long)t;
178    int todd = lt & 1;
179
180    p = w * t;
181    pp = fma(w, t, -p);
182    v = dx - p;
183    dx = v + (((dx - v) - p) - pp);
184    i = dx < 0.0;
185    todd ^= i;
186    dx += i ? w : 0.0;
187
188    lt -= i;
189
190    // At this point, dx lies in the range [0,dy)
191
192    // For the remainder function, we need to adjust dx
193    // so that it lies in the range (-y/2, y/2] by carefully
194    // subtracting w (== dy == y) if necessary. The rigmarole
195    // with todd is to get the correct sign of the result
196    // when x/y lies exactly half way between two integers,
197    // when we need to choose the even integer.
198
199    int al = (2.0*dx > w) | (todd & (2.0*dx == w));
200    double dxl = dx - (al ? w : 0.0);
201
202    int ag = (dx > 0.5*w) | (todd & (dx == 0.5*w));
203    double dxg = dx - (ag ? w : 0.0);
204
205    dx = dy < 0x1.0p+1022 ? dxl : dxg;
206    lt += dy < 0x1.0p+1022 ? al : ag;
207    int quo = ((int)lt & 0x7f) * qsgn;
208
209    double ret = as_double(xsgn ^ as_ulong(dx));
210    dx = as_double(ax);
211
212    // Now handle |x| == |y|
213    int c = dx == dy;
214    t = as_double(xsgn);
215    quo = c ? qsgn : quo;
216    ret = c ? t : ret;
217
218    // Next, handle |x| < |y|
219    c = dx < dy;
220    quo = c ? 0 : quo;
221    ret = c ? x : ret;
222
223    c &= (yexp < 1023 & 2.0*dx > dy) | (dx > 0.5*dy);
224    quo = c ? qsgn : quo;
225    // we could use a conversion here instead since qsgn = +-1
226    p = qsgn == 1 ? -1.0 : 1.0;
227    t = fma(y, p, x);
228    ret = c ? t : ret;
229
230    // We don't need anything special for |x| == 0
231
232    // |y| is 0
233    c = dy == 0.0;
234    quo = c ? 0 : quo;
235    ret = c ? as_double(QNANBITPATT_DP64) : ret;
236
237    // y is +-Inf, NaN
238    c = yexp > BIASEDEMAX_DP64;
239    quo = c ? 0 : quo;
240    t = y == y ? x : y;
241    ret = c ? t : ret;
242
243    // x is +=Inf, NaN
244    c = xexp > BIASEDEMAX_DP64;
245    quo = c ? 0 : quo;
246    ret = c ? as_double(QNANBITPATT_DP64) : ret;
247
248    *pquo = quo;
249    return ret;
250}
251__VEC_REMQUO(double, 2,)
252__VEC_REMQUO(double, 3, 2)
253__VEC_REMQUO(double, 4, 2)
254__VEC_REMQUO(double, 8, 4)
255__VEC_REMQUO(double, 16, 8)
256#endif
257