Приложение В: Коды на Си для алгоритмов из "Основные разложения".
#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
}
}
Приложение Г: Коды на Си для алгоритмов из "Еще алгоритмы".
|
#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
}
}
8 8 8
| 
*
Обратная связь