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 
24 namespace android {
25 
ABS(T x)26 template<class T> static inline T ABS(T x) { return (x < 0) ? -x : x; }
27 
28 // Static math methods involving linear transformations
scale_u64_to_u64(uint64_t val,uint32_t N,uint32_t D,uint64_t * res,bool round_up_not_down)29 static 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 
linear_transform_s64_to_s64(int64_t val,int64_t basis1,int32_t N,uint32_t D,bool invert_frac,int64_t basis2,int64_t * out)112 static 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 
doForwardTransform(int64_t a_in,int64_t * b_out) const187 bool 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 
doReverseTransform(int64_t b_in,int64_t * a_out) const200 bool 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 
reduce(T * N,T * D)213 template <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 
250 template void LinearTransform::reduce<uint64_t>(uint64_t* N, uint64_t* D);
251 template void LinearTransform::reduce<uint32_t>(uint32_t* N, uint32_t* D);
252 
reduce(int32_t * N,uint32_t * D)253 void 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