1 2 /* 3 * Copyright 2006 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 "SkCordic.h" 11 #include "SkMath.h" 12 #include "Sk64.h" 13 14 // 0x20000000 equals pi / 4 15 const int32_t kATanDegrees[] = { 0x20000000, 16 0x12E4051D, 0x9FB385B, 0x51111D4, 0x28B0D43, 0x145D7E1, 0xA2F61E, 0x517C55, 17 0x28BE53, 0x145F2E, 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C, 18 0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14, 19 0xA, 0x5, 0x2, 0x1 }; 20 21 const int32_t kFixedInvGain1 = 0x18bde0bb; // 0.607252935 22 23 static void SkCircularRotation(int32_t* x0, int32_t* y0, int32_t* z0) 24 { 25 int32_t t = 0; 26 int32_t x = *x0; 27 int32_t y = *y0; 28 int32_t z = *z0; 29 const int32_t* tanPtr = kATanDegrees; 30 do { 31 int32_t x1 = y >> t; 32 int32_t y1 = x >> t; 33 int32_t tan = *tanPtr++; 34 if (z >= 0) { 35 x -= x1; 36 y += y1; 37 z -= tan; 38 } else { 39 x += x1; 40 y -= y1; 41 z += tan; 42 } 43 } while (++t < 16); // 30); 44 *x0 = x; 45 *y0 = y; 46 *z0 = z; 47 } 48 49 SkFixed SkCordicSinCos(SkFixed radians, SkFixed* cosp) 50 { 51 int32_t scaledRadians = radians * 0x28be; // scale radians to 65536 / PI() 52 int quadrant = scaledRadians >> 30; 53 quadrant += 1; 54 if (quadrant & 2) 55 scaledRadians = -scaledRadians + 0x80000000; 56 /* |a| <= 90 degrees as a 1.31 number */ 57 SkFixed sin = 0; 58 SkFixed cos = kFixedInvGain1; 59 SkCircularRotation(&cos, &sin, &scaledRadians); 60 Sk64 scaled; 61 scaled.setMul(sin, 0x6488d); 62 sin = scaled.fHi; 63 scaled.setMul(cos, 0x6488d); 64 if (quadrant & 2) 65 scaled.fHi = - scaled.fHi; 66 *cosp = scaled.fHi; 67 return sin; 68 } 69 70 SkFixed SkCordicTan(SkFixed a) 71 { 72 int32_t cos; 73 int32_t sin = SkCordicSinCos(a, &cos); 74 return SkFixedDiv(sin, cos); 75 } 76 77 static int32_t SkCircularVector(int32_t* y0, int32_t* x0, int32_t vecMode) 78 { 79 int32_t x = *x0; 80 int32_t y = *y0; 81 int32_t z = 0; 82 int32_t t = 0; 83 const int32_t* tanPtr = kATanDegrees; 84 do { 85 int32_t x1 = y >> t; 86 int32_t y1 = x >> t; 87 int32_t tan = *tanPtr++; 88 if (y < vecMode) { 89 x -= x1; 90 y += y1; 91 z -= tan; 92 } else { 93 x += x1; 94 y -= y1; 95 z += tan; 96 } 97 } while (++t < 16); // 30 98 Sk64 scaled; 99 scaled.setMul(z, 0x6488d); // scale back into the SkScalar space (0x100000000/0x28be) 100 return scaled.fHi; 101 } 102 103 SkFixed SkCordicASin(SkFixed a) { 104 int32_t sign = SkExtractSign(a); 105 int32_t z = SkFixedAbs(a); 106 if (z >= SK_Fixed1) 107 return SkApplySign(SK_FixedPI>>1, sign); 108 int32_t x = kFixedInvGain1; 109 int32_t y = 0; 110 z *= 0x28be; 111 z = SkCircularVector(&y, &x, z); 112 z = SkApplySign(z, ~sign); 113 return z; 114 } 115 116 SkFixed SkCordicACos(SkFixed a) { 117 int32_t z = SkCordicASin(a); 118 z = (SK_FixedPI>>1) - z; 119 return z; 120 } 121 122 SkFixed SkCordicATan2(SkFixed y, SkFixed x) { 123 if ((x | y) == 0) 124 return 0; 125 int32_t xsign = SkExtractSign(x); 126 x = SkFixedAbs(x); 127 int32_t result = SkCircularVector(&y, &x, 0); 128 if (xsign) { 129 int32_t rsign = SkExtractSign(result); 130 if (y == 0) 131 rsign = 0; 132 SkFixed pi = SkApplySign(SK_FixedPI, rsign); 133 result = pi - result; 134 } 135 return result; 136 } 137 138 const int32_t kATanHDegrees[] = { 139 0x1661788D, 0xA680D61, 0x51EA6FC, 0x28CBFDD, 0x1460E34, 140 0xA2FCE8, 0x517D2E, 0x28BE6E, 0x145F32, 141 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C, 142 0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14, 143 0xA, 0x5, 0x2, 0x1 }; 144 145 const int32_t kFixedInvGain2 = 0x31330AAA; // 1.207534495 146 147 static void SkHyperbolic(int32_t* x0, int32_t* y0, int32_t* z0, int mode) 148 { 149 int32_t t = 1; 150 int32_t x = *x0; 151 int32_t y = *y0; 152 int32_t z = *z0; 153 const int32_t* tanPtr = kATanHDegrees; 154 int k = -3; 155 do { 156 int32_t x1 = y >> t; 157 int32_t y1 = x >> t; 158 int32_t tan = *tanPtr++; 159 int count = 2 + (k >> 31); 160 if (++k == 1) 161 k = -2; 162 do { 163 if (((y >> 31) & mode) | ~((z >> 31) | mode)) { 164 x += x1; 165 y += y1; 166 z -= tan; 167 } else { 168 x -= x1; 169 y -= y1; 170 z += tan; 171 } 172 } while (--count); 173 } while (++t < 30); 174 *x0 = x; 175 *y0 = y; 176 *z0 = z; 177 } 178 179 SkFixed SkCordicLog(SkFixed a) { 180 a *= 0x28be; 181 int32_t x = a + 0x28BE60DB; // 1.0 182 int32_t y = a - 0x28BE60DB; 183 int32_t z = 0; 184 SkHyperbolic(&x, &y, &z, -1); 185 Sk64 scaled; 186 scaled.setMul(z, 0x6488d); 187 z = scaled.fHi; 188 return z << 1; 189 } 190 191 SkFixed SkCordicExp(SkFixed a) { 192 int32_t cosh = kFixedInvGain2; 193 int32_t sinh = 0; 194 SkHyperbolic(&cosh, &sinh, &a, 0); 195 return cosh + sinh; 196 } 197 198 #ifdef SK_DEBUG 199 200 #ifdef SK_CAN_USE_FLOAT 201 #include "SkFloatingPoint.h" 202 #endif 203 204 void SkCordic_UnitTest() 205 { 206 #if defined(SK_SUPPORT_UNITTEST) && defined(SK_CAN_USE_FLOAT) 207 float val; 208 for (float angle = -720; angle < 720; angle += 30) { 209 float radian = angle * 3.1415925358f / 180.0f; 210 SkFixed f_angle = (int) (radian * 65536.0f); 211 // sincos 212 float sine = sinf(radian); 213 float cosine = cosf(radian); 214 SkFixed f_cosine; 215 SkFixed f_sine = SkCordicSinCos(f_angle, &f_cosine); 216 float sine2 = (float) f_sine / 65536.0f; 217 float cosine2 = (float) f_cosine / 65536.0f; 218 float error = fabsf(sine - sine2); 219 if (error > 0.001) 220 SkDebugf("sin error : angle = %g ; sin = %g ; cordic = %g\n", angle, sine, sine2); 221 error = fabsf(cosine - cosine2); 222 if (error > 0.001) 223 SkDebugf("cos error : angle = %g ; cos = %g ; cordic = %g\n", angle, cosine, cosine2); 224 // tan 225 float _tan = tanf(radian); 226 SkFixed f_tan = SkCordicTan(f_angle); 227 float tan2 = (float) f_tan / 65536.0f; 228 error = fabsf(_tan - tan2); 229 if (error > 0.05 && fabsf(_tan) < 1e6) 230 SkDebugf("tan error : angle = %g ; tan = %g ; cordic = %g\n", angle, _tan, tan2); 231 } 232 for (val = -1; val <= 1; val += .1f) { 233 SkFixed f_val = (int) (val * 65536.0f); 234 // asin 235 float arcsine = asinf(val); 236 SkFixed f_arcsine = SkCordicASin(f_val); 237 float arcsine2 = (float) f_arcsine / 65536.0f; 238 float error = fabsf(arcsine - arcsine2); 239 if (error > 0.001) 240 SkDebugf("asin error : val = %g ; asin = %g ; cordic = %g\n", val, arcsine, arcsine2); 241 } 242 #if 1 243 for (val = -1; val <= 1; val += .1f) { 244 #else 245 val = .5; { 246 #endif 247 SkFixed f_val = (int) (val * 65536.0f); 248 // acos 249 float arccos = acosf(val); 250 SkFixed f_arccos = SkCordicACos(f_val); 251 float arccos2 = (float) f_arccos / 65536.0f; 252 float error = fabsf(arccos - arccos2); 253 if (error > 0.001) 254 SkDebugf("acos error : val = %g ; acos = %g ; cordic = %g\n", val, arccos, arccos2); 255 } 256 // atan2 257 #if 1 258 for (val = -1000; val <= 1000; val += 500.f) { 259 for (float val2 = -1000; val2 <= 1000; val2 += 500.f) { 260 #else 261 val = 0; { 262 float val2 = -1000; { 263 #endif 264 SkFixed f_val = (int) (val * 65536.0f); 265 SkFixed f_val2 = (int) (val2 * 65536.0f); 266 float arctan = atan2f(val, val2); 267 SkFixed f_arctan = SkCordicATan2(f_val, f_val2); 268 float arctan2 = (float) f_arctan / 65536.0f; 269 float error = fabsf(arctan - arctan2); 270 if (error > 0.001) 271 SkDebugf("atan2 error : val = %g ; val2 = %g ; atan2 = %g ; cordic = %g\n", val, val2, arctan, arctan2); 272 } 273 } 274 // log 275 #if 1 276 for (val = 0.125f; val <= 8.f; val *= 2.0f) { 277 #else 278 val = .5; { 279 #endif 280 SkFixed f_val = (int) (val * 65536.0f); 281 // acos 282 float log = logf(val); 283 SkFixed f_log = SkCordicLog(f_val); 284 float log2 = (float) f_log / 65536.0f; 285 float error = fabsf(log - log2); 286 if (error > 0.001) 287 SkDebugf("log error : val = %g ; log = %g ; cordic = %g\n", val, log, log2); 288 } 289 // exp 290 #endif 291 } 292 293 #endif 294
