# Экспонента, синус и косинус

На данной странице представлен вариант вычисления экспоненты (функции exp()) с использованием последовательных умножений на табличные данные. Поскольку синус и косинус также представляют из себя экспоненту, только от мнимого аргумента, то аналогичный алгоритм (альтернативный тому, что описан здесь) можно привести и для их вычисления.

В отличие от других страниц элементарных функций, на этой странице проводится их вычисление с точностью 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.

    /* Sine/cosine/exponent using special tables. Double precision argument. */
    /* double precision implementation */

    /* functions return constants */
    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 - total (limited) 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 and determination of bits number */
    #include "sicoex.h"

    /* decompose a double precision number into 5 unsigned short numbers + 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)
    #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;
      /* divide x by 2*pi */
      x*=M_R2PI;
      /* 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);
    }

Скачать программу sicoex.c и включаемые таблицы sicoex.h

Элементарные функции
Индекс