pelhomot.cc

Gambit 是一个游戏库理论软件

    /* Adjustment of maximum number of steps for smaller step size */
    if (tweak == 1)
      limit = 10000;
    else if (tweak == 2)
      limit = 100000;
    else if (tweak == 3)
      limit = 1000000;

/* *****  END OF INITIALIZATION BLOCK.  ***** */


/* *****  MAIN LOOP.  ***** */

L120:
    i__1 = limit;
    for (iter = 1; iter <= i__1; ++iter) {
	if (y[1] < (float)0.) {
	    *arclen = s;
	    *iflag = 5;
	    return 0;
	}

/* SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH 
*/
/* CURVATURE COMPONENTS. */
/* L140: */
	curtol = hold * 10.;
	relerr = *arcre;
	abserr = *arcae;
	i__2 = np1;
	for (jw = 1; jw <= i__2; ++jw) {
	    if ((d__1 = yp[jw] - ypold[jw], abs((int)d__1)) > curtol) {
		relerr = *ansre;
		abserr = *ansae;
		goto L200;
	    }
/* L160: */
	}

/* TAKE A STEP ALONG THE CURVE. */
L200:
	stepnf_((integer *)&nc, (integer *)&nfec, (integer *)&iflagc, &start, 
		&crash, &hold, &h, &relerr, &abserr, 
		&s, &y[1], &yp[1], &yold[1], &ypold[1], &a[1], &qr[
		qr_offset], &alpha[1], &tz[1], (integer *)&pivot[1], 
		&w[1], &wp[1], &z0[
		1], &z1[1], &sspar[1], &par[1], (integer *)&ipar[1]);
/* PRINT LATEST POINT ON CURVE IF REQUESTED. */
	if (*trace > 0) { 
     print_homog(y+2,&coord_r, &coord_i);
#ifdef HOM_PRINT    
 fprintf(Hom_LogFile,"C %g %g",coord_r,coord_i)
#endif
; 
     for(jw=2; jw<=np1; jw++){
      #ifdef HOM_PRINT
   fprintf(Hom_LogFile," %d = %g",jw,y[jw])
#endif
;
     }
#ifdef HOM_PRINT  
   fprintf(Hom_LogFile," %g",y[1])
#endif
;
     } 
#ifdef HOM_PRINT
     fprintf(Hom_LogFile," 4 0 %d %d %f\n",iter,nfec,s)
#endif
;
	*nfe = nfec;
/* CHECK IF THE STEP WAS SUCCESSFUL. */
	if (iflagc > 0) {
	    *arclen = s;
	    *iflag = iflagc;
	    return 0;
	}
	if (crash) {
/* RETURN CODE FOR ERROR TOLERANCE TOO SMALL. */
	    *iflag = 2;
/* CHANGE ERROR TOLERANCES. */
	    if (*arcre < relerr) {
		*arcre = relerr;
	    }
	    if (*ansre < relerr) {
		*ansre = relerr;
	    }
	    if (*arcae < abserr) {
		*arcae = abserr;
	    }
	    if (*ansae < abserr) {
		*ansae = abserr;
	    }
/* CHANGE LIMIT ON NUMBER OF ITERATIONS. */
	    limit -= iter;
	    return 0;
	}

	if (y[1] >= (float)1.) {

/* USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE
 */
/* ANSWER AT LAMBDA = 1.0 . */

/* SAVE  YOLD  FOR ARC LENGTH CALCULATION LATER. */
	    i__2 = np1;
	    for (jw = 1; jw <= i__2; ++jw) {
		z0[jw] = yold[jw];
/* L260: */
	    }
	    rootnf_(&nc, &nfec, &iflagc, ansre, ansae, &y[1], &yp[1], &yold[1],
		    &ypold[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1],
		    &pivot[1], &w[1], &wp[1], &par[1], &ipar[1]);

	    *nfe = nfec;
	    *iflag = 1;
/* SET ERROR FLAG IF  ROOTNF  COULD NOT GET THE POINT ON THE ZERO 
*/
/* CURVE AT  LAMBDA = 1.0  . */
	    if (iflagc > 0) {
		*iflag = iflagc;
	    }
/* CALCULATE FINAL ARC LENGTH. */
	    i__2 = np1;
	    for (jw = 1; jw <= i__2; ++jw) {
		w[jw] = y[jw] - z0[jw];
/* L290: */
	    }
	    *arclen = s - hold + dnrm2_((integer *)&np1, &w[1], &c__1);
	    return 0;
	}

/* FOR POLYNOMIAL SYSTEMS AND THE  POLSYS  HOMOTOPY MAP, */
/* D LAMBDA/DS .GE. 0 NECESSARILY.  THIS CONDITION IS FORCED HERE IF 
*/
/* THE ENTRY POINT WAS  POLYNF . */

/*	if (polsys) {   */
	    if (yp[1] < (float)0.) {
/* REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 . */
		i__2 = np1;
		for (jw = 1; jw <= i__2; ++jw) {
		    yp[jw] = -yp[jw];
		    ypold[jw] = yp[jw];
/* L310: */
/*		} */
/* FORCE  STEPNF  TO USE THE LINEAR PREDICTOR FOR THE NEXT STE
P ONLY. */
		start = TRUE_;
	    }
	}

/* L400: */
    }

/* *****  END OF MAIN LOOP.  ***** */

/* LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. */
    *iflag = 3;
    *arclen = s;
    return 0;

} /* fixpnf_ */

/* Subroutine */ int fixpnf_(int    *n, 
			     double *y, 
			     int    *iflag, 
			     double *arcre, 
			     double *arcae, 
			     double *ansre, 
			     double *ansae, 
			     int    *trace, 
			     double *a, 
			     int    *nfe, 
			     double *arclen, 
			     double *yp, 
			     double *yold, 
			     double *ypold, 
			     double *qr, 
			     double *alpha, 
			     double *tz, 
			     int    *pivot, 
			     double *w, 
			     double *wp, 
			     double *z0, 
			     double *z1, 
			     double *sspar, 
			     double *par, 
			     int    *ipar,

			     int    tweak)
{
    return fixpnf_0_(0, n, y, iflag, arcre, arcae, ansre, ansae, trace, a, 
	    nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1,
	     sspar, par, ipar, tweak);
    }

/* Subroutine */ int polynf_(int    *n, 
			     double *y, 
			     int    *iflag, 
			     double *arcre, 
			     double *arcae, 
			     double *ansre, 
			     double *ansae, 
			     int    *trace, 
			     double *a, 
			     int    *nfe, 
			     double *arclen, 
			     double *yp, 
			     double *yold, 
			     double *ypold, 
			     double *qr, 
			     double *alpha, 
			     double *tz, 
			     int    *pivot, 
			     double *w, 
			     double *wp, 
			     double *z0, 
			     double *z1, 
			     double *sspar, 
			     double *par, 
			     int    *ipar,

			     int    tweak)
{
    return fixpnf_0_(1, n, y, iflag, arcre, arcae, ansre, ansae, trace, a, 
	    nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1,
	     sspar, par, ipar, tweak);
    }



/* end fixpnf.c */

/**********************************************************************/
/***************** implementation code from Hom_Mem.c *****************/
/**********************************************************************/

/*--------------------------------------------------------------------
 Workspace and Workspace management macros
-------------------------------------------------------------------*/
#define TOPN 20       /*Assume no more than 20 variables*/
#define TOPM 500      /* Assume no more than 500 monomials total*/
static int Istore[TOPM*(TOPN+1)+3*TOPN+2+TOPM];
static double Dstore[TOPM*2+12*(TOPN+1)+8+TOPN*(TOPN+2)+1];
static int didx=0,iidx=0;

/* initialization for storage  currently static should be
   made dynamic*/

/* access funtions to double storage */
double *Dres(int sz){ int v=didx; didx+=sz; return Dstore+v;}
int Dtop(){return didx;}
int Dfree(int ntop){ didx=ntop; return 0;}

/* access funtions to int storage */
int *Ires(int sz){ int v=iidx; iidx+=sz; return Istore+v;}
int Itop(){return iidx;}
int Ifree(int ntop){ iidx=ntop; return 0;}

/* end Hom_Mem.c */

/**********************************************************************/
/***************** implementation code from d1mach.c ******************/
/**********************************************************************/

/* Table of constant values */

doublereal d1mach_(integer *i)
{
    /* Initialized data */

    static integer sc = 987;
    static struct {
	integer e_1[10];
	doublereal e_2;
	} equiv_4 = { 0, 1048576, -1, 2146435071, 0, 1017118720, 0, 
		1018167296, 1352628735, 1070810131, 0.0 };



    /* System generated locals */
    doublereal ret_val = 0.0;

    /* Builtin functions */
    /* Subroutine */ int s_stop();
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
#define log10 ((integer *)&equiv_4 + 8)
#define dmach ((doublereal *)&equiv_4)
#define large ((integer *)&equiv_4 + 2)
#define small ((integer *)&equiv_4)
#define diver ((integer *)&equiv_4 + 6)
#define right ((integer *)&equiv_4 + 4)




/*  DOUBLE-PRECISION MACHINE CONSTANTS */

/*  D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. */

/*  D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. */

/*  D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. */

/*  D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. */

/*  D1MACH( 5) = LOG10(B) */

/*  TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, */
/*  THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY */
/*  REMOVING THE C FROM COLUMN 1. */
/*  ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. */
/*  (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) */

/*  FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST */
/*  TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE.  IF YOU DO NOT */
/*  KNOW WHICH SET TO USE, TRY BOTH AND SEE WHICH GIVES PLAUSIBLE */
/*  VALUES. */

/*  WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED */
/*  TO SPECIFY THE CONSTANTS EXACTLY.  SOMETIMES THIS REQUIRES USING */
/*  EQUIVALENT INTEGER ARRAYS.  IF YOUR COMPILER USES HALF-WORD */
/*  INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO */
/*  CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER */
/*  TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. */

/*  COMMENTS JUST BEFORE THE END STATEMENT (LINES STARTING WITH *) */
/*  GIVE C SOURCE FOR D1MACH. */




/*     MACHINE CONSTANTS FOR BIG-ENDIAN IEEE ARITHMETIC (BINARY FORMAT) */
/*     MACHINES IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST, */
/*     SUCH AS THE AT&T 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. */
/*     SUN 3), AND MACHINES THAT USE SPARC, HP, OR IBM RISC CHIPS. */

/*      DATA SMALL(1),SMALL(2) /    1048576,          0 / */
/*      DATA LARGE(1),LARGE(2) / 2146435071,         -1 / */
/*      DATA RIGHT(1),RIGHT(2) / 1017118720,          0 / */
/*      DATA DIVER(1),DIVER(2) / 1018167296,          0 / */
/*      DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ */

/*     MACHINE CONSTANTS FOR LITTLE-ENDIAN (BINARY) IEEE ARITHMETIC */
/*     MACHINES IN WHICH THE LEAST SIGNIFICANT BYTE IS STORED FIRST, */
/*     E.G. IBM PCS AND OTHER MACHINES THAT USE INTEL 80X87 OR DEC */
/*     ALPHA CHIPS. */


/*     MACHINE CONSTANTS FOR AMDAHL MACHINES. */

/*      DATA SMALL(1),SMALL(2) /    1048576,          0 / */
/*      DATA LARGE(1),LARGE(2) / 2147483647,         -1 / */
/*      DATA RIGHT(1),RIGHT(2) /  856686592,          0 / */
/*      DATA DIVER(1),DIVER(2) /  873463808,          0 / */
/*      DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ */

/*     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. */

/*      DATA SMALL(1) / ZC00800000 / */
/*      DATA SMALL(2) / Z000000000 / */

/*      DATA LARGE(1) / ZDFFFFFFFF / */
/*      DATA LARGE(2) / ZFFFFFFFFF / */

/*      DATA RIGHT(1) / ZCC5800000 / */
/*      DATA RIGHT(2) / Z000000000 / */

/*      DATA DIVER(1) / ZCC6800000 / */
/*      DATA DIVER(2) / Z000000000 / */

/*      DATA LOG10(1) / ZD00E730E7 / */
/*      DATA LOG10(2) / ZC77800DC0 /, SC/987/ */

/*     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. */

/*      DATA SMALL(1) / O1771000000000000 / */
/*      DATA SMALL(2) / O0000000000000000 / */

/*      DATA LARGE(1) / O0777777777777777 / */
/*      DATA LARGE(2) / O0007777777777777 / */

/*      DATA RIGHT(1) / O1461000000000000 / */
/*      DATA RIGHT(2) / O0000000000000000 / */

/*      DATA DIVER(1) / O1451000000000000 / */
/*      DATA DIVER(2) / O0000000000000000 / */

/*      DATA LOG10(1) / O1157163034761674 / */
/*      DATA LOG10(2) / O0006677466732724 /, SC/987/ */

/*     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. */

/*      DATA SMALL(1) / O1771000000000000 / */
/*      DATA SMALL(2) / O7770000000000000 / */

/*      DATA LARGE(1) / O0777777777777777 / */
/*      DATA LARGE(2) / O7777777777777777 / */

/*      DATA RIGHT(1) / O1461000000000000 / */
/*      DATA RIGHT(2) / O0000000000000000 / */

/*      DATA DIVER(1) / O1451000000000000 / */
/*      DATA DIVER(2) / O0000000000000000 / */

/*      DATA LOG10(1) / O1157163034761674 / */
/*      DATA LOG10(2) / O0006677466732724 /, SC/987/ */

/*     MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. */

/*      DATA SMALL(1) / 00564000000000000000B / */
/*      DATA SMALL(2) / 00000000000000000000B / */

/*      DATA LARGE(1) / 37757777777777777777B / */
/*      DATA LARGE(2) / 37157777777777777774B / */

/*      DATA RIGHT(1) / 15624000000000000000B / */
/*      DATA RIGHT(2) / 00000000000000000000B / */

/*      DATA DIVER(1) / 15634000000000000000B / */
/*      DATA DIVER(2) / 00000000000000000000B / */

/*      DATA LOG10(1) / 17164642023241175717B / */
/*      DATA LOG10(2) / 16367571421742254654B /, SC/987/ */

/*     MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. */

/*      DATA SMALL(1) / O"00564000000000000000" / */
/*      DATA SMALL(2) / O"00000000000000000000" / */

/*      DATA LARGE(1) / O"37757777777777777777" / */
/*      DATA LARGE(2) / O"37157777777777777774" / */

/*      DATA RIGHT(1) / O"15624000000000000000" / */
/*      DATA RIGHT(2) / O"00000000000000000000" / */

/*      DATA DIVER(1) / O"15634000000000000000" / */
/*      DATA DIVER(2) / O"00000000000000000000" / */

/*      DATA LOG10(1) / O"17164642023241175717" / */
/*      DATA LOG10(2) / O"16367571421742254654" /, SC/987/ */

/*     MACHINE CONSTANTS FOR CONVEX C-1 */

/*      DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / */
/*      DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / */
/*      DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / */
/*      DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / */
/*      DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ */

/*     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. */

/*      DATA SMALL(1) / 201354000000000000000B / */
/*      DATA SMALL(2) / 000000000000000000000B / */

/*      DATA LARGE(1) / 577767777777777777777B / */
/*      DATA LARGE(2) / 000007777777777777776B / */

/*      DATA RIGHT(1) / 376434000000000000000B / */
/*      DATA RIGHT(2) / 000000000000000000000B / */

/*      DATA DIVER(1) / 376444000000000000000B / */
/*      DATA DIVER(2) / 000000000000000000000B / */

/*      DATA LOG10(1) / 377774642023241175717B / */
/*      DATA LOG10(2) / 000007571421742254654B /, SC/987/ */

/*     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 */

/*     SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED */
/*     INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) */

/*     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - */
/*     STATIC DMACH(5) */

/*      DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ */
/*      DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ */
/*      DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ */