|
Страница 5 из 5
Программы вычисления экспоненты, синуса и косинуса
Ниже находятся программы разложения числа на целую и дробную часть (дробная часть в специальном формате), вычисления экспоненты, вычисления синуса с косинусом. Кроме этого, на отдельной странице помещены необходимые табличные данные. В начале программы обратите внимание на литералы (перечислимые константы), характеризующие возвращаемые ошибки. Заметим также, что на аппаратном уровне (при знании формата представления действительного числа) программа может быть сильно оптимизирована. Для одинарной точности следует уменьшить число обрабатываемых битов до интервала [-24,4] (вместо [-48,7]) и объявить табличные данные как float вместо double. /* Синус/косинус/экспонента с использованием специальных таблиц. Аргумент двойной точности. */
/* функция, возвращающая константы */
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);
} |