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