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

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

В отличие от других страниц элементарных функций, на этой странице проводится их вычисление с точностью double precision. Кроме того, вообще не используется операция деления, что может быть важно для процессоров с затрудненной его реализацией.

Область значений аргумента изначально неограничена, но в результате вычисления для некоторых его значений могут возникнуть следующие ошибки: OVERFLOW, ISNAN, TLOSS (объяснение ниже).

Все три функции представляют из себя экспоненту от комплексного числа. По формуле Эйлера:

Если число a положительно, то его можно представить в виде суммы целой и дробной части: a=[a]+{a}. Целая часть может быть записана в беззнаковое длинное значение:

Более того, поскольку для положительных k уже значение exp(2⁸) дает число, близкое к предельно возможному для чисел двойной точности, то имеет смысл перемножать не более 8 значений для положительных k, а при обнаружении ненулевого бита a_k при k>7 возвращать признак переполнения (OVERFLOW).

Для отрицательных значений a=-c можно поступить двумя путями: либо вычислить экспоненту для c, а затем обратить, либо провести аналогичное разложение: a=-[c]-{c}, но при этом использовать табличные данные для обратных величин exp(-2^k). Вместо переполнения при обнаружении ненулевого бита для k>7 следует возвращать нуль. На данной странице приведен вариант с использованием табличных данных для обратных величин; это вообще исключает деление из процесса вычисления экспоненты.

Экспонента от мнимого значения (действительной частью является косинус, а мнимой - синус) может быть вычислена аналогично, но в данном случае вначале удобно разделить аргумент на 2*pi (умножив на 1/(2*pi)). Либо считать сразу, что вычисляется exp(2*pi*i*x) (для некоторых приложений важно именно это).

Для положительного значения аргумента x, он представляется в виде суммы целой и дробной части, однако, в этом случае целая часть нам не требуется, а дробную представляем в виде двоичной дроби:

Для расчета экспоненты мнимого аргумента следует аналогично предыдущему перемножить табличные значения exp(2*pi*i*2^k) (k<0; комплекснозначные, то есть пары косинус-синус) для тех k, при которых биты x_k ненулевые. При этом используются формулы комплексного умножения.

В случае отрицательного аргумента, нужно вычислить значения синуса и косинуса для соответствующего положительного, а затем обратить знак синуса.

Еще один вид ошибки ISNAN следует генерировать тогда, когда подставляемый для обработки аргумент уже является испорченным (nan, испорченное действительное число).

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

Для одинарной точности следует уменьшить число обрабатываемых битов до интервала [-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);
}