src/common/tcos.c

#ifndef lint
static const char       RCSid[] = "$Id: tcos.c,v 3.9 2013/02/08 16:39:01 greg Exp $";
#endif
/*
 * Table-based cosine approximation.
 *
 * Use doubles in table even though we're not nearly that accurate just
 * to avoid conversion and guarantee that tsin(x)^2 + tcos(x)^2 == 1.
 *
 * No interpolation in this version.
 *
 * External symbols declared in rtmath.h
 */

#include "copyright.h"

#include <math.h>

#include "rtmath.h"

#ifndef __FAST_MATH__

#ifndef NCOSENTRY
#define NCOSENTRY       1024
#endif

double
tcos(double x)                          /* approximate cosine */
{
        static double   costab[NCOSENTRY+1];
        register int    i;

        if (costab[0] < 0.5)            /* initialize table */
                for (i = 0; i <= NCOSENTRY; i++)
                        costab[i] = cos((PI/2./NCOSENTRY)*i);
                                        /* normalize angle */
        if (x < 0.)
                x = -x;
        i = (NCOSENTRY*2./PI) * x  +  0.5;
        while (i >= 4*NCOSENTRY)
                i -= 4*NCOSENTRY;
        switch (i / NCOSENTRY) {
        case 0:
                return(costab[i]);
        case 1:
                return(-costab[(2*NCOSENTRY)-i]);
        case 2:
                return(-costab[i-(2*NCOSENTRY)]);
        case 3:
                return(costab[(4*NCOSENTRY)-i]);
        }
        return(0.);             /* should never be reached */
}

#endif

/* Fast arctangent approximation due to Rajan et al. 2006 */
double
atan2a(double y, double x)
{
        double  ratio, aratio, val;

        if (x == 0)
                return (y > 0) ? PI/2. : 3./2.*PI;

        aratio = (ratio = y/x) >= 0 ? ratio : -ratio;

        if (aratio > 1.01)
                return PI/2. - atan2a(x, y);

        val = PI/4.*ratio - ratio*(aratio - 1.)*(0.2447 + 0.0663*aratio);

        return val + PI*(x < 0);
}
Revision:	3.10
Committed:	Fri Jul 12 05:16:02 2013 UTC (10 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R2P1, rad5R3, HEAD
Changes since 3.9:	+4 -2 lines
Log Message:	Took out definiton of tcos() when __FAST_MATH__ is defined
#	User	Rev	Content
1	gwlarson	3.1	#ifndef lint
2	greg	3.10	static const char RCSid[] = "$Id: tcos.c,v 3.9 2013/02/08 16:39:01 greg Exp $";
3	gwlarson	3.1	#endif
4			/*
5			* Table-based cosine approximation.
6			*
7			* Use doubles in table even though we're not nearly that accurate just
8			* to avoid conversion and guarantee that tsin(x)^2 + tcos(x)^2 == 1.
9			*
10			* No interpolation in this version.
11	greg	3.2	*
12	greg	3.5	* External symbols declared in rtmath.h
13	greg	3.2	*/
14
15	greg	3.3	#include "copyright.h"
16	gwlarson	3.1
17			#include <math.h>
18	schorsch	3.4
19			#include "rtmath.h"
20	gwlarson	3.1
21	greg	3.10	#ifndef __FAST_MATH__
22
23	gwlarson	3.1	#ifndef NCOSENTRY
24	greg	3.6	#define NCOSENTRY 1024
25	gwlarson	3.1	#endif
26
27			double
28	greg	3.5	tcos(double x) /* approximate cosine */
29	gwlarson	3.1	{
30			static double costab[NCOSENTRY+1];
31			register int i;
32
33			if (costab[0] < 0.5) /* initialize table */
34			for (i = 0; i <= NCOSENTRY; i++)
35			costab[i] = cos((PI/2./NCOSENTRY)*i);
36			/* normalize angle */
37			if (x < 0.)
38			x = -x;
39			i = (NCOSENTRY2./PI) x + 0.5;
40	greg	3.5	while (i >= 4*NCOSENTRY)
41			i -= 4*NCOSENTRY;
42	gwlarson	3.1	switch (i / NCOSENTRY) {
43			case 0:
44			return(costab[i]);
45			case 1:
46			return(-costab[(2*NCOSENTRY)-i]);
47			case 2:
48			return(-costab[i-(2*NCOSENTRY)]);
49			case 3:
50			return(costab[(4*NCOSENTRY)-i]);
51			}
52			return(0.); /* should never be reached */
53			}
54	greg	3.7
55	greg	3.10	#endif
56	greg	3.7
57			/* Fast arctangent approximation due to Rajan et al. 2006 */
58			double
59	greg	3.8	atan2a(double y, double x)
60	greg	3.7	{
61			double ratio, aratio, val;
62
63			if (x == 0)
64			return (y > 0) ? PI/2. : 3./2.*PI;
65
66			aratio = (ratio = y/x) >= 0 ? ratio : -ratio;
67
68			if (aratio > 1.01)
69	greg	3.9	return PI/2. - atan2a(x, y);
70	greg	3.7
71			val = PI/4.ratio - ratio(aratio - 1.)(0.2447 + 0.0663aratio);
72
73			return val + PI*(x < 0);
74			}