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