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