Вычисление экспоненты, синуса и косинуса
ОГЛАВЛЕНИЕ
В отличие от других страниц элементарных функций, на этой странице проводится их вычисление с точностью double precision. Кроме того, вообще не используется операция деления, что может быть важно для процессоров с затрудненной его реализацией.
Область значений аргумента изначально неограничена, но в результате вычисления для некоторых его значений могут возникнуть следующие ошибки: OVERFLOW, ISNAN, TLOSS (см. ниже).
Все три функции представляют из себя экспоненту от комплексного числа. По формуле Эйлера:
exp(a+ib)=ea(cos(b)+i*sin(b)).
Рассмотрим теперь последовательно действительную часть (экспоненту) и мнимую (синус и косинус).
Экспонента
Если число a положительно, то его можно представить в виде суммы целой и дробной части: a=[a]+{a}. Целая часть может быть записана в беззнаковое длинное значение:
[a]=a020+a121+...+a31231.
Дробная же часть может быть (для двойной точности) записана в формате двоичной дроби с фиксированной точкой, в которой имеется не более 48 битов:
{a}=a-12-1+a-22-2+...+a-482-48.
Поскольку, по свойству действительной экспоненты, ex+y=exey, то для ее вычисления необходимо просто перемножить те значения exp(2k), где k меняется от -48 до 31 (включая нуль), при которых соответствующие биты ak ненулевые.
Более того, поскольку для положительных k уже значение exp(28) дает число, близкое к предельно возможному для чисел двойной точности, то имеет смысл перемножать не более 8 значений для положительных k, а при обнаружении ненулевого бита ak при k>7 возвращать признак переполнения (OVERFLOW).
Для отрицательных значений a=-c можно поступить двумя путями: либо вычислить экспоненту для c, а затем обратить, либо провести аналогичное разложение: a=-[c]-{c}, но при этом использовать табличные данные для обратных величин exp(-2k). Вместо переполнения при обнаружении ненулевого бита для k>7 следует возвращать нуль. На данной странице приведен вариант с использованием табличных данных для обратных величин; это вообще исключает деление из процесса вычисления экспоненты.
Синус и косинус
Экспонента от мнимого значения (действительной частью является косинус, а мнимой - синус) может быть вычислена аналогично, но в данном случае вначале удобно разделить аргумент на 2*pi (умножив на 1/(2*pi)). Либо считать сразу, что вычисляется exp(2*pi*i*x) (для некоторых приложений важно именно это).
Для положительного значения аргумента x, он представляется в виде суммы целой и дробной части, однако, в этом случае целая часть нам не требуется, а дробную представляем в виде двоичной дроби:
{x}=x-12-1+x-22-2+...+x-482-48.
Заметим, однако, что в случае, когда целая часть достаточно велика (в нашей дальнейшей программе ограничение [x]<232), то точность дробной части числа снижена, либо вообще не определена. В этом случае следует генерировать ошибку TLOSS (полная потеря точности, при ограниченном результате).
Для расчета экспоненты мнимого аргумента следует аналогично предыдущему перемножить табличные значения exp(2*pi*i*2k) (k<0; комплекснозначные, то есть пары косинус-синус) для тех k, при которых биты xk ненулевые. При этом используются формулы комплексного умножения.
В случае отрицательного аргумента, нужно вычислить значения синуса и косинуса для соответствующего положительного, а затем обратить знак синуса.
Еще один вид ошибки ISNAN следует генерировать тогда, когда подставляемый для обработки аргумент уже является испорченным (nan, испорченное действительное число).
Комплексная экспонента
Для вычисления комплексной экспоненты используется вышеприведенная формула Эйлера. Единственная сложность -- классифицировать исключительные случаи и сгенерировать код ошибки. В нашей программе сделано следующее:
| Результат для Expo | Результат для CosSin | Результат для CExpo |
| ISNAN | любой | ISNAN |
| любой | ISNAN | ISNAN |
| OVERFLOW | любой | ISNAN |
| NORMAL =0 | TLOSS | NORMAL =0 |
| NORMAL =0 | NORMAL | NORMAL =0 |
| NORMAL !=0 | TLOSS | TLOSS |
| NORMAL !=0 | NORMAL | NORMAL |
Программы вычисления экспоненты, синуса и косинуса
Ниже находятся программы разложения числа на целую и дробную часть (дробная часть в специальном формате), вычисления экспоненты, вычисления синуса с косинусом. Кроме этого, на отдельной странице помещены необходимые табличные данные. В начале программы обратите внимание на литералы (перечислимые константы), характеризующие возвращаемые ошибки. Заметим также, что на аппаратном уровне (при знании формата представления действительного числа) программа может быть сильно оптимизирована.
Для одинарной точности следует уменьшить число обрабатываемых битов до интервала [-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);
}