INLINE P YgooSmathYPfsinh(P x) {
  INTFLO iz, ix; ix.i = (PINT)x; 
  iz.f = (float)sinh((double)ix.f);
  return (P)iz.i;
}

inline double _sinh(double arg)            { return  sinh(arg); }

float
sinhf (float x)
{
  return (float) sinh ((double) x);
}

double Formulaeditor::factor(qint32& nPosition, QString& strCharacter)
{
    qreal f = 0.0;
    qint32 wI = 0, wL = 0, wBeginn = 0, wError = 0;

    if	(strCharacter == str_char(0)) return 0.0;
    // read digit and save as float in f
    if (((strCharacter >= "0") && (strCharacter <= "9")) || (strCharacter == "."))
    {
        wBeginn = nPosition;

        do
        {
            char_n(nPosition, strCharacter);
        }
        while ((((strCharacter >= "0") && (strCharacter <= "9")) || (strCharacter == ".")));

        if (strCharacter == ".")
        {
            do
            {
                char_n(nPosition, strCharacter);
            }
            while (!(((qint8)strCharacter.at(0).digitValue() >= 0) && ((qint8)strCharacter.at(0).digitValue() <=  9))  || (strCharacter.at(0) == '.'));
        }

        QString g_strF = m_strFunction.mid(wBeginn - 1, nPosition - wBeginn);
        f = g_strF.toFloat();
    }
    else
    {
        QString strCharacterUpper = strCharacter.toUpper();
        if (strCharacter == "(")
        {
            char_n(nPosition, strCharacter);
            f = expression(nPosition, strCharacter);
            if (strCharacter == ")")
                char_n(nPosition, strCharacter);
        }
        else if (strCharacterUpper == "X")
        {
            char_n(nPosition, strCharacter);
            f = m_dFktValue;
        }
        else
        {
            bool gefunden = false;
            qint32 AnzStdFunctions = m_strStandardFunction.length() - 1;
            for (wI = 1; wI <= AnzStdFunctions; wI++)
            {
                wL = m_strStandardFunction.at(wI).length();
                QString strFunktionUpper = m_strFunction.mid(nPosition - 1, wL);
                strFunktionUpper = strFunktionUpper.toUpper();
                QString strDummy(m_strStandardFunction.at(wI));
                strDummy = strDummy.toUpper();
                if (strFunktionUpper == strDummy)
                {
                    gefunden = true;
                    nPosition = nPosition + wL - 1;
                    char_n(nPosition, strCharacter);
                    // ! recursion !!!!!!!!!!!!!!!!!!!!!!
                    f = factor(nPosition, strCharacter);
                    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    if (strFunktionUpper == "ABS")
                        f = fabs(f);
                    else if (strFunktionUpper == "SQRT")
                        if (f >= 0)
                            f = sqrt(f);
                        else
                            wError = -1;
                    else if (strFunktionUpper == "SINH")
                        f = sinh(f);
                    else if (strFunktionUpper == "COSH")
                        f = cosh(f);
                    else if (strFunktionUpper == "TANH")
                        f = tanh(f);
                    else if (strFunktionUpper == "ARCTAN")
                        f = atan(f);
                    else if (strFunktionUpper == "LN")
                    {
                        if (f >= 0)
                            f = log(f);
                        else
                            wError = -1;
                    }
                    else if (strFunktionUpper == "LOG")
                    {
                        if (f >= 0)
                            f = log10(f);
                        else
                            wError = -1;
                    }
                    else if (strFunktionUpper == "EXP")
                    {
                        //if (f <= 41)
                            f = exp(f);
                        //else
                            //wError = -1;
                    }
                    else if (strFunktionUpper == "SIN")
//.........这里部分代码省略.........

void c_qd_sinh(const double *a, double *b) {
  qd_real bb;
  bb = sinh(qd_real(a));
  TO_DOUBLE_PTR(bb, b);
}

static VALUE
math_sinh(VALUE obj, VALUE x)
{
    Need_Float(x);
    return DBL2NUM(sinh(RFLOAT_VALUE(x)));
}

static int math_sinh (lua_State *L) {
  lua_pushnumber(L, sinh(luaL_checknumber(L, 1)));
  return 1;
}

//.........这里部分代码省略.........
        //            yParticle = fRand->Uniform( fYmin+yStringShift, fYmax+yStringShift ); //use it to "shuffle" particles in y!
        //        yParticle += fRand->Gaus(0,1);
        //        yParticle = fRand->Uniform( -yStringSize/2, yStringSize/2 );
        yParticle = fRand->Uniform( -yStringSize/2+yStringShift, yStringSize/2+yStringShift );



        // to get pions = 0.8 when half of them goes from rho decays:
        // probabilities for rho, pions, kaons+protons should be 0.25, 0.5, 0.25
        //        k=1 //ratio of pions from rho-s to pions from string
        //        a=0.8 //ratio of final state pions to all charged
        //        x=(a*k)/(2*k+2-a*k)
        //        y=(a*k)/(2*k+2-a*k)*2/k
        //        z=1-x-y
        //        check: (2*x+y)/(2*x+y+z)

        // TMP: just assign mass for particle. Todo?: separate mechanisms for mesons and proton (?)
        double particleMass = mPion;
        // ##### tune particle ratios!
        if (0) // only rho mesons!!!
        {
            while ( fabs(particleMass-mRho) > mRhoWidth/2 )
                particleMass = fRand->Gaus(mRho,mRhoWidth/2);//( fRand->Uniform(0,1) > 0.5 ? fRand->Gaus(mRho,mRhoWidth/2) : mPion );
        }
        //        else //if (0)
        //        {
        //            double probPID = fRand->Uniform(0,1);
        //            if ( probPID < 0.0 )
        //                particleMass = fRand->Gaus(mRho,mRhoWidth/2);//( fRand->Uniform(0,1) > 0.5 ? fRand->Gaus(mRho,mRhoWidth/2) : mPion );
        //            else if ( probPID < 0.75 )
        //                particleMass = mPion;
        //            else //kaon or proton (13% and 4%)
        //            {
        //                if ( fRand->Uniform(0,1) < 13./(13+4) )
        //                    particleMass = mKaon;
        //                else
        //                    particleMass = mProton;
        //            }
        //        }
        else
        {
            short q1 = breakPointType[iBreak];
            short q2 = breakPointType[iBreak+1];
            if ( q1 == 0 && q2 == 0 ) // u/d quarks => pions
                particleMass = mPion;
            else if ( (q1 == 0 && q2 == 1)
                      || (q1 == 1 && q2 == 0) ) // u/d and s quarks => kaons
                particleMass = mKaon;
            else if ( q1 == 1 && q2 == 1 ) // two s quarks => phi
                particleMass = mPhi;
            else if ( (q1 == 0 && q2 == 2)
                      || (q1 == 2 && q2 == 0) ) // u/d quarks and diquark  => protons/neutrons (!!! also neutrons!)
                particleMass = mProton;
            else if ( (q1 == 1 && q2 == 2)
                      || (q1 == 2 && q2 == 1) ) // s quark and diquark  => Lambda
                particleMass = mLambda;
            else if ( q1 == 3 || q2 == 3 ) // KOSTYL': if at least one of the two string ends is a c-quark
                //(q1 == 0 && q2 == 3)
//                      || (q1 == 3 && q2 == 0) ) // u/d quarks and c quark  => D-meson
                particleMass = mD0;
            else
            {
                cout << "breakPointTypes: impossible configuration! "
                     << q1 << " and " <<  q2 << endl;
                particleMass = mLambda;
            }
        }



        //        cout << particleMass << endl;

        //        ptParticle = fRand->Exp(0.45);

        // prepare lorentz vector
        if (1) //use direct sampling of pt "boltzman" distr
        {
            if ( fabs( particleMass-mPion) < 0.001 )
                ptParticle = funcPtBoltzmanLikePion->GetRandom();
            else if ( fabs( particleMass-mKaon) < 0.001 )
                ptParticle = funcPtBoltzmanLikeKaon->GetRandom();
            else if ( fabs( particleMass-mProton) < 0.001 )
                ptParticle = funcPtBoltzmanLikeProton->GetRandom();
            else if ( fabs( particleMass-mD0) < 0.001 )
                ptParticle = funcPtBoltzmanLikeDmeson->GetRandom();
        }
        if (0)
            phiParticle = fRand->Uniform( 0, TMath::TwoPi() );

        double mT = sqrt( ptParticle*ptParticle + particleMass*particleMass );
        double pX = ptParticle * cos(phiParticle);
        double pY = ptParticle * sin(phiParticle);
        double pZ = mT*sinh(yParticle);
        vArr[iBreak].SetXYZM( pX, pY, pZ, particleMass );

    }

    return nParticlesInString;

}

//function to fit
double fun_fit(double Z2,double M,int t)
{
  if(parity==1) return Z2*exp(-M*TH)*cosh(M*(TH-t))/sinh(M);
  else return Z2*exp(-M*TH)*sin(M*(TH-t))/sinh(M);
}

Constante& Rationnel::sinush()const
{
    Reel* res = new Reel(sinh((float)num/den));
    return *res;
}

double HHVM_FUNCTION(sinh, double arg) { return sinh(arg); }

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

printf("%f\n", sin(0.12));
printf("%f\n", cos(0.12));
printf("%f\n", tan(0.12));
printf("%f\n", asin(0.12));
printf("%f\n", acos(0.12));
printf("%f\n", atan(0.12));
printf("%f\n", sinh(0.12));
printf("%f\n", cosh(0.12));
printf("%f\n", tanh(0.12));
printf("%f\n", exp(0.12));
printf("%f\n", fabs(-0.12));
printf("%f\n", log(0.12));
printf("%f\n", log10(0.12));
printf("%f\n", pow(0.12, 0.12));
printf("%f\n", sqrt(0.12));
printf("%f\n", round(12.34));
printf("%f\n", ceil(12.34));
printf("%f\n", floor(12.34));

void main() {}

double r_sinh(real *x)
#endif
{
return( sinh(*x) );
}

//.........这里部分代码省略.........
    double nCNa = ((nConIntraNa/ConExtraNa)*exp(nCurlyPhi)- 1)/
      (exp(nCurlyPhi)- 1);
    double ngNaNMDA = ((nPK*F*F)/(RT))*ConExtraNa; 
    double nINaNMDA = 1.e-3 * ngNaNMDA * nA2 * nB2 * nCNa * nVm; 
          
    double nCCa = ((nConIntraCa/ConExtraCa)*exp(2*nCurlyPhi)- 1)/
      (exp(2*nCurlyPhi)- 1);
    double ngCaNMDA = ((4*6*nPK*F*F)/(RT))*ConExtraCa; 
    double nICaNMDA = 1.e-3 * ngCaNMDA * nA2 * nB2 * nCCa * nVm; 
    //nIKNMDA*=bNMDA;
    //nINaNMDA*=bNMDA;
    //nICaNMDA*=bNMDA;
     
    double aConGLU=ConExtraglu;
          
    double aA2 = 1100*aConGLU/(1100*aConGLU + 190) ; 
    double aCK = ((aConIntraK/ConExtraK)*exp(aCurlyPhi)- 1)/
      (exp(aCurlyPhi)- 1);

    double agKAMPA = ((aPK*F*F)/(RT))*ConExtraK ; 
    double aIKAMPA = 1.e-3 * agKAMPA * aA2 * aCK * aVm; 
              
    double aCNa = ((aConIntraNa/ConExtraNa)*exp(aCurlyPhi)- 1)/
      (exp(aCurlyPhi)- 1);
    double agNaAMPA = ((aPK*F*F)/(RT))*ConExtraNa; 
    double aINaAMPA = 1.e-3 * agNaAMPA * aA2 * aCNa * aVm; 
    //aIKAMPA*=bAMPA;
    //aINaAMPA*= bAMPA;
    //courants de la pompe Na+/K+  
          
    double nCapitalPhi = FRT*(nVm + 176.5); 
    double nA3 = pow(ConExtraK/(ConExtraK + 3.7),2) ; 
    double nB3 = pow(nConIntraNa/(nConIntraNa + 0.6),3) ; 
    double nC3 = (0.052*sinh(nCapitalPhi))/(0.026*exp(nCapitalPhi) + 
  					  22.5*exp(-nCapitalPhi)); 
    double nIKpompeKNa = -0.01 * nrNaK * depol * nA3 * nB3 * nC3 ; 
          
    //double nINapompeKNa = (3/2)*(-10)^(-2) * nrNaK * depol * nA3 * nB3 * nC3; 
    double nINapompeKNa = (1.5)*0.01 * nrNaK * depol * nA3 * nB3 * nC3;

    double aCapitalPhi = FRT*(aVm + 176.5); 
    double aA3 = pow(ConExtraK/(ConExtraK + 3.7),2) ; 
    double aB3 = pow(aConIntraNa/(aConIntraNa + 0.6),3) ; 
    double aC3 = (0.052*sinh(aCapitalPhi))/(0.026*exp(aCapitalPhi) +
  					    22.5*exp(-aCapitalPhi)); 
    double aIKpompeKNa = -0.01 * arNaK * depol * aA3 * aB3 * aC3 ; 
  
    double aINapompeKNa = (1.5)*0.01  * arNaK * depol * aA3 * aB3 * aC3; 

    //courant de la pompe Ca2+
          
    double nICapompe = depol*ngCapompe*0.01*nConIntraCa/(nConIntraCa+0.0002);
    double aICapompe = depol*agCapompe*0.01*aConIntraCa/(aConIntraCa+0.0002);

    //courant de la pompe Cl-
           
    double nIClpompe=-0.01 * depol *ngClpompe * nConIntraCl/(nConIntraCl+25);
    double aIClpompe=-0.01 * depol *agClpompe * aConIntraCl/(aConIntraCl+25);  


    //courants de l'antiport Na+/Ca2+   
          
    double aa=F/(2*RT);
 
    double nA4 = pow(nConIntraNa,3)*ConExtraCa*exp(aa*nVm)
      - pow(ConExtraNa,3)*nConIntraCa*exp(-aa*nVm); 

void ctss_calculate_biquad_coeff(CTSS_DSPNode *node,
                                 CTSS_BiquadType type,
                                 float freq,
                                 float dbGain,
                                 float bandwidth) {
  CTSS_BiquadState *state = node->state;
  float a0, a1, a2, b0, b1, b2;

  float A     = powf(10.0f, dbGain / 40.0f);
  float omega = HZ_TO_RAD(freq);
  float sn    = sinf(omega);
  float cs    = cosf(omega);
  float alpha = sn * sinh(CT_LN2 / 2.0f * bandwidth * omega / sn);
  float beta  = sqrtf(A + A);

  switch (type) {
    case LPF:
    default:
      b0 = (1.0f - cs) / 2.0f;
      b1 = 1.0f - cs;
      b2 = (1.0f - cs) / 2.0f;
      a0 = 1.0f + alpha;
      a1 = -2.0f * cs;
      a2 = 1.0f - alpha;
      break;
    case HPF:
      b0 = (1.0f + cs) / 2.0f;
      b1 = -(1.0f + cs);
      b2 = (1.0f + cs) / 2.0f;
      a0 = 1.0f + alpha;
      a1 = -2.0f * cs;
      a2 = 1.0f - alpha;
      break;
    case BPF:
      b0 = alpha;
      b1 = 0;
      b2 = -alpha;
      a0 = 1.0f + alpha;
      a1 = -2.0f * cs;
      a2 = 1.0f - alpha;
      break;
    case NOTCH:
      b0 = 1.0f;
      b1 = -2.0f * cs;
      b2 = 1.0f;
      a0 = 1.0f + alpha;
      a1 = -2.0f * cs;
      a2 = 1.0f - alpha;
      break;
    case PEQ:
      b0 = 1.0f + (alpha * A);
      b1 = -2.0f * cs;
      b2 = 1.0f - (alpha * A);
      a0 = 1.0f + (alpha / A);
      a1 = -2.0f * cs;
      a2 = 1.0f - (alpha / A);
      break;
    case LSH:
      b0 = A * ((A + 1.0f) - (A - 1.0f) * cs + beta * sn);
      b1 = 2.0f * A * ((A - 1.0f) - (A + 1.0f) * cs);
      b2 = A * ((A + 1.0f) - (A - 1.0f) * cs - beta * sn);
      a0 = (A + 1.0f) + (A - 1.0f) * cs + beta * sn;
      a1 = -2.0f * ((A - 1.0f) + (A + 1.0f) * cs);
      a2 = (A + 1.0f) + (A - 1.0f) * cs - beta * sn;
      break;
    case HSH:
      b0 = A * ((A + 1.0f) + (A - 1.0f) * cs + beta * sn);
      b1 = -2.0f * A * ((A - 1.0f) + (A + 1.0f) * cs);
      b2 = A * ((A + 1.0f) + (A - 1.0f) * cs - beta * sn);
      a0 = (A + 1.0f) - (A - 1.0f) * cs + beta * sn;
      a1 = 2.0f * ((A - 1.0f) - (A + 1.0f) * cs);
      a2 = (A + 1.0f) - (A - 1.0f) * cs - beta * sn;
      break;
  }

  a0 = 1.0f / a0;

  state->f[0] = b0 * a0;
  state->f[1] = b1 * a0;
  state->f[2] = b2 * a0;
  state->f[3] = a1 * a0;
  state->f[4] = a2 * a0;

  state->f[5] = state->f[6] = state->f[7] = state->f[8] = 0.0f;
}

//.........这里部分代码省略.........
  
  //fit each combo
  int ic=0;
  for(int ims=0;ims<nmass;ims++)
    for(int imc=ims;imc<nmass;imc++)
      {
	//take into account corr
	jvec corr1(T,njack),corr2(T,njack),corr;
	int ic1=0+2*(ims+nmass*(0+2*imc));
	int ic2=1+2*(ims+nmass*(1+2*imc));
	cout<<ims<<" "<<imc<<" "<<ic1<<" "<<ic2<<endl;
	corr1.put(buf+ic1*T*(njack+1));
	corr2.put(buf+ic2*T*(njack+1));
	corr=(corr1+corr2)/2;
	cout<<corr1[0]<<" "<<corr2[0]<<endl;
	
	//choose the index of the fitting interval
	if(ims>=nlights) ifit_int=2;
	else
	  if(imc>=nlights) ifit_int=1;
	  else ifit_int=0;
	
	//simmetrize
	corr=corr.simmetrized(parity);
	int ttmin=tmin[ifit_int];
	int ttmax=tmax[ifit_int];
	jvec Mcor=effective_mass(corr),Z2cor(TH+1,njack);
	jack Meff=constant_fit(Mcor,ttmin,ttmax);
	for(int t=0;t<=TH;t++)
	  for(int ijack=0;ijack<=njack;ijack++)
	    Z2cor[t].data[ijack]=corr[t].data[ijack]/fun_fit(1,Meff[ijack],t);
	jack Z2eff=constant_fit(Z2cor,ttmin,ttmax);
	
	if(!isnan(Z2eff[0])) minu.DefineParameter(0,"Z2",Z2eff[0],Z2eff.err(),0,2*Z2eff[0]);
	if(!isnan(Meff[0])) minu.DefineParameter(1,"M",Meff[0],Meff.err(),0,2*Meff[0]);
	for(int t=tmin[ifit_int];t<=tmax[ifit_int];t++) corr_err[t]=corr.data[t].err();
	
	//jacknife analysis
	for(int ijack=0;ijack<njack+1;ijack++)
	  {
	    //copy data so that glob function may access it
	    for(int t=tmin[ifit_int];t<=tmax[ifit_int];t++) corr_fit[t]=corr.data[t].data[ijack];
	  
	    //fit
	    double dum;
	    minu.Migrad();	    
	    minu.GetParameter(0,Z2.data[ic].data[ijack],dum);
	    minu.GetParameter(1,M.data[ic].data[ijack],dum);
	  }
	
	//if((ims==iml_un||ims==nlights-1||ims==nlights||ims==nmass-1)&&
	//(imc==iml_un||imc==nlights-1||imc==nlights||imc==nmass-1))
	  {
	    //plot eff mass
	    {
	      ofstream out(combine("eff_mass_plot_%02d_%02d.xmg",ims,imc).c_str());
	      out<<"@type xydy"<<endl;
	      out<<"@s0 line type 0"<<endl;
	      out<<Mcor<<endl;
	      out<<"&"<<endl;
	      out<<"@type xy"<<endl;
	      double av_mass=M[ic].med();
	      double er_mass=M[ic].err();
	      out<<tmin[ifit_int]<<" "<<av_mass-er_mass<<endl;
	      out<<tmax[ifit_int]<<" "<<av_mass-er_mass<<endl;
	      out<<tmax[ifit_int]<<" "<<av_mass+er_mass<<endl;
	      out<<tmin[ifit_int]<<" "<<av_mass+er_mass<<endl;
	      out<<tmin[ifit_int]<<" "<<av_mass-er_mass<<endl;
	    }
	    //plot fun
	    {
	      ofstream out(combine("fun_plot_%02d_%02d.xmg",ims,imc).c_str());
	      out<<"@type xydy"<<endl;
	      out<<"@s0 line type 0"<<endl;
	      out<<corr<<endl;
	      out<<"&"<<endl;
	      out<<"@type xy"<<endl;
	      for(int t=tmin[ifit_int];t<tmax[ifit_int];t++)
		out<<t<<" "<<fun_fit(Z2[ic][njack],M[ic][njack],t)<<endl;
	    }
	  }
	
	  cout<<mass[ims]<<" "<<mass[imc]<<"  "<<M[ic]<<" "<<Z2[ic]<<" f:"<<sqrt(Z2[ic])/(sinh(M[ic])*M[ic])*(mass[ims]+mass[imc])<<" fV:"<<Za_med[ibeta]*sqrt(Z2[ic])/M[ic]/lat[ibeta].med()<<endl;
	  ic++;
      }
  
  ofstream out("fitted_mass.xmg");
  out<<"@type xydy"<<endl;
  for(int ims=0;ims<nmass;ims++)
    {
      //out<<"s0 line type 0"<<endl;
      for(int imc=0;imc<nmass;imc++) out<<mass[imc]<<" "<<M[icombo(ims,imc,nmass,nlights,0)]<<endl;
      out<<"&"<<endl;
    }

  M.write_to_binfile(out_file);
  Z2.append_to_binfile(out_file);
  
  return 0;
}

void MathSinh(struct ParseState *Parser, struct Value *ReturnValue, struct Value **Param, int NumArgs)
{
    ReturnValue->Val->FP = sinh(Param[0]->Val->FP);
}

void bessjy(double  x, double xnu, double *rj, double *ry)
{
    void beschb(double x, double *gam1, double *gam2, double *gampl, double
*gammi);

int i,isign,l,nl;
double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2,
              fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl,
              rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1,
              temp,w,x2,xi,xi2,xmu,xmu2;
if (x <= 0.0 || xnu<0.0) nrerror("bad argumnets in bessjy");
nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu+1.5))); 
xmu=xnu-nl;
xmu2=xmu*xmu;
xi=1.0/x;
xi2=2.0*xi;
w=xi2/M_PI;
isign=1;
h=xnu*xi;
if(h<FPMIN) h=FPMIN;
b=xi2*xnu;
d=0.0;
c=h;
for (i=1;i<=MAXIT;i++)
{
 b+=xi2;
 d=b-d;
 if(fabs(d)< FPMIN) d=FPMIN;
 c=b-1.0/c;
 if(fabs(c)<FPMIN) c=FPMIN;
 d=1.0/d;
 del=c*d;
 h=del*h;
 if(d<0.0)  isign=-isign;
 if(fabs(del-1.0)<EPS) break;
}
if(i>MAXIT) nrerror("x too large in bessjy; try asymptotic expansion");
rjl=isign*FPMIN;
rjpl=h*rjl;
rjl1=rjl;
rjp1=rjpl;
fact=xnu*xi;
for (l=nl;l>=1;l--){
rjtemp=fact*rjl+rjpl;
fact-=xi;
rjpl=fact*rjtemp-rjl;
rjl=rjtemp;
}
if(rjl ==0.0) rjl=EPS;
f=rjpl/rjl;
if (x<XMIN) {
x2=0.5*x;
pimu=M_PI*xmu;
fact=(fabs(pimu)<EPS ? 1.0 : pimu/sin(pimu));
d=-log(x2);
e=xmu*d;
fact2=(fabs(e)< EPS ? 1.0 : sinh(e)/e);
beschb(xmu,&gam1,&gam2,&gampl,&gammi);
ff=2.0/M_PI*fact*(gam1*cosh(e)+gam2*fact2*d);
e=exp(e);
p=e/(gampl*M_PI);
q=1.0/(e*M_PI*gammi);
pimu2=0.5*pimu;
fact3=(fabs(pimu2) <EPS ? 1.0 : sin(pimu2)/pimu2);
r=M_PI*pimu2*fact3*fact3;
c=1.0;
d=-x2*x2;
sum=ff+r*q;
sum1=p;
for (i=1;i<=MAXIT;i++)
{ff=(i*ff+p+q)/(i*i-xmu2);
c *=(d/i);
p/= (i-xmu);
q/=(i+xmu);
del=c*(ff+r*q);
sum+=del;
del1=c*p-i*del;
sum1 +=del1;
if(fabs(del) < (1.0+fabs(sum))*EPS) break;
}
if(i > MAXIT) nrerror("bessy series failed to converge");
rymu=-sum;
ry1=-sum1*xi2;
rymup=xmu*xi*rymu-ry1;
rjmu=w/(rymup-f*rymu);
} else {
a=0.25-xmu2;
p=-0.5*xi;
q=1.0;
br=2.0*x;
bi=2.0;
fact=a*xi/(p*p+q*q);
cr=br+q*fact;
ci=bi+p*fact;
den=br*br+bi*bi;
dr=br/den;
di=-bi/den;
dlr=cr*dr-ci*di;
dli=cr*di+ci*dr;
temp=p*dlr-q*dli;
//.........这里部分代码省略.........

double
tanh(double x)
{
    return sinh(x) / cosh(x);
}

int decode_function(char *ptr,int argnum,double *argptr,
                    double *value,void *data) {


  int i; 
 
  for (i=0;fnlist[i] !=0;i++) if (strcmp(ptr,fnlist[i])==0) break;

  switch(i) {
    case 0:
      *value=fabs(argptr[0]);
      break;
    case 1: 
      *value=acos(argptr[0]);
      break;
    case 2: 
      *value=acosh(argptr[0]);
      break;
    case 3: 
      *value=asin(argptr[0]);
      break;
    case 4: 
      *value=asinh(argptr[0]);
      break;
    case 5: 
      *value=atan(argptr[0]);
      break;
    case 6: 
      *value=atan2(argptr[0],argptr[1]);
      break;
    case 7: 
      *value=atanh(argptr[0]);
      break;
   case 8: 
      *value=ceil(argptr[0]);
      break;
   case 9: 
      *value=cos(argptr[0]);
      break;
  case 10: 
      *value=cosh(argptr[0]);
      break;
  case 11: 
      *value=exp(argptr[0]);
      break;
  case 12: 
      *value=floor(argptr[0]);
      break;
  case 13: 
      *value=(int) argptr[0];
      break;
  case 14: 
      *value=log(argptr[0]);
      break;
  case 15: 
      *value=log10(argptr[0]);
      break;
  case 16: 
      *value=sin(argptr[0]);
      break;
  case 17: 
      *value=sinh(argptr[0]);
      break;
  case 18: 
      *value=sqrt(argptr[0]);
      break;
  case 19: 
      *value=tan(argptr[0]);
      break;
  case 20: 
      *value=tanh(argptr[0]);
      break;
      
  }
  return 0;
}