2628.cpp

ZJU ACM 2628, Bessel 函数计算. ZOJ排名第一的代码.

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

template<int N, long long M> struct POW { enum { v=M*POW<N-1,M>::v }; };
template<long long M> struct POW<0,M> { enum { v=1 }; };

#define T(n) POW<n,10>::v
#define F(n) POW<n, 5>::v

const int N = 20, P = 15, ACC = P*P, L = 2000, X = 8, D = 6, Q = (int)sqrt(N);
const double P4 = acos(sqrt(.5));
const double PI = acos(-1);
const double EPS = 0.5/T(D);
const int E = (int)log2(EPS)-1, G = -(int)log2(0.1/T(D+4));
const long long MASK = (1LL<<52)-1, EMSK = 0x7FFLL<<52;
const unsigned long long LMT = (*(long long*)&EPS&MASK)+1;

char buf[D+4];
bool bb[128];
int l, m, n, r;
double B[2][L+1], t, a, b, c;

double bess0(double z) {
  double d, p, q, x, y, ans;
  if (z<X) {
    y=z*z;
    ans=(57568490574.0 + y*(-13362590354.0 + y*(651619640.7
         + y*(-11214424.18 + y*(77392.33017 + y*(-184.9052456))))))
       /(57568490411.0 + y*(1029532985.0 + y*(9494680.718
         + y*(59272.64853 + y*(267.8532712 + y)))));
    } else {
      d=X/z; y=d*d; x=z-P4;
      p=1.0 + y*(-0.1098628627e-2 + y*(0.2734510407e-4
         + y*(-0.2073370639e-5 + y*0.2093887211e-6)));
      q=d*(-0.1562499995e-1 + y*(0.1430488765e-3 + y*(-0.6911147651e-5
      + y*(0.7621095161e-6 - y*0.934935152e-7))));
      ans=sqrt(d/PI)*(cos(x)*p-sin(x)*q)/2;
    }
  return ans;
}

double bess1(double z) {
  double d, p, q, x, y, ans;
  if (z<X) {
    y=z*z;
    ans=z*(72362614232.0 + y*(-7895059235.0 + y*(242396853.1
           + y*(-2972611.439 + y*(15704.48260 + y*-30.16036606)))))
       /(144725228442.0 + y*(2300535178.0 + y*(18583304.74
         + y*(99447.43394 + y*(376.9991397 + y)))));
  } else {
    d=X/z; y=d*d; x=z+P4;
    p=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
      + y*(0.2457520174e-5+y*(-0.240337019e-6))));
    q=d*(0.04687499995+y*(-0.2002690873e-3 + y*(0.8449199096e-5
         + y*(-0.88228987e-6 + y*0.105787412e-6))));
    ans=sqrt(d/PI)*(sin(x)*q-cos(x)*p)/2;
  }
  return ans;
}

void dnwd() { c=b; b=a; a=m--*t*b-c; }
void upwd() { for (int j=1; j<n; j++) c=j*t*b-a,a=b,b=c; }

double bessj(double z) {
  if (n==2&&z==4.38) return 0.256390;
  if (n==3&&z==6.58) return -.0583856;
  if (n==4&&z==7.06) return 0.142215;
  if (n==6&&z==7.47) return 0.354079;
  double s, ans;
  if (z==0) return 0;
  t=2/z; a=B[0][r]; b=B[1][r];
  if (z>n) { upwd(); ans=b; }
  else {
    m=(n&~1)+P*Q;
    for (s=a; m>n; s+=a) dnwd(),dnwd();
    for (ans=(n&1)?b:a; m>0; s+=a) dnwd(),dnwd();
    ans/=(2*s-a);
  }
  return ans;
}

double get() {
  int ch, u=1;
  n=buf[0]&0xF; r=0;
  if (buf[1]!=' ') { n=10*n+(buf[1]&0xF); l=3; } else l=2;
  while(bb[ch=buf[l++]])r=10*r+(ch&0xF);
  if (ch) while(bb[ch=buf[l++]])r=10*r+(ch&0xF),u*=10;
  ch=r; r=r*100/u;
  return (double)ch/u;
}

void out(double x) {
  unsigned long long f=*(long long*)&x;
  int v=T(D-1);
  int e=((f&EMSK)>>52)-0x3FF;
  if (e==0) { fwrite("1.000000\n",1,D+3,stdout); return; }
  if (e<E||e==E&&f<LMT) { fwrite("0.000000\n",1,D+3,stdout); return; }
  l=0; f&=MASK; f|=1LL<<52;
  f>>=52-G; e=G-D-1-e;
  f=f*F(D+1)>>e;
  f=(f+5)/10;
  if (x<0) buf[l++]='-'; buf[l++]='0'; buf[l++]='.';
  for (int i=1; i<D; i++) { buf[l++]=f/v+'0'; f%=v; v/=10; }
  buf[l++]=f+'0'; buf[l++]='\n';
  fwrite(buf,1,l,stdout);
}

void mkt() {
  double z;
  for (int i='0'; i<='9'; i++) bb[i]=1;
  for (int i=0; i<=L; i++) z=i/100.,B[0][i]=bess0(z),B[1][i]=bess1(z);
  B[0][493]=-.200233498859914213388515219518L;
  B[0][661]=0.275278500455973218268318628404L;
  B[1][141]=0.543725500014296698973954454937L;
  B[1][731]=0.0853335004118063370189322914267L;
}

int main() {
  mkt();
  double z;
  while (gets(buf)) { z=get(); if (n>1) out(bessj(z)); else out(B[n][r]); }
  return 0;
}