📄 cbess3_cpp.txt
字号:
*NZ = 0;
FMR = 1.0*MR;
SGN = -SIGN(PI,FMR);
/*----------------------------------------------------------------------
! CSPN AND CSGN ARE COEFF OF K AND I FUNCIONS RESP.
!---------------------------------------------------------------------*/
CSGNI = SGN;
if (YY <= 0.0) CSGNI = -CSGNI;
IFN = INU + N - 1;
ANG = FNF*SGN;
CSPNR = COS(ANG);
CSPNI = SIN(ANG);
if ((IFN % 2) == 0) goto e190;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
/*----------------------------------------------------------------------
! CS=COEFF OF THE J FUNCTION TO GET THE I FUNCTION. I(FNU,Z) IS
! COMPUTED FROM EXP(I*FNU*HPI)*J(FNU,-I*Z) WHERE Z IS IN THE FIRST
! QUADRANT. FOURTH QUADRANT VALUES (YY.LE.0.0E0) ARE COMPUTED BY
! CONJUGATION SINCE THE I FUNCTION IS REAL ON THE POSITIVE REAL AXIS
!---------------------------------------------------------------------*/
e190: CSR = SAR*CSGNI;
CSI = CAR*CSGNI;
IN0 = (IFN % 4) + 1;
C2R = CIPR[IN0];
C2I = CIPI[IN0];
STR = CSR*C2R + CSI*C2I;
CSI = -CSR*C2I + CSI*C2R;
CSR = STR;
ASC = BRY[1];
IUF = 0;
KK = N;
KDFLG = 1;
IB = IB - 1;
IC = IB - 1;
for (K=1; K<=N; K++) {
FN = FNU + 1.0*(KK-1);
/*----------------------------------------------------------------------
! LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
! FUNCTION ABOVE
!---------------------------------------------------------------------*/
if (N > 2) goto e175;
e172: PHIDR = PHIR[J];
PHIDI = PHII[J];
ARGDR = ARGR[J];
ARGDI = ARGI[J];
ZET1DR = ZETA1R[J];
ZET1DI = ZETA1I[J];
ZET2DR = ZETA2R[J];
ZET2DI = ZETA2I[J];
ASUMDR = ASUMR[J];
ASUMDI = ASUMI[J];
BSUMDR = BSUMR[J];
BSUMDI = BSUMI[J];
J = 3 - J;
goto e210;
e175: if (KK == N && IB < N) goto e210;
if (KK == IB || KK == IC) goto e172;
ZUNHJ(ZNR, ZNI, FN, 0, TOL, &PHIDR, &PHIDI, &ARGDR,
&ARGDI, &ZET1DR, &ZET1DI, &ZET2DR, &ZET2DI, &ASUMDR,
&ASUMDI, &BSUMDR, &BSUMDI);
e210: if (KODE == 1) goto e220;
STR = ZBR + ZET2DR;
STI = ZBI + ZET2DI;
RAST = FN/ZABS(STR,STI);
STR = STR*RAST*RAST;
STI = -STI*RAST*RAST;
S1R = -ZET1DR + STR;
S1I = -ZET1DI + STI;
goto e230;
e220: S1R = -ZET1DR + ZET2DR;
S1I = -ZET1DI + ZET2DI;
/*----------------------------------------------------------------------
! TEST FOR UNDERFLOW AND OVERFLOW
!---------------------------------------------------------------------*/
e230: RS1 = S1R;
if (ABS(RS1) > ELIM) goto e280;
if (KDFLG == 1) IFLAG = 2;
if (ABS(RS1) < ALIM) goto e240;
/*----------------------------------------------------------------------
! REFINE TEST AND SCALE
!---------------------------------------------------------------------*/
APHI = ZABS(PHIDR,PHIDI);
AARG = ZABS(ARGDR,ARGDI);
RS1 = RS1 + log(APHI) - 0.25*log(AARG) - AIC;
if (ABS(RS1) > ELIM) goto e280;
if (KDFLG == 1) IFLAG = 1;
if (RS1 < 0.0) goto e240;
if (KDFLG == 1) IFLAG = 3;
e240: ZAIRY(ARGDR, ARGDI, 0, 2, &AIR, &AII, &NAI, &IDUM);
ZAIRY(ARGDR, ARGDI, 1, 2, &DAIR, &DAII, &NDAI, &IDUM);
STR = DAIR*BSUMDR - DAII*BSUMDI;
STI = DAIR*BSUMDI + DAII*BSUMDR;
STR = STR + (AIR*ASUMDR-AII*ASUMDI);
STI = STI + (AIR*ASUMDI+AII*ASUMDR);
PTR = STR*PHIDR - STI*PHIDI;
PTI = STR*PHIDI + STI*PHIDR;
S2R = PTR*CSR - PTI*CSI;
S2I = PTR*CSI + PTI*CSR;
STR = EXP(S1R)*CSSR[IFLAG];
S1R = STR*COS(S1I);
S1I = STR*SIN(S1I);
STR = S2R*S1R - S2I*S1I;
S2I = S2R*S1I + S2I*S1R;
S2R = STR;
if (IFLAG != 1) goto e250;
ZUCHK(S2R, S2I, &NW, BRY[1], TOL);
if (NW == 0) goto e250;
S2R = ZEROR;
S2I = ZEROI;
e250: if (YY <= 0.0) S2I = -S2I;
CYR[KDFLG] = S2R;
CYI[KDFLG] = S2I;
C2R = S2R;
C2I = S2I;
S2R = S2R*CSRR[IFLAG];
S2I = S2I*CSRR[IFLAG];
/*----------------------------------------------------------------------
! ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
!---------------------------------------------------------------------*/
S1R = YR[KK];
S1I = YI[KK];
if (KODE == 1) goto e270;
ZS1S2(&ZRR, &ZRI, &S1R, &S1I, &S2R, &S2I, &NW, ASC, ALIM, &IUF);
*NZ = *NZ + NW;
e270: YR[KK] = S1R*CSPNR - S1I*CSPNI + S2R;
YI[KK] = S1R*CSPNI + S1I*CSPNR + S2I;
KK = KK - 1;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
STR = CSI;
CSI = -CSR;
CSR = STR;
if (C2R != 0.0 || C2I != 0.0) goto e255;
KDFLG = 1;
goto e290;
e255: if (KDFLG == 2) goto e295;
KDFLG = 2;
goto e290;
e280: if (RS1 > 0.0) goto e320;
S2R = ZEROR;
S2I = ZEROI;
goto e250;
e290:;} // K loop
K = N;
e295: IL = N - K;
if (IL == 0) {
vmfree(vmblock);
return;
}
/*----------------------------------------------------------------------
! RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
! K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
! INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
!---------------------------------------------------------------------*/
S1R = CYR[1];
S1I = CYI[1];
S2R = CYR[2];
S2I = CYI[2];
CSR = CSRR[IFLAG];
ASCLE = BRY[IFLAG];
FN = 1.0*(INU+IL);
for (I=1; I<=IL; I++) {
C2R = S2R;
C2I = S2I;
S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I);
S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R);
S1R = C2R;
S1I = C2I;
FN = FN - 1.0;
C2R = S2R*CSR;
C2I = S2I*CSR;
CKR = C2R;
CKI = C2I;
C1R = YR[KK];
C1I = YI[KK];
if (KODE == 1) goto e300;
ZS1S2(&ZRR, &ZRI, &C1R, &C1I, &C2R, &C2I, &NW, ASC, ALIM, &IUF);
*NZ = *NZ + NW;
e300: YR[KK] = C1R*CSPNR - C1I*CSPNI + C2R;
YI[KK] = C1R*CSPNI + C1I*CSPNR + C2I;
KK--;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
if (IFLAG >= 3) goto e310;
C2R = ABS(CKR);
C2I = ABS(CKI);
C2M = DMAX(C2R,C2I);
if (C2M <= ASCLE) goto e310;
IFLAG++;
ASCLE = BRY[IFLAG];
S1R = S1R*CSR;
S1I = S1I*CSR;
S2R = CKR;
S2I = CKI;
S1R = S1R*CSSR[IFLAG];
S1I = S1I*CSSR[IFLAG];
S2R = S2R*CSSR[IFLAG];
S2I = S2I*CSSR[IFLAG];
CSR = CSRR[IFLAG];
e310:;} // I loop
vmfree(vmblock);
return;
e320: *NZ = -1;
} //ZUNK2()
void ZACON(REAL ZR, REAL ZI, REAL FNU, int KODE, int MR, int N, REAL *YR, REAL *YI, int *NZ,
REAL RL, REAL FNUL, REAL TOL, REAL ELIM, REAL ALIM) {
/***BEGIN PROLOGUE ZACON
!***REFER TO ZBESK,ZBESH
!
! ZACON APPLIES THE ANALYTIC CONTINUATION FORMULA
!
! K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
! MP=PI*MR*CMPLX(0.0,1.0)
!
! TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
! HALF Z PLANE
!
!***ROUTINES CALLED ZBINU,ZBKNU,ZS1S2,D1MACH,ZABS,ZMLT
!***END PROLOGUE ZACON
! COMPLEX CK,CONE,CSCL,CSCR,CSGN,CSPN,CY,CZERO,C1,C2,RZ,SC1,SC2,ST,
! S1,S2,Y,Z,ZN */
//Labels: e10,e20,e30,e40,e50,e60,e70,e80,e90
REAL ARG, ASCLE, AS2, AZN, BSCLE, CKI, CKR, CONEI, CONER, CPN, CSCL, CSCR,
CSGNI, CSGNR, CSPNI, CSPNR, CSR, C1I, C1M, C1R, C2I, C2R, FMR, FN,
PTI, PTR, RAZN, RZI, RZR, SC1I, SC1R, SC2I, SC2R, SGN, SPN, STI, STR,
S1I, S1R, S2I, S2R, YY, ZEROI, ZEROR, ZNI, ZNR;
int I, INU, IUF, KFLAG, NN, NW;
REAL *CYR, *CYI, *CSSR, *CSRR, *BRY;
void *vmblock = NULL;
//*** First executable statement ZACON
//initialize pointers to vectors
vmblock = vminit();
BRY = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0); //index 0 not used
CYR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
CYI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0);
CSSR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
CSRR = (REAL *) vmalloc(vmblock, VEKTOR, 4, 0);
if (! vmcomplete(vmblock)) {
LogError ("No Memory", 0, __FILE__, __LINE__);
return;
}
ZEROR=0.0; ZEROI=0.0; CONER=1.0; CONEI=0.0;
NZ = 0;
ZNR = -ZR;
ZNI = -ZI;
NN = N;
ZBINU(ZNR, ZNI, FNU, KODE, NN, YR, YI, &NW, RL, FNUL, TOL, ELIM, ALIM);
if (NW < 0) goto e90;
/*----------------------------------------------------------------------
! ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
!---------------------------------------------------------------------*/
NN = IMIN(2,N);
ZBKNU(ZNR, ZNI, FNU, KODE, NN, CYR, CYI, &NW, TOL, ELIM, ALIM);
if (NW != 0) goto e90;
S1R = CYR[1];
S1I = CYI[1];
FMR = 1.0*MR;
SGN = -SIGN(PI,FMR);
CSGNR = ZEROR;
CSGNI = SGN;
if (KODE == 1) goto e10;
YY = -ZNI;
CPN = COS(YY);
SPN = SIN(YY);
ZMLT(CSGNR, CSGNI, CPN, SPN, &CSGNR, &CSGNI);
/*----------------------------------------------------------------------
! CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
! WHEN FNU IS LARGE
!---------------------------------------------------------------------*/
e10: INU = (int) floor(FNU);
ARG = (FNU-1.0*INU)*SGN;
CPN = COS(ARG);
SPN = SIN(ARG);
CSPNR = CPN;
CSPNI = SPN;
if ((INU % 2) == 0) goto e20;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
e20: IUF = 0;
C1R = S1R;
C1I = S1I;
C2R = YR[1];
C2I = YI[1];
ASCLE = 1000.0*D1MACH(1)/TOL;
if (KODE == 1) goto e30;
ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
*NZ = *NZ + NW;
SC1R = C1R;
SC1I = C1I;
e30: ZMLT(CSPNR, CSPNI, C1R, C1I, &STR, &STI);
ZMLT(CSGNR, CSGNI, C2R, C2I, &PTR, &PTI);
YR[1] = STR + PTR;
YI[1] = STI + PTI;
if (N == 1) {
vmfree(vmblock);
return;
}
CSPNR = -CSPNR;
CSPNI = -CSPNI;
S2R = CYR[2];
S2I = CYI[2];
C1R = S2R;
C1I = S2I;
C2R = YR[2];
C2I = YI[2];
if (KODE == 1) goto e40;
ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
*NZ = *NZ + NW;
SC2R = C1R;
SC2I = C1I;
e40: ZMLT(CSPNR, CSPNI, C1R, C1I, &STR, &STI);
ZMLT(CSGNR, CSGNI, C2R, C2I, &PTR, &PTI);
YR[2] = STR + PTR;
YI[2] = STI + PTI;
if (N == 2) {
vmfree(vmblock);
return;
}
CSPNR = -CSPNR;
CSPNI = -CSPNI;
AZN = ZABS(ZNR,ZNI);
RAZN = 1.0/AZN;
STR = ZNR*RAZN;
STI = -ZNI*RAZN;
RZR = (STR+STR)*RAZN;
RZI = (STI+STI)*RAZN;
FN = FNU + 1.0;
CKR = FN*RZR;
CKI = FN*RZI;
/*---------------------------------------------------------------------
! SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
!---------------------------------------------------------------------*/
CSCL = 1.0/TOL;
CSCR = TOL;
CSSR[1] = CSCL;
CSSR[2] = CONER;
CSSR[3] = CSCR;
CSRR[1] = CSCR;
CSRR[2] = CONER;
CSRR[3] = CSCL;
BRY[1] = ASCLE;
BRY[2] = 1.0/ASCLE;
BRY[3] = D1MACH(2);
AS2 = ZABS(S2R,S2I);
KFLAG = 2;
if (AS2 > BRY[1]) goto e50;
KFLAG = 1;
goto e60;
e50: if (AS2 < BRY[2]) goto e60;
KFLAG = 3;
e60: BSCLE = BRY[KFLAG];
S1R = S1R*CSSR[KFLAG];
S1I = S1I*CSSR[KFLAG];
S2R = S2R*CSSR[KFLAG];
S2I = S2I*CSSR[KFLAG];
CSR = CSRR[KFLAG];
for (I=3; I<=N; I++) {
STR = S2R;
STI = S2I;
S2R = CKR*STR - CKI*STI + S1R;
S2I = CKR*STI + CKI*STR + S1I;
S1R = STR;
S1I = STI;
C1R = S2R*CSR;
C1I = S2I*CSR;
STR = C1R;
STI = C1I;
C2R = YR[I];
C2I = YI[I];
if (KODE == 1) goto e70;
if (IUF < 0) goto e70;
ZS1S2(&ZNR, &ZNI, &C1R, &C1I, &C2R, &C2I, &NW, ASCLE, ALIM, &IUF);
*NZ = *NZ + NW;
SC1R = SC2R;
SC1I = SC2I;
SC2R = C1R;
SC2I = C1I;
if (IUF != 3) goto e70;
IUF = -4;
S1R = SC1R*CSSR[KFLAG];
S1I = SC1I*CSSR[KFLAG];
S2R = SC2R*CSSR[KFLAG];
S2I = SC2I*CSSR[KFLAG];
STR = SC2R;
STI = SC2I;
e70: PTR = CSPNR*C1R - CSPNI*C1I;
PTI = CSPNR*C1I + CSPNI*C1R;
YR[I] = PTR + CSGNR*C2R - CSGNI*C2I;
YI[I] = PTI + CSGNR*C2I + CSGNI*C2R;
CKR = CKR + RZR;
CKI = CKI + RZI;
CSPNR = -CSPNR;
CSPNI = -CSPNI;
if (KFLAG >= 3) goto e80;
PTR = ABS(C1R);
PTI = ABS(C1I);
C1M = DMAX(PTR,PTI);
if (C1M <= BSCLE) goto e80;
KFLAG++;
BSCLE = BRY[KFLAG];
S1R = S1R*CSR;
S1I = S1I*CSR;
S2R = STR;
S2I = STI;
S1R = S1R*CSSR[KFLAG];
S1I = S1I*CSSR[KFLAG];
S2R = S2R*CSSR[KFLAG];
S2I = S2I*CSSR[KFLAG];
CSR = CSRR[KFLAG];
e80: ;} // I loop
vmfree(vmblock);
return;
e90: *NZ = -1;
if (NW = -2) *NZ=-2;
vmfree(vmblock);
} //ZACON()
//end of file Cbess3.cpp
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -