# include # include # include # include # include # include # include "subpak.h" /******************************************************************************/ double angle_shift ( double alpha, double beta ) /******************************************************************************/ /* Purpose: ANGLE_SHIFT shifts angle ALPHA to lie between BETA and BETA+2PI. Discussion: The input angle ALPHA is shifted by multiples of 2 * PI to lie between BETA and BETA+2*PI. The resulting angle GAMMA has all the same trigonometric function values as ALPHA. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 March 2008 Author: John Burkardt Parameters: Input, double ALPHA, the angle to be shifted. Input, double BETA, defines the lower endpoint of the angle range. Output, double ANGLE_SHIFT, the shifted angle. */ { double gamma; double pi = 3.141592653589793; if ( alpha < beta ) { gamma = beta - fmod ( beta - alpha, 2.0 * pi ) + 2.0 * pi; } else { gamma = beta + fmod ( alpha - beta, 2.0 * pi ); } return gamma; } /******************************************************************************/ double angle_shift_deg ( double alpha, double beta ) /******************************************************************************/ /* Purpose: ANGLE_SHIFT_DEG shifts angle ALPHA to lie between BETA and BETA+360. Discussion: The input angle ALPHA is shifted by multiples of 360 to lie between BETA and BETA+360. The resulting angle GAMMA has all the same trigonometric function values as ALPHA. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 June 2007 Author: John Burkardt Parameters: Input, double ALPHA, the angle to be shifted. Input, double BETA, defines the lower endpoint of the angle range. Output, double ANGLE_SHIFT, the shifted angle. */ { double gamma; if ( alpha < beta ) { gamma = beta - fmod ( beta - alpha, 360.0 ) + 360.0; } else { gamma = beta + fmod ( alpha - beta, 360.0 ); } return gamma; } /******************************************************************************/ double *angle_to_rgb ( double angle ) /******************************************************************************/ /* Purpose: ANGLE_TO_RGB returns a color on the perimeter of the color hexagon. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, double ANGLE, the angle in the color hexagon. The sextants are defined by the following points: 0 degrees, 1, 0, 0, red; 60 degrees, 1, 1, 0, yellow; 120 degrees, 0, 1, 0, green; 180 degrees, 0, 1, 1, cyan; 240 degrees, 0, 0, 1, blue; 300 degrees, 1, 0, 1, magenta. Output, double ANGLE_TO_RGB[3], the RGB specifications for the color that lies at the given angle, on the perimeter of the color hexagon. One value will be 1, and one value will be 0. */ { # define DEGREES_TO_RADIANS ( 3.141592653589793 / 180.0 ) double angle2; double *rgb; rgb = ( double * ) malloc ( 3 * sizeof ( double ) ); angle = r8_modp ( angle, 360.0 ); if ( angle <= 60.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * angle / 4.0; rgb[0] = 1.0; rgb[1] = tan ( angle2 ); rgb[2] = 0.0; } else if ( angle <= 120.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * angle / 4.0; rgb[0] = cos ( angle2 ) / sin ( angle2 ); rgb[1] = 1.0; rgb[2] = 0.0; } else if ( angle <= 180.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * ( angle - 120.0 ) / 4.0; rgb[0] = 0.0; rgb[1] = 1.0; rgb[2] = tan ( angle2 ); } else if ( angle <= 240.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * ( angle - 120.0 ) / 4.0; rgb[0] = 0.0; rgb[1] = cos ( angle2 ) / sin ( angle2 ); rgb[2] = 1.0; } else if ( angle <= 300.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * ( angle - 240.0 ) / 4.0; rgb[0] = tan ( angle2 ); rgb[1] = 0.0; rgb[2] = 1.0; } else if ( angle <= 360.0 ) { angle2 = DEGREES_TO_RADIANS * 3.0 * ( angle - 240.0 ) / 4.0; rgb[0] = 1.0; rgb[1] = 0.0; rgb[2] = cos ( angle2 ) / sin ( angle2 ); } return rgb; # undef DEGREES_TO_RADIANS } /******************************************************************************/ void axis_limits ( double xmin, double xmax, int ndivs, double *pxmin, double *pxmax, double *pxdiv, int *nticks ) /******************************************************************************/ /* Purpose: AXIS_LIMITS returns "nice" axis limits for a plot. Discussion: The routine is given information about the range of a variable, and the number of divisions desired. It returns suggestions for labeling a plotting axis for the variable, including the starting and ending points, the length of a single division, and a suggested tick marking for the axis. Licensing: This code is distributed under the GNU LGPL license. Modified: 22 March 2004 Author: John Burkardt Parameters: Input, double XMIN, XMAX, the lower and upper values that must be included on the axis. XMIN must be less than XMAX. Input, int NDIVS, the number of divisions desired along the axis. Output, double *PXMIN, *PXMAX, the recommended lower and upper axis bounds. It will be the case that PXMIN <= XMIN < XMAX <= PXMAX. Output, double *PXDIV, the recommended size of a single division. Output, int *NTICKS, a suggested number of ticks to use, if subdividing each of the NDIVS divisions of the axis. */ { # define NSTEPS 5 double best; double good; int i; int ihi; int ilo; int intlog; int iticks[NSTEPS] = { 5, 4, 4, 5, 5 }; int ival; int j; double pxmax2; double pxmin2; double pxdiv2; double reldif; double steps[NSTEPS] = { 1.0, 2.0, 4.0, 5.0, 10.0 }; double temp; if ( xmin == xmax ) { xmin = xmin - 0.5; xmax = xmax + 0.5; } else if ( xmax < xmin ) { temp = xmin; xmin = xmax; xmax = temp; } if ( ndivs <= 0 ) { ndivs = 5; } /* Set RELDIF, the size of the X interval divided by the largest X. */ if ( xmax != xmin ) { reldif = ( xmax - xmin ) / r8_max ( fabs ( xmax ), fabs ( xmin ) ); } else { reldif = 0.0; } /* If RELDIF tells us that XMIN and XMAX are extremely close, do some simple things. */ if ( reldif < 0.00001 ) { if ( xmax == 0.0 ) { *pxdiv = 1.0; } else { intlog = ( int ) ( r8_log_10 ( xmax ) ); if ( intlog < 0 ) { intlog = intlog - 1; } *pxdiv = pow ( 10.0, intlog ); if ( 1.0 < *pxdiv ) { *pxdiv = 1.0; } } *nticks = 5; *pxmin = xmax - ( double ) ( ndivs / 2 ) * (*pxdiv); *pxmax = xmax + ( double ) ( ndivs - ( ndivs / 2 ) ) * (*pxdiv); } /* But now handle the more general case, when XMIN and XMAX are relatively far apart. */ else { best = -999.0; /* On second loop, increase INTLOG by 1. */ for ( j = 1; j <= 2; j++ ) { /* Compute INTLOG, roughly the logarithm base 10 of the range divided by the number of divisions. */ intlog = ( int ) ( r8_log_10 ( ( xmax - xmin ) / ( double ) ( ndivs ) ) ) + ( j - 1 ); if ( xmax - xmin < ( double ) ( ndivs ) ) { intlog = intlog - 1; } /* Now consider taking 1, 2, 4, 5 or 10 steps of size 10**INTLOG: */ for ( i = 1; i <= NSTEPS; i++ ) { /* Compute the size of each step. */ pxdiv2 = steps[i-1] * pow ( 10.0, intlog ); /* Make sure NDIVS steps can reach from XMIN to XMAX, at least. */ if ( xmax <= xmin + ndivs * pxdiv2 ) { /* Now decide where to start the axis. Start the axis at PXMIN2, to the left of XMIN, and representing a whole number of steps of size PXDIV2. */ if ( 0.0 <= xmin ) { ival = ( int ) ( xmin / pxdiv2 ); } else { ival = ( int ) ( xmin / pxdiv2 ) - 1; } pxmin2 = ival * pxdiv2; /* PXMAX2 is, of course, NDIVS steps above PXMIN2. */ pxmax2 = pxmin2 + ndivs * pxdiv2; /* Only consider going on if PXMAX2 is at least XMAX. */ if ( xmax <= pxmax2 ) { /* Now judge this grid by the relative amount of wasted axis length. */ good = ( xmax - xmin ) / ( pxmax2 - pxmin2 ); if ( best < good ) { best = good; *pxmax = pxmax2; *pxmin = pxmin2; *pxdiv = pxdiv2; *nticks = iticks[i-1]; } } } } } } /* If necessary, adjust the locations of PXMIN and PXMAX so that the interval is more symmetric in containing XMIN through XMAX. */ for ( ; ; ) { ilo = ( int ) ( ( xmin - *pxmin ) / (*pxdiv) ); ihi = ( int ) ( ( *pxmax - xmax ) / (*pxdiv) ); if ( ihi < ilo + 2 ) { break; } *pxmin = *pxmin - (*pxdiv); *pxmax = *pxmax - (*pxdiv); } return; # undef NSTEPS } /******************************************************************************/ int bar_check ( int digit[12] ) /******************************************************************************/ /* Purpose: BAR_CHECK computes the check digit for a barcode. Formula: CHECK = SUM ( I = 1, 11, by 2's ) DIGIT(I) + 3 * SUM ( I = 2, 10, by 2's ) DIGIT(I) CHECK = MOD ( 10 - MOD ( CHECK, 10 ), 10 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 23 March 2004 Author: John Burkardt Parameters: Input, int DIGIT[12], entries 1 through 11 of DIGIT contain the digits of the bar code. Each entry must be between 0 and 9. The 12th digit should be the check digit. Output, int BAR_CHECK, the correct check digit. If the bar code is correct, then DIGIT(12) should equal BAR_CHECK. */ { int check; check = ( digit[0] + digit[2] + digit[4] + digit[6] + digit[8] + digit[10] ) + 3 * ( digit[1] + digit[3] + digit[5] + digit[7] + digit[9] ); check = ( 10 - ( check % 10 ) ) % 10; return check; } /******************************************************************************/ char *bar_code ( int digit[] ) /******************************************************************************/ /* Purpose: BAR_CODE constructs the 113 character barcode from 11 digits. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 March 2004 Author: John Burkardt Parameters: Input/output, int DIGIT(12). On input, the first 11 entries of DIGIT contain a code to be turned into a barcode. On output, the 12-th entry of DIGIT is a check digit. Output, char BAR_CODE[114], the 113 character bar code corresponding to the digit information. */ { char *bar; int check; char *codel; char *coder; int i; bar = ( char * ) malloc ( 114 * sizeof ( char ) ); /* 9 character quiet zone. */ strcpy ( bar, "000000000" ); /* 3 character guard pattern. */ strcpy ( bar+9, "101" ); /* 7 character product category. */ codel = bar_digit_code_left ( digit[0] ); strcpy ( bar+12, codel ); free ( codel ); /* 35 characters contain the 5 digit manufacturer code. */ for ( i = 0; i < 5; i++ ) { codel = bar_digit_code_left ( digit[i+1] ); strcpy ( bar+19+i*7, codel ); free ( codel ); } /* Center guard pattern. */ strcpy ( bar+54, "01010" ); /* 35 characters contain the 5 digit product code. */ for ( i = 0; i < 5; i++ ) { coder = bar_digit_code_right ( digit[i+6] ); strcpy ( bar+59+i*7, coder ); free ( coder ); } /* Compute the check digit. */ check = bar_check ( digit ); digit[11] = check; coder = bar_digit_code_right ( check ); strcpy ( bar+94, coder ); /* Guard pattern. */ strcpy ( bar+101, "101" ); /* Quiet zone. */ strcpy ( bar+104, "000000000" ); bar[113] = '\0'; return bar; } /******************************************************************************/ char *bar_digit_code_left ( int digit ) /******************************************************************************/ /* Purpose: BAR_DIGIT_CODE_LEFT returns the 7 character left bar code for a digit. Example: DIGIT = 3 CODEL = '0111101' Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2004 Author: John Burkardt Parameters: Input, int DIGIT, the digit, between 0 and 9. Output, char BAR_CODE_DIGIT_LEFT[8], the 7 character left code for the digit. */ { char *codel; codel = ( char * ) malloc ( 8 * sizeof ( char ) ); if ( digit == 0 ) { strcpy ( codel, "0001101" ); } else if ( digit == 1 ) { strcpy ( codel, "0011001" ); } else if ( digit == 2 ) { strcpy ( codel, "0010011" ); } else if ( digit == 3 ) { strcpy ( codel, "0111101" ); } else if ( digit == 4 ) { strcpy ( codel, "0100011" ); } else if ( digit == 5 ) { strcpy ( codel, "0110001" ); } else if ( digit == 6 ) { strcpy ( codel, "0101111" ); } else if ( digit == 7 ) { strcpy ( codel, "0111011" ); } else if ( digit == 8 ) { strcpy ( codel, "0110111" ); } else if ( digit == 9 ) { strcpy ( codel, "0001011" ); } else { strcpy ( codel, "???????" ); } return codel; } /******************************************************************************/ char *bar_digit_code_right ( int digit ) /******************************************************************************/ /* Purpose: BAR_DIGIT_CODE_RIGHT returns the 7 character right bar code for a digit. Example: DIGIT = 3 CODER = "1000010" Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2004 Author: John Burkardt Parameters: Input, int DIGIT, the digit, between 0 and 9. Output, char BAR_DIGIT_CODE_RIGHT[8], the 7 character right code. */ { char *coder; coder = ( char * ) malloc ( 8 * sizeof ( char ) ); if ( digit == 0 ) { strcpy ( coder, "1110010" ); } else if ( digit == 1 ) { strcpy ( coder, "1100110" ); } else if ( digit == 2 ) { strcpy ( coder, "1101100" ); } else if ( digit == 3 ) { strcpy ( coder, "1000010" ); } else if ( digit == 4 ) { strcpy ( coder, "1011100" ); } else if ( digit == 5 ) { strcpy ( coder, "1001110" ); } else if ( digit == 6 ) { strcpy ( coder, "1010000" ); } else if ( digit == 7 ) { strcpy ( coder, "1000100" ); } else if ( digit == 8 ) { strcpy ( coder, "1001000" ); } else if ( digit == 9 ) { strcpy ( coder, "1110100" ); } else { strcpy ( coder, "???????" ); } return coder; } /******************************************************************************/ double bmi_english ( double w_lb, double h_ft, double h_in ) /******************************************************************************/ /* Purpose: BMI_ENGLISH computes the body mass index given English measurements. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2004 Author: John Burkardt Parameters: Input, double W_LB, the body weight in pounds. Input, double H_FT, H_IN, the body height in feet and inches Output, double BMI_ENGLISH, the body mass index. */ { double h_m; double value; double w_kg; w_kg = pounds_to_kilograms ( w_lb ); h_m = feet_to_meters ( h_ft + ( h_in / 12.0 ) ); value = bmi_metric ( w_kg, h_m ); return value; } /******************************************************************************/ double bmi_metric ( double w_kg, double h_m ) /******************************************************************************/ /* Purpose: BMI_METRIC computes the body mass index given metric measurements. Licensing: This code is distributed under the GNU LGPL license. Modified: 17 April 2004 Author: John Burkardt Parameters: Input, double W_KG, the body weight in kilograms. Input, double H_M, the body height in meters. Output, double BMI_METRIC, the body mass index. */ { double value; value = w_kg / ( h_m * h_m ); return value; } /******************************************************************************/ char ch_cap ( char ch ) /******************************************************************************/ /* Purpose: CH_CAP capitalizes a single character. Discussion: This routine should be equivalent to the library "toupper" function. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 July 1998 Author: John Burkardt Parameters: Input, char CH, the character to capitalize. Output, char CH_CAP, the capitalized character. */ { if ( 97 <= ch && ch <= 122 ) { ch = ch - 32; } return ch; } /******************************************************************************/ int ch_eqi ( char ch1, char ch2 ) /******************************************************************************/ /* Purpose: CH_EQI is TRUE (1) if two characters are equal, disregarding case. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH1, CH2, the characters to compare. Output, int CH_EQI, is TRUE (1) if the two characters are equal, disregarding case and FALSE (0) otherwise. */ { int value; if ( 97 <= ch1 && ch1 <= 122 ) { ch1 = ch1 - 32; } if ( 97 <= ch2 && ch2 <= 122 ) { ch2 = ch2 - 32; } if ( ch1 == ch2 ) { value = 1; } else { value = 0; } return value; } /******************************************************************************/ int ch_is_digit ( char c ) /******************************************************************************/ /* Purpose: CH_IS_DIGIT returns TRUE if a character is a decimal digit. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 December 2003 Author: John Burkardt Parameters: Input, char C, the character to be analyzed. Output, bool CH_IS_DIGIT, is TRUE if C is a digit. */ { if ( '0' <= c && c <= '9' ) { return 1; } else { return 0; } } /******************************************************************************/ int ch_to_digit ( char ch ) /******************************************************************************/ /* Purpose: CH_TO_DIGIT returns the integer value of a base 10 digit. Example: CH DIGIT --- ----- '0' 0 '1' 1 ... ... '9' 9 ' ' 0 'X' -1 Licensing: This code is distributed under the GNU LGPL license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH, the decimal digit, '0' through '9' or blank are legal. Output, int CH_TO_DIGIT, the corresponding integer value. If the character was 'illegal', then DIGIT is -1. */ { int digit; if ( '0' <= ch && ch <= '9' ) { digit = ch - '0'; } else if ( ch == ' ' ) { digit = 0; } else { digit = -1; } return digit; } /******************************************************************************/ double degrees_to_radians ( double degrees ) /******************************************************************************/ /* Purpose: DEGREES_TO_RADIANS converts an angle measure from degrees to radians. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 April 2008 Author: John Burkardt Parameters: Input, double DEGREES, the angle measure in degrees. Output, double DEGREES_TO_RADIANS, the angle measure in radians. */ { double value; value = ( degrees / 180.0 ) * 3.141592653589793; return value; } /******************************************************************************/ double e_constant ( void ) /******************************************************************************/ /* Purpose: E_CONSTANT returns the value of E. Discussion: "E" was named in honor of Euler, but is known as Napier's constant. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2005 Author: John Burkardt Parameters: Output, double E_CONSTANT, the base of the natural logarithm system. */ { double value = 2.718281828459045; return value; } /******************************************************************************/ double euler_constant ( void ) /******************************************************************************/ /* Purpose: EULER_CONSTANT returns the value of the Euler-Mascheroni constant. Discussion: The Euler-Mascheroni constant is often denoted by a lower-case Gamma. Gamma is defined as Gamma = limit ( M -> Infinity ) ( Sum ( 1 <= N <= M ) 1 / N ) - Log ( M ) Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2005 Author: John Burkardt Parameters: Output, double EULER_CONSTANT, the value of the Euler-Mascheroni constant. */ { double value = 0.5772156649015328; return value; } /******************************************************************************/ void fac_div ( int prime_num, int npower1[], int npower2[], int npower3[] ) /******************************************************************************/ /* Purpose: FAC_DIV divides two quantities represented as prime factors. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER1[PRIME_NUM], the powers of primes in the representation of the first quantity. Input, int NPOWER2[PRIME_NUM], the powers of primes in the representation of the second quantity. Output, int NPOWER3[PRIME_NUM], the powers of primes in the representation of the quotient. */ { int i; for ( i = 0; i < prime_num; i++ ) { npower3[i] = npower1[i] - npower2[i]; } return; } /******************************************************************************/ void fac_gcd ( int prime_num, int npower1[], int npower2[], int npower3[] ) /******************************************************************************/ /* Purpose: FAC_GCD finds the GCD of two products of prime factors. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER1[PRIME_NUM], the powers of primes in the representation of the first quantity. All the powers must be nonnegative. Input, int NPOWER2[PRIME_NUM], the powers of primes in the representation of the second quantity. All the powers must be nonnegative. Output, int NPOWER3[PRIME_NUM], the powers of primes in the representation of the GCD. */ { int i; for ( i = 0; i < prime_num; i++ ) { if ( npower1[i] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FAC_GCD - Fatal error!\n" ); fprintf ( stderr, " One of the powers is negative.\n" ); exit ( 1 ); } if ( npower2[i] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FAC_GCD - Fatal error!\n" ); fprintf ( stderr, " One of the powers is negative.\n" ); exit ( 1 ); } npower3[i] = i4_min ( npower1[i], npower2[i] ); } return; } /******************************************************************************/ void fac_lcm ( int prime_num, int npower1[], int npower2[], int npower3[] ) /******************************************************************************/ /* Purpose: FAC_LCM finds the LCM of two products of prime factors. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER1[PRIME_NUM], the powers of primes in the representation of the first quantity. Input, int NPOWER2[PRIME_NUM], the powers of primes in the representation of the second quantity. Output, int NPOWER3[PRIME_NUM], the powers of primes in the representation of the LCM. */ { int i; for ( i = 0; i < prime_num; i++ ) { if ( npower1[i] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FAC_LCM - Fatal error!\n" ); fprintf ( stderr, " One of the powers is negative.\n" ); exit ( 1 ); } if ( npower2[i] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FAC_LCM - Fatal error!\n" ); fprintf ( stderr, " One of the powers is negative.\n" ); exit ( 1 ); } npower3[i] = i4_max ( npower1[i], npower2[i] ); } return; } /******************************************************************************/ void fac_mul ( int prime_num, int npower1[], int npower2[], int npower3[] ) /******************************************************************************/ /* Purpose: FAC_MUL multiplies two quantities represented as prime factors. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER1[PRIME_NUM], the powers of primes in the representation of the first quantity. Input, int NPOWER2[PRIME_NUM], the powers of primes in the representation of the second quantity. Output, int NPOWER3[PRIME_NUM], the powers of primes in the representation of the product. */ { int i; for ( i = 0; i < prime_num; i++ ) { npower3[i] = npower1[i] + npower2[i]; } return; } /******************************************************************************/ void fac_print ( int prime_num, int npower[] ) /******************************************************************************/ /* Purpose: FAC_PRINT prints a product of prime factors. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 July 2010 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER[PRIME_NUM], the powers of primes in the representation of the quantity. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, " Prime Power\n" ); fprintf ( stdout, "\n" ); for ( i = 0; i < prime_num; i++ ) { if ( npower[i] != 0 ) { fprintf ( stdout, " %8d %8d\n", prime ( i+1 ), npower[i] ); } } return; } /******************************************************************************/ int fac_to_i4 ( int prime_num, int npower[] ) /******************************************************************************/ /* Purpose: FAC_TO_I4 converts a product of prime factors into an I4. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER[PRIME_NUM], the powers of primes in the representation of the quantity. If any of these powers are negative, then INTVAL will be set to 0. Output, int FAC_TO_I4, the integer represented by the product of the prime factors. */ { int factor; int i; int intval; int j; intval = 1; for ( i = 0; i < prime_num; i++ ) { if ( npower[i] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FAC_TO_I4 - Fatal error!\n" ); fprintf ( stderr, " One of the powers is negative.\n" ); exit ( 1 ); } factor = prime ( i+1 ); for ( j = 1; j <= npower[i]; j++ ) { intval = intval * factor; } } return intval; } /******************************************************************************/ void fac_to_rat ( int prime_num, int npower[], int *top, int *bot ) /******************************************************************************/ /* Purpose: FAC_TO_RAT converts a prime factorization into a rational value. Example: Start with the prime factorization representation: 40/9 = 2**3 * 3**(-2) * 5 Input: NPOWER = ( 3, -2, 1 ) Output: TOP = 40 ( = 2**3 * 5**1 = PRIME(1)**3 * PRIME(3)**1 ) BOT = 9 ( = 3**2 = PRIME(2)**2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int PRIME_NUM, the index of the highest prime number used in the representations. Input, int NPOWER[PRIME_NUM]. NPOWER(I) is the power of the I-th prime in the prime factorization. NPOWER(I) may be positive or negative. Output, int TOP, BOT, the top and bottom of a rational value. */ { int i; *top = 1; *bot = 1; for ( i = 0; i < prime_num; i++ ) { if ( 0 < npower[i] ) { *top = *top * i4_power ( prime ( i + 1 ), npower[i] ); } else if ( npower[i] < 0 ) { *bot = *bot * i4_power ( prime ( i + 1 ), -npower[i] ); } } return; } /******************************************************************************/ double feet_to_meters ( double ft ) /******************************************************************************/ /* Purpose: FEET_TO_METERS converts a measurement in feet to meters. Licensing: This code is distributed under the GNU LGPL license. Modified: 17 April 2004 Author: John Burkardt Parameters: Input, double FT, the length in feet. Output, double FEET_TO_METERS, the corresponding length in meters. */ { double value; value = 0.0254 * 12.0 * ft; return value; } /******************************************************************************/ double gauss_sum ( int ndim, int n, double amplitude[], double center[], double width[], double x[] ) /******************************************************************************/ /* Purpose: GAUSS_SUM evaluates a function that is the sum of Gaussians. Discussion: Gauss_Sum(X) = Sum ( 1 <= J <= Ngauss ) Amplitude(I) * exp ( -Arg ) where Arg = sum ( 1 <= I <= NDIM ) ( ( ( X(I) - Center(I,J) ) / Width(J) )^2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int NDIM, the spatial dimension. Input, int N, the number of component Gaussian functions. Input, double AMPLITUDE[N], CENTER[NDIM*N], WIDTH[N], the amplitude, center and width for the component Gaussian functions. Input, double X[NDIM], the point at which the function is to be evaluated. Output, double GAUSS_SUM, the value of the function. */ { double arg; int i; int j; double value; value = 0.0; for ( j = 0; j < n; j++ ) { arg = 0.0; for ( i = 0; i < ndim; i++ ) { arg = arg + pow ( ( x[i] - center[i+j*ndim] ) / width[j], 2 ); } value = value + amplitude[j] * exp ( -arg ); } return value; } /******************************************************************************/ int get_seed ( ) /******************************************************************************/ /* Purpose: GET_SEED returns a random seed for the random number generator. Licensing: This code is distributed under the GNU LGPL license. Modified: 17 November 2004 Author: John Burkardt Parameters: Output, int GET_SEED, a random seed value. */ { time_t clock; int i4_huge = 2147483647; int ihour; int imin; int isec; int seed; struct tm *lt; time_t tloc; /* If the internal seed is 0, generate a value based on the time. */ clock = time ( &tloc ); lt = localtime ( &clock ); /* Hours is 1, 2, ..., 12. */ ihour = lt->tm_hour; if ( 12 < ihour ) { ihour = ihour - 12; } /* Move Hours to 0, 1, ..., 11 */ ihour = ihour - 1; imin = lt->tm_min; isec = lt->tm_sec; seed = isec + 60 * ( imin + 60 * ihour ); /* We want values in [1,43200], not [0,43199]. */ seed = seed + 1; /* Remap SEED from [1,43200] to [1,HUGE]. */ seed = ( int ) ( ( ( double ) seed ) * ( ( double ) i4_huge ) / ( 60.0 * 60.0 * 12.0 ) ); /* Never use a seed of 0. */ if ( seed == 0 ) { seed = 1; } return seed; } /******************************************************************************/ double *grid1 ( int dim_num, int nstep, double x1[], double x2[] ) /******************************************************************************/ /* Purpose: GRID1 finds grid points between X1 and X2 in N dimensions. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2008 Author: John Burkardt Parameters: Input, int DIM_NUM, the dimension of the points X1 and X2. Input, int NSTEP, the number of points to be generated. NSTEP must be at least 2. Input, double X1[DIM_NUM], X2[DIM_NUM], the first and last points, between which the equally spaced points are to be computed. Output, double X[DIM_NUM*NSTEP], the set of equally spaced points. Each column of X represents one point, with X[*,1] = X1 and X[*,NSTEP] = X2. */ { int i; int j; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID1 - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID1 - Fatal error!\n" ); fprintf ( stderr, " NSTEP < 2.\n" ); fprintf ( stderr, " NSTEP = %d\n", nstep ); exit ( 1 ); } x = ( double * ) malloc ( dim_num * nstep * sizeof ( double ) ); for ( j = 1; j <= nstep; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i+(j-1)*dim_num] = ( ( double ) ( nstep - j ) * x1[i] + ( double ) ( j - 1 ) * x2[i] ) / ( double ) ( nstep - 1 ); } } return x; } /******************************************************************************/ double *grid1n ( int j, int dim_num, int nstep, double x1[], double x2[] ) /******************************************************************************/ /* Purpose: GRID1N finds the I-th grid point between X1 and X2 in N dimensions. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 July 2010 Author: John Burkardt Parameters: Input, int J, the number of the desired point. Normally J would be between 1 and NSTEP, but that is not necessary. Note that J = 1 returns X1 and J = NSTEP returns X2. Input, int DIM_NUM, the dimension of the points X, X1 and X2. Input, int NSTEP, this is the number of equally spaced points that are between X1 and X2. NSTEP must be at least 2, because X1 and X2 are always included in the set of points. Input, double X1[DIM_NUM], X2[DIM_NUM], the first and last points, between which the equally spaced points lie. Output, double GRID1N[DIM_NUM], the J-th grid point between X1 and X2. */ { int i; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID1N - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID1N - Fatal error!\n" ); fprintf ( stderr, " NSTEP < 2.\n" ); fprintf ( stderr, " NSTEP = %d\n", nstep ); exit ( 1 ); } x = ( double * ) malloc ( dim_num * sizeof ( double ) ); for ( i = 0; i < dim_num; i++ ) { x[i] = ( ( double ) ( nstep - j ) * x1[i] + ( double ) ( j - 1 ) * x2[i] ) / ( double ) ( nstep - 1 ); } return x; } /******************************************************************************/ double *grid2 ( int j1, int j2, int dim_num, int nstep, double x1[], double x2[] ) /******************************************************************************/ /* // Purpose: // // GRID2 computes grid points between X1 and X2 in N dimensions. // // Discussion: // // GRID2 computes grid points between X1 and X2 in N dimensions. // // However, X1 need not be the first point computed, nor X2 the last. // The user must specify the steps on which X1 and X2 are passed // through. These steps may even be outside the range of 1 through NSTEP. // // We assume that a set of equally spaced points have // been drawn on the line through X1 and X2, and that // they have been numbered, with X1 labeled J1 and X2 // labeled J2. J1 or J2 may be between 1 and NSTEP, // in which case X1 or X2 will actually be returned in the // X array, but there is no requirement that J1 or J2 // satisfy this condition. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int J1, J2. J1 specifies the step on which // X1 would be computed, and similarly for J2. // J1 and J2 must be distinct. // // Input, int DIM_NUM, the dimension of the points X1 and X2. // // Input, int NSTEP, this is the number of equally // spaced points that are to be generated. // NSTEP should be at least 1. // // Input, double X1[DIM_NUM], X2[DIM_NUM], the points that define // the line along which the equally spaced points are generated, and // which may or may not be included in the set of computed points. // // Output, double GRID2[DIM_NUM*NSTEP], the set of equally spaced // points. Each column of X represents one point. // If 1 <= J1 <= NSTEP, then X(*,J1) = X1, and similarly for J2. */ { int i; int j; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID2 - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( j1 == j2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID2 - Fatal error!\n" ); fprintf ( stderr, " J1 = J2, leading to zero denominator.\n" ); fprintf ( stderr, " J1 = %d\n", j1 ); fprintf ( stderr, " J2 = %d\n", j2 ); exit ( 1 ); } x = ( double * ) malloc ( nstep * dim_num * sizeof ( double ) ); for ( j = 1; j <= nstep; j++ ) { for ( i = 0; i < dim_num; i++ ) { x[i+(j-1)*dim_num] = ( ( double ) ( j2 - j ) * x1[i] + ( double ) ( j - j1 ) * x2[i] ) / ( double ) ( j2 - j1 ); } } return x; } /******************************************************************************/ double *grid2n ( int j, int j1, int j2, int dim_num, double x1[], double x2[] ) /******************************************************************************/ /* // Purpose: // // GRID2N computes one grid point between X1 and X2 in N dimensions. // // Discussion: // // However, X1 need not be the first point computed, nor X2 the last. // The user must specify the steps on which X1 and X2 are passed through. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int J, the J coordinate of the desired point. // Note that if J = J1, X will be returned as X1, and if // J = J2, X will be returned as X2. // // Input, int J1, J2. J1 specifies the step on which // X1 would be computed, and similarly for J2. That is, // we assume that a set of equally spaced points have // been drawn on the line through X1 and X2, and that // they have been numbered, with X1 labeled J1 and X2 // labeled J2. J1 and J2 must be distinct. // // Input, int DIM_NUM, the dimension of the points X1 and X2. // // Input, double X1[DIM_NUM], X2[DIM_NUM], the points that define // the line along which the equally spaced points are // generated, and which may or may not be included in the // set of computed points. // // Output, double GRID_2N[DIM_NUM]. X(I) is the J-th point from the // set of equally spaced points. */ { int i; double *x; x = ( double * ) malloc ( dim_num * sizeof ( double ) ); if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID2N - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( j1 == j2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID2N - Fatal error!\n" ); fprintf ( stderr, " J1 = J2, leading to zero denominator.\n" ); fprintf ( stderr, " J1 = %d\n", j1 ); fprintf ( stderr, " J2 = %d\n", j2 ); exit ( 1 ); } for ( i = 0; i < dim_num; i++ ) { x[i] = ( ( double ) ( j2 - j ) * x1[i] + ( double ) ( j - j1 ) * x2[i] ) / ( double ) ( j2 - j1 ); } return x; } /******************************************************************************/ double *grid3 ( int dim_num, int nstep1, int nstep2, double x1[], double x2[], double x3[] ) /******************************************************************************/ /* // Purpose: // // GRID3 computes a grid on the parallelogram set by X1, X2 and X3 in N space. // // Discussion: // // The line between X1 and X2 will have NSTEP1 points generated along // it, and the line between X1 and X3 will have NSTEP2 points generated // along it. // // Fixing the second and third indices of X represents one point, with // the following special values: // // X(*,1,1) = X1 // X(*,NSTEP1,1) = X2 // X(*,1,NSTEP2) = X3. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int DIM_NUM, the dimension of the points X1, X2 and X3. // // Input, int NSTEP1, NSTEP2. These are the number of // equally spaced points to generate in the first and second // directions. NSTEP1 and NSTEP2 must be at least 2, because X1, X2 and // X3 are always included in the set of points. // // Input, double X1[DIM_NUM], X2[DIM_NUM], X3[DIM_NUM], the points // which define three corners of the parallelogram on // which the grid will be generated. // // Output, double GRID3[DIM_NUM*NSTEP1*NSTEP2], the set of equally // spaced points. */ { int i; int j; int k; double psi1; double psi2; double psi3; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3 - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep1 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3 - Fatal error!\n" ); fprintf ( stderr, " NSTEP1 < 2.\n" ); fprintf ( stderr, " NSTEP1 = %d\n", nstep1 ); exit ( 1 ); } if ( nstep2 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3 - Fatal error!\n" ); fprintf ( stderr, " NSTEP2 < 2.\n" ); fprintf ( stderr, " NSTEP2 = %d\n", nstep2 ); exit ( 1 ); } x = ( double * ) malloc ( nstep1 * nstep2 * dim_num * sizeof ( double ) ); for ( j = 1; j <= nstep1; j++ ) { psi2 = ( double ) ( j - 1 ) / ( double ) ( nstep1 - 1 ); for ( k = 1; k <= nstep2; k++ ) { psi3 = ( double ) ( k - 1 ) / ( double ) ( nstep2 - 1 ); psi1 = 1.0 - psi2 - psi3; for ( i = 0; i < dim_num; i++ ) { x[i+(j-1)*dim_num+(k-1)*dim_num*nstep1] = psi1 * x1[i] + psi2 * x2[i] + psi3 * x3[i]; } } } return x; } /******************************************************************************/ double *grid3n ( int j, int k, int dim_num, int nstep1, int nstep2, double x1[], double x2[], double x3[] ) /******************************************************************************/ /* // Purpose: // // GRID3N computes a parallelogram grid on 3 points in N dimensions. // // Discussion: // // The line between X1 and X2 will have NSTEP1 // points generated along it, and the line between X1 and // X3 will have NSTEP2 points generated along it. // // The following special values are: // // J K X // // 1 1 X1 // NSTEP1 1 X2 // 1 NSTEP2 X3 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int J, K, the parallelogram coordinates // of the point. J measures steps from X1 to X2, and // K measures steps from X1 to X3. Normally, J would // be between 1 and NSTEP1, K between 1 and NSTEP2, // but this is not necessary. // // Input, int DIM_NUM, the dimension of the points X1, X2 and X3. // // Input, int NSTEP1, NSTEP2. These are the number of // equally spaced points to generate in the first and second // directions. NSTEP1 and NSTEP2 must be at least 2, because X1, X2 and // X3 are always included in the set of points. // // Input, double X1[DIM_NUM], X2[DIM_NUM], X3[DIM_NUM], the points // which define three corners of the parallelogram on // which the grid will be generated. // // Output, double GRID3N[DIM_NUM], the point with coordinates (J,K) // from the the set of equally spaced points. */ { int i; double psi1; double psi2; double psi3; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3N - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep1 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3N - Fatal error!\n" ); fprintf ( stderr, " NSTEP1 < 2.\n" ); fprintf ( stderr, " NSTEP1 = %d\n", nstep1 ); exit ( 1 ); } if ( nstep2 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID3N - Fatal error!\n" ); fprintf ( stderr, " NSTEP2 < 2.\n" ); fprintf ( stderr, " NSTEP2 = %d\n", nstep2 ); exit ( 1 ); } x = ( double * ) malloc ( dim_num * sizeof ( double ) ); psi2 = ( double ) ( j - 1 ) / ( double ) ( nstep1 - 1 ); psi3 = ( double ) ( k - 1 ) / ( double ) ( nstep2 - 1 ); psi1 = 1.0 - psi2 - psi3; for ( i = 0; i < dim_num; i++ ) { x[i] = psi1 * x1[i] + psi2 * x2[i] + psi3 * x3[i]; } return x; } /******************************************************************************/ double *grid4 ( int j1, int j2, int k1, int k2, int dim_num, int nstep1, int nstep2, double x1[], double x2[], double x3[] ) /******************************************************************************/ /* // Purpose: // // GRID4 computes a grid on the parallelogram set by X1, X2 and X3 in N space. // // Discussion: // // Unlike GRID3, GRID4 does not necessarily place X1 at the // "origin" of the parallelogram, with X2 and X3 set at the // extreme J and K coordinates. Instead, the user is free // to specify the J and K coordinates of the points, although // they are required to lie on a subparallelogram of the // larger one. // // The line through X1 and X2 will have NSTEP1 // points generated along it, and the line through X1 and // X3 will have NSTEP2 points generated along it. // // If we imagine that the // main parallelogram is drawn first, with coordinate // ranges 1 <= J <= NSTEP1 and 1 <= K <= NSTEP2, then // these indices determine the (J,K) coordinates of the // three points, namely: // // X1 : (J1,K1) // X2 : (J2,K1) // X3 : (J1,K2) // // Of course, we actually start with the points X1, X2, // and X3, and they define a parallelogram and a (J,K) // coordinate system over the plane containing them. We // then are free to consider the parallelogram defined // by the three points (1,1), (NSTEP1,1) and (1,NSTEP2), // which may or may not contain any of the points X1, X2 // and X3. // // Assuming that the indices J1, J2, K1 and K2 are "within // bounds", the following special values will be computed: // // X(*,J1,K1) = X1 // X(*,J2,K1) = X2 // X(*,J1,K2) = X3. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int J1, J2, K1, K2, the indices. // // Input, int DIM_NUM, the dimension of the points X1, X2 and X3. // // Input, int NSTEP1, NSTEP2. These are the number of // equally spaced points to generate in the first and second // directions. NSTEP1 and NSTEP2 should be at least 1. // // Input, double X1[DIM_NUM], X2[DIM_NUM], X3[DIM_NUM], the points // which define three corners of the parallelogram on // which the grid will be generated. // // Output, double X[DIM_NUM*NSTEP1*NSTEP2], the set of equally // spaced points. Fixing the second and third indices // of X represents one point. */ { int i; int j; int k; double psi1; double psi2; double psi3; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4 - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep1 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4 - Fatal error!\n" ); fprintf ( stderr, " NSTEP1 < 2.\n" ); fprintf ( stderr, " NSTEP1 = %d\n", nstep1 ); exit ( 1 ); } if ( nstep2 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4 - Fatal error!\n" ); fprintf ( stderr, " NSTEP2 < 2.\n" ); fprintf ( stderr, " NSTEP2 = %d\n", nstep2 ); exit ( 1 ); } if ( j1 == j2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4 - Fatal error!\n" ); fprintf ( stderr, " J1 = J2, leading to zero denominator.\n" ); fprintf ( stderr, " J1 = %d\n", j1 ); fprintf ( stderr, " J2 = %d\n", j2 ); exit ( 1 ); } if ( k1 == k2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4 - Fatal error!\n" ); fprintf ( stderr, " K1 = K2, leading to zero denominator.\n" ); fprintf ( stderr, " K1 = %d\n", k1 ); fprintf ( stderr, " K2 = %d\n", k2 ); exit ( 1 ); } x = ( double * ) malloc ( nstep1 * nstep2 * dim_num * sizeof ( double ) ); for ( j = 1; j <= nstep1; j++ ) { psi2 = ( double ) ( j - j1 ) / ( double ) ( j2 - j1 ); for ( k = 1; k <= nstep2; k++ ) { psi3 = ( double ) ( k - k1 ) / ( double ) ( k2 - k1 ); psi1 = 1.0 - psi2 - psi3; for ( i = 0; i < dim_num; i++ ) { x[i+(j-1)*dim_num+(k-1)*dim_num*nstep1] = psi1 * x1[i] + psi2 * x2[i] + psi3 * x3[i]; } } } return x; } /******************************************************************************/ double *grid4n ( int j, int j1, int j2, int k, int k1, int k2, int dim_num, int nstep1, int nstep2, double x1[], double x2[], double x3[] ) /******************************************************************************/ /* // Purpose: // // GRID4N computes a single point on a parallelogram grid in N space. // // Discussion: // // The computation is identical to that of GRID4, except that // only one point at a time is computed. // // The line through X1 and X2 will have NSTEP1 // points generated along it, and the line through X1 and // X3 will have NSTEP2 points generated along it. // // The following special values will be computed: // // J K X // // J1 K1 X1 // J2 K2 X2 // J1 K2 X3 // // If we imagine that the main parallelogram is drawn first, with // coordinate ranges 1 <= J <= NSTEP1 and 1 <= K <= NSTEP2, then // the indices J and K determine the (J,K) coordinates of the // three points X1, X2, and X3, namely: // // X1 : (J1,K1) // X2 : (J2,K1) // X3 : (J1,K2) // // Of course, we actually start with the points X1, X2, // and X3, and they define a parallelogram and an (J,K) // coordinate system over the plane containing them. We // then are free to consider the parallelogram defined // by the three points (1,1), (NSTEP1,1) and (1,NSTEP2), // which may or may not contain any of the points X1, X2 // and X3. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 July 2010 // // Author: // // John Burkardt // // Parameters: // // Input, int J, the J coordinate of the point X. // // Input, int J1, J2. See discussion. // // Input, int K, the K coordinate of the point X. // // Input, int K1, K2. See discussion. // // Input, int DIM_NUM, the dimension of the points X, X1, X2 and X3. // // Input, int NSTEP1, NSTEP2. These are the number of // equally spaced points generated in the first and second // directions. // NSTEP1 and NSTEP2 should be at least 1. // // Input, double X1[DIM_NUM], X2[DIM_NUM], X3[DIM_NUM], the points // which define three corners of the parallelogram on // which the grid will be generated. // // Output, double GRID4N[DIM_NUM], the point whose parallelogram // coordinates are (J,K). */ { int i; double psi1; double psi2; double psi3; double *x; if ( dim_num < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4N - Fatal error!\n" ); fprintf ( stderr, " DIM_NUM < 1.\n" ); fprintf ( stderr, " DIM_NUM = %d\n", dim_num ); exit ( 1 ); } if ( nstep1 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4N - Fatal error!\n" ); fprintf ( stderr, " NSTEP1 < 2.\n" ); fprintf ( stderr, " NSTEP1 = %d\n", nstep1 ); exit ( 1 ); } if ( nstep2 < 2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4N - Fatal error!\n" ); fprintf ( stderr, " NSTEP2 < 2.\n" ); fprintf ( stderr, " NSTEP2 = %d\n", nstep2 ); exit ( 1 ); } if ( j1 == j2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4N - Fatal error!\n" ); fprintf ( stderr, " J1 = J2, leading to zero denominator.\n" ); fprintf ( stderr, " J1 = %d\n", j1 ); fprintf ( stderr, " J2 = %d\n", j2 ); exit ( 1 ); } if ( k1 == k2 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GRID4N - Fatal error!\n" ); fprintf ( stderr, " K1 = K2, leading to zero denominator.\n" ); fprintf ( stderr, " K1 = %d\n", k1 ); fprintf ( stderr, " K2 = %d\n", k2 ); exit ( 1 ); } psi2 = ( double ) ( j - j1 ) / ( double ) ( j2 - j1 ); psi3 = ( double ) ( k - k1 ) / ( double ) ( k2 - k1 ); psi1 = 1.0 - psi2 - psi3; x = ( double * ) malloc ( dim_num * sizeof ( double ) ); for ( i = 0; i < dim_num; i++ ) { x[i] = psi1 * x1[i] + psi2 * x2[i] + psi3 * x3[i]; } return x; } /******************************************************************************/ short int i2_reverse_bytes ( short int x ) /******************************************************************************/ /* Purpose: I2_REVERSE_BYTES reverses the two bytes in an I2. Licensing: This code is distributed under the GNU LGPL license. Modified: 12 May 2007 Author: John Burkardt Parameters: Input, short int X, a value whose bytes are to be reversed. Output, short int I2_REVERSE_BYTES, a value with bytes in reverse order. */ { char c; union { short int yshortint; char ychar[2]; } y; y.yshortint = x; c = y.ychar[0]; y.ychar[0] = y.ychar[1]; y.ychar[1] = c; return ( y.yshortint ); } /******************************************************************************/ int i4_log_10 ( int i ) /******************************************************************************/ /* Purpose: I4_LOG_10 returns the integer part of the logarithm base 10 of an I4. Example: I I4_LOG_10 ----- -------- 0 0 1 0 2 0 9 0 10 1 11 1 99 1 100 2 101 2 999 2 1000 3 1001 3 9999 3 10000 4 Discussion: I4_LOG_10 ( I ) + 1 is the number of decimal digits in I. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, the number whose logarithm base 10 is desired. Output, int I4_LOG_10, the integer part of the logarithm base 10 of the absolute value of X. */ { int i_abs; int ten_pow; int value; if ( i == 0 ) { value = 0; } else { value = 0; ten_pow = 10; i_abs = abs ( i ); while ( ten_pow <= i_abs ) { value = value + 1; ten_pow = ten_pow * 10; } } return value; } /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the smaller of two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 29 August 2006 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_modp ( int i, int j ) /******************************************************************************/ /* Purpose: I4_MODP returns the nonnegative remainder of I4 division. Discussion: If NREM = I4_MODP ( I, J ) NMULT = ( I - NREM ) / J then I = J * NMULT + NREM where NREM is always nonnegative. The MOD function computes a result with the same sign as the quantity being divided. Thus, suppose you had an angle A, and you wanted to ensure that it was between 0 and 360. Then mod(A,360) would do, if A was positive, but if A was negative, your result would be between -360 and 0. On the other hand, I4_MODP(A,360) is between 0 and 360, always. Example: I J MOD I4_MODP I4_MODP Factorization 107 50 7 7 107 = 2 * 50 + 7 107 -50 7 7 107 = -2 * -50 + 7 -107 50 -7 43 -107 = -3 * 50 + 43 -107 -50 -7 43 -107 = 3 * -50 + 43 Licensing: This code is distributed under the GNU LGPL license. Modified: 12 January 2007 Author: John Burkardt Parameters: Input, int I, the number to be divided. Input, int J, the number that divides I. Output, int I4_MODP, the nonnegative remainder when I is divided by J. */ { int value; if ( j == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_MODP - Fatal error!\n" ); fprintf ( stderr, " I4_MODP ( I, J ) called with J = %d\n", j ); exit ( 1 ); } value = i % j; if ( value < 0 ) { value = value + abs ( j ); } return value; } /******************************************************************************/ int i4_power ( int i, int j ) /******************************************************************************/ /* Purpose: I4_POWER returns the value of I^J. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, J, the base and the power. J should be nonnegative. Output, int I4_POWER, the value of I^J. */ { int k; int value; if ( j < 0 ) { if ( i == 1 ) { value = 1; } else if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J negative.\n" ); exit ( 1 ); } else { value = 0; } } else if ( j == 0 ) { if ( i == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_POWER - Fatal error!\n" ); fprintf ( stderr, " I^J requested, with I = 0 and J = 0.\n" ); exit ( 1 ); } else { value = 1; } } else if ( j == 1 ) { value = i; } else { value = 1; for ( k = 1; k <= j; k++ ) { value = value * i; } } return value; } /******************************************************************************/ int i4_sign ( int i ) /******************************************************************************/ /* Purpose: I4_SIGN returns the sign of an I4. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input, int I, the integer whose sign is desired. Output, int I4_SIGN, the sign of I. */ { int value; if ( i < 0 ) { value = -1; } else { value = 1; } return value; } /******************************************************************************/ void i4_swap ( int *i, int *j ) /******************************************************************************/ /* Purpose: I4_SWAP switches two I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 October 2007 Author: John Burkardt Parameters: Input/output, int *I, *J. On output, the values of I and J have been interchanged. */ { int k; k = *i; *i = *j; *j = k; return; } /******************************************************************************/ int *i4_to_digits_decimal ( int i, int n ) /******************************************************************************/ /* Purpose: I4_TO_DIGITS_DECIMAL determines the last N decimal digits of an I4. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int I, the integer to be analyzed. Input, int N, the number of digits to determine. Output, int I4_TO_DIGITS_DECIMAL[N], the last N decimal digits of I. DIGIT[I-1] is the "coefficient" of 10**(I-1). */ { int *digit; int j; digit = ( int * ) malloc ( n * sizeof ( n ) ); i = abs ( i ); for ( j = 1; j <= n; j++ ) { digit[j-1] = i % 10; i = ( i - digit[j-1] ) / 10; } return digit; } /******************************************************************************/ int *i4_to_fac ( int i, int prime_num ) /******************************************************************************/ /* Purpose: I4_TO_FAC converts an I4 into a product of prime factors. Discussion: This routine will fail if the input integer is not positive, or if PRIME_NUM is too small to account for the factors of the integer. The formula is: I = Product ( 1 <= J <= PRIME_NUM ) PRIME(J)**NPOWER(J). Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int I, the integer to be factored. Input, int PRIME_NUM, the number of prime factors for which storage has been allocated. Output, int I4_TO_FAC[PRIME_NUM], the powers of the primes. */ { int j; int *npower; int p; if ( i <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_TO_FAC - Fatal error!\n" ); fprintf ( stderr, " Input integer I is not positive.\n" ); exit ( 1 ); } npower = ( int * ) malloc ( prime_num * sizeof ( int ) ); /* Try dividing the remainder by each prime. */ for ( j = 1; j <= prime_num; j++ ) { npower[j-1] = 0; p = prime ( j ); while ( ( i % p ) == 0 ) { npower[j-1] = npower[j-1] + 1; i = i / p; } } return npower; } /******************************************************************************/ int i4_uniform_ab ( int a, int b, int *seed ) /******************************************************************************/ /* Purpose: I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 May 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int A, B, the limits of the interval. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, int I4_UNIFORM_AB, a number between A and B. */ { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( float ) ( *seed ) * 4.656612875E-10; /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ double i4int_to_r8int ( int imin, int imax, int i, double rmin, double rmax ) /******************************************************************************/ /* Purpose: I4INT_TO_R8INT maps an I4 interval to an R8 interval. Discussion: The formula is R := RMIN + ( RMAX - RMIN ) * ( I - IMIN ) / ( IMAX - IMIN ) Licensing: This code is distributed under the GNU LGPL license. Modified: 05 June 2005 Author: John Burkardt Parameters: Input, int IMIN, IMAX, the range. Input, int I, the integer to be converted. Input, double RMIN, RMAX, the range. Output, double R, the corresponding value in [RMIN,RMAX]. */ { double r; if ( imax == imin ) { r = 0.5 * ( rmin + rmax ); } else { r = ( ( double ) ( imax - i ) * rmin + ( double ) ( i - imin ) * rmax ) / ( double ) ( imax - imin ); } return r; } /******************************************************************************/ void i4vec_copy ( int n, int a1[], int a2[] ) /******************************************************************************/ /* Purpose: I4VEC_COPY copies an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 25 April 2007 Author: John Burkardt Parameters: Input, int N, the number of entries in the vectors. Input, int A1[N], the vector to be copied. Input, int A2[N], the copy of A1. */ { int i; for ( i = 0; i < n; i++ ) { a2[i] = a1[i]; } return; } /******************************************************************************/ int *i4vec_indicator_new ( int n ) /******************************************************************************/ /* Purpose: I4VEC_INDICATOR_NEW sets an I4VEC to the indicator vector. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of elements of A. Output, int I4VEC_INDICATOR_NEW[N], the array. */ { int *a; int i; a = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a[i] = i + 1; } return a; } /******************************************************************************/ int i4vec_min ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_MIN returns the minimum element in an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 17 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, int A[N], the array to be checked. Output, int I4VEC_MIN, the value of the minimum element. This is set to 0 if N <= 0. */ { int i; int value; if ( n <= 0 ) { return 0; } value = a[0]; for ( i = 1; i < n; i++ ) { if ( a[i] < value ) { value = a[i]; } } return value; } /******************************************************************************/ void i4vec_permute ( int n, int p[], int base, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_PERMUTE permutes an I4VEC in place. Discussion: An I4VEC is a vector of I4's. This routine permutes an array of integer "objects", but the same logic can be used to permute an array of objects of any arithmetic type, or an array of objects of any complexity. The only temporary storage required is enough to store a single object. The number of data movements made is N + the number of cycles of order 2 or more, which is never more than N + N/2. Example: Input: N = 5 P = ( 1, 3, 4, 0, 2 ) A = ( 1, 2, 3, 4, 5 ) Output: A = ( 2, 4, 5, 1, 3 ). Licensing: This code is distributed under the GNU LGPL license. Modified: 05 June 2010 Author: John Burkardt Parameters: Input, int N, the number of objects. Input, int P[N], the permutation. P(I) = J means that the I-th element of the output array should be the J-th element of the input array. Input, int BASE, is 0 for a 0-based permutation and 1 for a 1-based permutation. Input/output, int A[N], the array to be permuted. */ { int a_temp; int i; int iget; int iput; int istart; if ( !perm_check ( n, p, base ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_PERMUTE - Fatal error!\n" ); fprintf ( stderr, " PERM_CHECK rejects this permutation.\n" ); exit ( 1 ); } /* In order for the sign negation trick to work, we need to assume that the entries of P are strictly positive. Presumably, the lowest number is BASE. So temporarily add 1-BASE to each entry to force positivity. */ for ( i = 0; i < n; i++ ) { p[i] = p[i] + 1 - base; } /* Search for the next element of the permutation that has not been used. */ for ( istart = 1; istart <= n; istart++ ) { if ( p[istart-1] < 0 ) { continue; } else if ( p[istart-1] == istart ) { p[istart-1] = - p[istart-1]; continue; } else { a_temp = a[istart-1]; iget = istart; /* Copy the new value into the vacated entry. */ for ( ; ; ) { iput = iget; iget = p[iget-1]; p[iput-1] = - p[iput-1]; if ( iget < 1 || n < iget ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4VEC_PERMUTE - Fatal error!\n" ); fprintf ( stderr, " Entry IPUT = %d of the permutation has\n", iput ); fprintf ( stderr, " an illegal value IGET = %d.\n", iget ); exit ( 1 ); } if ( iget == istart ) { a[iput-1] = a_temp; break; } a[iput-1] = a[iget-1]; } } } /* Restore the signs of the entries. */ for ( i = 0; i < n; i++ ) { p[i] = - p[i]; } /* Restore the base of the entries. */ for ( i = 0; i < n; i++ ) { p[i] = p[i] - 1 + base; } return; } /******************************************************************************/ void i4vec_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_PRINT prints an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 14 November 2003 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %6d: %8d\n", i, a[i] ); } return; } /******************************************************************************/ int i4vec_sum ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SUM sums the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_SUM = 10 Licensing: This code is distributed under the GNU LGPL license. Modified: 29 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector to be summed. Output, int I4VEC_SUM, the sum of the entries of A. */ { int i; int sum; sum = 0; for ( i = 0; i < n; i++ ) { sum = sum + a[i]; } return sum; } /******************************************************************************/ void i4vec_zero ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_ZERO zeroes an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, int A[N], a vector of zeroes. */ { int i; for ( i = 0; i < n; i++ ) { a[i] = 0; } return; } /******************************************************************************/ int *i4vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: I4VEC_ZERO_NEW creates and zeroes an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 September 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, int I4VEC_ZERO_NEW[N], a vector of zeroes. */ { int *a; int i; a = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a[i] = 0; } return a; } /******************************************************************************/ void ij_next ( int *i, int *j, int n ) /******************************************************************************/ /* Purpose: IJ_NEXT returns the next matrix index. Discussion: For N = 3, the sequence of indices returned is: (1,1), (1,2), (1,3), (2,1), (2,2), (2,3), (3,1), (3,2), (3,3), (0,0). Note that once the value (N,N) is returned, the next value returned will be (0,0). Licensing: This code is distributed under the GNU LGPL license. Modified: 24 July 2010 Author: John Burkardt Parameters: Input/output, int *I, *J. On input, the current pair of indices. On output, the next pair of indices. If either index is illegal on input, the output value of (I,J) will be (1,1). Input, int N, the maximum value for I and J. */ { if ( n < 1 ) { *i = 0; *j = 0; return; } if ( *i < 1 || n < *i || *j < 1 || n < *j ) { *i = 1; *j = 1; return; } if ( *j < n ) { *j = *j + 1; } else if ( *i < n ) { *i = *i + 1; *j = 1; } else { *i = 0; *j = 0; } return; } /******************************************************************************/ void ij_next_gt ( int *i, int *j, int n ) /******************************************************************************/ /* Purpose: IJ_NEXT_GT returns the next matrix index, with the constraint that I < J. Discussion: For N = 5, the sequence of indices returned is: (1,2), (1,3), (1,4), (1,5), (2,3), (2,4), (2,5), (3,4), (3,5), (4,5). Licensing: This code is distributed under the GNU LGPL license. Modified: 24 July 2010 Author: John Burkardt Parameters: Input/output, int *I, *J. On input, the current pair of indices. On output, the next pair of indices. If either index is illegal on input, the output value of (I,J) will be (1,2). Input, int N, the maximum value for I and J. A value of N less than 2 is nonsense. */ { if ( n < 2 ) { *i = 0; *j = 0; return; } if ( *i < 1 || n < *i || *j < 1 || n < *j || *j <= *i ) { *i = 1; *j = 2; return; } if ( *j < n ) { *j = *j + 1; } else if ( *i < n - 1 ) { *i = *i + 1; *j = *i + 1; } else { *i = 0; *j = 0; } return; } /******************************************************************************/ int index1_col ( int i_min, int i, int i_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX1_COL indexes a 1D vector by columns. Discussion: This 1D routine is provided primarily for analogy. Moreover, in 1D there is no difference between row and column indexing. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for the first index, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of element I_MIN. Typically, this is 0 or 1. Output, int INDEX1_COL, the index of element I. */ { int value; value = index_min + ( i - i_min ); return value; } /******************************************************************************/ int index1_row ( int i_min, int i, int i_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX1_ROW indexes a 1D vector by rows. Discussion: This 1D routine is provided primarily for analogy. Moreover, in 1D there is no difference between row and column indexing. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for the first index, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of element I_MIN. Typically, this is 0 or 1. Output, int INDEX1_ROW, the index of element I. */ { int value; value = index_min + ( i - i_min ); return value; } /******************************************************************************/ int index2_col ( int i_min, int i, int i_max, int j_min, int j, int j_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX2_COL indexes a 2D array by columns. Discussion: Entries of the array are indexed starting at entry (I_MIN,J_MIN), and increasing the row index first. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for row indices, the minimum, the index, and the maximum. Input, int J_MIN, J, J_MAX, for column indices, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of element (I_MIN,J_MIN). Typically, this is 0 or 1. Output, int INDEX2_COL, the index of element (I,J). */ { int value; value = index_min + ( i - i_min ) + ( j - j_min ) * ( i_max + 1 - i_min ); return value; } /******************************************************************************/ int index2_row ( int i_min, int i, int i_max, int j_min, int j, int j_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX2_ROW indexes a 2D array by row. Discussion: Entries of the array are indexed starting at entry (I_MIN,J_MIN), and increasing the column index first. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for row indices, the minimum, the index, and the maximum. Input, int J_MIN, J, J_MAX, for column indices, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of element (I_MIN,J_MIN). Typically, this is 0 or 1. Output, int INDEX2_ROW, the index of element (I,J). */ { int value; value = index_min + ( j - j_min ) + ( i - i_min ) * ( j_max + 1 - j_min ); return value; } /******************************************************************************/ int index3_col ( int i_min, int i, int i_max, int j_min, int j, int j_max, int k_min, int k, int k_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX3_COL indexes a 3D array by columns. Discussion: Entries of the array are indexed starting at entry (I_MIN,J_MIN,K_MIN), and increasing the row index first, then the column index. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for row indices, the minimum, the index, and the maximum. Input, int J_MIN, J, J_MAX, for column indices, the minimum, the index, and the maximum. Input, int K_MIN, K, K_MAX, for plane indices, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of (I_MIN,J_MIN,K_MIN). Typically, this is 0 or 1. Output, int INDEX3_COL, the index of element (I,J,K). */ { int value; value = index_min + ( i - i_min ) + ( j - j_min ) * ( i_max + 1 - i_min ) + ( k - k_min ) * ( j_max + 1 - j_min ) * ( i_max + 1 - i_min ); return value; } /******************************************************************************/ int index3_row ( int i_min, int i, int i_max, int j_min, int j, int j_max, int k_min, int k, int k_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX3_ROW indexes a 3D array by rows. Discussion: When we say "by rows", we really just mean that entries of the array are indexed starting at entry (I_MIN,J_MIN,K_MIN), and the increasing the LAST index first, then the next-to-the-last, and so on. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I_MIN, I, I_MAX, for row indices, the minimum, the index, and the maximum. Input, int J_MIN, J, J_MAX, for column indices, the minimum, the index, and the maximum. Input, int K_MIN, K, K_MAX, for plane indices, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of (I_MIN,J_MIN,K_MIN). Typically, this is 0 or 1. Output, int INDEX3_ROW, the index of element (I,J,K). */ { int value; value = index_min + ( k - k_min ) + ( j - j_min ) * ( k_max + 1 - k_min ) + ( i - i_min ) * ( j_max + 1 - j_min ) * ( k_max + 1 - k_min ); return value; } /******************************************************************************/ int index4_col ( int i1_min, int i1, int i1_max, int i2_min, int i2, int i2_max, int i3_min, int i3, int i3_max, int i4_min, int i4, int i4_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX4_COL indexes a 4D array by columns. Discussion: Entries of the array are indexed starting at (I1_MIN,I2_MIN,I3_MIN,I4_MIN), and increasing the initial index first, then the second, third and so on. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I1_MIN, I1, I1_MAX, for index 1, the minimum, the index, and the maximum. Input, int I2_MIN, I2, I2_MAX, for index 2, the minimum, the index, and the maximum. Input, int I3_MIN, I3, I3_MAX, for index 3, the minimum, the index, and the maximum. Input, int I4_MIN, I4, I4_MAX, for index 4, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of (I1_MIN,I2_MIN,I3_MIN,I4_MIN). Typically, this is 0 or 1. Output, int INDEX4_COL, the index of element (I1,I2,I3,I4). */ { int value; value = index_min + ( i1 - i1_min ) + ( i2 - i2_min ) * ( i1_max + 1 - i1_min ) + ( i3 - i3_min ) * ( i2_max + 1 - i2_min ) * ( i1_max + 1 - i1_min ) + ( i4 - i4_min ) * ( i3_max + 1 - i3_min ) * ( i2_max + 1 - i2_min ) * ( i1_max + 1 - i1_min ); return value; } /******************************************************************************/ int index4_row ( int i1_min, int i1, int i1_max, int i2_min, int i2, int i2_max, int i3_min, int i3, int i3_max, int i4_min, int i4, int i4_max, int index_min ) /******************************************************************************/ /* Purpose: INDEX4_ROW indexes a 4D array by rows. Discussion: Entries of the array are indexed starting at (I1_MIN,I2_MIN,I3_MIN,I4_MIN), and increasing the last index, then the next to last, and so on. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int I1_MIN, I1, I1_MAX, for index 1, the minimum, the index, and the maximum. Input, int I2_MIN, I2, I2_MAX, for index 2, the minimum, the index, and the maximum. Input, int I3_MIN, I3, I3_MAX, for index 3, the minimum, the index, and the maximum. Input, int I4_MIN, I4, I4_MAX, for index 4, the minimum, the index, and the maximum. Input, int INDEX_MIN, the index of element (I1_MIN,I2_MIN,I3_MIN,I4_MIN). Typically, this is 0 or 1. Output, int INDEX4_ROW, the index of element (I1,I2,I3,I4). */ { int value; value = index_min + ( i4 - i4_min ) + ( i3 - i3_min ) * ( i4_max + 1 - i4_min ) + ( i2 - i2_min ) * ( i3_max + 1 - i3_min ) * ( i4_max + 1 - i4_min ) + ( i1 - i1_min ) * ( i2_max + 1 - i2_min ) * ( i3_max + 1 - i3_min ) * ( i4_max + 1 - i4_min ); return value; } /******************************************************************************/ int indexn_col ( int n, int i_min[], int i[], int i_max[], int index_min ) /******************************************************************************/ /* Purpose: INDEXN_COL indexes an ND array by columns. Discussion: Entries of the array are indexed starting at entry ( I_MIN(1), I_MIN(2),...,I_MIN(N) ), and increasing the first index up to I_MAX(1), then the second and so on. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int N, the number of indices. Input, int I_MIN[N], the minimum indices. Input, int I[N], the indices. Input, int I_MAX[N], for maximum indices. Input, int INDEX_MIN, the index of ( I_MIN[0], I_MIN[1],...,I_MIN[N-1] ). Typically, this is 0 or 1. Output, int INDEXN_COL, the index of element I. */ { int j; int value; value = ( i[n-1] - i_min[n-1] ); for ( j = n - 2; 0 <= j; j-- ) { value = value * ( i_max[j] + 1 - i_min[j] ) + ( i[j] - i_min[j] ); } value = value + index_min; return value; } /******************************************************************************/ int indexn_row ( int n, int i_min[], int i[], int i_max[], int index_min ) /******************************************************************************/ /* Purpose: INDEXN_ROW indexes an ND array by rows. Discussion: Entries of the array are indexed starting at entry ( I_MIN(1), I_MIN(2),...,I_MIN(N) ), and increasing the last index up to I_MAX(N), then the next-to-last and so on. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 April 2010 Author: John Burkardt Parameters: Input, int N, the number of indices. Input, int I_MIN[N], the minimum indices. Input, int I[N], the indices. Input, int I_MAX[N], for maximum indices. Input, int INDEX_MIN, the index of ( I_MIN[0], I_MIN[1],...,I_MIN[N-1] ). Typically, this is 0 or 1. Output, int INDEXN_ROW, the index of element I. */ { int j; int value; value = ( i[0] - i_min[0] ); for ( j = 1; j < n; j++ ) { value = value * ( i_max[j] + 1 - i_min[j] ) + ( i[j] - i_min[j] ); } value = value + index_min; return value; } /******************************************************************************/ int lcm_12n ( int n ) /******************************************************************************/ /* Purpose: LCM_12N computes the least common multiple of the integers 1 through N. Example: N LCM_12N 1 1 2 2 3 3 4 12 5 60 6 60 7 420 8 840 9 2520 10 2520 Licensing: This code is distributed under the GNU LGPL license. Modified: 20 July 2010 Author: John Burkardt Parameters: Input, int N, the value of N. Output, int LCM_12N, the least common multiple of the integers 1 to N. */ { int i; int imult; int j; int value; value = 1; for ( i = 2; i <= n; i++ ) { imult = i; for ( j = 1; j < i; j++ ) { if ( ( imult % ( i - j ) ) == 0 ) { imult = imult / ( i - j ); } } value = value * imult; } return value; } /******************************************************************************/ int pause_input ( ) /******************************************************************************/ /* Purpose: PAUSE_INPUT waits until an input character is entered. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Output, int PAUSE_INPUT, the character read from STDIN. */ { int value; printf ( "Press RETURN to continue.\n" ); value = getc ( stdin ); return value; } /******************************************************************************/ void pause_seconds ( int seconds ) /******************************************************************************/ /* Purpose: PAUSE_SECONDS waits a specified number of seconds. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Input, int SECONDS, the number of seconds to pause. */ { time_t t1; time_t t2; t1 = time ( NULL ); t2 = t1; while ( t2 - t1 < seconds ) { t2 = time ( NULL ); } return; } /******************************************************************************/ int perm_check ( int n, int p[], int base ) /******************************************************************************/ /* Purpose: PERM_CHECK checks that a vector represents a permutation. Discussion: The routine verifies that each of the integers from BASE to to BASE+N-1 occurs among the N entries of the permutation. Set the input quantity BASE to 0, if P is a 0-based permutation, or to 1 if P is a 1-based permutation. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 June 2009 Author: John Burkardt Parameters: Input, int N, the number of entries. Input, int P[N], the array to check. Input, int BASE, the index base. Output, int PERM_CHECK, is TRUE if the permutation is OK. */ { int found; int i; int seek; for ( seek = base; seek < base + n; seek++ ) { found = 0; for ( i = 0; i < n; i++ ) { if ( p[i] == seek ) { found = 1; break; } } if ( found == 0 ) { return 0; } } return 1; } /******************************************************************************/ void perm_cycle ( int n, int p[], int *isgn, int *ncycle, int iopt ) /******************************************************************************/ /* Purpose: PERM_CYCLE analyzes a permutation. Discussion: The routine will count cycles, find the sign of a permutation, and tag a permutation. Example: Input: N = 9 IOPT = 1 P = 2, 3, 9, 6, 7, 8, 5, 4, 1 Output: NCYCLE = 3 ISGN = +1 P = -2, 3, 9, -6, -7, 8, 5, 4, 1 Licensing: This code is distributed under the GNU LGPL license. Modified: 07 July 2010 Author: FORTRAN77 original version by Albert Nijenhuis, Herbert Wilf. C version by John Burkardt. Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the number of objects being permuted. Input/output, int P[N]. On input, P describes a permutation, in the sense that entry I is to be moved to P[I]. If IOPT = 0, then P will not be changed by this routine. If IOPT = 1, then on output, P will be "tagged". That is, one element of every cycle in P will be negated. In this way, a user can traverse a cycle by starting at any entry I1 of P which is negative, moving to I2 = ABS(P[I1]), then to P[I2], and so on, until returning to I1. Output, int *ISGN, the "sign" of the permutation, which is +1 if the permutation is even, -1 if odd. Every permutation may be produced by a certain number of pairwise switches. If the number of switches is even, the permutation itself is called even. Output, int *NCYCLE, the number of cycles in the permutation. Input, int IOPT, requests tagging. 0, the permutation will not be tagged. 1, the permutation will be tagged. */ { int base = 1; int i; int i1; int i2; int is; if ( !perm_check ( n, p, base ) ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_CYCLE - Fatal error!\n" ); fprintf ( stderr, " PERM_CHECK rejects this permutation.\n" ); exit ( 1 ); } is = 1; *ncycle = n; for ( i = 1; i <= n; i++ ) { i1 = p[i-1]; while ( i < i1 ) { *ncycle = *ncycle - 1; i2 = p[i1-1]; p[i1-1] = -i2; i1 = i2; } if ( iopt != 0 ) { is = - i4_sign ( p[i-1] ); } p[i-1] = abs ( p[i-1] ) * i4_sign ( is ); } *isgn = 1 - 2 * ( ( n - *ncycle ) % 2 ); return; } /******************************************************************************/ int *perm_free ( int npart, int ipart[], int nfree ) /******************************************************************************/ /* Purpose: PERM_FREE reports the number of unused items in a partial permutation. Discussion: It is assumed that the N objects being permuted are the integers from 1 to N, and that IPART contains a "partial" permutation, that is, the NPART entries of IPART represent the beginning of a permutation of all N items. The routine returns in IFREE the items that have not been used yet. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2005 Author: John Burkardt Parameters: Input, int NPART, the number of entries in IPART. NPART may be 0. Input, int IPART[NPART], the partial permutation, which should contain, at most once, some of the integers between 1 and NPART+NFREE. Input, int NFREE, the number of integers that have not been used in IPART. This is simply N - NPART. NFREE may be zero. Output, int PERM_FREE[NFREE], the integers between 1 and NPART+NFREE that were not used in IPART. */ { int i; int *ifree; int j; int k; int match; int n; n = npart + nfree; if ( npart < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_FREE - Fatal error!\n" ); fprintf ( stderr, " NPART < 0.\n" ); fprintf ( stderr, " NPART = %d\n", npart ); exit ( 1 ); } else if ( npart == 0 ) { ifree = i4vec_indicator_new ( n ); return ifree; } else if ( nfree < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_FREE - Fatal error!\n" ); fprintf ( stderr, " NFREE < 0.\n" ); fprintf ( stderr, " NFREE = %d\n", nfree ); exit ( 1 ); } else if ( nfree == 0 ) { return NULL; } else { ifree = ( int * ) malloc ( nfree * sizeof ( int ) ); k = 0; for ( i = 1; i <= n; i++ ) { match = 0; for ( j = 1; j <= npart; j++ ) { if ( ipart[j-1] == i ) { match = j; break; } } if ( match == 0 ) { k = k + 1; if ( nfree < k ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_FREE - Fatal error!\n" ); fprintf ( stderr, " The partial permutation is illegal.\n" ); fprintf ( stderr, " It should contain, at most once, some of\n" ); fprintf ( stderr, " the integers between 1 and N = %d\n", n ); fprintf ( stderr, " The number of integers that have not\n" ); fprintf ( stderr, " been used is at least K = %d\n", k ); fprintf ( stderr, " This should be exactly NFREE = %d\n", nfree ); i4vec_print ( npart, ipart, " The partial permutation:" ); exit ( 1 ); } ifree[k-1] = i; } } } return ifree; } /******************************************************************************/ void perm_inverse ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_INVERSE inverts a permutation "in place". Discussion: This algorithm assumes that the entries in the permutation vector are strictly positive. In particular, the value 0 must not occur. When necessary, this function shifts the data temporarily so that this requirement is satisfied. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 June 2009 Author: John Burkardt Parameters: Input, int N, the number of objects being permuted. Input/output, int P[N], the permutation, in standard index form. On output, P describes the inverse permutation */ { int base; int i; int i0; int i1; int i2; int is; int p_min; if ( n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_INVERSE - Fatal error!\n" ); fprintf ( stderr, " Input value of N = %d\n", n ); exit ( 1 ); } /* Find the least value, and shift data so it begins at 1. */ p_min = i4vec_min ( n, p ); base = 1; for ( i = 0; i < n; i++ ) { p[i] = p[i] - p_min + base; } /* Now we can safely check the permutation. */ if ( perm_check ( n, p, base ) == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_INVERSE - Fatal error!\n" ); fprintf ( stderr, " PERM_CHECK rejects this permutation.\n" ); exit ( 1 ); } /* Now we can invert the permutation. */ is = 1; for ( i = 1; i <= n; i++ ) { i1 = p[i-1]; while ( i < i1 ) { i2 = p[i1-1]; p[i1-1] = -i2; i1 = i2; } is = - i4_sign ( p[i-1] ); p[i-1] = i4_sign ( is ) * abs ( p[i-1] ); } for ( i = 1; i <= n; i++ ) { i1 = - p[i-1]; if ( 0 <= i1 ) { i0 = i; for ( ; ; ) { i2 = p[i1-1]; p[i1-1] = i0; if ( i2 < 0 ) { break; } i0 = i1; i1 = i2; } } } /* Now we can restore the permutation. */ for ( i = 0; i < n; i++ ) { p[i] = p[i] + p_min - base; } return; } /******************************************************************************/ void perm_print ( int n, int p[], char *title ) /******************************************************************************/ /* Purpose: PERM_PRINT prints a permutation. Example: Input: P = 7 2 4 1 5 3 6 Printed output: "This is the permutation:" 1 2 3 4 5 6 7 7 2 4 1 5 3 6 Licensing: This code is distributed under the GNU LGPL license. Modified: 29 April 2003 Author: John Burkardt Parameters: Input, int N, the number of objects permuted. Input, int P[N], the permutation, in standard index form. Input, char *TITLE, a title. */ { int i; int ihi; int ilo; int inc = 20; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); for ( ilo = 0; ilo < n; ilo = ilo + inc ) { ihi = ilo + inc; if ( n < ihi ) { ihi = n; } fprintf ( stdout, "\n" ); for ( i = ilo; i < ihi; i++ ) { fprintf ( stdout, "%4d", i + 1 ); } fprintf ( stdout, "\n" ); for ( i = ilo; i < ihi; i++ ) { fprintf ( stdout, "%4d", p[i] ); } fprintf ( stdout, "\n" ); } return; } /******************************************************************************/ int *perm_uniform_new ( int n, int base, int *seed ) /******************************************************************************/ /* Purpose: PERM_UNIFORM_NEW selects a random permutation of N objects. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 June 2010 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the number of objects to be permuted. Input, int BASE, is 0 for a 0-based permutation and 1 for a 1-based permutation. Input/output, int *SEED, a seed for the random number generator. Output, int PERM_UNIFORM_NEW[N], a permutation of (BASE, BASE+1, ..., BASE+N-1). */ { int i; int j; int k; int *p; p = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { p[i] = i + base; } for ( i = 0; i < n; i++ ) { j = i4_uniform_ab ( i, n - 1, seed ); k = p[i]; p[i] = p[j]; p[j] = k; } return p; } /******************************************************************************/ double pounds_to_kilograms ( double lb ) /******************************************************************************/ /* Purpose: POUNDS_TO_KILOGRAMS converts a measurement in pounds to kilograms. Licensing: This code is distributed under the GNU LGPL license. Modified: 24 March 2004 Author: John Burkardt Parameters: Input, double LB, the weight in pounds. Output, double POUNDS_TO_KILOGRAMS, the corresponding weight in kilograms. */ { double value; value = 0.4535924 * lb; return value; } /******************************************************************************/ int prime ( int n ) /******************************************************************************/ /* Purpose: PRIME returns any of the first PRIME_MAX prime numbers. Discussion: PRIME_MAX is 1600, and the largest prime stored is 13499. Thanks to Bart Vandewoestyne for pointing out a typo, 18 February 2005. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 March 2009 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, US Department of Commerce, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Daniel Zwillinger, CRC Standard Mathematical Tables and Formulae, 30th Edition, CRC Press, 1996, pages 95-98. Parameters: Input, int N, the index of the desired prime number. In general, is should be true that 0 <= N <= PRIME_MAX. N = -1 returns PRIME_MAX, the index of the largest prime available. N = 0 is legal, returning PRIME = 1. Output, int PRIME, the N-th prime. If N is out of range, PRIME is returned as -1. */ { # define PRIME_MAX 1600 int npvec[PRIME_MAX] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219, 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291, 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387, 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501, 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597, 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677, 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741, 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831, 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929, 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011, 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109, 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199, 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283, 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377, 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439, 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533, 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631, 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733, 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811, 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887, 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007, 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099, 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177, 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271, 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343, 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459, 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567, 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657, 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739, 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859, 10861,10867,10883,10889,10891,10903,10909,10937,10939,10949, 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059, 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149, 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251, 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329, 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443, 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527, 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657, 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777, 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833, 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933, 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011, 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109, 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211, 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289, 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401, 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487, 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553, 12569,12577,12583,12589,12601,12611,12613,12619,12637,12641, 12647,12653,12659,12671,12689,12697,12703,12713,12721,12739, 12743,12757,12763,12781,12791,12799,12809,12821,12823,12829, 12841,12853,12889,12893,12899,12907,12911,12917,12919,12923, 12941,12953,12959,12967,12973,12979,12983,13001,13003,13007, 13009,13033,13037,13043,13049,13063,13093,13099,13103,13109, 13121,13127,13147,13151,13159,13163,13171,13177,13183,13187, 13217,13219,13229,13241,13249,13259,13267,13291,13297,13309, 13313,13327,13331,13337,13339,13367,13381,13397,13399,13411, 13417,13421,13441,13451,13457,13463,13469,13477,13487,13499 }; if ( n == -1 ) { return PRIME_MAX; } else if ( n == 0 ) { return 1; } else if ( n <= PRIME_MAX ) { return npvec[n-1]; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "PRIME - Fatal error!\n" ); fprintf ( stderr, " Unexpected input value of n = %d\n", n ); exit ( 1 ); } return 0; # undef PRIME_MAX } /******************************************************************************/ int prime_ge ( int n ) /******************************************************************************/ /* Purpose: PRIME_GE returns the smallest prime greater than or equal to N. Example: N PRIME_GE -10 2 1 2 2 2 3 3 4 5 5 5 6 7 7 7 8 11 9 11 10 11 Licensing: This code is distributed under the GNU LGPL license. Modified: 25 March 2009 Author: John Burkardt Parameters: Input, int N, the number to be bounded. Output, int PRIME_GE, the smallest prime number that is greater than or equal to N. However, if N is larger than the largest prime stored, then PRIME_GE is returned as -1. */ { int i_hi; int i_lo; int i_mid; int p; int p_hi; int p_mid; if ( n <= 2 ) { p = 2; } else { i_lo = 1; i_hi = prime(-1); p_hi = prime(i_hi); if ( p_hi < n ) { p = - p_hi; } else { for ( ; ; ) { if ( i_lo + 1 == i_hi ) { p = p_hi; break; } i_mid = ( i_lo + i_hi ) / 2; p_mid = prime(i_mid); if ( p_mid < n ) { i_lo = i_mid; } else if ( n <= p_mid ) { i_hi = i_mid; p_hi = p_mid; } } } } return p; } /******************************************************************************/ int r4_nint ( float x ) /******************************************************************************/ /* Purpose: R4_NINT returns the nearest integer to an R4. Example: X R4_NINT 1.3 1 1.4 1 1.5 1 or 2 1.6 2 0.0 0 -0.7 -1 -1.1 -1 -1.6 -2 Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, float X, the value. Output, int R4_NINT, the nearest integer to X. */ { int s; int value; if ( x < 0.0 ) { s = - 1; } else { s = + 1; } value = s * ( int ) ( fabs ( x ) + 0.5 ); return value; } /******************************************************************************/ double r8_abs ( double x ) /******************************************************************************/ /* Purpose: R8_ABS returns the absolute value of an R8. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, the quantity whose absolute value is desired. Output, double R8_ABS, the absolute value of X. */ { double value; if ( 0.0 <= x ) { value = x; } else { value = - x; } return value; } /******************************************************************************/ double r8_huge ( ) /******************************************************************************/ /* Purpose: R8_HUGE returns a "huge" R8. Discussion: The value returned by this function is NOT required to be the maximum representable R8. This value varies from machine to machine, from compiler to compiler, and may cause problems when being printed. We simply want a "very large" but non-infinite number. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 October 2007 Author: John Burkardt Parameters: Output, double R8_HUGE, a "huge" R8 value. */ { double value; value = 1.0E+30; return value; } /******************************************************************************/ double r8_log_10 ( double x ) /******************************************************************************/ /* Purpose: R8_LOG_10 returns the logarithm base 10 of an R8. Discussion: R8_LOG_10 ( X ) = Log10 ( |X| ) Licensing: This code is distributed under the GNU LGPL license. Modified: 11 May 2006 Author: John Burkardt Parameters: Input, double X, the number whose base 2 logarithm is desired. X should not be 0. Output, double R8_LOG_10, the logarithm base 10 of the absolute value of X. It should be true that |X| = 10**R_LOG_10. */ { double value; if ( x == 0.0 ) { value = - r8_huge ( ); } else { value = log10 ( fabs ( x ) ); } return value; } /******************************************************************************/ double r8_max ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MAX returns the maximum of two R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 07 May 2006 Author: John Burkardt Parameters: Input, double X, Y, the quantities to compare. Output, double R8_MAX, the maximum of X and Y. */ { double value; if ( y < x ) { value = x; } else { value = y; } return value; } /******************************************************************************/ double r8_modp ( double x, double y ) /******************************************************************************/ /* Purpose: R8_MODP returns the nonnegative remainder of R8 division. Formula: If REM = R8_MODP ( X, Y ) RMULT = ( X - REM ) / Y then X = Y * RMULT + REM where REM is always nonnegative. Discussion: The MOD function computes a result with the same sign as the quantity being divided. Thus, suppose you had an angle A, and you wanted to ensure that it was between 0 and 360. Then mod(A,360.0) would do, if A was positive, but if A was negative, your result would be between -360 and 0. On the other hand, R8_MODP(A,360.0) is between 0 and 360, always. Example: I J MOD R8_MODP R8_MODP Factorization 107 50 7 7 107 = 2 * 50 + 7 107 -50 7 7 107 = -2 * -50 + 7 -107 50 -7 43 -107 = -3 * 50 + 43 -107 -50 -7 43 -107 = 3 * -50 + 43 Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, double X, the number to be divided. Input, double Y, the number that divides X. Output, double R8_MODP, the nonnegative remainder when X is divided by Y. */ { double value; if ( y == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8_MODP - Fatal error!\n" ); fprintf ( stderr, " R8_MODP ( X, Y ) called with Y = %f\n", y ); exit ( 1 ); } value = x - ( ( double ) ( ( int ) ( x / y ) ) ) * y; if ( value < 0.0 ) { value = value + fabs ( y ); } return value; } /******************************************************************************/ double r8_uniform ( double b, double c, int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM returns a scaled pseudorandom R8. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the GNU LGPL license. Modified: 21 November 2004 Author: John Burkardt Parameters: Input, double B, C, the limits of the interval. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, double R8_UNIFORM, a number strictly between A and B. */ { double value; value = b + ( c - b ) * r8_uniform_01 ( seed ); return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a unit pseudorandom R8. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the GNU LGPL license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } /* Although SEED can be represented exactly as a 32 bit integer, it generally cannot be represented exactly as a 32 bit real number! */ r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8's, which may be stored as a vector in column-major order. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; j2hi = i4_min ( j2hi, n ); j2hi = i4_min ( j2hi, jhi ); fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ i2lo = i4_max ( ilo, 1 ); i2hi = i4_min ( ihi, m ); for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14f", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Licensing: This code is distributed under the GNU LGPL license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14f", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ int r8poly_degree ( int na, double a[] ) /******************************************************************************/ /* Purpose: R8POLY_DEGREE returns the degree of a polynomial. Discussion: The degree of a polynomial is the index of the highest power of X with a nonzero coefficient. The degree of a constant polynomial is 0. The degree of the zero polynomial is debatable, but this routine returns the degree as 0. Licensing: This code is distributed under the GNU LGPL license. Modified: 15 May 2010 Author: John Burkardt Parameters: Input, int NA, the dimension of A. Input, double A[NA+1], the coefficients of the polynomials. Output, int R8POLY_DEGREE, the degree of A. */ { int degree; degree = na; while ( 0 < degree ) { if ( a[degree] != 0.0 ) { return degree; } degree = degree - 1; } return degree; } /******************************************************************************/ void r8poly_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8POLY_PRINT prints out a polynomial. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Input, int N, the dimension of A. Input, double A[N+1], the polynomial coefficients. A(0) is the constant term and A(N) is the coefficient of X**N. Input, char *TITLE, a title. */ { int i; double mag; int n2; char plus_minus; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); n2 = r8poly_degree ( n, a ); if ( n2 <= 0 ) { printf ( " p(x) = 0\n" ); return; } if ( a[n2] < 0.0 ) { plus_minus = '-'; } else { plus_minus = ' '; } mag = r8_abs ( a[n2] ); if ( 2 <= n2 ) { printf ( " p(x) = %c%f * x^%d\n", plus_minus, mag, n2 ); } else if ( n2 == 1 ) { printf ( " p(x) = %c%f * x\n", plus_minus, mag ); } else if ( n2 == 0 ) { printf ( " p(x) = %c%f\n", plus_minus, mag ); } for ( i = n2 - 1; 0 <= i; i-- ) { if ( a[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = r8_abs ( a[i] ); if ( mag != 0.0 ) { if ( 2 <= i ) { printf ( " %c%f * x^%d\n", plus_minus, mag, i ); } else if ( i == 1 ) { printf ( " %c%f * x\n", plus_minus, mag ); } else if ( i == 0 ) { printf ( " %c%f\n", plus_minus, mag ); } } } return; } /******************************************************************************/ double *r8vec_indicator_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_INDICATOR_NEW sets an R8VEC to the indicator vector {1,2,3...}. Licensing: This code is distributed under the GNU LGPL license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of elements of A. Output, double R8VEC_INDICATOR_NEW[N], the array. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i <= n-1; i++ ) { a[i] = ( double ) ( i + 1 ); } return a; } /******************************************************************************/ double r8vec_max ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MAX returns the value of the maximum element in a R8VEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], a pointer to the first entry of the array. Output, double R8VEC_MAX, the value of the maximum element. This is set to 0.0 if N <= 0. */ { int i; double value; value = -r8_huge ( ); if ( n <= 0 ) { return value; } for ( i = 0; i < n; i++ ) { if ( value < r8vec[i] ) { value = r8vec[i]; } } return value; } /******************************************************************************/ double r8vec_mean ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_MEAN returns the mean of a R8VEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector whose mean is desired. Output, double R8VEC_MEAN, the mean, or average, of the vector entries. */ { int i; double mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + x[i]; } mean = mean / ( double ) n; return mean; } /******************************************************************************/ double r8vec_min ( int n, double r8vec[] ) /******************************************************************************/ /* Purpose: R8VEC_MIN returns the value of the minimum element in a R8VEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double R8VEC[N], the array to be checked. Output, double R8VEC_MIN, the value of the minimum element. */ { int i; double value; value = r8_huge ( ); if ( n <= 0 ) { return value; } for ( i = 0; i < n; i++ ) { if ( r8vec[i] < value ) { value = r8vec[i]; } } return value; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14f\n", i, a[i] ); } return; } /******************************************************************************/ double r8vec_variance ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_VARIANCE returns the variance of an R8VEC. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 May 2006 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector whose variance is desired. Output, double R8VEC_VARIANCE, the variance of the vector entries. */ { int i; double mean; double variance; mean = r8vec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( x[i] - mean ) * ( x[i] - mean ); } if ( 1 < n ) { variance = variance / ( double ) ( n - 1 ); } else { variance = 0.0; } return variance; } /******************************************************************************/ double *r8vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_ZERO_NEW creates and zeroes an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the GNU LGPL license. Modified: 25 March 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double R8VEC_ZERO_NEW[N], a vector of zeroes. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return a; } /******************************************************************************/ double radians_to_degrees ( double angle ) /******************************************************************************/ /* Purpose: RADIANS_TO_DEGREES converts an angle from radians to degrees. Licensing: This code is distributed under the GNU LGPL license. Modified: 27 August 2003 Author: John Burkardt Parameters: Input, double ANGLE, an angle in radians. Output, double RADIANS_TO_DEGREES, the equivalent angle in degrees. */ { double pi = 3.141592653589793; double value; value = ( angle / pi ) * 180.0; return value; } /******************************************************************************/ unsigned long rand_initialize ( unsigned long seed ) /******************************************************************************/ /* Purpose: RAND_INITIALIZE initializes the random number generator. Discussion: If you don't initialize RAND, the random number generator, it will behave as though it were seeded with value 1. This routine will either take a user-specified seed, or (if the user passes a 0) make up a "random" one. In either case, the seed is passed to SRAND (the appropriate routine to call when setting the seed for RAND). The seed is also returned to the user as the value of the function. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Input, unsigned long SEED, is either 0, which means that the user wants this routine to come up with a seed, or nonzero, in which case the user has supplied the seed. Output, unsigned long RAND_INITIALIZE, is the value of the seed passed to SRAND, which is either the user's input value, or if that was zero, the value selected by this routine. */ { if ( seed != 0 ) { printf ( "\n" ); printf ( "RAND_INITIALIZE:\n" ); printf ( " Initialize RAND with user SEED = %lud\n", seed ); } else { seed = get_seed ( ); printf ( "\n" ); printf ( "RAND_INITIALIZE:\n" ); printf ( " Initialize RAND with arbitrary SEED = %lud\n", seed ); } /* Now set the seed. */ srand ( seed ); return seed; } /******************************************************************************/ unsigned long random_initialize ( unsigned long seed ) /******************************************************************************/ /* Purpose: RANDOM_INITIALIZE initializes the RANDOM random number generator. Discussion: If you don't initialize RANDOM, the random number generator, it will behave as though it were seeded with value 1. This routine will either take a user-specified seed, or (if the user passes a 0) make up a "random" one. In either case, the seed is passed to SRANDOM (the appropriate routine to call when setting the seed for RANDOM). The seed is also returned to the user as the value of the function. Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Input, unsigned long SEED, is either 0, which means that the user wants this routine to come up with a seed, or nonzero, in which case the user has supplied the seed. Output, unsigned long RANDOM_INITIALIZE, is the value of the seed passed to SRANDOM, which is either the user's input value, or if that was zero, the value selected by this routine. */ { # define DEBUG 0 if ( seed != 0 ) { if ( DEBUG ) { printf ( "\n" ); printf ( "RANDOM_INITIALIZE:\n" ); printf ( " Initialize RANDOM with user SEED = %lud\n", seed ); } } else { seed = get_seed ( ); if ( DEBUG ) { printf ( "\n" ); printf ( "RANDOM_INITIALIZE:\n" ); printf ( " Initialize RANDOM with arbitrary SEED = %lud\n", seed ); } } /* Now set the seed. */ srandom ( seed ); return seed; # undef DEBUG } /******************************************************************************/ void rat_factor ( int m, int n, int factor_max, int *factor_num, int factor[], int power[], int *mleft, int *nleft ) /******************************************************************************/ /* Purpose: RAT_FACTOR factors a rational value into a product of prime factors. Discussion: ( M / N ) = ( MLEFT / NLEFT ) * Product ( I = 1 to FACTOR_NUM ) FACTOR(I)**POWER(I). Licensing: This code is distributed under the GNU LGPL license. Modified: 23 July 2010 Author: John Burkardt Parameters: Input, int M, N, the top and bottom of a rational value. The ratio of M and N must be positive. Input, int FACTOR_MAX, the maximum number of factors for which storage has been allocated. Output, int *FACTOR_NUM, the number of prime factors of M/N. Output, int FACTOR[FACTOR_MAX], the prime factors of M/N. Output, int POWER[FACTOR_MAX]. POWER(I) is the power of the FACTOR(I) in the representation of M/N. Output, int *MLEFT, *NLEFT, the top and bottom of the factor of M / N that remains. If ABS ( MLEFT / NLEFT ) is not 1, then the rational value was not completely factored. */ { int i; int p; int prime_max; *factor_num = 0; *mleft = m; *nleft = n; /* NLEFT should be nonnegative. */ if ( *nleft < 0 ) { *mleft = -(*mleft); *nleft = -(*nleft); } if ( m == 0 || n == 0 ) { return; } if ( m == n ) { *factor_num = 1; factor[0] = 1; power[0] = 1; return; } /* Find out how many primes we stored. */ prime_max = prime ( -1 ); for ( i = 1; i <= prime_max; i++ ) { p = prime ( i ); if ( ( *nleft % p ) == 0 || ( abs ( *mleft ) % p ) == 0 ) { if ( *factor_num < factor_max ) { *factor_num = *factor_num + 1; factor[*factor_num-1] = p; power[*factor_num-1] = 0; /* Divide MLEFT by PRIME(I) as often as you can. */ if ( ( abs ( *mleft ) % p ) == 0 ) { for ( ; ; ) { power[*factor_num-1] = power[*factor_num-1] + 1; *mleft = *mleft / p; if ( ( abs ( *mleft ) % p ) != 0 ) { break; } } } /* Divide NLEFT by PRIME(I) as often as you can. */ if ( ( *nleft % p ) == 0 ) { for ( ; ; ) { power[*factor_num-1] = power[*factor_num-1] - 1; *nleft = *nleft / p; if ( ( *nleft % p ) != 0 ) { break; } } } if ( power[*factor_num-1] == 0 ) { *factor_num = *factor_num - 1; } } } } return; } /******************************************************************************/ double rickey ( int ab, int bb, int er, double f, int h, int hb, int hp, int r, int so, int tb ) /******************************************************************************/ /* Purpose: RICKEY evaluates Branch Rickey's baseball index. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 October 2005 Author: John Burkardt Reference: Alan Schwarz, Looking Beyond the Batting Average, The New York Times, Sunday, 1 August 2004. Branch Rickey, Goodby to Some Old Baseball Ideas, Life Magazine, 2 August 1954. Parameters: Input, int AB, number of at-bats. Input, int BB, base on balls. Input, int ER, earned runs. Input, double F, the fielding rating. Input, int H, number of hits. Input, int HB, hit batsmen. Input, int HP, hit by pitcher. Input, int R, runs. Input, int SO, strike outs. Input, int TB, total bases. Output, double RICKEY, the Branch Rickey index, an estimate for the expected winning percentage of a team with the given statistics. (0.5 has already been subtracted from this value.) */ { double g;; double hitting; double pitching; hitting = ( double ) ( h + bb + hp ) / ( double ) ( ab + bb + hp ) + ( double ) ( 3 * ( tb - h ) ) / ( double ) ( 4 * ab ) + ( double ) ( r ) / ( double ) ( h + bb + hp ); pitching = ( double ) ( h ) / ( double ) ( ab ) + ( double ) ( bb + hb ) / ( double ) ( ab + bb + hb ) + ( double ) ( er ) / ( double ) ( h + bb + hb ) - ( double ) ( so ) / ( double ) ( 8 * ( ab + bb + hb ) ); g = hitting - pitching - f; return g; } /******************************************************************************/ int *roots_to_i4poly ( int n, int x[] ) /******************************************************************************/ /* Purpose: ROOTS_TO_I4POLY converts polynomial roots to polynomial coefficients. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 July 2010 Author: John Burkardt Parameters: Input, int N, the number of roots specified. Input, int X[N], the roots. Output, int ROOTS_TO_I4POLY[N+1], the coefficients of the polynomial. */ { int *c; int i; int j; c = i4vec_zero_new ( n + 1 ); /* Initialize C to (0, 0, ..., 0, 1). Essentially, we are setting up a divided difference table. */ c[n] = 1; /* Convert to standard polynomial form by shifting the abscissas of the divided difference table to 0. */ for ( j = 1; j <= n; j++ ) { for ( i = 1; i <= n+1-j; i++ ) { c[n-i] = c[n-i] - x[n+1-i-j] * c[n-i+1]; } } return c; } /******************************************************************************/ double *roots_to_r8poly ( int n, double x[] ) /******************************************************************************/ /* Purpose: ROOTS_TO_R8POLY converts polynomial roots to polynomial coefficients. Licensing: This code is distributed under the GNU LGPL license. Modified: 20 July 2010 Author: John Burkardt Parameters: Input, int N, the number of roots specified. Input, double X[N], the roots. Output, double ROOTS_TO_R8POLY[N+1], the coefficients of the polynomial. */ { double *c; int i; int j; c = r8vec_zero_new ( n + 1 ); /* Initialize C to (0, 0, ..., 0, 1). Essentially, we are setting up a divided difference table. */ c[n] = 1.0; /* Convert to standard polynomial form by shifting the abscissas of the divided difference table to 0. */ for ( j = 1; j <= n; j++ ) { for ( i = 1; i <= n+1-j; i++ ) { c[n-i] = c[n-i] - x[n+1-i-j] * c[n-i+1]; } } return c; } /******************************************************************************/ int s_eqi ( char *s1, char *s2 ) /******************************************************************************/ /* Purpose: S_EQI reports whether two strings are equal, ignoring case. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 June 2008 Author: John Burkardt Parameters: Input, char *S1, char *S2, pointers to two strings. Output, int S_EQI, is true if the strings are equal. */ { int i; int nchar; int nchar1; int nchar2; nchar1 = strlen ( s1 ); nchar2 = strlen ( s2 ); if ( nchar1 < nchar2 ) { nchar = nchar1; } else { nchar = nchar2; } /* The strings are not equal if they differ over their common length. */ for ( i = 0; i < nchar; i++ ) { if ( ch_cap ( s1[i] ) != ch_cap ( s2[i] ) ) { return 0; } } /* The strings are not equal if the longer one includes nonblanks in the tail. */ if ( nchar < nchar1 ) { for ( i = nchar; i < nchar1; i++ ) { if ( s1[i] != ' ' ) { return 0; } } } else if ( nchar < nchar2 ) { for ( i = nchar; i < nchar2; i++ ) { if ( s2[i] != ' ' ) { return 0; } } } return 1; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Discussion: It turns out that I also want to ignore the '\n' character! Licensing: This code is distributed under the GNU LGPL license. Modified: 05 October 2014 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' && *t != '\n' ) { return n; } t--; n--; } return n; } /******************************************************************************/ int s_to_i4 ( char *s, int *last, int *error ) /******************************************************************************/ /* Purpose: S_TO_I4 reads an I4 from a string. Licensing: This code is distributed under the GNU LGPL license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *S, a string to be examined. Output, int *LAST, the last character of S used to make IVAL. Output, int *ERROR is TRUE (1) if an error occurred and FALSE (0) otherwise. Output, int *S_TO_I4, the integer value read from the string. If the string is blank, then IVAL will be returned 0. */ { char c; int i; int isgn; int istate; int ival; *error = 0; istate = 0; isgn = 1; i = 0; ival = 0; while ( *s ) { c = s[i]; i = i + 1; /* Haven't read anything. */ if ( istate == 0 ) { if ( c == ' ' ) { } else if ( c == '-' ) { istate = 1; isgn = -1; } else if ( c == '+' ) { istate = 1; isgn = + 1; } else if ( '0' <= c && c <= '9' ) { istate = 2; ival = c - '0'; } else { *error = 1; return ival; } } /* Have read the sign, expecting digits. */ else if ( istate == 1 ) { if ( c == ' ' ) { } else if ( '0' <= c && c <= '9' ) { istate = 2; ival = c - '0'; } else { *error = 1; return ival; } } /* Have read at least one digit, expecting more. */ else if ( istate == 2 ) { if ( '0' <= c && c <= '9' ) { ival = 10 * (ival) + c - '0'; } else { ival = isgn * ival; *last = i - 1; return ival; } } } /* If we read all the characters in the string, see if we're OK. */ if ( istate == 2 ) { ival = isgn * ival; *last = s_len_trim ( s ); } else { *error = 1; *last = 0; } return ival; } /******************************************************************************/ int s_to_i4vec ( char *s, int n, int ivec[] ) /******************************************************************************/ /* Purpose: S_TO_I4VEC reads an I4VEC from a string. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 February 2001 Author: John Burkardt Parameters: Input, char *S, the string to be read. Input, int N, the number of values expected. Output, int IVEC[N], the values read from the string. Output, int S_TO_I4VEC, is TRUE (1) if an error occurred and FALSE (0) otherwise. */ { int error; int i; int lchar; error = 0; for ( i = 0; i < n; i++ ) { ivec[i] = s_to_i4 ( s, &lchar, &error ); if ( error ) { return error; } s = s + lchar; } return error; } /******************************************************************************/ double s_to_r8 ( char *s, int *lchar, int *error ) /******************************************************************************/ /* Purpose: S_TO_R8 reads an R8 value from a string. Discussion: We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. This routine will read as many characters as possible until it reaches the end of the string, or encounters a character which cannot be part of the real number. Legal input is: 1 blanks, 2 '+' or '-' sign, 2.5 spaces 3 integer part, 4 decimal point, 5 fraction part, 6 'E' or 'e' or 'D' or 'd', exponent marker, 7 exponent sign, 8 exponent integer part, 9 exponent decimal point, 10 exponent fraction part, 11 blanks, 12 final comma or semicolon. with most quantities optional. Example: S R '1' 1.0 ' 1 ' 1.0 '1A' 1.0 '12,34,56' 12.0 ' 34 7' 34.0 '-1E2ABCD' -100.0 '-1X2ABCD' -1.0 ' 2E-1' 0.2 '23.45' 23.45 '-4.2E+2' -420.0 '17d2' 1700.0 '-14e-2' -0.14 'e2' 100.0 '-12.73e-9.23' -12.73 * 10.0**(-9.23) Licensing: This code is distributed under the GNU LGPL license. Modified: 24 June 2005 Author: John Burkardt Parameters: Input, char *S, the string containing the data to be read. Reading will begin at position 1 and terminate at the end of the string, or when no more characters can be read to form a legal real. Blanks, commas, or other nonnumeric data will, in particular, cause the conversion to halt. Output, int *LCHAR, the number of characters read from the string to form the number, including any terminating characters such as a trailing comma or blanks. Output, int *ERROR, is TRUE (1) if an error occurred and FALSE (0) otherwise. Output, double S_TO_R8, the value that was read from the string. */ { char c; int ihave; int isgn; int iterm; int jbot; int jsgn; int jtop; int nchar; int ndig; double r; double rbot; double rexp; double rtop; char TAB = 9; nchar = s_len_trim ( s ); *error = 0; r = 0.0; *lchar = -1; isgn = 1; rtop = 0.0; rbot = 1.0; jsgn = 1; jtop = 0; jbot = 1; ihave = 1; iterm = 0; for ( ; ; ) { c = s[*lchar+1]; *lchar = *lchar + 1; /* Blank or TAB character. */ if ( c == ' ' || c == TAB ) { if ( ihave == 2 ) { } else if ( ihave == 6 || ihave == 7 ) { iterm = 1; } else if ( 1 < ihave ) { ihave = 11; } } /* Comma. */ else if ( c == ',' || c == ';' ) { if ( ihave != 1 ) { iterm = 1; ihave = 12; *lchar = *lchar + 1; } } /* Minus sign. */ else if ( c == '-' ) { if ( ihave == 1 ) { ihave = 2; isgn = -1; } else if ( ihave == 6 ) { ihave = 7; jsgn = -1; } else { iterm = 1; } } /* Plus sign. */ else if ( c == '+' ) { if ( ihave == 1 ) { ihave = 2; } else if ( ihave == 6 ) { ihave = 7; } else { iterm = 1; } } /* Decimal point. */ else if ( c == '.' ) { if ( ihave < 4 ) { ihave = 4; } else if ( 6 <= ihave && ihave <= 8 ) { ihave = 9; } else { iterm = 1; } } /* Exponent marker. */ else if ( ch_eqi ( c, 'E' ) || ch_eqi ( c, 'D' ) ) { if ( ihave < 6 ) { ihave = 6; } else { iterm = 1; } } /* Digit. */ else if ( ihave < 11 && '0' <= c && c <= '9' ) { if ( ihave <= 2 ) { ihave = 3; } else if ( ihave == 4 ) { ihave = 5; } else if ( ihave == 6 || ihave == 7 ) { ihave = 8; } else if ( ihave == 9 ) { ihave = 10; } ndig = ch_to_digit ( c ); if ( ihave == 3 ) { rtop = 10.0 * rtop + ( double ) ndig; } else if ( ihave == 5 ) { rtop = 10.0 * rtop + ( double ) ndig; rbot = 10.0 * rbot; } else if ( ihave == 8 ) { jtop = 10 * jtop + ndig; } else if ( ihave == 10 ) { jtop = 10 * jtop + ndig; jbot = 10 * jbot; } } /* Anything else is regarded as a terminator. */ else { iterm = 1; } /* If we haven't seen a terminator, and we haven't examined the entire string, go get the next character. */ if ( iterm == 1 || nchar <= *lchar + 1 ) { break; } } /* If we haven't seen a terminator, and we have examined the entire string, then we're done, and LCHAR is equal to NCHAR. */ if ( iterm != 1 && (*lchar) + 1 == nchar ) { *lchar = nchar; } /* Number seems to have terminated. Have we got a legal number? Not if we terminated in states 1, 2, 6 or 7! */ if ( ihave == 1 || ihave == 2 || ihave == 6 || ihave == 7 ) { *error = 1; return r; } /* Number seems OK. Form it. We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. */ if ( jtop == 0 ) { rexp = 1.0; } else { if ( jbot == 1 ) { if ( jsgn * jtop < -20 ) { rexp = 0.0; } else { rexp = pow ( ( double ) 10.0, ( double ) ( jsgn * jtop ) ); } } else { if ( jsgn * jtop < -20 * jbot ) { rexp = 0.0; } else { rexp = jsgn * jtop; rexp = rexp / jbot; rexp = pow ( ( double ) 10.0, ( double ) rexp ); } } } r = isgn * rexp * rtop / rbot; return r; } /******************************************************************************/ int s_to_r8vec ( char *s, int n, double rvec[] ) /******************************************************************************/ /* Purpose: S_TO_R8VEC reads an R8VEC from a string. Licensing: This code is distributed under the GNU LGPL license. Modified: 19 February 2001 Author: John Burkardt Parameters: Input, char *S, the string to be read. Input, int N, the number of values expected. Output, double RVEC[N], the values read from the string. Output, int S_TO_R8VEC, is TRUE (1) if an error occurred and FALSE (0) otherwise. */ { int error; int i; int lchar; error = 0; for ( i = 0; i < n; i++ ) { rvec[i] = s_to_r8 ( s, &lchar, &error ); if ( error ) { return error; } s = s + lchar; } return error; } /******************************************************************************/ void sort_heap_external ( int n, int *indx, int *i, int *j, int isgn ) /******************************************************************************/ /* Purpose: SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. Discussion: The actual list is not passed to the routine. Hence it may consist of integers, reals, numbers, names, etc. The user, after each return from the routine, will be asked to compare or interchange two items. The current version of this code mimics the FORTRAN version, so the values of I and J, in particular, are FORTRAN indices. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 February 2004 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C version by John Burkardt. Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the length of the input list. Input/output, int *INDX. The user must set INDX to 0 before the first call. On return, if INDX is greater than 0, the user must interchange items I and J and recall the routine. If INDX is less than 0, the user is to compare items I and J and return in ISGN a negative value if I is to precede J, and a positive value otherwise. If INDX is 0, the sorting is done. Output, int *I, *J. On return with INDX positive, elements I and J of the user's list should be interchanged. On return with INDX negative, elements I and J are to be compared by the user. Input, int ISGN. On return with INDX negative, the user should compare elements I and J of the list. If item I is to precede item J, set ISGN negative, otherwise set ISGN positive. */ { static int i_save = 0; static int j_save = 0; static int k = 0; static int k1 = 0; static int n1 = 0; /* INDX = 0: This is the first call. */ if ( *indx == 0 ) { i_save = 0; j_save = 0; k = n / 2; k1 = k; n1 = n; } /* INDX < 0: The user is returning the results of a comparison. */ else if ( *indx < 0 ) { if ( *indx == -2 ) { if ( isgn < 0 ) { i_save = i_save + 1; } j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } if ( 0 < isgn ) { *indx = 2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; } *i = i_save; *j = j_save; return; } k = k - 1; k1 = k; } /* 0 < INDX: the user was asked to make an interchange. */ else if ( *indx == 1 ) { k1 = k; } for ( ; ; ) { i_save = 2 * k1; if ( i_save == n1 ) { j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } else if ( i_save <= n1 ) { j_save = i_save + 1; *indx = -2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { break; } k = k - 1; k1 = k; } if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; *i = i_save; *j = j_save; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; *i = i_save; *j = j_save; } return; } /******************************************************************************/ void timestamp ( void ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ void tuple_next2 ( int n, int xmin[], int xmax[], int x[], int *rank ) /******************************************************************************/ /* Purpose: TUPLE_NEXT2 computes the next element of an integer tuple space. Discussion: The elements X are N vectors. Each entry X(I) is constrained to lie between XMIN(I) and XMAX(I). The elements are produced one at a time. The first element is (XMIN(1), XMIN(2), ..., XMIN(N)), the second is (probably) (XMIN(1), XMIN(2), ..., XMIN(N)+1), and the last element is (XMAX(1), XMAX(2), ..., XMAX(N)) Intermediate elements are produced in a lexicographic order, with the first index more important than the last, and the ordering of values at a fixed index implicitly defined by the sign of XMAX(I) - XMIN(I). Examples: N = 2, XMIN = (/ 1, 10 /) XMAX = (/ 3, 8 /) RANK X ---- ----- 1 1 10 2 1 9 3 1 8 4 2 10 5 2 9 6 2 8 7 3 10 8 3 9 9 3 8 Licensing: This code is distributed under the GNU LGPL license. Modified: 07 June 2008 Author: John Burkardt Parameters: Input, int N, the number of components. Input, int XMIN[N], XMAX[N], the "minimum" and "maximum" entry values. These values are minimum and maximum only in the sense of the lexicographic ordering. In fact, XMIN(I) may be less than, equal to, or greater than XMAX(I). Input/output, int X[N], on input the previous tuple. On output, the next tuple. Input/output, int *RANK, the rank of the item. On first call, set RANK to 0 to start up the sequence. On return, if RANK is zero, there are no more items in the sequence. */ { int i; int test; if ( *rank < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TUPLE_NEXT2 - Fatal error!\n" ); fprintf ( stderr, " Illegal value of RANK = %d\n", *rank ); exit ( 1 ); } test = 1; for ( i = 0; i < n; i++ ) { test = test * ( 1 + abs ( xmax[i] - xmin[i] ) ); } if ( test < *rank ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TUPLE_NEXT2 - Fatal error!\n" ); fprintf ( stderr, " Illegal value of RANK = %d\n", *rank ); exit ( 1 ); } if ( *rank == 0 ) { for ( i = 0; i < n; i++ ) { x[i] = xmin[i]; } *rank = 1; return; } *rank = *rank + 1; i = n - 1; for ( ; ; ) { if ( x[i] != xmax[i] ) { if ( xmin[i] < xmax[i] ) { x[i] = x[i] + 1; } else { x[i] = x[i] - 1; } break; } x[i] = xmin[i]; if ( i == 0 ) { *rank = 0; break; } i = i - 1; } return; } /******************************************************************************/ double *tvec_even ( int nt ) /******************************************************************************/ /* Purpose: TVEC_EVEN computes an evenly spaced set of angles between 0 and 2*PI. Discussion: The computation realizes that 0 = 2 * PI. Example: NT = 4 T = ( 0, PI/2, PI, 3*PI/2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Output, double TVEC[NT], the evenly spaced angles, in radians. */ { int i; double pi = 3.141592653589793; double *t; if ( nt < 1 ) { return NULL; } t = ( double * ) malloc ( nt * sizeof ( double ) ); for ( i = 0; i < nt; i++ ) { t[i] = ( double ) ( 2 * i ) * pi / ( double ) ( nt ); } return t; } /******************************************************************************/ double *tvec_even2 ( int nt ) /******************************************************************************/ /* Purpose: TVEC_EVEN2 computes evenly spaced angles between 0 and 2*PI. Discussion: The computation realizes that 0 = 2 * PI. The values are equally spaced in the circle, do not include 0, and are symmetric about 0. Example: NT = 4 T = ( PI/4, 3*PI/4, 5*PI/4, 7*PI/4 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Output, double TVEC[NT], the evenly spaced angles, in radians. */ { int i; double pi = 3.141592653589793; double *t; if ( nt < 1 ) { return NULL; } t = ( double * ) malloc ( nt * sizeof ( double ) ); for ( i = 0; i < nt; i++ ) { t[i] = ( double ) ( 2 * i + 1 ) * pi / ( double ) ( nt ); } return t; } /******************************************************************************/ double *tvec_even3 ( int nt ) /******************************************************************************/ /* Purpose: TVEC_EVEN3 computes an evenly spaced set of angles between 0 and 2*PI. Discussion: The angles begin with 0 and end with 2*PI. Example: NT = 4 T = ( 0, 2*PI/3, 4*PI/3 2*PI ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Output, double TVEC[NT], the evenly spaced angles, in radians. */ { int i; double pi = 3.141592653589793; double *t; if ( nt < 1 ) { return NULL; } t = ( double * ) malloc ( nt * sizeof ( double) ); if ( nt == 1 ) { t[0] = pi; } else { for ( i = 0; i < nt; i++ ) { t[i] = ( double ) ( 2 * i ) * pi / ( double ) ( nt - 1 ); } } return t; } /******************************************************************************/ double *tvec_even_bracket ( int nt, double theta1, double theta2 ) /******************************************************************************/ /* Purpose: TVEC_EVEN_BRACKET computes evenly spaced angles between THETA1 and THETA2. Example: NT = 4 THETA1 = 30 THETA2 = 90 T = ( 30, 50, 70, 90 ) Discussion: The interval between THETA1 and THETA2 is divided into NT-1 subintervals. The angles returned are the breakpoints of these subintervals, including THETA1 and THETA2. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Input, double THETA1, THETA2, the limiting angles. Output, double TVEC_EVEN_BRACKET[NT], the evenly spaced angles. */ { int i; double *t; if ( nt < 1 ) { return NULL; } t = ( double * ) malloc ( nt * sizeof ( double ) ); if ( nt == 1 ) { t[0] = ( theta1 + theta2 ) / 2.0; } else { for ( i = 1; i <= nt; i++ ) { t[i-1] = ( ( double ) ( nt - i ) * theta1 + ( double ) ( i - 1 ) * theta2 ) / ( double ) ( nt - 1 ); } } return t; } /******************************************************************************/ double *tvec_even_bracket2 ( int nt, double theta1, double theta2 ) /******************************************************************************/ /* Purpose: TVEC_EVEN_BRACKET2 computes evenly spaced angles from THETA1 to THETA2. Discussion: The interval between THETA1 and THETA2 is divided into NT+1 subintervals. The angles returned are the internal NT breakpoints of the subintervals. Example: NT = 5 THETA1 = 30 THETA2 = 90 T = ( 40, 50, 60, 70, 80 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Input, double THETA1, THETA2, the limiting angles. Output, double TVEC_EVEN_BRACKET2[NT], the evenly spaced angles. */ { int i; double *t; if ( nt < 1 ) { return NULL; } t = ( double * ) malloc ( nt * sizeof ( double ) ); for ( i = 1; i <= nt; i++ ) { t[i-1] = ( ( double ) ( nt + 1 - i ) * theta1 + ( double ) ( i ) * theta2 ) / ( double ) ( nt + 1 ); } return t; } /******************************************************************************/ double *tvec_even_bracket3 ( int nt, double theta1, double theta2 ) /******************************************************************************/ /* Purpose: TVEC_EVEN_BRACKET3 computes evenly spaced angles from THETA1 to THETA2. Discussion: The interval between THETA1 and THETA2 is divided into NT subintervals. The angles returned are the midpoints of each subinterval. Example: NT = 3 THETA1 = 30 THETA2 = 90 T = ( 40, 60, 80 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 July 2008 Author: John Burkardt Parameters: Input, int NT, the number of values to compute. Input, double THETA1, THETA2, the limiting angles. Output, double TVEC_EVEN_BRACKET3[NT], the evenly spaced angles. */ { int i; double *t; t = ( double * ) malloc ( nt * sizeof ( double ) ); for ( i = 1; i <= nt; i++ ) { t[i-1] = ( ( double ) ( 2 * nt - 2 * i + 1 ) * theta1 + ( double ) ( 2 * i - 1 ) * theta2 ) / ( double ) ( 2 * nt ); } return t; } /******************************************************************************/ double versine_pulse ( double t, double ta, double tb, double v1, double amp ) /******************************************************************************/ /* Purpose: VERSINE_PULSE adds a versine pulse to a constant. Licensing: This code is distributed under the GNU LGPL license. Modified: 01 March 2010 Author: John Burkardt Parameters: Input, double T, the current time. Input, double TA, the time at which the pulse begins. Input, double TB, the time at which the pulse finishes. Input, double V1, the constant value. Input, double AMP, the amplitude of the pulse. Output, double VERSINE_PULSE, the value of the signal at time T. */ { double pi = 3.141592653589793; double v; v = v1; if ( ta <= t && t <= tb ) { v = v + ( 0.5 * amp * ( 1.0 - cos ( 2.0 * pi * ( t - ta ) / ( tb - ta ) ) ) ); } return v; }