1/*
2 * Copyright (C) 2011 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17#define __STDC_LIMIT_MACROS
18
19#include <assert.h>
20#include <stdint.h>
21
22#include <utils/LinearTransform.h>
23
24namespace android {
25
26template<class T> static inline T ABS(T x) { return (x < 0) ? -x : x; }
27
28// Static math methods involving linear transformations
29static bool scale_u64_to_u64(
30        uint64_t val,
31        uint32_t N,
32        uint32_t D,
33        uint64_t* res,
34        bool round_up_not_down) {
35    uint64_t tmp1, tmp2;
36    uint32_t r;
37
38    assert(res);
39    assert(D);
40
41    // Let U32(X) denote a uint32_t containing the upper 32 bits of a 64 bit
42    // integer X.
43    // Let L32(X) denote a uint32_t containing the lower 32 bits of a 64 bit
44    // integer X.
45    // Let X[A, B] with A <= B denote bits A through B of the integer X.
46    // Let (A | B) denote the concatination of two 32 bit ints, A and B.
47    // IOW X = (A | B) => U32(X) == A && L32(X) == B
48    //
49    // compute M = val * N (a 96 bit int)
50    // ---------------------------------
51    // tmp2 = U32(val) * N (a 64 bit int)
52    // tmp1 = L32(val) * N (a 64 bit int)
53    // which means
54    // M = val * N = (tmp2 << 32) + tmp1
55    tmp2 = (val >> 32) * N;
56    tmp1 = (val & UINT32_MAX) * N;
57
58    // compute M[32, 95]
59    // tmp2 = tmp2 + U32(tmp1)
60    //      = (U32(val) * N) + U32(L32(val) * N)
61    //      = M[32, 95]
62    tmp2 += tmp1 >> 32;
63
64    // if M[64, 95] >= D, then M/D has bits > 63 set and we have
65    // an overflow.
66    if ((tmp2 >> 32) >= D) {
67        *res = UINT64_MAX;
68        return false;
69    }
70
71    // Divide.  Going in we know
72    // tmp2 = M[32, 95]
73    // U32(tmp2) < D
74    r = tmp2 % D;
75    tmp2 /= D;
76
77    // At this point
78    // tmp1      = L32(val) * N
79    // tmp2      = M[32, 95] / D
80    //           = (M / D)[32, 95]
81    // r         = M[32, 95] % D
82    // U32(tmp2) = 0
83    //
84    // compute tmp1 = (r | M[0, 31])
85    tmp1 = (tmp1 & UINT32_MAX) | ((uint64_t)r << 32);
86
87    // Divide again.  Keep the remainder around in order to round properly.
88    r = tmp1 % D;
89    tmp1 /= D;
90
91    // At this point
92    // tmp2      = (M / D)[32, 95]
93    // tmp1      = (M / D)[ 0, 31]
94    // r         =  M % D
95    // U32(tmp1) = 0
96    // U32(tmp2) = 0
97
98    // Pack the result and deal with the round-up case (As well as the
99    // remote possiblility over overflow in such a case).
100    *res = (tmp2 << 32) | tmp1;
101    if (r && round_up_not_down) {
102        ++(*res);
103        if (!(*res)) {
104            *res = UINT64_MAX;
105            return false;
106        }
107    }
108
109    return true;
110}
111
112static bool linear_transform_s64_to_s64(
113        int64_t  val,
114        int64_t  basis1,
115        int32_t  N,
116        uint32_t D,
117        bool     invert_frac,
118        int64_t  basis2,
119        int64_t* out) {
120    uint64_t scaled, res;
121    uint64_t abs_val;
122    bool is_neg;
123
124    if (!out)
125        return false;
126
127    // Compute abs(val - basis_64). Keep track of whether or not this delta
128    // will be negative after the scale opertaion.
129    if (val < basis1) {
130        is_neg = true;
131        abs_val = basis1 - val;
132    } else {
133        is_neg = false;
134        abs_val = val - basis1;
135    }
136
137    if (N < 0)
138        is_neg = !is_neg;
139
140    if (!scale_u64_to_u64(abs_val,
141                          invert_frac ? D : ABS(N),
142                          invert_frac ? ABS(N) : D,
143                          &scaled,
144                          is_neg))
145        return false; // overflow/undeflow
146
147    // if scaled is >= 0x8000<etc>, then we are going to overflow or
148    // underflow unless ABS(basis2) is large enough to pull us back into the
149    // non-overflow/underflow region.
150    if (scaled & INT64_MIN) {
151        if (is_neg && (basis2 < 0))
152            return false; // certain underflow
153
154        if (!is_neg && (basis2 >= 0))
155            return false; // certain overflow
156
157        if (ABS(basis2) <= static_cast<int64_t>(scaled & INT64_MAX))
158            return false; // not enough
159
160        // Looks like we are OK
161        *out = (is_neg ? (-scaled) : scaled) + basis2;
162    } else {
163        // Scaled fits within signed bounds, so we just need to check for
164        // over/underflow for two signed integers.  Basically, if both scaled
165        // and basis2 have the same sign bit, and the result has a different
166        // sign bit, then we have under/overflow.  An easy way to compute this
167        // is
168        // (scaled_signbit XNOR basis_signbit) &&
169        // (scaled_signbit XOR res_signbit)
170        // ==
171        // (scaled_signbit XOR basis_signbit XOR 1) &&
172        // (scaled_signbit XOR res_signbit)
173
174        if (is_neg)
175            scaled = -scaled;
176        res = scaled + basis2;
177
178        if ((scaled ^ basis2 ^ INT64_MIN) & (scaled ^ res) & INT64_MIN)
179            return false;
180
181        *out = res;
182    }
183
184    return true;
185}
186
187bool LinearTransform::doForwardTransform(int64_t a_in, int64_t* b_out) const {
188    if (0 == a_to_b_denom)
189        return false;
190
191    return linear_transform_s64_to_s64(a_in,
192                                       a_zero,
193                                       a_to_b_numer,
194                                       a_to_b_denom,
195                                       false,
196                                       b_zero,
197                                       b_out);
198}
199
200bool LinearTransform::doReverseTransform(int64_t b_in, int64_t* a_out) const {
201    if (0 == a_to_b_numer)
202        return false;
203
204    return linear_transform_s64_to_s64(b_in,
205                                       b_zero,
206                                       a_to_b_numer,
207                                       a_to_b_denom,
208                                       true,
209                                       a_zero,
210                                       a_out);
211}
212
213template <class T> void LinearTransform::reduce(T* N, T* D) {
214    T a, b;
215    if (!N || !D || !(*D)) {
216        assert(false);
217        return;
218    }
219
220    a = *N;
221    b = *D;
222
223    if (a == 0) {
224        *D = 1;
225        return;
226    }
227
228    // This implements Euclid's method to find GCD.
229    if (a < b) {
230        T tmp = a;
231        a = b;
232        b = tmp;
233    }
234
235    while (1) {
236        // a is now the greater of the two.
237        const T remainder = a % b;
238        if (remainder == 0) {
239            *N /= b;
240            *D /= b;
241            return;
242        }
243        // by swapping remainder and b, we are guaranteeing that a is
244        // still the greater of the two upon entrance to the loop.
245        a = b;
246        b = remainder;
247    }
248};
249
250template void LinearTransform::reduce<uint64_t>(uint64_t* N, uint64_t* D);
251template void LinearTransform::reduce<uint32_t>(uint32_t* N, uint32_t* D);
252
253void LinearTransform::reduce(int32_t* N, uint32_t* D) {
254    if (N && D && *D) {
255        if (*N < 0) {
256            *N = -(*N);
257            reduce(reinterpret_cast<uint32_t*>(N), D);
258            *N = -(*N);
259        } else {
260            reduce(reinterpret_cast<uint32_t*>(N), D);
261        }
262    }
263}
264
265}  // namespace android
266