К разделу 'Математика' сайта AlgoList.

Назад Оглавление Вперед 

Приложение В:
Коды на Си для алгоритмов из §1

 
#include <math.h>
#include <stdio.h>
#define sign(x) (x == 0 ? 0 : (x < 0 ? -1 : 1))
const double pi_const=0.3989422804014;
int N;

double
I(double x, double eps)
{
   double t=x, sum=t, x2=x*x, s=0;
   int n;
   N=0;
   for(n=2; fabs(s-sum) > eps; n+=2) {
      t*=-x2/n; s=sum; sum+=t/(n+1);
      N++; /* Внешняя переменная! */
   }
   return 0.5+sum*pi_const;
  /* 8+11*N операций */
}
double
T(double x, double eps)
{
   double t=x, sum=t, x2=x*x, s=0;
   int n;
   N=0;
   for(n=3; fabs(s-sum) > eps; n+=2) {
      t*=x2/n; s=sum; sum+=t;
      N++; /* Внешняя переменная! */
   }
   return 0.5+sign(x)*sum*exp(-x2/2.)*pi_const;
  /* 14+8*N операций */
}
double
S(double x, double eps)
{
   double x2=x*x,
          rho=1./(3.-x2),
          sum=fabs(x)*rho, term=x2*sum, t=0, s;
   int n;
   N=0;
   rho*=3.; s=(sum*=3.);
   for(n=2; (x2<0) || ((fabs(s-t) > eps) ||
                       (fabs(s-sum) > eps)); n++) {
      rho = 1./(1.+rho*x2*n/(4.*n*n-1));
      term*=rho-1; t=s; s=sum; sum+=term;
      x2=-x2;
      N++; /* Внешняя переменная! */
   }
   return 0.5+sign(x)*sum*exp(-x2/2.)*pi_const;
   /* 23+25*N операций */
}
double
L(double x, double eps)
{
   double x2=x*x, y=exp(-x2/2.)*pi_const,
          sum, term, r, rho, t, s;
   N=0;
   if (x == 0) return 0.5;
   x2=1./x2; term=sum=y/fabs(x); r=s=0; rho=1;
   do {
      r+=x2; rho=1./(1+r*rho); term*=rho-1;
      t=s; s=sum; sum+=term;
      N++; /* Внешняя переменная! */
   } while(fabs(s-t) > eps || fabs(s-sum) > eps);
   return (x > 0 ? 1-sum : sum);
  /* 16+10*N операций */
}
void main(void)
{
   double x;
   double f;
   for(x=0.25; x<8; x+=0.25) {
#if TABLE
      printf("\n  %4.2f", x);
      f = I(x,0);
      printf("\t%d\t%d", N, 8+11*N);
      f = T(x,0);
      printf("\t%d\t%d", N, 14+8*N);
      f = S(x,0);
      printf("\t%d\t%d", N, 23+25*N);
      f = L(x,0);
      printf("\t%d\t%d", N, 16+10*N);
#else
      f = I(x,0);
      printf("\nI(%4.2f)=%10.8g\t%d\t%d", x, f,N,8+11*N);
      f = T(x,0);
      printf("\nT(%4.2f)=%10.8g\t%d\t%d", x, f,N,14+8*N);
      f = S(x,0);
      printf("\nS(%4.2f)=%10.8g\t%d\t%d", x, f,N,23+25*N);
      f = L(x,0);
      printf("\nL(%4.2f)=%10.8g\t%d\t%d", x, f,N,16+10*N);
#endif
   }
}
 
Вверх 
Назад Оглавление Вперед 

Назад Оглавление Вперед 

Приложение Г:
Коды на Си для алгоритмов из §3

 
#include <math.h>
#include <stdio.h>

#define sign(x) (x == 0 ? 0 : (x < 0 ? -1 : 1))

const double pi_const=0.3989422804014;
int N;

double
D(double x, double eps)
{
   double x1, xn,g0, g1, s=0, sum, t, z;
   static const double PHI[] = {
      0.5, /* x=0 */
      0.841344746068543, /* x=1 */
      0.977249868051821, /* x=2 */
      0.998650101968370, /* x=3 */
      0.999968328758167, /* x=4 */
      0.999999713348428, /* x=5 */
      0.999999999013412, /* x=6 */
      0.999999999998720, /* x=7 */
      0.999999999999999 /* x=8 */
   },
   H=1.0;
   int i, n;
   x1=fabs(x); i=(int)(x1/H+0.5); z=i*H;
   g0=t=PHI[i]; g1=exp(-z*z/2.)*pi_const;
   xn=(x1-=z); sum=t+g1*xn;
   N=0;
   for(n=2;fabs(s-t) > eps || fabs(t-sum) > eps; n++) {
      s=-z*g1-(n-2)*g0; g0=g1; g1=s; xn*=x1/n;
      s=t; t=sum; sum+=g1*xn;
      N++; /* Внешняя переменная! */
   }
   return (x>0 ? sum : 1-sum);
}
double
P(double x, double eps)
{
   double x1, xn,c0, c1, rz, s, t=0, z;
   int i, n;
   static const double R[] = {
      0.,
      1.4106861306,
      8.8394391760,
      112.51515173,
      3735.8400745,
      336310.71435,
      82292564.750,
      54736210525.,
      98965389434000.
   },
   H=1.0;
   x1=fabs(x); i=(int)(x1/H+0.5); z=i*H;
   c0=s=R[i]; c1=1.+z*c0;
   xn=(x1-=z); rz=s+c1*xn;
   N=0;
   for(n=2; fabs(s-t) > eps || fabs(t-rz) > eps; n++) {
      s=(c0+c1*z)/n; c0=c1; c1=s; xn*=x1;
      s=t; t=rz; rz+=c1*xn;
      N++; /* Внешняя переменная! */
   }
   return 0.5+sign(x)*rz*exp(-x*x/2)*pi_const;
}

double
F(double x, double eps)
{
   double x1=fabs(x), a, z=0, fz=0;

   double t, s, sum, xn, g0, g1;
   int n;
   double h=0.25;
   N=0;
   while((a=x1-z)>0) {
      if (a > h) a = h;
      /* Подготовка к вычислению ряда. */
      xn = a; g0 = g1 = exp(-z*z/2)*pi_const;
      t = sum = g1*xn;
      for(n=2, s=0;(fabs(s-t)>eps) || (fabs(t-sum)>eps); n++) {
         N++; /* Внешняя переменная! */
         /* Вычислим очередную частичную сумму */
         xn *= a/n; s = -z*g1-(n-2)*g0;
         g0 = g1; g1 = s; s = t; t = sum; sum += g1*xn;
      };
      /* В т. z разложение посчитано - пошли дальше */
      fz += sum; z += a;
   }
   return 0.5 + sign(x)*fz;
}

double
C(double x, double eps)
{
   double x1, ra, rz, z, d;
   double t, s, xa, xn, c0, c1;
   int n, gotta=0;
   double h=0.5;
   x1=fabs(x); rz=t=z=0;
   N=0;
   while(!gotta) {
      xa=z; ra=rz; z+=h;
      if(z >= x1) {
         z = x1; gotta = 1;
      }
      c0=s=ra; c1=1.+xa*ra;
      xn=d=z-xa; rz=ra+c1*xn;
      for(n=2;(fabs(s-t)>eps) || (fabs(s-rz)>eps); n++) {
         N++; /* Внешняя переменная! */
         /* Вычислим очередную частичную сумму */
         t=(c0+xa*c1)/n; c0=c1; c1=t; xn*=d;
         t=s; s=rz; rz=s+c1*xn;
      }
   }
   return 0.5 + sign(x)*rz*exp(-x*x/2)*pi_const;
}

double
B(double x, double eps)
{
   double x1=fabs(x), fz=0;
   double t, s, sum, xn, g0, g1;
   int n;
   double a=0.5;

   N=0;
   while(x1>0) {
      if (x1 <= a) a = x1;
      /* Подготовка к вычислению ряда. */
      xn = a; g0 = g1 = exp(-x1*x1/2)*pi_const;
      t = sum = g1*xn;
      for(n=2, s=0;(fabs(s-t)>eps) || (fabs(t-sum)>eps); n++) {
         N++; /* Внешняя переменная! */
         /* Вычислим очередную частичную сумму */
         xn *= -a/n; s = -x1*g1-(n-2)*g0;
         g0 = g1; g1 = s; s = t; t = sum; sum += g1*xn;
      }
      /* В т. z разложение посчитано - пошли дальше */
      fz += sum; x1 -= a;
   }
   return 0.5 + sign(x)*fz;
}

double
E(double x, double eps)
{
   double x1=fabs(x)*0.707106781186548, z=0;
   double y, s, xn, u, v, w=0;
   int n;
   double a=0.5;
   N=0;
   while(x1>0) {
      if (x1 <= a) a = x1;
      /* Подготовка к вычислению ряда. */
      xn = a; u = v = exp(-x1*x1)*0.564189583547756;
      s = y = v*a;
      for(n=2;(fabs(w-s)>eps) || (fabs(s-y)>eps); n++) {
         N++; /* Внешняя переменная! */
         /* Вычислим очередную частичную сумму */
         w=-2.*(x1*v+u*(n-2.));
         u=v; v=w; xn*=-a/n; w=s; s=y; y+= v*xn;
      }
      /* В т. z разложение посчитано - пошли дальше */
      z+=y; x1-=a;
   }
   return 0.5 + sign(x)*z;
}

void main(void)
{
   double x;
   double f;

   for(x=0.25; x<8; x+=0.25) {
      #ifdef TABLE
      printf("\n %4.2f", x);
      f = D(x,0);
      printf("\t%d", N);
      f = P(x,0);
      printf("\t%d", N);
      f = C(x,0);
      printf("\t%d", N);
      f = B(x,0);
      printf("\t%d", N);
      f = E(x,0);
      printf("\t%d", N);
      #else
      f = D(x,0);
      printf("\nD(%4.2f)=%10.8g\t%d", x, f, N);
      f = P(x,0);
      printf("\nP(%4.2f)=%10.8g\t%d", x, f, N);
      f = C(x,0);
      printf("\nC(%4.2f)=%10.8g\t%d", x, f, N);
      f = B(x,0);
      printf("\nB(%4.2f)=%10.8g\t%d", x, f, N);
      f = E(x,0);
      printf("\nE(%4.2f)=%10.8g\t%d", x, f, N);
      #endif
   }
}
Вверх
Назад Оглавление Вперед