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 8 #include "SkCordic.h" 9 #include "SkMathPriv.h" 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 198 #include "SkFloatingPoint.h" 199 200 void SkCordic_UnitTest() 201 { 202 #if defined(SK_SUPPORT_UNITTEST) 203 float val; 204 for (float angle = -720; angle < 720; angle += 30) { 205 float radian = angle * 3.1415925358f / 180.0f; 206 SkFixed f_angle = SkFloatToFixed(radian); 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 = SkFloatToFixed(val); 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 = SkFloatToFixed(val); 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 = SkFloatToFixed(val); 261 SkFixed f_val2 = SkFloatToFixed(val2); 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 = SkFloatToFixed(val); 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 290
