bs.cpp

代码

double bs(int n, double a)
{
	int m;
    double b[20],temp,del,r,eps,s,ak,a2,omeg,da,ajp,ajm,q1,q2;
    if (a == 0.0)
	{
       if ((n / 2) == (n / 2.0))
	   {
         b[n] = 0.0;
         goto yw;
       }
       b[n] = 2.0 / n;
       goto yw;
    }
    if (fabs(a) < 8.0)
	{
       del = 0.00000001;
       if ((n / 2) != (n / 2.0))
	   {
         r = 2.0 / n;
         eps = r * del;
         s = r;
         ak = 0.0;
         a2 = a * a;
w:        ak = ak + 2.0;
         r = r * a2 * (n + ak - 2.0) / (ak * (ak - 1.0) * (n + ak));
         s = s + r;
         if (r > eps) 
			 goto w;
         b[n] = s + r;
         goto yw;
       }
       r = 2.0 * a / (n + 1.0);
       omeg = fabs(r * del);
       s = r;
       ak = 1.0;
       a2 = a * a;
l:      ak = ak + 2.0;
       r = r * a2 * (n + ak - 2.0) / (ak * (ak - 1.0) * (n + ak));
       s = s + r;
       if (fabs(r) > omeg)
		   goto l;
       b[n] = -(s + r);
       goto yw;
    }
    da = 1.0 / a;
    ajp = da * exp(a);
    ajm = (da * da) / ajp;
    b[1] = ajp - ajm;
    if (n == 1)
		goto yw;
    q1 = -1.0;
    q2 = 1.0;
    for (m = 2;m<=n;m++)
	{
        b[m] = q1 * ajp - ajm + q2 * da * b[m - 1];
        q1 = -q1;
        q2 = q2 + 1.0;
    }
yw:  temp = b[n];
	 return temp;
}