libcorecg/SkCordic.cpp - platform/external/skia - Git at Google

 /* libs/corecg/SkCordic.cpp
 **
 ** Copyright 2006, The Android Open Source Project
 **
 ** Licensed under the Apache License, Version 2.0 (the "License");
 ** you may not use this file except in compliance with the License.
 ** You may obtain a copy of the License at
 **
 **     http://www.apache.org/licenses/LICENSE-2.0
 **
 ** Unless required by applicable law or agreed to in writing, software
 ** distributed under the License is distributed on an "AS IS" BASIS,
 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 ** See the License for the specific language governing permissions and
 ** limitations under the License.
 */

 #include "SkCordic.h"
 #include "SkMath.h"
 #include "Sk64.h"

 // 0x20000000 equals pi / 4
 const int32_t kATanDegrees[] = { 0x20000000,
     0x12E4051D, 0x9FB385B, 0x51111D4, 0x28B0D43, 0x145D7E1, 0xA2F61E, 0x517C55,
     0x28BE53, 0x145F2E, 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
     0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
     0xA, 0x5, 0x2, 0x1 };

 const int32_t kFixedInvGain1 = 0x18bde0bb;  // 0.607252935

 static void SkCircularRotation(int32_t* x0, int32_t* y0, int32_t* z0)
 {
     int32_t t = 0;
     int32_t x = *x0;
     int32_t y = *y0;
     int32_t z = *z0;
     const int32_t* tanPtr = kATanDegrees;
    do {
         int32_t x1 = y >> t;
         int32_t y1 = x >> t;
         int32_t tan = *tanPtr++;
         if (z >= 0) {
             x -= x1;
             y += y1;
             z -= tan;
         } else {
             x += x1;
             y -= y1;
             z += tan;
         }
    } while (++t < 16); // 30);
     *x0 = x;
     *y0 = y;
     *z0 = z;
 }

 SkFixed SkCordicSinCos(SkFixed radians, SkFixed* cosp)
 {
     int32_t scaledRadians = radians * 0x28be;   // scale radians to 65536 / PI()
     int quadrant = scaledRadians >> 30;
     quadrant += 1;
     if (quadrant & 2)
         scaledRadians = -scaledRadians + 0x80000000;
     /* |a| <= 90 degrees as a 1.31 number */
     SkFixed sin = 0;
     SkFixed cos = kFixedInvGain1;
     SkCircularRotation(&cos, &sin, &scaledRadians);
     Sk64 scaled;
     scaled.setMul(sin, 0x6488d);
     sin = scaled.fHi;
     scaled.setMul(cos, 0x6488d);
     if (quadrant & 2)
         scaled.fHi = - scaled.fHi;
     *cosp = scaled.fHi;
     return sin;
 }

 SkFixed SkCordicTan(SkFixed a)
 {
     int32_t cos;
     int32_t sin = SkCordicSinCos(a, &cos);
     return SkFixedDiv(sin, cos);
 }

 static int32_t SkCircularVector(int32_t* y0, int32_t* x0, int32_t vecMode)
 {
     int32_t x = *x0;
     int32_t y = *y0;
     int32_t z = 0;
     int32_t t = 0;
     const int32_t* tanPtr = kATanDegrees;
    do {
         int32_t x1 = y >> t;
         int32_t y1 = x >> t;
         int32_t tan = *tanPtr++;
         if (y < vecMode) {
             x -= x1;
             y += y1;
             z -= tan;
         } else {
             x += x1;
             y -= y1;
             z += tan;
         }
    } while (++t < 16);  // 30
     Sk64 scaled;
     scaled.setMul(z, 0x6488d); // scale back into the SkScalar space (0x100000000/0x28be)
    return scaled.fHi;
 }

 SkFixed SkCordicASin(SkFixed a) {
     int32_t sign = SkExtractSign(a);
     int32_t z = SkFixedAbs(a);
     if (z >= SK_Fixed1)
         return SkApplySign(SK_FixedPI>>1, sign);
     int32_t x = kFixedInvGain1;
     int32_t y = 0;
     z *= 0x28be;
     z = SkCircularVector(&y, &x, z);
     z = SkApplySign(z, ~sign);
     return z;
 }

 SkFixed SkCordicACos(SkFixed a) {
     int32_t z = SkCordicASin(a);
     z = (SK_FixedPI>>1) - z;
     return z;
 }

 SkFixed SkCordicATan2(SkFixed y, SkFixed x) {
     if ((x | y) == 0)
         return 0;
     int32_t xsign = SkExtractSign(x);
     x = SkFixedAbs(x);
     int32_t result = SkCircularVector(&y, &x, 0);
     if (xsign) {
         int32_t rsign = SkExtractSign(result);
         if (y == 0)
             rsign = 0;
         SkFixed pi = SkApplySign(SK_FixedPI, rsign);
         result = pi - result;
     }
     return result;
 }

 const int32_t kATanHDegrees[] = {
     0x1661788D, 0xA680D61, 0x51EA6FC, 0x28CBFDD, 0x1460E34,
     0xA2FCE8, 0x517D2E, 0x28BE6E, 0x145F32,
     0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
     0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
     0xA, 0x5, 0x2, 0x1 };

 const int32_t kFixedInvGain2 = 0x31330AAA;  // 1.207534495

 static void SkHyperbolic(int32_t* x0, int32_t* y0, int32_t* z0, int mode)
 {
     int32_t t = 1;
     int32_t x = *x0;
     int32_t y = *y0;
     int32_t z = *z0;
     const int32_t* tanPtr = kATanHDegrees;
     int k = -3;
     do {
         int32_t x1 = y >> t;
         int32_t y1 = x >> t;
         int32_t tan = *tanPtr++;
         int count = 2 + (k >> 31);
         if (++k == 1)
             k = -2;
         do {
             if (((y >> 31) & mode) | ~((z >> 31) | mode)) {
                 x += x1;
                 y += y1;
                 z -= tan;
             } else {
                 x -= x1;
                 y -= y1;
                 z += tan;
             }
         } while (--count);
     } while (++t < 30);
     *x0 = x;
     *y0 = y;
     *z0 = z;
 }

 SkFixed SkCordicLog(SkFixed a) {
     a *= 0x28be;
     int32_t x = a + 0x28BE60DB; // 1.0
     int32_t y = a - 0x28BE60DB;
     int32_t z = 0;
     SkHyperbolic(&x, &y, &z, -1);
     Sk64 scaled;
     scaled.setMul(z, 0x6488d);
     z = scaled.fHi;
     return z << 1;
 }

 SkFixed SkCordicExp(SkFixed a) {
     int32_t cosh = kFixedInvGain2;
     int32_t sinh = 0;
     SkHyperbolic(&cosh, &sinh, &a, 0);
     return cosh + sinh;
 }

 #ifdef SK_DEBUG

 #ifdef SK_CAN_USE_FLOAT
     #include "SkFloatingPoint.h"
 #endif

 void SkCordic_UnitTest()
 {
 #if defined(SK_SUPPORT_UNITTEST) && defined(SK_CAN_USE_FLOAT)
     float val;
     for (float angle = -720; angle < 720; angle += 30) {
         float radian = angle * 3.1415925358f / 180.0f;
         SkFixed f_angle = (int) (radian * 65536.0f);
     // sincos
         float sine = sinf(radian);
         float cosine = cosf(radian);
         SkFixed f_cosine;
         SkFixed f_sine = SkCordicSinCos(f_angle, &f_cosine);
         float sine2 = (float) f_sine / 65536.0f;
         float cosine2 = (float) f_cosine / 65536.0f;
         float error = fabsf(sine - sine2);
         if (error > 0.001)
             SkDebugf("sin error : angle = %g ; sin = %g ; cordic = %g\n", angle, sine, sine2);
         error = fabsf(cosine - cosine2);
         if (error > 0.001)
             SkDebugf("cos error : angle = %g ; cos = %g ; cordic = %g\n", angle, cosine, cosine2);
     // tan
         float _tan = tanf(radian);
         SkFixed f_tan = SkCordicTan(f_angle);
         float tan2 = (float) f_tan / 65536.0f;
         error = fabsf(_tan - tan2);
         if (error > 0.05 && fabsf(_tan) < 1e6)
             SkDebugf("tan error : angle = %g ; tan = %g ; cordic = %g\n", angle, _tan, tan2);
     }
     for (val = -1; val <= 1; val += .1f) {
         SkFixed f_val = (int) (val * 65536.0f);
     // asin
         float arcsine = asinf(val);
         SkFixed f_arcsine = SkCordicASin(f_val);
         float arcsine2 = (float) f_arcsine / 65536.0f;
         float error = fabsf(arcsine - arcsine2);
         if (error > 0.001)
             SkDebugf("asin error : val = %g ; asin = %g ; cordic = %g\n", val, arcsine, arcsine2);
     }
 #if 1
     for (val = -1; val <= 1; val += .1f) {
 #else
     val = .5; {
 #endif
         SkFixed f_val = (int) (val * 65536.0f);
     // acos
         float arccos = acosf(val);
         SkFixed f_arccos = SkCordicACos(f_val);
         float arccos2 = (float) f_arccos / 65536.0f;
         float error = fabsf(arccos - arccos2);
         if (error > 0.001)
             SkDebugf("acos error : val = %g ; acos = %g ; cordic = %g\n", val, arccos, arccos2);
     }
     // atan2
 #if 1
     for (val = -1000; val <= 1000; val += 500.f) {
         for (float val2 = -1000; val2 <= 1000; val2 += 500.f) {
 #else
             val = 0; {
             float val2 = -1000; {
 #endif
             SkFixed f_val = (int) (val * 65536.0f);
             SkFixed f_val2 = (int) (val2 * 65536.0f);
             float arctan = atan2f(val, val2);
             SkFixed f_arctan = SkCordicATan2(f_val, f_val2);
             float arctan2 = (float) f_arctan / 65536.0f;
             float error = fabsf(arctan - arctan2);
             if (error > 0.001)
                 SkDebugf("atan2 error : val = %g ; val2 = %g ; atan2 = %g ; cordic = %g\n", val, val2, arctan, arctan2);
         }
     }
     // log
 #if 1
     for (val = 0.125f; val <= 8.f; val *= 2.0f) {
 #else
     val = .5; {
 #endif
         SkFixed f_val = (int) (val * 65536.0f);
     // acos
         float log = logf(val);
         SkFixed f_log = SkCordicLog(f_val);
         float log2 = (float) f_log / 65536.0f;
         float error = fabsf(log - log2);
         if (error > 0.001)
             SkDebugf("log error : val = %g ; log = %g ; cordic = %g\n", val, log, log2);
     }
     // exp
 #endif
 }

 #endif
	/* libs/corecg/SkCordic.cpp
	**
	** Copyright 2006, The Android Open Source Project
	**
	** Licensed under the Apache License, Version 2.0 (the "License");
	** you may not use this file except in compliance with the License.
	** You may obtain a copy of the License at
	**
	** http://www.apache.org/licenses/LICENSE-2.0
	**
	** Unless required by applicable law or agreed to in writing, software
	** distributed under the License is distributed on an "AS IS" BASIS,
	** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	** See the License for the specific language governing permissions and
	** limitations under the License.
	*/

	#include "SkCordic.h"
	#include "SkMath.h"
	#include "Sk64.h"

	// 0x20000000 equals pi / 4
	const int32_t kATanDegrees[] = { 0x20000000,
	0x12E4051D, 0x9FB385B, 0x51111D4, 0x28B0D43, 0x145D7E1, 0xA2F61E, 0x517C55,
	0x28BE53, 0x145F2E, 0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
	0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
	0xA, 0x5, 0x2, 0x1 };

	const int32_t kFixedInvGain1 = 0x18bde0bb; // 0.607252935

	static void SkCircularRotation(int32_t* x0, int32_t* y0, int32_t* z0)
	{
	int32_t t = 0;
	int32_t x = *x0;
	int32_t y = *y0;
	int32_t z = *z0;
	const int32_t* tanPtr = kATanDegrees;
	do {
	int32_t x1 = y >> t;
	int32_t y1 = x >> t;
	int32_t tan = *tanPtr++;
	if (z >= 0) {
	x -= x1;
	y += y1;
	z -= tan;
	} else {
	x += x1;
	y -= y1;
	z += tan;
	}
	} while (++t < 16); // 30);
	*x0 = x;
	*y0 = y;
	*z0 = z;
	}

	SkFixed SkCordicSinCos(SkFixed radians, SkFixed* cosp)
	{
	int32_t scaledRadians = radians * 0x28be; // scale radians to 65536 / PI()
	int quadrant = scaledRadians >> 30;
	quadrant += 1;
	if (quadrant & 2)
	scaledRadians = -scaledRadians + 0x80000000;
	/* \|a\| <= 90 degrees as a 1.31 number */
	SkFixed sin = 0;
	SkFixed cos = kFixedInvGain1;
	SkCircularRotation(&cos, &sin, &scaledRadians);
	Sk64 scaled;
	scaled.setMul(sin, 0x6488d);
	sin = scaled.fHi;
	scaled.setMul(cos, 0x6488d);
	if (quadrant & 2)
	scaled.fHi = - scaled.fHi;
	*cosp = scaled.fHi;
	return sin;
	}

	SkFixed SkCordicTan(SkFixed a)
	{
	int32_t cos;
	int32_t sin = SkCordicSinCos(a, &cos);
	return SkFixedDiv(sin, cos);
	}

	static int32_t SkCircularVector(int32_t* y0, int32_t* x0, int32_t vecMode)
	{
	int32_t x = *x0;
	int32_t y = *y0;
	int32_t z = 0;
	int32_t t = 0;
	const int32_t* tanPtr = kATanDegrees;
	do {
	int32_t x1 = y >> t;
	int32_t y1 = x >> t;
	int32_t tan = *tanPtr++;
	if (y < vecMode) {
	x -= x1;
	y += y1;
	z -= tan;
	} else {
	x += x1;
	y -= y1;
	z += tan;
	}
	} while (++t < 16); // 30
	Sk64 scaled;
	scaled.setMul(z, 0x6488d); // scale back into the SkScalar space (0x100000000/0x28be)
	return scaled.fHi;
	}

	SkFixed SkCordicASin(SkFixed a) {
	int32_t sign = SkExtractSign(a);
	int32_t z = SkFixedAbs(a);
	if (z >= SK_Fixed1)
	return SkApplySign(SK_FixedPI>>1, sign);
	int32_t x = kFixedInvGain1;
	int32_t y = 0;
	z *= 0x28be;
	z = SkCircularVector(&y, &x, z);
	z = SkApplySign(z, ~sign);
	return z;
	}

	SkFixed SkCordicACos(SkFixed a) {
	int32_t z = SkCordicASin(a);
	z = (SK_FixedPI>>1) - z;
	return z;
	}

	SkFixed SkCordicATan2(SkFixed y, SkFixed x) {
	if ((x \| y) == 0)
	return 0;
	int32_t xsign = SkExtractSign(x);
	x = SkFixedAbs(x);
	int32_t result = SkCircularVector(&y, &x, 0);
	if (xsign) {
	int32_t rsign = SkExtractSign(result);
	if (y == 0)
	rsign = 0;
	SkFixed pi = SkApplySign(SK_FixedPI, rsign);
	result = pi - result;
	}
	return result;
	}

	const int32_t kATanHDegrees[] = {
	0x1661788D, 0xA680D61, 0x51EA6FC, 0x28CBFDD, 0x1460E34,
	0xA2FCE8, 0x517D2E, 0x28BE6E, 0x145F32,
	0xA2F98, 0x517CC, 0x28BE6, 0x145F3, 0xA2F9, 0x517C,
	0x28BE, 0x145F, 0xA2F, 0x517, 0x28B, 0x145, 0xA2, 0x51, 0x28, 0x14,
	0xA, 0x5, 0x2, 0x1 };

	const int32_t kFixedInvGain2 = 0x31330AAA; // 1.207534495

	static void SkHyperbolic(int32_t* x0, int32_t* y0, int32_t* z0, int mode)
	{
	int32_t t = 1;
	int32_t x = *x0;
	int32_t y = *y0;
	int32_t z = *z0;
	const int32_t* tanPtr = kATanHDegrees;
	int k = -3;
	do {
	int32_t x1 = y >> t;
	int32_t y1 = x >> t;
	int32_t tan = *tanPtr++;
	int count = 2 + (k >> 31);
	if (++k == 1)
	k = -2;
	do {
	if (((y >> 31) & mode) \| ~((z >> 31) \| mode)) {
	x += x1;
	y += y1;
	z -= tan;
	} else {
	x -= x1;
	y -= y1;
	z += tan;
	}
	} while (--count);
	} while (++t < 30);
	*x0 = x;
	*y0 = y;
	*z0 = z;
	}

	SkFixed SkCordicLog(SkFixed a) {
	a *= 0x28be;
	int32_t x = a + 0x28BE60DB; // 1.0
	int32_t y = a - 0x28BE60DB;
	int32_t z = 0;
	SkHyperbolic(&x, &y, &z, -1);
	Sk64 scaled;
	scaled.setMul(z, 0x6488d);
	z = scaled.fHi;
	return z << 1;
	}

	SkFixed SkCordicExp(SkFixed a) {
	int32_t cosh = kFixedInvGain2;
	int32_t sinh = 0;
	SkHyperbolic(&cosh, &sinh, &a, 0);
	return cosh + sinh;
	}

	#ifdef SK_DEBUG

	#ifdef SK_CAN_USE_FLOAT
	#include "SkFloatingPoint.h"
	#endif

	void SkCordic_UnitTest()
	{
	#if defined(SK_SUPPORT_UNITTEST) && defined(SK_CAN_USE_FLOAT)
	float val;
	for (float angle = -720; angle < 720; angle += 30) {
	float radian = angle * 3.1415925358f / 180.0f;
	SkFixed f_angle = (int) (radian * 65536.0f);
	// sincos
	float sine = sinf(radian);
	float cosine = cosf(radian);
	SkFixed f_cosine;
	SkFixed f_sine = SkCordicSinCos(f_angle, &f_cosine);
	float sine2 = (float) f_sine / 65536.0f;
	float cosine2 = (float) f_cosine / 65536.0f;
	float error = fabsf(sine - sine2);
	if (error > 0.001)
	SkDebugf("sin error : angle = %g ; sin = %g ; cordic = %g\n", angle, sine, sine2);
	error = fabsf(cosine - cosine2);
	if (error > 0.001)
	SkDebugf("cos error : angle = %g ; cos = %g ; cordic = %g\n", angle, cosine, cosine2);
	// tan
	float _tan = tanf(radian);
	SkFixed f_tan = SkCordicTan(f_angle);
	float tan2 = (float) f_tan / 65536.0f;
	error = fabsf(_tan - tan2);
	if (error > 0.05 && fabsf(_tan) < 1e6)
	SkDebugf("tan error : angle = %g ; tan = %g ; cordic = %g\n", angle, _tan, tan2);
	}
	for (val = -1; val <= 1; val += .1f) {
	SkFixed f_val = (int) (val * 65536.0f);
	// asin
	float arcsine = asinf(val);
	SkFixed f_arcsine = SkCordicASin(f_val);
	float arcsine2 = (float) f_arcsine / 65536.0f;
	float error = fabsf(arcsine - arcsine2);
	if (error > 0.001)
	SkDebugf("asin error : val = %g ; asin = %g ; cordic = %g\n", val, arcsine, arcsine2);
	}
	#if 1
	for (val = -1; val <= 1; val += .1f) {
	#else
	val = .5; {
	#endif
	SkFixed f_val = (int) (val * 65536.0f);
	// acos
	float arccos = acosf(val);
	SkFixed f_arccos = SkCordicACos(f_val);
	float arccos2 = (float) f_arccos / 65536.0f;
	float error = fabsf(arccos - arccos2);
	if (error > 0.001)
	SkDebugf("acos error : val = %g ; acos = %g ; cordic = %g\n", val, arccos, arccos2);
	}
	// atan2
	#if 1
	for (val = -1000; val <= 1000; val += 500.f) {
	for (float val2 = -1000; val2 <= 1000; val2 += 500.f) {
	#else
	val = 0; {
	float val2 = -1000; {
	#endif
	SkFixed f_val = (int) (val * 65536.0f);
	SkFixed f_val2 = (int) (val2 * 65536.0f);
	float arctan = atan2f(val, val2);
	SkFixed f_arctan = SkCordicATan2(f_val, f_val2);
	float arctan2 = (float) f_arctan / 65536.0f;
	float error = fabsf(arctan - arctan2);
	if (error > 0.001)
	SkDebugf("atan2 error : val = %g ; val2 = %g ; atan2 = %g ; cordic = %g\n", val, val2, arctan, arctan2);
	}
	}
	// log
	#if 1
	for (val = 0.125f; val <= 8.f; val *= 2.0f) {
	#else
	val = .5; {
	#endif
	SkFixed f_val = (int) (val * 65536.0f);
	// acos
	float log = logf(val);
	SkFixed f_log = SkCordicLog(f_val);
	float log2 = (float) f_log / 65536.0f;
	float error = fabsf(log - log2);
	if (error > 0.001)
	SkDebugf("log error : val = %g ; log = %g ; cordic = %g\n", val, log, log2);
	}
	// exp
	#endif
	}

	#endif