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,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 int64_t basis2,
118 int64_t* out) {
119 uint64_t scaled, res;
120 uint64_t abs_val;
121 bool is_neg;
122
123 if (!out)
124 return false;
125
126 // Compute abs(val - basis_64). Keep track of whether or not this delta
127 // will be negative after the scale opertaion.
128 if (val < basis1) {
129 is_neg = true;
130 abs_val = basis1 - val;
131 } else {
132 is_neg = false;
133 abs_val = val - basis1;
134 }
135
136 if (N < 0)
137 is_neg = !is_neg;
138
139 if (!scale_u64_to_u64(abs_val,
140 ABS(N),
141 D,
142 &scaled,
143 is_neg))
144 return false; // overflow/undeflow
145
146 // if scaled is >= 0x8000<etc>, then we are going to overflow or
147 // underflow unless ABS(basis2) is large enough to pull us back into the
148 // non-overflow/underflow region.
149 if (scaled & INT64_MIN) {
150 if (is_neg && (basis2 < 0))
151 return false; // certain underflow
152
153 if (!is_neg && (basis2 >= 0))
154 return false; // certain overflow
155
156 if (ABS(basis2) <= static_cast<int64_t>(scaled & INT64_MAX))
157 return false; // not enough
158
159 // Looks like we are OK
160 *out = (is_neg ? (-scaled) : scaled) + basis2;
161 } else {
162 // Scaled fits within signed bounds, so we just need to check for
163 // over/underflow for two signed integers. Basically, if both scaled
164 // and basis2 have the same sign bit, and the result has a different
165 // sign bit, then we have under/overflow. An easy way to compute this
166 // is
167 // (scaled_signbit XNOR basis_signbit) &&
168 // (scaled_signbit XOR res_signbit)
169 // ==
170 // (scaled_signbit XOR basis_signbit XOR 1) &&
171 // (scaled_signbit XOR res_signbit)
172
173 if (is_neg)
174 scaled = -scaled;
175 res = scaled + basis2;
176
177 if ((scaled ^ basis2 ^ INT64_MIN) & (scaled ^ res) & INT64_MIN)
178 return false;
179
180 *out = res;
181 }
182
183 return true;
184 }
185
doForwardTransform(int64_t a_in,int64_t * b_out) const186 bool LinearTransform::doForwardTransform(int64_t a_in, int64_t* b_out) const {
187 if (0 == a_to_b_denom)
188 return false;
189
190 return linear_transform_s64_to_s64(a_in,
191 a_zero,
192 a_to_b_numer,
193 a_to_b_denom,
194 b_zero,
195 b_out);
196 }
197
doReverseTransform(int64_t b_in,int64_t * a_out) const198 bool LinearTransform::doReverseTransform(int64_t b_in, int64_t* a_out) const {
199 if (0 == a_to_b_numer)
200 return false;
201
202 return linear_transform_s64_to_s64(b_in,
203 b_zero,
204 a_to_b_denom,
205 a_to_b_numer,
206 a_zero,
207 a_out);
208 }
209
reduce(T * N,T * D)210 template <class T> void LinearTransform::reduce(T* N, T* D) {
211 T a, b;
212 if (!N || !D || !(*D)) {
213 assert(false);
214 return;
215 }
216
217 a = *N;
218 b = *D;
219
220 if (a == 0) {
221 *D = 1;
222 return;
223 }
224
225 // This implements Euclid's method to find GCD.
226 if (a < b) {
227 T tmp = a;
228 a = b;
229 b = tmp;
230 }
231
232 while (1) {
233 // a is now the greater of the two.
234 const T remainder = a % b;
235 if (remainder == 0) {
236 *N /= b;
237 *D /= b;
238 return;
239 }
240 // by swapping remainder and b, we are guaranteeing that a is
241 // still the greater of the two upon entrance to the loop.
242 a = b;
243 b = remainder;
244 }
245 };
246
247 template void LinearTransform::reduce<uint64_t>(uint64_t* N, uint64_t* D);
248 template void LinearTransform::reduce<uint32_t>(uint32_t* N, uint32_t* D);
249
reduce(int32_t * N,uint32_t * D)250 void LinearTransform::reduce(int32_t* N, uint32_t* D) {
251 if (N && D && *D) {
252 if (*N < 0) {
253 *N = -(*N);
254 reduce(reinterpret_cast<uint32_t*>(N), D);
255 *N = -(*N);
256 } else {
257 reduce(reinterpret_cast<uint32_t*>(N), D);
258 }
259 }
260 }
261
262 } // namespace android
263