#include <math.h>

void sincosint(float x, float *si, float *ci)
{
	void sincosfg(float, float *, float *);
	float chepolsum(int, float, float []);
	float absx,z,f,g;

	absx=fabs(x);
	if (absx <= 4.0) {
		float a[11],z2;
		a[0] =2.7368706803630e0;   a[1] = -1.1106314107894e0;
		a[2] =1.4176562194666e-1;  a[3] = -1.0252652579174e-2;
		a[4] =4.6494615619880e-4;  a[5] = -1.4361730896642e-5;
		a[6] =3.2093684948229e-7;  a[7] = -5.4251990770162e-9;
		a[8] =7.1776288639895e-11; a[9] = -7.6335493723482e-13;
		a[10]=6.6679958346983e-15;
		z=x/4.0;
		z2=z*z;
		g=z2+z2-1.0;
		*si = z*chepolsum(10,g,a);
		a[0] =2.9659610400727e0;   a[1] = -9.4297198341830e-1;
		a[2] =8.6110342738169e-2;  a[3] = -4.7776084547139e-3;
		a[4] =1.7529161205146e-4;  a[5] = -4.5448727803752e-6;
		a[6] =8.7515839180060e-8;  a[7] = -1.2998699938109e-9;
		a[8] =1.5338974898831e-11; a[9] = -1.4724256070277e-13;
		a[10]=1.1721420798429e-15;
		*ci = 0.577215664901533+log(absx)-z2*chepolsum(10,g,a);
	} else {
		float cx,sx;
		sincosfg(x,&f,&g);
		cx=cos(x);
		sx=sin(x);
		*si = 1.570796326794897;
		if (x < 0.0) *si = -(*si);
		*si -= f*cx+g*sx;
		*ci = f*sx-g*cx;
	}
}