/* Синус/косинус/экспонента с использованием специальных таблиц. Аргумент двойной точности. */

/* функция, возвращающая константы */
enum{NORMAL,OVERFLOW,UNDERFLOW,ISNAN,TLOSS};

/* Представлены следующие фукнции:
Разложение числа на целую и дробную части
int int_fract(int *sign,unsigned short *dest,double src);
sign -- знак результата; -1, 0 или 1 на выходе,
dest -- массив из 5 элементов: dest[0] и dest[1] содержат  
32-bit беззнаковую целую часть, dest[2]-dest[4]  
содержат дробную часть.  
src -- исходное число.
Результат:  
NORMAL - OK,
OVERFLOW - целая часть превышает 32 bits с положительным знаком,
UNDERFLOW - целая часть превышает 32 bits с отрицательным знаком,
ISNAN - некорректное исходное число.

Экпонент:
int Expo(double *e,double x);
e - результат,
x - аругмент,
Возвращает:
NORMAL - OK,
OVERFLOW - слишком большой аргумент,
ISNAN - некорректный аргумент.

Синус и косинус:
int CosSin(double *c,double *s,double x);
c - косинус,
s - синус,
x - аргумент,
Возвращет:
NORMAL - OK,
TLOSS - общая потеря точности, в этом случае *c=1., *s=0,
ISNAN - некорректный аргумент.

Комплексный экспонент:
int CExpo(souble *er,double *ei,double xr,double xi);
er - действительная часть результата,
ei - мнимания часть резултата,
xr - действительная часть аргумента,
xi - мнимая часть аргумента.
Возвращает:
NORMAL - OK
ISNAN - некорректный аргумент или общая потеря точности
TLOSS - общая потеря точности.
*/

/* Необходимые таблицы и определения разрядов числа */
#include "sicoex.h"

/* разложение числа двойной точности на 5 чисел без знака + знак.
*/
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);
}

/* вычисление экспоненты */
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;}
/* получение целой части */
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++;
}
/* получение дробной части */
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);
}

/* вычисление синуса и косинуса */
#define M_PI 3.141592653589793
#define M_2PI (2.*M_PI)
#define M_R2PI (1./M_2PI)
int CosSin(double *co,double *si,double x) {
int sign;
unsigned short xx[5],*xt;
int i,k,as;
double st,ct,*sit,*cot;
/* деление x на 2*pi */
x*=M_R2PI;
/* разложение 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(sign<0) *si=-(*si);
return(NORMAL);
}

/* комплексная экспонента */
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);
}