1 /*
2 * Copyright 2008 The Android Open Source Project
3 *
4 * Use of this source code is governed by a BSD-style license that can be
5 * found in the LICENSE file.
6 */
7
8 #include "SkMathPriv.h"
9 #include "SkCordic.h"
10 #include "SkFloatBits.h"
11 #include "SkFloatingPoint.h"
12 #include "Sk64.h"
13 #include "SkScalar.h"
14
15 #ifdef SK_SCALAR_IS_FLOAT
16 const uint32_t gIEEENotANumber = 0x7FFFFFFF;
17 const uint32_t gIEEEInfinity = 0x7F800000;
18 const uint32_t gIEEENegativeInfinity = 0xFF800000;
19 #endif
20
21 #define sub_shift(zeros, x, n) \
22 zeros -= n; \
23 x >>= n
24
SkCLZ_portable(uint32_t x)25 int SkCLZ_portable(uint32_t x) {
26 if (x == 0) {
27 return 32;
28 }
29
30 #ifdef SK_CPU_HAS_CONDITIONAL_INSTR
31 int zeros = 31;
32 if (x & 0xFFFF0000) {
33 sub_shift(zeros, x, 16);
34 }
35 if (x & 0xFF00) {
36 sub_shift(zeros, x, 8);
37 }
38 if (x & 0xF0) {
39 sub_shift(zeros, x, 4);
40 }
41 if (x & 0xC) {
42 sub_shift(zeros, x, 2);
43 }
44 if (x & 0x2) {
45 sub_shift(zeros, x, 1);
46 }
47 #else
48 int zeros = ((x >> 16) - 1) >> 31 << 4;
49 x <<= zeros;
50
51 int nonzero = ((x >> 24) - 1) >> 31 << 3;
52 zeros += nonzero;
53 x <<= nonzero;
54
55 nonzero = ((x >> 28) - 1) >> 31 << 2;
56 zeros += nonzero;
57 x <<= nonzero;
58
59 nonzero = ((x >> 30) - 1) >> 31 << 1;
60 zeros += nonzero;
61 x <<= nonzero;
62
63 zeros += (~x) >> 31;
64 #endif
65
66 return zeros;
67 }
68
SkMulDiv(int32_t numer1,int32_t numer2,int32_t denom)69 int32_t SkMulDiv(int32_t numer1, int32_t numer2, int32_t denom) {
70 SkASSERT(denom);
71
72 Sk64 tmp;
73 tmp.setMul(numer1, numer2);
74 tmp.div(denom, Sk64::kTrunc_DivOption);
75 return tmp.get32();
76 }
77
SkMulShift(int32_t a,int32_t b,unsigned shift)78 int32_t SkMulShift(int32_t a, int32_t b, unsigned shift) {
79 int sign = SkExtractSign(a ^ b);
80
81 if (shift > 63) {
82 return sign;
83 }
84
85 a = SkAbs32(a);
86 b = SkAbs32(b);
87
88 uint32_t ah = a >> 16;
89 uint32_t al = a & 0xFFFF;
90 uint32_t bh = b >> 16;
91 uint32_t bl = b & 0xFFFF;
92
93 uint32_t A = ah * bh;
94 uint32_t B = ah * bl + al * bh;
95 uint32_t C = al * bl;
96
97 /* [ A ]
98 [ B ]
99 [ C ]
100 */
101 uint32_t lo = C + (B << 16);
102 int32_t hi = A + (B >> 16) + (lo < C);
103
104 if (sign < 0) {
105 hi = -hi - Sk32ToBool(lo);
106 lo = 0 - lo;
107 }
108
109 if (shift == 0) {
110 #ifdef SK_DEBUGx
111 SkASSERT(((int32_t)lo >> 31) == hi);
112 #endif
113 return lo;
114 } else if (shift >= 32) {
115 return hi >> (shift - 32);
116 } else {
117 #ifdef SK_DEBUGx
118 int32_t tmp = hi >> shift;
119 SkASSERT(tmp == 0 || tmp == -1);
120 #endif
121 // we want (hi << (32 - shift)) | (lo >> shift) but rounded
122 int roundBit = (lo >> (shift - 1)) & 1;
123 return ((hi << (32 - shift)) | (lo >> shift)) + roundBit;
124 }
125 }
126
SkFixedMul_portable(SkFixed a,SkFixed b)127 SkFixed SkFixedMul_portable(SkFixed a, SkFixed b) {
128 #if 0
129 Sk64 tmp;
130
131 tmp.setMul(a, b);
132 tmp.shiftRight(16);
133 return tmp.fLo;
134 #elif defined(SkLONGLONG)
135 return static_cast<SkFixed>((SkLONGLONG)a * b >> 16);
136 #else
137 int sa = SkExtractSign(a);
138 int sb = SkExtractSign(b);
139 // now make them positive
140 a = SkApplySign(a, sa);
141 b = SkApplySign(b, sb);
142
143 uint32_t ah = a >> 16;
144 uint32_t al = a & 0xFFFF;
145 uint32_t bh = b >> 16;
146 uint32_t bl = b & 0xFFFF;
147
148 uint32_t R = ah * b + al * bh + (al * bl >> 16);
149
150 return SkApplySign(R, sa ^ sb);
151 #endif
152 }
153
SkFractMul_portable(SkFract a,SkFract b)154 SkFract SkFractMul_portable(SkFract a, SkFract b) {
155 #if 0
156 Sk64 tmp;
157 tmp.setMul(a, b);
158 return tmp.getFract();
159 #elif defined(SkLONGLONG)
160 return static_cast<SkFract>((SkLONGLONG)a * b >> 30);
161 #else
162 int sa = SkExtractSign(a);
163 int sb = SkExtractSign(b);
164 // now make them positive
165 a = SkApplySign(a, sa);
166 b = SkApplySign(b, sb);
167
168 uint32_t ah = a >> 16;
169 uint32_t al = a & 0xFFFF;
170 uint32_t bh = b >> 16;
171 uint32_t bl = b & 0xFFFF;
172
173 uint32_t A = ah * bh;
174 uint32_t B = ah * bl + al * bh;
175 uint32_t C = al * bl;
176
177 /* [ A ]
178 [ B ]
179 [ C ]
180 */
181 uint32_t Lo = C + (B << 16);
182 uint32_t Hi = A + (B >>16) + (Lo < C);
183
184 SkASSERT((Hi >> 29) == 0); // else overflow
185
186 int32_t R = (Hi << 2) + (Lo >> 30);
187
188 return SkApplySign(R, sa ^ sb);
189 #endif
190 }
191
SkFixedMulCommon(SkFixed a,int b,int bias)192 int SkFixedMulCommon(SkFixed a, int b, int bias) {
193 // this function only works if b is 16bits
194 SkASSERT(b == (int16_t)b);
195 SkASSERT(b >= 0);
196
197 int sa = SkExtractSign(a);
198 a = SkApplySign(a, sa);
199 uint32_t ah = a >> 16;
200 uint32_t al = a & 0xFFFF;
201 uint32_t R = ah * b + ((al * b + bias) >> 16);
202 return SkApplySign(R, sa);
203 }
204
205 #ifdef SK_DEBUGx
206 #define TEST_FASTINVERT
207 #endif
208
SkFixedFastInvert(SkFixed x)209 SkFixed SkFixedFastInvert(SkFixed x) {
210 /* Adapted (stolen) from gglRecip()
211 */
212
213 if (x == SK_Fixed1) {
214 return SK_Fixed1;
215 }
216
217 int sign = SkExtractSign(x);
218 uint32_t a = SkApplySign(x, sign);
219
220 if (a <= 2) {
221 return SkApplySign(SK_MaxS32, sign);
222 }
223
224 #ifdef TEST_FASTINVERT
225 SkFixed orig = a;
226 uint32_t slow = SkFixedDiv(SK_Fixed1, a);
227 #endif
228
229 // normalize a
230 int lz = SkCLZ(a);
231 a = a << lz >> 16;
232
233 // compute 1/a approximation (0.5 <= a < 1.0)
234 uint32_t r = 0x17400 - a; // (2.90625 (~2.914) - 2*a) >> 1
235
236 // Newton-Raphson iteration:
237 // x = r*(2 - a*r) = ((r/2)*(1 - a*r/2))*4
238 r = ( (0x10000 - ((a*r)>>16)) * r ) >> 15;
239 r = ( (0x10000 - ((a*r)>>16)) * r ) >> (30 - lz);
240
241 #ifdef TEST_FASTINVERT
242 SkDebugf("SkFixedFastInvert(%x %g) = %x %g Slow[%x %g]\n",
243 orig, orig/65536.,
244 r, r/65536.,
245 slow, slow/65536.);
246 #endif
247
248 return SkApplySign(r, sign);
249 }
250
251 ///////////////////////////////////////////////////////////////////////////////
252
253 #define DIVBITS_ITER(n) \
254 case n: \
255 if ((numer = (numer << 1) - denom) >= 0) \
256 result |= 1 << (n - 1); else numer += denom
257
SkDivBits(int32_t numer,int32_t denom,int shift_bias)258 int32_t SkDivBits(int32_t numer, int32_t denom, int shift_bias) {
259 SkASSERT(denom != 0);
260 if (numer == 0) {
261 return 0;
262 }
263
264 // make numer and denom positive, and sign hold the resulting sign
265 int32_t sign = SkExtractSign(numer ^ denom);
266 numer = SkAbs32(numer);
267 denom = SkAbs32(denom);
268
269 int nbits = SkCLZ(numer) - 1;
270 int dbits = SkCLZ(denom) - 1;
271 int bits = shift_bias - nbits + dbits;
272
273 if (bits < 0) { // answer will underflow
274 return 0;
275 }
276 if (bits > 31) { // answer will overflow
277 return SkApplySign(SK_MaxS32, sign);
278 }
279
280 denom <<= dbits;
281 numer <<= nbits;
282
283 SkFixed result = 0;
284
285 // do the first one
286 if ((numer -= denom) >= 0) {
287 result = 1;
288 } else {
289 numer += denom;
290 }
291
292 // Now fall into our switch statement if there are more bits to compute
293 if (bits > 0) {
294 // make room for the rest of the answer bits
295 result <<= bits;
296 switch (bits) {
297 DIVBITS_ITER(31); DIVBITS_ITER(30); DIVBITS_ITER(29);
298 DIVBITS_ITER(28); DIVBITS_ITER(27); DIVBITS_ITER(26);
299 DIVBITS_ITER(25); DIVBITS_ITER(24); DIVBITS_ITER(23);
300 DIVBITS_ITER(22); DIVBITS_ITER(21); DIVBITS_ITER(20);
301 DIVBITS_ITER(19); DIVBITS_ITER(18); DIVBITS_ITER(17);
302 DIVBITS_ITER(16); DIVBITS_ITER(15); DIVBITS_ITER(14);
303 DIVBITS_ITER(13); DIVBITS_ITER(12); DIVBITS_ITER(11);
304 DIVBITS_ITER(10); DIVBITS_ITER( 9); DIVBITS_ITER( 8);
305 DIVBITS_ITER( 7); DIVBITS_ITER( 6); DIVBITS_ITER( 5);
306 DIVBITS_ITER( 4); DIVBITS_ITER( 3); DIVBITS_ITER( 2);
307 // we merge these last two together, makes GCC make better ARM
308 default:
309 DIVBITS_ITER( 1);
310 }
311 }
312
313 if (result < 0) {
314 result = SK_MaxS32;
315 }
316 return SkApplySign(result, sign);
317 }
318
319 /* mod(float numer, float denom) seems to always return the sign
320 of the numer, so that's what we do too
321 */
SkFixedMod(SkFixed numer,SkFixed denom)322 SkFixed SkFixedMod(SkFixed numer, SkFixed denom) {
323 int sn = SkExtractSign(numer);
324 int sd = SkExtractSign(denom);
325
326 numer = SkApplySign(numer, sn);
327 denom = SkApplySign(denom, sd);
328
329 if (numer < denom) {
330 return SkApplySign(numer, sn);
331 } else if (numer == denom) {
332 return 0;
333 } else {
334 SkFixed div = SkFixedDiv(numer, denom);
335 return SkApplySign(SkFixedMul(denom, div & 0xFFFF), sn);
336 }
337 }
338
339 /* www.worldserver.com/turk/computergraphics/FixedSqrt.pdf
340 */
SkSqrtBits(int32_t x,int count)341 int32_t SkSqrtBits(int32_t x, int count) {
342 SkASSERT(x >= 0 && count > 0 && (unsigned)count <= 30);
343
344 uint32_t root = 0;
345 uint32_t remHi = 0;
346 uint32_t remLo = x;
347
348 do {
349 root <<= 1;
350
351 remHi = (remHi<<2) | (remLo>>30);
352 remLo <<= 2;
353
354 uint32_t testDiv = (root << 1) + 1;
355 if (remHi >= testDiv) {
356 remHi -= testDiv;
357 root++;
358 }
359 } while (--count >= 0);
360
361 return root;
362 }
363
SkCubeRootBits(int32_t value,int bits)364 int32_t SkCubeRootBits(int32_t value, int bits) {
365 SkASSERT(bits > 0);
366
367 int sign = SkExtractSign(value);
368 value = SkApplySign(value, sign);
369
370 uint32_t root = 0;
371 uint32_t curr = (uint32_t)value >> 30;
372 value <<= 2;
373
374 do {
375 root <<= 1;
376 uint32_t guess = root * root + root;
377 guess = (guess << 1) + guess; // guess *= 3
378 if (guess < curr) {
379 curr -= guess + 1;
380 root |= 1;
381 }
382 curr = (curr << 3) | ((uint32_t)value >> 29);
383 value <<= 3;
384 } while (--bits);
385
386 return SkApplySign(root, sign);
387 }
388
SkFixedMean(SkFixed a,SkFixed b)389 SkFixed SkFixedMean(SkFixed a, SkFixed b) {
390 Sk64 tmp;
391
392 tmp.setMul(a, b);
393 return tmp.getSqrt();
394 }
395
396 ///////////////////////////////////////////////////////////////////////////////
397
398 #ifdef SK_SCALAR_IS_FLOAT
SkScalarSinCos(float radians,float * cosValue)399 float SkScalarSinCos(float radians, float* cosValue) {
400 float sinValue = sk_float_sin(radians);
401
402 if (cosValue) {
403 *cosValue = sk_float_cos(radians);
404 if (SkScalarNearlyZero(*cosValue)) {
405 *cosValue = 0;
406 }
407 }
408
409 if (SkScalarNearlyZero(sinValue)) {
410 sinValue = 0;
411 }
412 return sinValue;
413 }
414 #endif
415
416 #define INTERP_SINTABLE
417 #define BUILD_TABLE_AT_RUNTIMEx
418
419 #define kTableSize 256
420
421 #ifdef BUILD_TABLE_AT_RUNTIME
422 static uint16_t gSkSinTable[kTableSize];
423
build_sintable(uint16_t table[])424 static void build_sintable(uint16_t table[]) {
425 for (int i = 0; i < kTableSize; i++) {
426 double rad = i * 3.141592653589793 / (2*kTableSize);
427 double val = sin(rad);
428 int ival = (int)(val * SK_Fixed1);
429 table[i] = SkToU16(ival);
430 }
431 }
432 #else
433 #include "SkSinTable.h"
434 #endif
435
436 #define SK_Fract1024SizeOver2PI 0x28BE60 /* floatToFract(1024 / 2PI) */
437
438 #ifdef INTERP_SINTABLE
interp_table(const uint16_t table[],int index,int partial255)439 static SkFixed interp_table(const uint16_t table[], int index, int partial255) {
440 SkASSERT((unsigned)index < kTableSize);
441 SkASSERT((unsigned)partial255 <= 255);
442
443 SkFixed lower = table[index];
444 SkFixed upper = (index == kTableSize - 1) ? SK_Fixed1 : table[index + 1];
445
446 SkASSERT(lower < upper);
447 SkASSERT(lower >= 0);
448 SkASSERT(upper <= SK_Fixed1);
449
450 partial255 += (partial255 >> 7);
451 return lower + ((upper - lower) * partial255 >> 8);
452 }
453 #endif
454
SkFixedSinCos(SkFixed radians,SkFixed * cosValuePtr)455 SkFixed SkFixedSinCos(SkFixed radians, SkFixed* cosValuePtr) {
456 SkASSERT(SK_ARRAY_COUNT(gSkSinTable) == kTableSize);
457
458 #ifdef BUILD_TABLE_AT_RUNTIME
459 static bool gFirstTime = true;
460 if (gFirstTime) {
461 build_sintable(gSinTable);
462 gFirstTime = false;
463 }
464 #endif
465
466 // make radians positive
467 SkFixed sinValue, cosValue;
468 int32_t cosSign = 0;
469 int32_t sinSign = SkExtractSign(radians);
470 radians = SkApplySign(radians, sinSign);
471 // scale it to 0...1023 ...
472
473 #ifdef INTERP_SINTABLE
474 radians = SkMulDiv(radians, 2 * kTableSize * 256, SK_FixedPI);
475 int findex = radians & (kTableSize * 256 - 1);
476 int index = findex >> 8;
477 int partial = findex & 255;
478 sinValue = interp_table(gSkSinTable, index, partial);
479
480 findex = kTableSize * 256 - findex - 1;
481 index = findex >> 8;
482 partial = findex & 255;
483 cosValue = interp_table(gSkSinTable, index, partial);
484
485 int quad = ((unsigned)radians / (kTableSize * 256)) & 3;
486 #else
487 radians = SkMulDiv(radians, 2 * kTableSize, SK_FixedPI);
488 int index = radians & (kTableSize - 1);
489
490 if (index == 0) {
491 sinValue = 0;
492 cosValue = SK_Fixed1;
493 } else {
494 sinValue = gSkSinTable[index];
495 cosValue = gSkSinTable[kTableSize - index];
496 }
497 int quad = ((unsigned)radians / kTableSize) & 3;
498 #endif
499
500 if (quad & 1) {
501 SkTSwap<SkFixed>(sinValue, cosValue);
502 }
503 if (quad & 2) {
504 sinSign = ~sinSign;
505 }
506 if (((quad - 1) & 2) == 0) {
507 cosSign = ~cosSign;
508 }
509
510 // restore the sign for negative angles
511 sinValue = SkApplySign(sinValue, sinSign);
512 cosValue = SkApplySign(cosValue, cosSign);
513
514 #ifdef SK_DEBUG
515 if (1) {
516 SkFixed sin2 = SkFixedMul(sinValue, sinValue);
517 SkFixed cos2 = SkFixedMul(cosValue, cosValue);
518 int diff = cos2 + sin2 - SK_Fixed1;
519 SkASSERT(SkAbs32(diff) <= 7);
520 }
521 #endif
522
523 if (cosValuePtr) {
524 *cosValuePtr = cosValue;
525 }
526 return sinValue;
527 }
528
529 ///////////////////////////////////////////////////////////////////////////////
530
SkFixedTan(SkFixed radians)531 SkFixed SkFixedTan(SkFixed radians) { return SkCordicTan(radians); }
SkFixedASin(SkFixed x)532 SkFixed SkFixedASin(SkFixed x) { return SkCordicASin(x); }
SkFixedACos(SkFixed x)533 SkFixed SkFixedACos(SkFixed x) { return SkCordicACos(x); }
SkFixedATan2(SkFixed y,SkFixed x)534 SkFixed SkFixedATan2(SkFixed y, SkFixed x) { return SkCordicATan2(y, x); }
SkFixedExp(SkFixed x)535 SkFixed SkFixedExp(SkFixed x) { return SkCordicExp(x); }
SkFixedLog(SkFixed x)536 SkFixed SkFixedLog(SkFixed x) { return SkCordicLog(x); }
537