/* Sine/cosine/exponent using the special tables */
/* double precision implementation */

/* functions return */
enum{NORMAL,OVERFLOW,UNDERFLOW,ISNAN,TLOSS};

/* The following functions are present:
   Number decomposition into integer and fractial parts
   int int_fract(int *sign,unsigned short *dest,double src);
       sign -- sign of result; -1, 0 or 1 on output,
       dest -- array of 5 elements: dest[0] and dest[1] contain 
               32-bit unsigned integer part, dest[2]-dest[4] 
               contain the fractial part. 
       src -- the source number.
   Returns: 
       NORMAL - OK,
       OVERFLOW - integer part exceed 32 bits and positive sign,
       UNDERFLOW - integer part exceed 32 bits and negative sign,
       ISNAN - bad source.

   Exponent:
   int Expo(double *e,double x);
       e - result on output,
       x - argument,
   Returns:
       NORMAL - OK,
       OVERFLOW - too big argument,
       ISNAN - bad argument.

   Sine and cosine:
   int CosSin(double *c,double *s,double x);
       c - cosine result,
       s - sine result,
       x - argument,
   Returns:
       NORMAL - OK,
       TLOSS - loss of precision, in this case *c=1., *s=0. is stored,
       ISNAN - bad argument.

   Complex exponent:
   int CExpo(souble *er,double *ei,double xr,double xi);
       er - real part of result,
       ei - imaginary part of result,
       xr - real part of the argument,
       xi - imaginary part of the argument.
   Returns:
       NORMAL - OK
       ISNAN - bad argument or total unlimited loss of precision
       TLOSS - total limited loss of precision.
*/

/* The necessary tables, determining also K0=8, K1=48 */
#include "sicoex.h"

/* decompose a double precision number into 5 unsigned short numbers + sign.
   Two numbers: 2: integer part + 3 numbers fractal part + sign 
*/
int int_fract(int *sign,unsigned short *dest,double src) {
 unsigned long intp;
 unsigned short *dst,as;
 double a,*fr;
 int i,k,ks;
 if(src<0.) {*sign=-1;src=-src;}
 else if(src>0.) *sign=1;
 else if(src==0.) {
  *sign=0;
  for(i=0;i<5;i++) dest[i]=0;
  return(NORMAL);
 }
 else return(ISNAN);
 if(src>(double)0xffffffff) return(*sign==1?OVERFLOW:UNDERFLOW);
 intp=(unsigned long)src;
 src-=(double)intp;
 dest[0]=(unsigned short)(intp&0xffff);
 dest[1]=(unsigned short)(intp>>16);
 dst=dest+2;fr=Frac;
 for(i=0;i<3;i++) {
  *dst=0;as=1;
  for(k=0;k<16;k++) {
   if((a=src-(*fr))>=0.) {src=a;*dst|=as;}
   fr++;as<<=1;
  }
  dst++;
 }
 return(NORMAL);
}

/* Calculate the exponent */
int Expo(double *e,double x) {
 int sign;
 unsigned short xx[5],as,*xt;
 int i,k,ks;
 double *ep,*en;
 switch(int_fract(&sign,xx,x)) {
  case OVERFLOW: return(OVERFLOW);
  case ISNAN: return(ISNAN);
  case UNDERFLOW: *e=0.;return(NORMAL);
 }
 *e=1.;
 if(sign==0) return(NORMAL);
 else if(sign>0) {ep=EpPos;en=EpNeg;}
 else {ep=EnPos;en=EnNeg;}
 /* compose the integer part */
 xt=xx;ks=0;
 for(i=0;i<2;i++) {
  as=1;
  for(k=0;k<16;k++) {
   if(*xt&as) {
    if(ks>=K0) {
     if(sign>0) return(OVERFLOW);
     else {*e=0.;return(NORMAL);}
    }
    (*e)*=(*ep);
   }
   ks++;ep++;as<<=1;
  }
  xt++;
 }
 /* compose the fractial part */
 xt=xx+2;
 for(i=0;i<3;i++) {
  as=1;
  for(k=0;k<16;k++) {
   if(*xt&as) (*e)*=(*en);
   en++;as<<=1;
  }
  xt++;
 }
 return(NORMAL);
}

/* calculate the sine and cosine */
#define M_PI 3.141592653589793
#define M_2PI (2.*M_PI)
int CosSin(double *co,double *si,double x) {
 int sign;
 unsigned short xx[5],*xt;
 int i,k,as;
 double st,ct,*sit,*cot;
 /* divide x by 2*pi */
 x/=M_2PI;
 /* decompose x */
 switch(int_fract(&sign,xx,x)) {
  case OVERFLOW:case UNDERFLOW: {*co=1.;(*si)=0.;return(TLOSS);}
  case ISNAN: return(ISNAN);
 }
 *co=1.;*si=0.;
 cot=Cosi;sit=Sine;
 xt=xx+2;
 for(i=0;i<3;i++) {
  as=1;
  for(k=0;k<16;k++) {
   if(*xt&as) {
    ct=(*co);st=(*si);
    *co=ct*(*cot)-st*(*sit);
    *si=ct*(*sit)+st*(*cot);
   }
   as<<=1;sit++;cot++;
  }
  xt++;
 }
 /* if the sign was changed, change sine result */
 if(sign<0) *si=-(*si);
 return(NORMAL);
}

/* complex exponent */
int CExpo(double *er,double *ei,double xr,double xi) {
  int re,rcs;
  double c,s;
  re=Expo(er,xr);
  if(re==ISNAN || re==OVERFLOW) return(ISNAN);
  rcs=CosSin(&c,&s,xi);
  if(rcs==ISNAN) return(ISNAN);
  if(*er==0.) {*ei=0.; return(NORMAL);}
  *ei=(*er)*s; (*er)*=c;
  if(rcs==TLOSS) return(TLOSS);
  else return(NORMAL);
}