📄 cbess0_cpp.txt
字号:
ZLOG(ARGR, ARGI, &STR, &STI, &IDUM);
CZR = CZR - 0.25*STR - AIC;
CZI = CZI - 0.25*STI;
e120: AX = EXP(RCZ)/TOL;
AY = CZI;
CZR = AX*COS(AY);
CZI = AX*SIN(AY);
ZUCHK(CZR, CZI, &NW, ASCLE, TOL);
if (NW != 0) goto e90;
e130: if (IKFLG == 2) {
vmfree(vmblock);
return;
}
if (N == 1) {
vmfree(vmblock);
return;
}
/*----------------------------------------------------------------------
! SET UNDERFLOWS ON I SEQUENCE
!---------------------------------------------------------------------*/
e140: GNU = FNU + 1.0*(NN-1);
if (IFORM == 2) goto e150;
INIT = 0;
ZUNIK(ZRR, ZRI, GNU, IKFLG, 1, TOL, INIT, &PHIR, &PHII, &ZETA1R, &ZETA1I,
&ZETA2R, &ZETA2I, &SUMR, &SUMI, CWRKR, CWRKI);
CZR = -ZETA1R + ZETA2R;
CZI = -ZETA1I + ZETA2I;
goto e160;
e150: ZUNHJ(ZNR, ZNI, GNU, 1, TOL, &PHIR, &PHII, &ARGR, &ARGI, &ZETA1R, &ZETA1I,
&ZETA2R, &ZETA2I, &ASUMR, &ASUMI, &BSUMR, &BSUMI);
CZR = -ZETA1R + ZETA2R;
CZI = -ZETA1I + ZETA2I;
AARG = ZABS(ARGR,ARGI);
e160: if (KODE == 1) goto e170;
CZR = CZR - ZBR;
CZI = CZI - ZBI;
e170: APHI = ZABS(PHIR,PHII);
RCZ = CZR;
if (RCZ < -ELIM) goto e180;
if (RCZ > -ALIM) {
vmfree(vmblock);
return;
}
RCZ += log(APHI);
if (IFORM == 2) RCZ -= (0.25*log(AARG) - AIC);
if (RCZ > -ELIM) goto e190;
e180: YR[NN] = ZEROR;
YI[NN] = ZEROI;
NN = NN - 1;
*NUF = *NUF + 1;
if (NN == 0) {
vmfree(vmblock);
return;
}
goto e140;
e190: ASCLE = 1000*D1MACH(1)/TOL;
ZLOG(PHIR, PHII, &STR, &STI, &IDUM);
CZR = CZR + STR;
CZI = CZI + STI;
if (IFORM == 1) goto e200;
ZLOG(ARGR, ARGI, &STR, &STI, &IDUM);
CZR = CZR - 0.25*STR - AIC;
CZI = CZI - 0.25*STI;
e200: AX = EXP(RCZ)/TOL;
AY = CZI;
CZR = AX*COS(AY);
CZI = AX*SIN(AY);
ZUCHK(CZR, CZI, &NW, ASCLE, TOL);
if (NW != 0) goto e180;
vmfree(vmblock);
return;
e210: *NUF = -1;
vmfree(vmblock);
} //ZUOIK()
void ZS1S2(REAL *ZRR, REAL *ZRI, REAL *S1R, REAL *S1I, REAL *S2R, REAL *S2I,
int *NZ, REAL ASCLE, REAL ALIM, int *IUF) {
/***BEGIN PROLOGUE ZS1S2
!***REFER TO ZBESK,ZAIRY
!
! ZS1S2 TESTS FOR A POSSIBLE UNDERFLOW RESULTING FROM THE
! ADDITION OF THE I AND K FUNCTIONS IN THE ANALYTIC CON-
! TINUATION FORMULA WHERE S1=K FUNCTION AND S2=I FUNCTION.
! ON KODE=1 THE I AND K FUNCTIONS ARE DIFFERENT ORDERS OF
! MAGNITUDE, BUT FOR KODE=2 THEY CAN BE OF THE SAME ORDER
! OF MAGNITUDE AND THE MAXIMUM MUST BE AT LEAST ONE
! PRECISION ABOVE THE UNDERFLOW LIMIT.
!
!***ROUTINES CALLED ZABS,ZEXP,ZLOG
!***END PROLOGUE ZS1S2
! COMPLEX CZERO,C1,S1,S1D,S2,ZR */
//Label e10
REAL AA, ALN, AS1, AS2, C1I, C1R, S1DI, S1DR, ZEROI, ZEROR;
int IDUM;
ZEROR=0.0; ZEROI=0.0;
*NZ = 0;
AS1 = ZABS(*S1R,*S1I);
AS2 = ZABS(*S2R,*S2I);
if (*S1R == 0.0 && *S1I == 0.0) goto e10;
if (AS1 == 0.0) goto e10;
ALN = -(*ZRR) -(*ZRR) + log(AS1);
S1DR = *S1R;
S1DI = *S1I;
*S1R = ZEROR;
*S1I = ZEROI;
AS1 = ZEROR;
if (ALN < -ALIM) goto e10;
ZLOG(S1DR, S1DI, &C1R, &C1I, &IDUM);
C1R = C1R - (*ZRR) -(*ZRR);
C1I = C1I - (*ZRI) -(*ZRI);
ZEXP(C1R, C1I, S1R, S1I);
AS1 = ZABS(*S1R,*S1I);
*IUF = *IUF + 1;
e10: AA = DMAX(AS1,AS2);
if (AA > ASCLE) return;
*S1R = ZEROR;
*S1I = ZEROI;
*S2R = ZEROR;
*S2I = ZEROI;
*NZ = 1;
*IUF = 0;
} //ZS1S2()
//for debug only
void test(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
int *NZ, REAL TOL, REAL ELIM, REAL ALIM) {
printf(" %f %f %f %d %d %f %f %d %f %f %f\n", ZR, ZI, FNU, KODE, N,
YR[1], YI[1], *NZ, TOL, ELIM, ALIM);
}
void ZSERI(REAL ZR, REAL ZI, REAL FNU, int KODE, int N, REAL *YR, REAL *YI,
int *NZ, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE ZSERI
!***REFER TO ZBESI,ZBESK
!
! ZSERI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z) >= 0 BY
! MEANS OF THE POWER SERIES FOR LARGE CABS(Z) IN THE
! REGION CABS(Z) <= 2*SQRT(FNU+1). NZ=0 IS A NORMAL RETURN.
! NZ > 0 MEANS THAT THE LAST NZ COMPONENTS WERE SET TO ZERO
! DUE TO UNDERFLOW. NZ < 0 MEANS UNDERFLOW OCCURRED, BUT THE
! CONDITION CABS(Z) <= 2*SQRT(FNU+1) WAS VIOLATED AND THE
! COMPUTATION MUST BE COMPLETED IN ANOTHER ROUTINE WITH N=N-ABS(NZ).
!
!***ROUTINES CALLED DGAMLN,D1MACH,ZUCHK,ZABS,ZDIV,ZLOG,ZMLT
!***END PROLOGUE ZSERI
! COMPLEX AK1,CK,COEF,CONE,CRSC,CSCL,CZ,CZERO,HZ,RZ,S1,S2,Y,Z */
//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90,e100,e120,e140,e150,e160,
// e170,e190
REAL AA, ACZ, AK, AK1I, AK1R, ARM, ASCLE, ATOL,AZ, CKI, CKR, COEFI,
COEFR, CONEI, CONER, CRSCR, CZI, CZR, DFNU, FNUP, HZI, HZR,
RAZ, RS, RTR1, RZI, RZR, S, SS, STI, STR, S1I, S1R, S2I, S2R,
ZEROI, ZEROR;
int I, IB, idebug,IDUM, IFLAG, IL, K, L, M, NN, NW;
REAL *WR, *WI;
void *vmblock = NULL;
FILE *fp;
*NZ=0;
// Initialize WR, WI
vmblock = vminit();
WR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); //index 0 not used
WI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
if (! vmcomplete(vmblock)) {
LogError ("No Memory", 0, __FILE__, __LINE__);
return;
}
ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
idebug=1;
// Note: In C++ version, results are not correct with idebug=0 (any idea?).
if (idebug) fp=fopen("cbess0.txt","w");
AZ = ZABS(ZR,ZI);
if (AZ == 0.0) goto e160;
ARM = 1000*D1MACH(1);
RTR1 = SQRT(ARM);
CRSCR = 1.0;
IFLAG = 0;
if (AZ < ARM) goto e150;
HZR = 0.5*ZR;
HZI = 0.5*ZI;
CZR = ZEROR;
CZI = ZEROI;
if (AZ <= RTR1) goto e10;
ZMLT(HZR, HZI, HZR, HZI, &CZR, &CZI);
e10: ACZ = ZABS(CZR,CZI);
NN = N;
ZLOG(HZR, HZI, &CKR, &CKI, &IDUM);
e20: DFNU = FNU + 1.0*(NN-1);
FNUP = DFNU + 1.0;
/*----------------------------------------------------------------------
! UNDERFLOW TEST
!---------------------------------------------------------------------*/
AK1R = CKR*DFNU;
AK1I = CKI*DFNU;
AK = DGAMLN(FNUP,&IDUM);
AK1R -= AK;
if (KODE == 2) AK1R -= ZR;
if (AK1R > -ELIM) goto e40;
e30: *NZ = *NZ + 1;
YR[NN] = ZEROR;
YI[NN] = ZEROI;
if (ACZ > DFNU) goto e190;
NN--;
if (NN == 0) {
vmfree(vmblock);
return;
}
goto e20;
e40: if (AK1R > -ALIM) goto e50;
IFLAG = 1;
SS = 1.0/TOL;
CRSCR = TOL;
ASCLE = ARM*SS;
e50: AA = EXP(AK1R);
if (IFLAG == 1) AA *= SS;
COEFR = AA*COS(AK1I);
COEFI = AA*SIN(AK1I);
ATOL = TOL*ACZ/FNUP;
IL = IMIN(2,NN);
for (I=1; I<=IL; I++) {
DFNU = FNU + 1.0*(NN-I);
FNUP = DFNU + 1.0;
S1R = CONER;
S1I = CONEI;
if (ACZ < TOL*FNUP) goto e70;
AK1R = CONER;
AK1I = CONEI;
AK = FNUP + 2.0;
S = FNUP;
AA = 2.0;
e60: RS = 1.0/S;
if (idebug) fprintf(fp," ** RS=%e\n", RS);
STR = AK1R*CZR - AK1I*CZI;
STI = AK1R*CZI + AK1I*CZR;
AK1R = STR*RS;
AK1I = STI*RS;
S1R = S1R + AK1R;
S1I = S1I + AK1I;
S = S + AK;
AK = AK + 2.0;
AA = AA*ACZ*RS;
if (AA > ATOL) goto e60;
e70: S2R = S1R*COEFR - S1I*COEFI;
S2I = S1R*COEFI + S1I*COEFR;
WR[I] = S2R;
WI[I] = S2I;
if (IFLAG == 0) goto e80;
ZUCHK(S2R, S2I, &NW, ASCLE, TOL);
if (NW != 0) goto e30;
e80: M = NN - I + 1;
YR[M] = S2R*CRSCR;
YI[M] = S2I*CRSCR;
if (I == IL) goto e90;
ZDIV(COEFR, COEFI, HZR, HZI, &STR, &STI);
COEFR = STR*DFNU;
COEFI = STI*DFNU;
e90: ;} // I loop
if (idebug) fclose(fp);
if (NN <= 2) {
vmfree(vmblock);
return;
}
K = NN - 2;
AK = 1.0*K;
RAZ = 1.0/AZ;
STR = ZR*RAZ;
STI = -ZI*RAZ;
RZR = (STR+STR)*RAZ;
RZI = (STI+STI)*RAZ;
if (IFLAG == 1) goto e120;
IB = 3;
e100: for (I=IB; I<=NN; I++) {
YR[K] = (AK+FNU)*(RZR*YR[K+1]-RZI*YI[K+1]) + YR[K+2];
YI[K] = (AK+FNU)*(RZR*YI[K+1]+RZI*YR[K+1]) + YI[K+2];
AK -= 1.0;
K--;
}
vmfree(vmblock);
return;
/*----------------------------------------------------------------------
! RECUR BACKWARD WITH SCALED VALUES
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
! EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION ABOVE THE
! UNDERFLOW LIMIT = ASCLE = D1MACH(1)*SS*1000.
!---------------------------------------------------------------------*/
e120: S1R = WR[1];
S1I = WI[1];
S2R = WR[2];
S2I = WI[2];
for (L=3; L<=NN; L++) {
CKR = S2R;
CKI = S2I;
S2R = S1R + (AK+FNU)*(RZR*CKR-RZI*CKI);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -