typedef struct {
float cos1, sin1, cos2, sin2;
} ARG_FF_2; /* vector 2 pairs */
ARG_FF_2 cosf_sinf_2(x1, x2)
union {
float flt;
int iflt;
} x1, x2;
{ /* 2 pair single precision
sin/cos function */
#include "float.h"
#if FLT_ROUNDS != 1
#error "rounding mode not nearest, fix code"
#endif
#include <errno.h>
#ifndef errno
extern int errno;
#endif
#include <math.h>
#define M_2_PI 0.63661977236758134308
#define T2PI 2*M_2_PI
#define TP2 T2PI*T2PI
#define TP3 TP2*T2PI
#define TP4 TP2*TP2
#define TP5 TP3*TP2
ARG_FF_2 res;
double tn, td, r;
/* integer compare with 0 same as float, for IEEE
** since arg comes in int register, this is faster
**
** doing everything in double, we won't lose accuracy
** by converting arg to multiple of PI/2
** this allows range reduction by subtracting an integer
**
** reduce to range +- 2, divide by 2 later
** when it cannot underflow */
tn = x1.flt * M_2_PI + (td = x1.iflt >= 0 ?
4 / LDBL_EPSILON : -4 / LDBL_EPSILON);
if (fabs(x1.flt * M_2_PI) >= 4 / LDBL_EPSILON |
fabs(x2.flt * M_2_PI) >= 4 / LDBL_EPSILON)
errno = ERANGE;
tn -= td;
tn = x1.flt * M_2_PI - tn;
td = tn * tn;
/* divide arg by 2 and rationalize numerator and denominator
** numerator of rational approx for tan(x1/2)
** Horner polynomials 1st time */
tn *= 886.77348 * TP4 + td * (-99.398954 * TP2 + td);
/* denominator */
td = 886.77346 * TP5 + td *
(-394.98971 * TP3 + td * 14.425694 * T2PI);
/* cos, sin half angle formulae, rationalized */
res.cos1 = (td * td - tn * tn) / (td * td + tn * tn);
res.sin1 = (tn * td + tn * td) / (td * td + tn * tn);
/* copy 2 */
tn = x2.flt * M_2_PI + (td = x2.iflt >= 0 ?
4 / LDBL_EPSILON : -4 / LDBL_EPSILON);
tn -= td;
tn = x2.flt * M_2_PI - tn;
td = tn * tn;
/* distribute terms to finish polynomials quicker */
tn *= 886.77348 * TP4 - td * 99.398954 * TP2 + td * td;
r = 886.77346 * TP5 - td * 394.98971 * TP3;
td = r + td * td * 14.425694 * T2PI;
res.cos2 = (td * td - tn * tn) / (td * td + tn * tn);
res.sin2 = (tn * td + tn * td) / (td * td + tn * tn);
return res;
}
/* End of File */