real ln_sqrtx2y2(const real& x, const real& y) 
                                throw(STD_FKT_OUT_OF_DEF)
// ln( sqrt(x^2+y^2) ) == 0.5*ln(x^2+y^2); Blomquist, 21.11.03;
// Relative error bound: 5.160563E-016;
// Absolute error bound: 2.225075E-308; if x=1 and 0<=y<=b0;
{
    int j,N;
    real a,b,r,r1;
    dotprecision dot;
   
    a = sign(x)<0 ? -x : x;  // a = |x| >= 0;
    b = sign(y)<0 ? -y : y;  // b = |y| >= 0;
    int exa=expo(a), exb=expo(b), ex;
    if (b > a)
    {
	r = a;   a = b;   b = r;
	ex = exa;   exa = exb;   exb = ex;
    }
    // It holds now:   0 <= b <= a 
    if (sign(a)==0)
        cxscthrow(STD_FKT_OUT_OF_DEF
                    ("real ln_sqrtx2y2(const real&, const real&)"));
    if (exa>20) // to avoid overflow by calculating a^2 + b^2
    {  // a>=2^(20):
	j = Interval_Nr(B_lnx2y2_1,21,exa); // j: No. of subinterval
	N = B_lnx2y2_N1[j];    // N: Optimal int value
	if (exb-exa > -25)
	{   // For (exb-exa>-25) we use the complete term:  
            // N*ln(2) + [ln(2^(-N)*a)+0.5*ln(1+(b/a)^2)]
	    b = b/a;  // a > 0
	    b = lnp1(b*b);
	    times2pown(b,-1);  // exact division by 2
	    times2pown(a,-N);
	    r = b + ln(a); // [ ... ] calculated!
	    r += B_lnx2y2_c1[j];
	}
	else { // For (exb-exa<=-25) only two summands!: 
	    times2pown(a,-N);
	    r = ln(a) + B_lnx2y2_c1[j];
	}
    }
    else  // exa<=20 or a<2^(20):
    {     // Now calculation of a^2+b^2 without overflow:
	if (exa<=-20) // to avoid underflow by calculating a^2+b^2
	    if (exa<=-1022) // a in the denormalized range
	    {
		r = b/a;
		r = lnp1(r*r);  times2pown(r,-1); // r: 0.5*ln(1+..)
		times2pown(a,1067);
		r += ln(a);     // [ .... ] ready
		r -= ln2_1067;  // rel. error = 2.459639e-16;
	    }
	    else // MinReal=2^(-1022) <= a < 2^(-20)
	    {    // Calculating the number j of the subinterval:
		j = 20 - Interval_Nr(B_lnx2y2_2,21,exa);
		r = a;  times2pown(r,B_lnx2y2_N1[j]);
		r = ln(r);  // r: ln(2^N*a);
		if (exb-exa > -25) { // calculating the complete term
		    b = b/a;
		    a = lnp1(b*b);
		    times2pown(a,-1);
		    r += a;  // [ ... ] ready now
		}
                // We now have: exb-exa<=-25,  ==>  b/a <= 2^(-24);
		r -= B_lnx2y2_c1[j]; // 0.5*ln(1+(b/a)^2) neglected!
                // relative error = 4.524090e-16 in both cases;
	    }
	else // calculation of a^2+b^2 without overflow or underflow:
	{   // exa>-20  respective  a>=2^(-20):
	    dot = 0;
	    accumulate(dot,a,a);
	    accumulate(dot,b,b);  // dot = a^2+b^2, exact!
	    real s = rnd(dot);    // s = a^2 + b^2, rounded!
	    if (s>=0.25 && s<=1.75)
		if (s>=0.828125 && s<=1.171875)
		{ // Series:
		    if (a==1 && exb<=-28)
		    {
			r = b;  times2pown(r,-1);
			r *= b;
		    }
		    else {
			dot -= 1;
			r = rnd(dot); // r = a^2+b^2-1 rounded!
			r = lnp1(r);
			times2pown(r,-1);
		    }
		}
		else { // Reading dot = a^2+b^2 twice:
		    r = rnd(dot);
		    dot -= r;
		    r1 = rnd(dot); // a^2+b^2 = r+r1, rounded!
		    r1 = lnp1(r1/r);
		    r = ln(r) + r1;
		    times2pown(r,-1); // exact division by 2
		}
	    else { // calculating straight from: 0.5*ln(x^2+y^2)
		r = ln(s);
		times2pown(r,-1);
	    } 
//.........这里部分代码省略.........

/* returns 0 on success, 1 on error with errno set */
static int calc_lon(privpath_t *pp) {
  point_t *pt = pp->pt;
  int n = pp->len;
  int i, j, k, k1;
  int ct[4], dir;
  point_t constraint[2];
  point_t cur;
  point_t off;
  int *pivk = NULL;  /* pivk[n] */
  int *nc = NULL;    /* nc[n]: next corner */
  point_t dk;  /* direction of k-k1 */
  int a, b, c, d;

  SAFE_MALLOC(pivk, n, int);
  SAFE_MALLOC(nc, n, int);

  /* initialize the nc data structure. Point from each point to the
     furthest future point to which it is connected by a vertical or
     horizontal segment. We take advantage of the fact that there is
     always a direction change at 0 (due to the path decomposition
     algorithm). But even if this were not so, there is no harm, as
     in practice, correctness does not depend on the word "furthest"
     above.  */
  k = 0;
  for (i=n-1; i>=0; i--) {
    if (pt[i].x != pt[k].x && pt[i].y != pt[k].y) {
      k = i+1;  /* necessarily i<n-1 in this case */
    }
    nc[i] = k;
  }

  SAFE_MALLOC(pp->lon, n, int);

  /* determine pivot points: for each i, let pivk[i] be the furthest k
     such that all j with i<j<k lie on a line connecting i,k. */
  
  for (i=n-1; i>=0; i--) {
    ct[0] = ct[1] = ct[2] = ct[3] = 0;

    /* keep track of "directions" that have occurred */
    dir = (3+3*(pt[mod(i+1,n)].x-pt[i].x)+(pt[mod(i+1,n)].y-pt[i].y))/2;
    ct[dir]++;

    constraint[0].x = 0;
    constraint[0].y = 0;
    constraint[1].x = 0;
    constraint[1].y = 0;

    /* find the next k such that no straight line from i to k */
    k = nc[i];
    k1 = i;
    while (1) {
      
      dir = (3+3*sign(pt[k].x-pt[k1].x)+sign(pt[k].y-pt[k1].y))/2;
      ct[dir]++;

      /* if all four "directions" have occurred, cut this path */
      if (ct[0] && ct[1] && ct[2] && ct[3]) {
	pivk[i] = k1;
	goto foundk;
      }

      cur.x = pt[k].x - pt[i].x;
      cur.y = pt[k].y - pt[i].y;

      /* see if current constraint is violated */
      if (xprod(constraint[0], cur) < 0 || xprod(constraint[1], cur) > 0) {
	goto constraint_viol;
      }

      /* else, update constraint */
      if (abs(cur.x) <= 1 && abs(cur.y) <= 1) {
	/* no constraint */
      } else {
	off.x = cur.x + ((cur.y>=0 && (cur.y>0 || cur.x<0)) ? 1 : -1);
	off.y = cur.y + ((cur.x<=0 && (cur.x<0 || cur.y<0)) ? 1 : -1);
	if (xprod(constraint[0], off) >= 0) {
	  constraint[0] = off;
	}
	off.x = cur.x + ((cur.y<=0 && (cur.y<0 || cur.x<0)) ? 1 : -1);
	off.y = cur.y + ((cur.x>=0 && (cur.x>0 || cur.y<0)) ? 1 : -1);
	if (xprod(constraint[1], off) <= 0) {
	  constraint[1] = off;
	}
      }	
      k1 = k;
      k = nc[k1];
      if (!cyclic(k,i,k1)) {
	break;
      }
    }
  constraint_viol:
    /* k1 was the last "corner" satisfying the current constraint, and
       k is the first one violating it. We now need to find the last
       point along k1..k which satisfied the constraint. */
    dk.x = sign(pt[k].x-pt[k1].x);
    dk.y = sign(pt[k].y-pt[k1].y);
    cur.x = pt[k1].x - pt[i].x;
    cur.y = pt[k1].y - pt[i].y;
//.........这里部分代码省略.........

int srTDriftSpace::PropagateRadiationSimple_AnalytTreatQuadPhaseTerm(srTSRWRadStructAccessData* pRadAccessData)
{// e in eV; Length in m !!!
	int result = 0;

	SetupPropBufVars_AnalytTreatQuadPhaseTerm(pRadAccessData);
	if(pRadAccessData->Pres != 0) if(result = SetRadRepres(pRadAccessData, 0)) return result;
		
	PropBufVars.PassNo = 1; //Remove quadratic term from the Phase in coord. repres.
	if(result = TraverseRadZXE(pRadAccessData)) return result;

		//testOC09302011
		//if(Length == 34.63) return 0;

	double xStartOld = pRadAccessData->xStart, zStartOld = pRadAccessData->zStart;
	pRadAccessData->xStart = -(pRadAccessData->nx >> 1)*pRadAccessData->xStep;
	pRadAccessData->zStart = -(pRadAccessData->nz >> 1)*pRadAccessData->zStep;
	double xShift = pRadAccessData->xStart - xStartOld, zShift = pRadAccessData->zStart - zStartOld;

	pRadAccessData->xWfrMin += xShift; pRadAccessData->xWfrMax += xShift;
	pRadAccessData->zWfrMin += zShift; pRadAccessData->zWfrMax += zShift;

		pRadAccessData->WfrEdgeCorrShouldBeDone = 0;

	if(result = SetRadRepres(pRadAccessData, 1)) return result; //To angular repres.

	PropBufVars.PassNo = 2; //Loop in angular repres.
	if(result = TraverseRadZXE(pRadAccessData)) return result;
		if(pRadAccessData->UseStartTrToShiftAtChangingRepresToCoord)
		{
			pRadAccessData->xStartTr += xShift;
			pRadAccessData->zStartTr += zShift;
		}

	if(result = SetRadRepres(pRadAccessData, 0)) return result; //Back to coord. repres.

	pRadAccessData->xStart = xStartOld; pRadAccessData->zStart = zStartOld;
		if(pRadAccessData->UseStartTrToShiftAtChangingRepresToCoord)
		{
			pRadAccessData->xStart = pRadAccessData->xStartTr - xShift;
			pRadAccessData->zStart = pRadAccessData->zStartTr - zShift;
		}

	//Change scale
	//double kx = (pRadAccessData->RobsX + Length)/pRadAccessData->RobsX;
	//pRadAccessData->xStart = kx*(pRadAccessData->xStart) - (Length/pRadAccessData->RobsX)*(pRadAccessData->xc);
	//pRadAccessData->xStep *= kx;
	pRadAccessData->xStart = PropBufVars.kx_AnalytTreatQuadPhaseTerm*(pRadAccessData->xStart) - PropBufVars.kxc_AnalytTreatQuadPhaseTerm*(pRadAccessData->xc);
	pRadAccessData->xStep *= PropBufVars.kx_AnalytTreatQuadPhaseTerm;

	//double kz = (pRadAccessData->RobsZ + Length)/pRadAccessData->RobsZ;
	//pRadAccessData->zStart = kz*(pRadAccessData->zStart) - (Length/pRadAccessData->RobsZ)*(pRadAccessData->zc);
	//pRadAccessData->zStep *= kz;
	pRadAccessData->zStart = PropBufVars.kz_AnalytTreatQuadPhaseTerm*(pRadAccessData->zStart) - PropBufVars.kzc_AnalytTreatQuadPhaseTerm*(pRadAccessData->zc);
	pRadAccessData->zStep *= PropBufVars.kz_AnalytTreatQuadPhaseTerm;

	PropBufVars.PassNo = 3; //Add new quadratic term to the Phase in coord. repres.
	if(result = TraverseRadZXE(pRadAccessData)) return result;

	//pRadAccessData->MirrorFieldData(sign(kx), sign(kz));
	pRadAccessData->MirrorFieldData((int)sign(PropBufVars.kx_AnalytTreatQuadPhaseTerm), (int)sign(PropBufVars.kz_AnalytTreatQuadPhaseTerm));
	//if(kx < 0)
	if(PropBufVars.kx_AnalytTreatQuadPhaseTerm < 0)
	{
		double xEnd = pRadAccessData->xStart + (pRadAccessData->nx - 1)*pRadAccessData->xStep; //or nx ???
		pRadAccessData->xStart = xEnd;
		pRadAccessData->xStep *= -1;
	}
	//if(kz < 0)
	if(PropBufVars.kz_AnalytTreatQuadPhaseTerm < 0)
	{
		double zEnd = pRadAccessData->zStart + (pRadAccessData->nz - 1)*pRadAccessData->zStep; //or nz ???
		pRadAccessData->zStart = zEnd;
		pRadAccessData->zStep *= -1;
	}

	pRadAccessData->SetNonZeroWavefrontLimitsToFullRange();
	return result;
}

EXPORT int ridders_method(			
	    double level_set_left, 			// left hand initial value for root finding
	    double level_set_right, 			// right hand initial value for root finding
	    double function_g_left, 			// left hand initial function value for root finding
	    double function_g_right,			// right hand initial function value for root findin
	    double volume_of_fluid, 			// volume of fluid value for which we need the 
							// the corresponding level-set value
	    double &level_set,				// the level-set value corresponding to given
							// volume of fluid value
	    double d_level_set_d_x1, 			// first partial derivative in x1 
							// direction of level-set
	    double d_level_set_d_x2, 			// first partial derivative in x2 
							// direction of level-set 
	    double d_level_set_d_x3, 			// first partial derivative in x3 
							// direction of level-set
	    double vof_2_level_set_tolerance,		// tolerance in the conversion from volume
							// of fluid value to level-set value
	    int number_iterations_ridder,		// maximum number of iterations allowed in the
							// nonlinear root finding algorithm
	    int lower_bound_derivatives			// lower bound for the first partial derivatives
							// to consider it a limiting case of vanishing
							// partial derivatives
			)

{
      /* function definitions */
      double new_mid_point_beta;			// new estimate of the root location based on
							// midpoint search interval
      double function_g_value_beta;			// function value in the new estimate of the
							// root location based on the midpoint of the
							// search interval
      double function_g_value=10E10;			// function value in the new estimate of the
							// root location
      double new_volume_of_fluid;			// volume of fluid at the new estimate of the
							// root location
      double ridders_denominator;			// denominator of Ridders expression for the
							// estimate of the root
      double ridders_enumerator;			// enumerator of Ridders expression for the
							// estimate of the root
      int iteration_index_ridder=1;			// index of the iteration in the root finding
							// algorithm
      
      /* first check if either of the two starting values are sufficiently */
      /* close to the root */
      
      std::cout<<"start with vof "<< volume_of_fluid << " \n";
      
      if( fabs(function_g_left)< vof_2_level_set_tolerance)
      {
	    /* left hand starting value is a root */
	    level_set=level_set_left;
	    return 0;
	 
      }
      else
      {
	  if(fabs(function_g_right)< vof_2_level_set_tolerance)
	  {
   /*	     right hand starting value is a root */
	      level_set=level_set_right;
	      return 0;
	  }
	  else
	  {
	      while( fabs(function_g_value)>vof_2_level_set_tolerance &&
		      iteration_index_ridder<number_iterations_ridder)
	      {
		  new_mid_point_beta=0.5*(level_set_left+level_set_right);
		  if(level_set_2_vof(new_mid_point_beta, 
			    d_level_set_d_x1, d_level_set_d_x2, d_level_set_d_x3, 
			      new_volume_of_fluid, lower_bound_derivatives))
		  {
		      std::cout<< "error in computation of volume of fluid in ridders_method \n";
		      std::cout<< "new mid point, aborting \n";
		      return 1;
		  }
		  function_g_value_beta=new_volume_of_fluid-volume_of_fluid;
		  ridders_enumerator=0.5*(level_set_right-level_set_left)*
				    sign(1.0,function_g_left)*function_g_value_beta;
		  ridders_denominator=sqrt(function_g_value_beta*function_g_value_beta-
					function_g_right*function_g_left);
	      
		  /* make a new estimate of the root */
	      
		  level_set=new_mid_point_beta+ridders_enumerator/ridders_denominator;
		  if(level_set_2_vof(level_set, 
			    d_level_set_d_x1, d_level_set_d_x2, d_level_set_d_x3, 
			      new_volume_of_fluid, lower_bound_derivatives))
		  {
		      std::cout<< "error in computation of volume of fluid in ridders_method \n";
		      std::cout<< "new mid point, aborting \n";
		      return 1;
		  }
	      
// 		  std::cerr<<" new_volume_of_fluid "<< new_volume_of_fluid << " \n";
	      /* compute the new residual */
	      
		  function_g_value=new_volume_of_fluid-volume_of_fluid;
		  std::cout<<"g_value "<< function_g_value << " \n";
	      
//.........这里部分代码省略.........

PHPAPI int
php_version_compare(const char *orig_ver1, const char *orig_ver2)
{
	char *ver1;
	char *ver2;
	char *p1, *p2, *n1, *n2;
	long l1, l2;
	int compare = 0;

	if (!*orig_ver1 || !*orig_ver2) {
		if (!*orig_ver1 && !*orig_ver2) {
			return 0;
		} else {
			return *orig_ver1 ? 1 : -1;
		}
	}
	if (orig_ver1[0] == '#') {
		ver1 = estrdup(orig_ver1);
	} else {
		ver1 = php_canonicalize_version(orig_ver1);
	}
	if (orig_ver2[0] == '#') {
		ver2 = estrdup(orig_ver2);
	} else {
		ver2 = php_canonicalize_version(orig_ver2);
	}
	p1 = n1 = ver1;
	p2 = n2 = ver2;
	while (*p1 && *p2 && n1 && n2) {
		if ((n1 = strchr(p1, '.')) != NULL) {
			*n1 = '\0';
		}
		if ((n2 = strchr(p2, '.')) != NULL) {
			*n2 = '\0';
		}
		if (isdigit(*p1) && isdigit(*p2)) {
			/* compare element numerically */
			l1 = strtol(p1, NULL, 10);
			l2 = strtol(p2, NULL, 10);
			compare = sign(l1 - l2);
		} else if (!isdigit(*p1) && !isdigit(*p2)) {
			/* compare element names */
			compare = compare_special_version_forms(p1, p2);
		} else {
			/* mix of names and numbers */
			if (isdigit(*p1)) {
				compare = compare_special_version_forms("#N#", p2);
			} else {
				compare = compare_special_version_forms(p1, "#N#");
			}
		}
		if (compare != 0) {
			break;
		}
		if (n1 != NULL) {
			p1 = n1 + 1;
		}
		if (n2 != NULL) {
			p2 = n2 + 1;
		}
	}
	if (compare == 0) {
		if (n1 != NULL) {
			if (isdigit(*p1)) {
				compare = 1;
			} else {
				compare = php_version_compare(p1, "#N#");
			}
		} else if (n2 != NULL) {
			if (isdigit(*p2)) {
				compare = -1;
			} else {
				compare = php_version_compare("#N#", p2);
			}
		}
	}
	efree(ver1);
	efree(ver2);
	return compare;
}

Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
:
    U_(A),
    V_(A.m(), A.m()),
    S_(A.m()),
    VSinvUt_(A.m(), A.n()),
    nZeros_(0)
{
    // SVDcomp to find U_, V_ and S_ - the singular values

    const label Um = U_.m();
    const label Un = U_.n();

    scalarList rv1(Um);
    scalar g = 0;
    scalar scale = 0;
    scalar s = 0;
    scalar anorm = 0;
    label l = 0;

    for (label i = 0; i < Um; i++)
    {
        l = i+2;
        rv1[i] = scale*g;
        g = s = scale = 0;

        if (i < Un)
        {
            for (label k = i; k < Un; k++)
            {
                scale += mag(U_[k][i]);
            }

            if (scale != 0)
            {
                for (label k = i; k < Un; k++)
                {
                    U_[k][i] /= scale;
                    s += U_[k][i]*U_[k][i];
                }

                scalar f = U_[i][i];
                g = -sign(Foam::sqrt(s), f);
                scalar h = f*g - s;
                U_[i][i] = f - g;

                for (label j = l-1; j < Um; j++)
                {
                    s = 0;
                    for (label k = i; k < Un; k++)
                    {
                        s += U_[k][i]*U_[k][j];
                    }

                    f = s/h;
                    for (label k = i; k < A.n(); k++)
                    {
                        U_[k][j] += f*U_[k][i];
                    }
                }

                for (label k = i; k < Un; k++)
                {
                    U_[k][i] *= scale;
                }
            }
        }

        S_[i] = scale*g;

        g = s = scale = 0;

        if (i+1 <= Un && i != Um)
        {
            for (label k = l-1; k < Um; k++)
            {
                scale += mag(U_[i][k]);
            }

            if (scale != 0)
            {
                for (label k=l-1; k < Um; k++)
                {
                    U_[i][k] /= scale;
                    s += U_[i][k]*U_[i][k];
                }

                scalar f = U_[i][l-1];
                g = -sign(Foam::sqrt(s),f);
                scalar h = f*g - s;
                U_[i][l-1] = f - g;

                for (label k = l-1; k < Um; k++)
                {
                    rv1[k] = U_[i][k]/h;
                }

                for (label j = l-1; j < Un; j++)
                {
                    s = 0;
//.........这里部分代码省略.........

void Table::addCoord(Coordinate c)
{
    switch(mType)
    {
        case PolyEdit::Polygon:
        {
            int row = rowCount()-1;
            setItem(row,0, new Cell(QString::number(c.rx())));
            setItem(row,1, new Cell(QString::number(c.ry())));
            break;
        }
        case PolyEdit::Circle:
        {
            Cell *item[3];
            item[0] = getItem(0,0);
            item[1] = getItem(1,0);
            item[2] = getItem(2,0);

            if (!item[2]->text().isEmpty() && mMask->size() == 1)
            {
                clearContents();
                mMask->clear();
                mMask->addc(c);
                mMask->setRadius(0);
                setItem(0,0, new Cell(QString::number(c.rx())));
                setItem(1,0, new Cell(QString::number(c.ry())));
            }
            else if (item[0]->text().isEmpty() || item[1]->text().isEmpty())
            {
                setItem(0,0, new Cell(QString::number(c.rx())));
                setItem(1,0, new Cell(QString::number(c.ry())));
                mMask->addc(c);
            }
            else
            {
                int x1 = item[0]->getInt();
                int y1 = item[1]->getInt();
                Coordinate c1(x1,y1);
                //int x2 = c.rx();
                //int y2 = c.ry();

                int radius = distance(c, c1);//sqrt(pow((x2-x1),2) + pow((y2-y1),2));

                setItem(2,0, new Cell(QString::number(radius)));
            }
            break;
        }
        case PolyEdit::Box:
        {
            bool item[4];
            item[0] = getItem(0,0)->text().isEmpty();
            item[1] = getItem(1,0)->text().isEmpty();
            item[2] = getItem(2,0)->text().isEmpty();
            item[3] = getItem(3,0)->text().isEmpty();

            if (!item[0] && !item[1] && !item[2] && !item[3])
            {
                clearContents();
                setItem(0,0, new Cell(QString::number(c.rx())));
                setItem(1,0, new Cell(QString::number(c.ry())));
                mMask->clear();
                mMask->addc(c);
            }
            if (item[0] || item[1])
            {
                setItem(0,0, new Cell(QString::number(c.rx())));
                setItem(1,0, new Cell(QString::number(c.ry())));
                mMask->addc(c);
            }
            else if (item[2] || item[3])
            {
                int width = c.rx() - getItem(0,0)->getInt();
                int height = getItem(1,0)->getInt() - c.ry();

                if (sign(width) == -1)
                {
                    setItem(0,0, new Cell(QString::number(c.rx())));
                    setItem(2,0, new Cell(QString::number(abs(width))));
                }
                else
                    setItem(2,0, new Cell(QString::number(width)));

                if (sign(height) == -1)
                {
                    setItem(1,0, new Cell(QString::number(c.ry())));
                    setItem(3,0, new Cell(QString::number(abs(height))));
                }
                else
                    setItem(3,0, new Cell(QString::number(height)));
            }

            break;
        }
        case PolyEdit::Invalid:
        {
            break;
        }
    }
}

//---------------------------------------------------------------------------
void __fastcall TQrySaleForm::btnQryYearClick(TObject *Sender)
{
        CString szSQL = "";
        /////////////////////////////////////////////////////////////////////////////////////////////
        //业绩查询
        szSQL = "select op_code_sl, isnull(busi_volume_cur,0) gross, \
                        dbo.gen_should_com(busi_volume_cur, dbo.gen_new_cli_prize_yearmonth(year_month,op_code_sl), busi_volume_base, ratio) should_com, \
                        cast(100*(1-isnull(busi_volume_cur,0)/isnull(busi_volume_base,100)) as decimal(9,2)) as pect_busi,  \
                        case sign(busi_volume_base-busi_volume_cur)           \
                                when 1 then (busi_volume_base-busi_volume_cur)\
                                else 0 end as last_busi,\
                        *\
                        from commition \
                        where 1=1 ";
        if (g_theOperator.op_class==E_OPERATOR_TYPE_SALER){
                szSQL += " and op_code_sl="; szSQL += Str2DBString(g_theOperator.op_code);
        }
        if (!cbbYear->Text.IsEmpty()){
                szSQL += " and year(year_month)="; szSQL += Str2DBString(cbbYear->Text.c_str());
        }
        if (!cbbMonth->Text.IsEmpty()){
                szSQL += " and month(year_month)="; szSQL += Str2DBString(cbbMonth->Text.c_str());
        }
        szSQL += "order by year_month, op_code_sl";
        Edit1->Text = AnsiString(szSQL);
        //run
        RunSQL(szSQL, true);
        if (dm1->Query1->Eof){
                ShowMessage("记录不存在");
        }

        lvTask->Clear();

        //show
        double sum_should_com = 0;
	while(!dm1->Query1->Eof)
	{
                TListItem *pItem = lvTask->Items->Add();
                double should_com = dm1->Query1->FieldByName("should_com")->AsFloat;
                sum_should_com += should_com;
                pItem->Caption = dm1->Query1->FieldByName("op_code_sl")->AsString;
                pItem->SubItems->Add(dm1->Query1->FieldByName("gross")->AsString);
                pItem->SubItems->Add(should_com);
                pItem->SubItems->Add(0.85*should_com);
                pItem->SubItems->Add(0.15*should_com);
                pItem->SubItems->Add(dm1->Query1->FieldByName("passed")->AsString);
                pItem->SubItems->Add(dm1->Query1->FieldByName("busi_volume_cur")->AsString);
                pItem->SubItems->Add(dm1->Query1->FieldByName("last_busi")->AsString);
                double last_busi = dm1->Query1->FieldByName("pect_busi")->AsFloat;
                last_busi = last_busi<0?0:last_busi;
                pItem->SubItems->Add(last_busi);

        	dm1->Query1->Next();
        }
        lbTotalPrize->Caption =  0.85*sum_should_com;
        /////////////////////////////////////////////////////////////////////////////////////////////
        //回扣查询
        szSQL = "select *, \
                CONVERT(varchar, acceptdate, 111) as date_accept, \
                case year(date) when 2099 then '-' else CONVERT(varchar, date, 111) end as date_kb \
                from detail_kb, operation \
                where dkoid=oid ";
        if (g_theOperator.op_class==E_OPERATOR_TYPE_SALER){   
                szSQL += " and op_code_sl="; szSQL += Str2DBString(g_theOperator.op_code);
        }
        if (!cbbClShortName->Text.IsEmpty()){
                szSQL += " and cl_shortname="; szSQL += Str2DBString(cbbClShortName->Text.c_str());
        }
        if (!cbbYear->Text.IsEmpty()){
                szSQL += " and year(acceptdate)="; szSQL += Str2DBString(cbbYear->Text.c_str());
        }
        if (!cbbMonth->Text.IsEmpty()){
                szSQL += " and month(acceptdate)="; szSQL += Str2DBString(cbbMonth->Text.c_str());
        }

        RunSQL(szSQL, true);
        if (dm1->Query1->Eof){
                ShowMessage("回扣记录不存在");
        }
        lvKB->Clear();
	while(!dm1->Query1->Eof)
	{
                TListItem *pItem = lvKB->Items->Add();
                pItem->Caption = dm1->Query1->FieldByName("oid")->AsString;
                pItem->SubItems->Add(dm1->Query1->FieldByName("date_accept")->AsString);
                pItem->SubItems->Add(dm1->Query1->FieldByName("cl_shortname")->AsString);
                pItem->SubItems->Add(dm1->Query1->FieldByName("kb_value")->AsFloat);
                pItem->SubItems->Add(dm1->Query1->FieldByName("date_kb")->AsString);

        	dm1->Query1->Next();
        }

}

/*
 * When doing (row,column) to (lat,lon) conversions it's sometimes useful
 * to have some auxillary values around to simplify the calculations.  This
 * function computes those auxillary values, if any, for a projection and
 * stores them in the projection's auxargs array.
 * Input:  proj - the projection
 */
static void compute_aux_proj_args( struct projection *proj )
{
   float lat1, lat2;

   switch (proj->Kind) {

      case PROJ_LAMBERT:   /* Lambert conformal */
#define Lat1 proj->Args[0]
#define Lat2 proj->Args[1]
#define PoleRow proj->Args[2]
#define PoleCol proj->Args[3]
#define CentralLon proj->Args[4]
#define ColInc proj->Args[5]
#define Hemisphere proj->AuxArgs[0]
#define ConeFactor proj->AuxArgs[1]
#define Cone proj->AuxArgs[2]
         proj->AuxArgs = (float *) MALLOC( 3 * sizeof(float) );
         if (Lat1==Lat2) {
            /* polar stereographic? */
            if (Lat1>0.0) {
               lat1 = (90.0 - Lat1) * DEG2RAD;
               /* Cone = 1.0; */
            }
            else {
               lat1 = (90.0 + Lat1) * DEG2RAD;
               /* Cone = -1.0; */
            }
            Cone = cos(lat1); /* WLH 9-27-96 */
            Hemisphere = 1.0;
         }
         else {
            /* general Lambert conformal */
            float a, b;
            if (sign(Lat1) != sign(Lat2)) {
               printf("Error: standard latitudes must have the same sign.\n");
               exit(1);
            }
            if (Lat1<Lat2) {
               printf("Error: Lat1 must be >= Lat2\n");
               exit(1);
            }
            Hemisphere = 1.0;
            lat1 = (90.0 - Lat1) * DEG2RAD;
            lat2 = (90.0 - Lat2) * DEG2RAD;
            a = log( sin(lat1) ) - log( sin(lat2) );
            b = log( tan(lat1/2.0) ) - log( tan(lat2/2.0) );
            Cone = a / b;
         }

         /* Cone is in [-1,1] */
         ConeFactor = EARTH_RADIUS * sin(lat1)
                          / (ColInc * Cone * pow(tan(lat1/2.0), Cone) );
#undef Lat1
#undef Lat2
#undef PoleRow
#undef PoleCol
#undef CentralLon
#undef ColInc
#undef Hemisphere
#undef ConeFactor
#undef Cone
         break;

      default:
         /* No auxillary args */
         proj->AuxArgs = NULL;
   }

}

static int point2DcmpX(const Point2D &p1, const Point2D &p2) {
  return sign(p1.x-p2.x);
}

//.........这里部分代码省略.........
          SP_fft3 = (double (*)[2]) fftw_malloc ( sizeof ( fftw_complex ) * (( (SP_fft_length) / 2 ) + 1) );
          
          SP_fft_plan1 = fftw_plan_dft_r2c_1d(FP_FFT_LENGTH, SP_adc1_samples, SP_fft1, FFTW_ESTIMATE);
          SP_fft_plan2 = fftw_plan_dft_r2c_1d(FP_FFT_LENGTH, SP_adc1_samples, SP_fft2, FFTW_ESTIMATE);
          SP_fft_plan3 = fftw_plan_dft_r2c_1d(FP_FFT_LENGTH, SP_adc1_samples, SP_fft3, FFTW_ESTIMATE);
          
          // Perform FFT on current input data
          fftw_execute(SP_fft_plan1);
          fftw_execute(SP_fft_plan2);
          fftw_execute(SP_fft_plan3);
          
          FFTBinAffinePrediction SP_f1_bin(bitrate,
                                           SP_fft_length,
                                           target_frequency1,
                                           BIN_PREDICTION_M,
                                           BIN_PREDICTION_B,
                                           3);
          /*FFTBinAffinePrediction SP_f2_bin(bitrate,
                                           SP_fft_length,
                                           target_frequency2,
                                           BIN_PREDICTION_M,
                                           BIN_PREDICTION_B,
                                           3);*/
          
          int SP_target_bin = SP_f1_bin.getTargetBin();
          
          fprintf(stderr, "SP Target Bin: %d\n", SP_target_bin);
          
          // Calculate Phase of each signal found on each ADC
          double phase1 = atan2((double)SP_fft1[(int)SP_target_bin][1], (double)SP_fft1[(int)SP_target_bin][0]);
          double phase2 = atan2((double)SP_fft2[(int)SP_target_bin][1], (double)SP_fft2[(int)SP_target_bin][0]);
          double phase3 = atan2((double)SP_fft3[(int)SP_target_bin][1], (double)SP_fft3[(int)SP_target_bin][0]);
          
          // Calculate usable phase differences
          double delta1 = phase2 - phase1;
          double delta2 = phase3 - phase2;
          double delta3 = phase1 - phase3;
          
          fprintf(stderr, "Phase: %5.5f  %5.5f  %5.5f\n", phase1, phase2, phase3);
          fprintf(stderr, "Deltas: %5.5f  %5.5f  %5.5f\n", delta1, delta2, delta3);
          
          // Free Data Array Memory
          fftw_free(SP_fft1);
          fftw_free(SP_fft2);
          fftw_free(SP_fft3);
          
          // Determine minimum phase difference for pair selection
          int min_index = 3;
          if (fabs(delta2) < fabs(delta1) && fabs(delta2) < fabs(delta3))
            min_index = 2;
          if (fabs(delta1) < fabs(delta2) && fabs(delta1) < fabs(delta3))
            min_index = 1;
          
          double delta_tmp;
          switch (min_index) {
            case 1:
              delta_tmp = delta1;
              if (delta3 > delta2)
                delta_tmp = -1.0 * delta1 + sign(delta_tmp) * 2.0 * M_PI;
              angle_estimate1 = delta_tmp * (double)(((double)target_wavelength1 / 2.0) / SENSOR_SPACING_2_1) * 180.0 / M_PI / 2.0;
              break;
            case 2:
              delta_tmp = delta2;
              if (delta1 > delta3)
                delta_tmp = -1.0 * delta2 + 2.0 * M_PI;
              angle_estimate1 = (delta_tmp - 4.0 / 3.0 * M_PI) * ((target_wavelength1 / 2.0) / SENSOR_SPACING_3_2) * 180.0 / M_PI / 2.0;
              break;
            case 3:
              delta_tmp = delta3;
              if (delta2 > delta1)
                delta_tmp = -1.0 * delta3 - 2.0 * M_PI;
              angle_estimate1 = (delta_tmp + 4.0 / 3.0 * M_PI ) * (((double)target_wavelength1 / 2.0) / SENSOR_SPACING_1_3) * 180.0 / M_PI / 2.0;
              break;
            default:
              fprintf(stderr, "Sonar: Invalid min_index for phase difference.");
          }
          fprintf(stderr, "Min Index: %d\n", min_index);
          fprintf(stderr, "Angle Estimate (1): %3.2f\n", angle_estimate1);
          
          // DEBUG
          for (int blah = 0; blah < SP_bin_count * FP_FFT_LENGTH; ++blah)
            printf("%5.0f", SP_adc1_samples[blah]);
          printf("\n");
          
        }
      }
  
      if (FP_bin_indicator.size() > FP_BIN_HISTORY_LENGTH)
        FP_bin_indicator.erase(FP_bin_indicator.begin());
    }
  }
  fftw_destroy_plan(FP_fft_plan1);
  fftw_destroy_plan(FP_fft_plan2);
  fftw_destroy_plan(FP_fft_plan3);
  fftw_destroy_plan(SP_fft_plan1);
  fftw_destroy_plan(SP_fft_plan2);
  fftw_destroy_plan(SP_fft_plan3);

  return 0;
}

int main(int argc, char * const argv[])
{
	int err = 0, r, c, long_optind = 0;
	int do_decipher = 0;
	int do_sign = 0;
	int action_count = 0;
        struct sc_pkcs15_object *key;
	sc_context_param_t ctx_param;

	while (1) {
		c = getopt_long(argc, argv, "sck:r:i:o:Rp:vw", options, &long_optind);
		if (c == -1)
			break;
		if (c == '?')
			util_print_usage_and_die(app_name, options, option_help, NULL);
		switch (c) {
		case 's':
			do_sign++;
			action_count++;
			break;
		case 'c':
			do_decipher++;
			action_count++;
			break;
		case 'k':
			opt_key_id = optarg;
			action_count++;
			break;
		case 'r':
			opt_reader = optarg;
			break;
		case 'i':
			opt_input = optarg;
			break;
		case 'o':
			opt_output = optarg;
			break;
		case 'R':
			opt_raw = 1;
			break;
		case OPT_SHA1:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_SHA1;
			break;
		case OPT_SHA256:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_SHA256;
			break;
		case OPT_SHA384:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_SHA384;
			break;
		case OPT_SHA512:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_SHA512;
			break;
		case OPT_SHA224:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_SHA224;
			break;
		case OPT_MD5:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_MD5;
			break;
		case OPT_HASH_NONE:
			opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_NONE;
			break;
		case OPT_PKCS1:
			opt_crypt_flags |= SC_ALGORITHM_RSA_PAD_PKCS1;
			break;
		case 'v':
			verbose++;
			break;
		case 'p':
			opt_pincode = optarg;
			break;
		case OPT_BIND_TO_AID:
			opt_bind_to_aid = optarg;
			break;
		case 'w':
			opt_wait = 1;
			break;
		}
	}
	if (action_count == 0)
		util_print_usage_and_die(app_name, options, option_help, NULL);

	if (!(opt_crypt_flags & SC_ALGORITHM_RSA_HASHES))
		opt_crypt_flags |= SC_ALGORITHM_RSA_HASH_NONE;

	memset(&ctx_param, 0, sizeof(ctx_param));
	ctx_param.ver      = 0;
	ctx_param.app_name = app_name;

	r = sc_context_create(&ctx, &ctx_param);
	if (r) {
		fprintf(stderr, "Failed to establish context: %s\n", sc_strerror(r));
		return 1;
	}

	if (verbose > 1) {
		ctx->debug = verbose;
		sc_ctx_log_to_file(ctx, "stderr");
	}

	err = util_connect_card(ctx, &card, opt_reader, opt_wait, verbose);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
				LOG_PRINT1("Extra shots: %d\n", o_ptr->pval);

				if (o_ptr->pval >= INHIBIT_SHOTS || o_ptr->pval < 0)
				{
					p += INHIBIT_POWER;	/* inhibit */
					LOG_PRINT("INHIBITING - too many extra shots\n");
				}
				else if (o_ptr->pval > 0)
				{
					p = (p * (1 + o_ptr->pval));
					LOG_PRINT2("Multiplying power by 1 + %d, total is %d\n",
						o_ptr->pval, p);
				}
			}

			/* Apply the correct slay multiplier */
			p = (p * slay_power(o_ptr, verbose, log_file, flags)) / tot_mon_power;
			LOG_PRINT1("Adjusted for slay power, total is %d\n", p);

			if (o_ptr->weight < k_ptr->weight)
			{
				p++;
				LOG_PRINT("Incrementing power by one for low weight\n");
			}

			/*
			 * Correction to match ratings to melee damage ratings.
			 * We multiply all missile weapons by 1.5 in order to compare damage.
			 * (CR 11/20/01 - changed this to 1.25).
			 * (CC 25/07/10 - changed this back to 1.5).
			 * Melee weapons assume MAX_BLOWS per turn, so we must
			 * also divide by MAX_BLOWS to get equal ratings.
			 */
			p = sign(p) * (ABS(p) * BOW_RESCALER / (10 * MAX_BLOWS));
			LOG_PRINT1("Rescaling bow power, total is %d\n", p);

			p += sign(o_ptr->to_h) * (ABS(o_ptr->to_h) * TO_HIT_POWER / 2);
			LOG_PRINT1("Adding power from to_hit, total is %d\n", p);

			break;
		}
		case TV_SHOT:
		case TV_ARROW:
		case TV_BOLT:
		case TV_DIGGING:
		case TV_HAFTED:
		case TV_POLEARM:
		case TV_SWORD:
		{
			p += (o_ptr->dd * (o_ptr->ds + 1) * DAMAGE_POWER / 4);
			LOG_PRINT1("Adding power for dam dice, total is %d\n", p);

			/* Apply the correct slay multiplier */
			p = (p * slay_power(o_ptr, verbose, log_file, flags)) / tot_mon_power;
			LOG_PRINT1("Adjusted for slay power, total is %d\n", p);

			p += (o_ptr->to_d * DAMAGE_POWER / 2);
			LOG_PRINT1("Adding power for to_dam, total is %d\n", p);

			if (of_has(flags, OF_BLOWS) &&
			    (known || object_pval_is_visible(o_ptr)))
			{
				LOG_PRINT1("Extra blows: %d\n", o_ptr->pval);
				if (o_ptr->pval >= INHIBIT_BLOWS || o_ptr->pval < 0)
				{
					p += INHIBIT_POWER;	/* inhibit */

//.........这里部分代码省略.........
            targetFunction(sigma_, m, mStar);
            secondaryQuadrature(mStar);

            return;
        }

        // Apply Ridder's algorithm to find sigma
        for (label iter = 0; iter < maxSigmaIter_; iter++)
        {
            scalar fMid, sigmaMid;

            sigmaMid = (sigmaLow + sigmaHigh)/2.0;
            fMid = targetFunction(sigmaMid, m, mStar);

            scalar s = sqrt(sqr(fMid) - fLow*fHigh);

            if (s == 0.0)
            {
                FatalErrorIn
                (
                    "Foam::extendedMomentInversion::invert\n"
                    "(\n"
                    "   const univariateMomentSet& moments\n"
                    ")"
                )<< "Singular value encountered while attempting to find root."
                 << "Moment set = " << m << endl
                 << "sigma = " << sigma_ << endl
                 << "fLow = " << fLow << endl
                 << "fMid = " << fMid << endl
                 << "fHigh = " << fHigh
                 << abort(FatalError);
            }

            sigma_ = sigmaMid + (sigmaMid - sigmaLow)*sign(fLow - fHigh)*fMid/s;

            momentsToMomentsStar(sigma_, m, mStar);

            scalar fNew = targetFunction(sigma_, m, mStar);
            scalar dSigma = (sigmaHigh - sigmaLow)/2.0;

            // Check for convergence
            if (mag(fNew) <= targetFunctionTol_ || mag(dSigma) <= sigmaTol_)
            {
                if (mag(sigma_) < sigmaTol_)
                {
                    m.invert();
                    sigma_ = 0.0;
                    nullSigma_ = true;
                    secondaryQuadrature(m);

                    return;
                }

                scalar momentError = normalizedMomentError(sigma_, m, mStar);

                if
                (
                    momentError < momentsTol_
                )
                {
                    // Found a value of sigma that preserves all the moments
                    secondaryQuadrature(mStar);

                    return;
                }
                else

//------------------------------------------------------------------------------
// getTriDiagonal
//------------------------------------------------------------------------------
bool Matrix::getTriDiagonal(Matrix* const pA) const
{
   //-------------------------------------------------------
   // initial compatibility and error checks
   //-------------------------------------------------------
   //if (!isSymmetric()) return 0;

   const int N = getRows();
   const auto pAI = new Matrix(*this);

   for (int k=0; k<N-2; k++) {
      double gama = (*pAI)(k+1,k);
      double alfa = (gama * gama);
      for (int j=k+2; j<N; j++) {
         double zeta = (*pAI)(j,k);
         alfa += (zeta * zeta);
      }

      alfa = -sign(gama) * std::sqrt(alfa);
      double beta = std::sqrt(0.5 * alfa * (alfa - gama));

      //----------------------------------------------------
      // construct column vector X
      //----------------------------------------------------
      const auto pX = new CVector(N);
      for (int p=0; p<k+1; p++) {
         (*pX)[p] = 0.0;
      }

      (*pX)[k+1] = (gama - alfa) / beta / 2.0;

      for (int q=k+2; q<N; q++) {
         (*pX)[q] = (*pAI)(q,k) / beta / 2.0;
      }

      //----------------------------------------------------
      // construct row vector Y = X'
      //----------------------------------------------------
      RVector* pY = pX->getTranspose();

      //----------------------------------------------------
      // M = 2*X*Y = 2*X*X'
      //----------------------------------------------------
      Matrix* pM = outerProduct(*pX, *pY);
      pM->multiply(2.0);

      //----------------------------------------------------
      // H = I - M = I - 2*X*X'
      //----------------------------------------------------
      const auto pH = new Matrix(N,N);
      pH->makeIdent();
      pH->subtract(*pM);

      //----------------------------------------------------
      // A = H*A*H
      //----------------------------------------------------
      Matrix* pAI0 = base::multiply(*pH,