#include <stdio.h>
#include <math.h>
typedef double R;
typedef struct{R cos; R sin;} A;

R pw(int j) {if(j<0) return 1./pw(-j); if(j>3) return 16.*pw(j-4);
  if (j>0) return 2.*pw(j-1); return 1;}
A cs(R z){A a = {0.6072529350088812561694, 0}; {int j; for(j=0; j<53; ++j) {
  R p = (z<0?-1:1)*pw(-j); z -= atan(p);
  {R t = a.sin; a.sin += p*a.cos; a.cos -= p*t;}
  if(0) printf("%e %2d %19.15f %19.15f %19.15f %e\n", p,
     j, z, a.cos, a.sin, atan(a.cos/a.sin));}}
return a;}

#include <stdlib.h>
void mx(R t){A a = cs(t); printf("%19.15f %19.15f %19.15f %19.15f %19.15f\n",
   t, a.sin, sin(t), a.cos, cos(t));}
int main(){mx(0); R pi=4*atan(1);
  mx(pi/6); mx(pi/2); mx(pi/4); mx(pi/3); mx(-1.1);

{int j=100000; R f = pw(-31)*pi; R me = 0, ev;
  while(j--){R v = (random()*f) - pi/2;
  void ue(R t){if (fabs(t)> me) {me = fabs(t); ev=v;}}
  A a = cs(v); ue(a.cos - cos(v)); ue(a.sin - sin(v));}
  printf("max error = %e at %20.16f\n", me, ev);}
return 0;}

/*
  0.000000000000000  -0.000000000000000   0.000000000000000   1.000000000000000   1.000000000000000
  0.523598775598299   0.500000000000000   0.500000000000000   0.866025403784438   0.866025403784439
  1.570796326794897   1.000000000000000   1.000000000000000  -0.000000000000000   0.000000000000000
  0.785398163397448   0.707106781186548   0.707106781186547   0.707106781186548   0.707106781186548
  1.047197551196598   0.866025403784438   0.866025403784439   0.500000000000000   0.500000000000000
 -1.100000000000000  -0.891207360061435  -0.891207360061435   0.453596121425577   0.453596121425577
max error = 1.110223e-15 at  -0.4460695575018367
*/