subroutine bvec_print ( n, bvec, title ) !*****************************************************************************80 ! !! bvec_print() prints a BVEC, with an optional title. ! ! Discussion: ! ! A BVEC is a vector of binary digits representing an integer. ! ! BVEC(1) is 0 for positive values and 1 for negative values, which ! are stored in 2's complement form. ! ! For positive values, BVEC(N) contains the units digit, BVEC(N-1) ! the coefficient of 2, BVEC(N-2) the coefficient of 4 and so on, ! so that printing the digits in order gives the binary form of the number. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 20 March 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, integer BVEC(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer n integer bvec(n) integer ihi integer ilo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do ilo = 1, n, 70 ihi = min ( ilo + 70 - 1, n ) write ( *, '(2x,70i1)' ) bvec(ilo:ihi) end do return end subroutine bvec_uniform ( n, seed, bvec ) !*****************************************************************************80 ! !! BVEC_UNIFORM returns a pseudorandom BVEC. ! ! Discussion: ! ! An BVEC is a vector of binary (0/1) values, representing an integer. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 December 2014 ! ! 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, integer N, the order of the vector. ! ! Input/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, integer BEC(N), a pseudorandom binary vector. ! implicit none integer n integer bvec(n) integer, parameter :: i4_huge = 2147483647 integer, parameter :: i4_huge_half = 1073741823 integer i integer k integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BVEC_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if if ( i4_huge_half < seed ) then bvec(i) = 0 else bvec(i) = 1 end if end do return end function c8_uniform_01 ( seed ) !*****************************************************************************80 ! !! C8_UNIFORM_01 returns a unit pseudorandom C8. ! ! Discussion: ! ! A C8 is a complex ( kind = ck ) value. ! ! For now, the input quantity SEED is an integer variable. ! ! The angle should be uniformly distributed between 0 and 2 * PI, ! the square root of the radius uniformly distributed between 0 and 1. ! ! This results in a uniform distribution of values in the unit circle. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, complex ( kind = ck ) C8_UNIFORM_01, a pseudorandom complex value. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) complex ( kind = ck ) c8_uniform_01 integer, parameter :: i4_huge = 2147483647 integer k real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) r integer seed real ( kind = rk ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r = sqrt ( real ( seed, kind = rk ) * 4.656612875D-10 ) k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if theta = 2.0D+00 * pi * ( real ( seed, kind = rk ) * 4.656612875D-10 ) c8_uniform_01 = r * cmplx ( cos ( theta ), sin ( theta ), kind = ck ) return end subroutine c8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! C8MAT_PRINT prints a C8MAT. ! ! Discussion: ! ! A C8MAT is a matrix of C8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the matrix. ! ! Input, complex ( kind = ck ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer m integer n complex ( kind = ck ) a(m,n) character ( len = * ) title call c8mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine c8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! C8MAT_PRINT_SOME prints some of a C8MAT. ! ! Discussion: ! ! A C8MAT is a matrix of C8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the matrix. ! ! Input, complex ( kind = ck ) A(M,N), the matrix. ! ! Input, integer ILO, JLO, IHI, JHI, the first row and ! column, and the last row and column to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: incx = 4 integer m integer n complex ( kind = ck ) a(m,n) character ( len = 20 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title complex ( kind = ck ) zero zero = cmplx ( 0.0D+00, 0.0D+00, kind = ck ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if ! ! Print the columns of the matrix, in strips of INCX. ! do j2lo = jlo, min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i10,10x)' ) j end do write ( *, '(a,4a20)' ) ' Col: ', ( ctemp(j2), j2 = 1, inc ) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ---' ! ! Determine the range of the rows in this strip. ! i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi ! ! Print out (up to) INCX entries in row I, that lie in the current strip. ! do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(2g10.3)' ) a(i,j) end do write ( *, '(i5,a1,4a20)' ) i, ':', ( ctemp(j2), j2 = 1, inc ) end do end do return end subroutine c8mat_uniform_01 ( m, n, seed, c ) !*****************************************************************************80 ! !! C8MAT_UNIFORM_01 returns a unit pseudorandom C8MAT. ! ! Discussion: ! ! A C8MAT is a matrix of C8's. ! ! For now, the input quantity SEED is an integer variable. ! ! The angles should be uniformly distributed between 0 and 2 * PI, ! the square roots of the radius uniformly distributed between 0 and 1. ! ! This results in a uniform distribution of values in the unit circle. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer M, N, the number of rows and columns ! in the matrix. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, complex ( kind = ck ) C(M,N), the pseudorandom complex matrix. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n complex ( kind = ck ) c(m,n) integer i integer, parameter :: i4_huge = 2147483647 integer j real ( kind = rk ) r integer k real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 integer seed real ( kind = rk ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8MAT_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r = sqrt ( real ( seed, kind = rk ) * 4.656612875D-10 ) k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if theta = 2.0D+00 * pi * ( real ( seed, kind = rk ) * 4.656612875D-10 ) c(i,j) = r * cmplx ( cos ( theta ), sin ( theta ), kind = ck ) end do end do return end subroutine c8vec_print ( n, a, title ) !*****************************************************************************80 ! !! C8VEC_PRINT prints a C8VEC, with an optional title. ! ! Discussion: ! ! A C8VEC is a vector of C8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, complex ( kind = ck ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer n complex ( kind = ck ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,2g14.6)' ) i, a(i) end do return end subroutine c8vec_uniform_01 ( n, seed, c ) !*****************************************************************************80 ! !! C8VEC_UNIFORM_01 returns a unit pseudorandom C8VEC. ! ! Discussion: ! ! A C8VEC is a vector of C8's. ! ! For now, the input quantity SEED is an integer variable. ! ! The angles should be uniformly distributed between 0 and 2 * PI, ! the square roots of the radius uniformly distributed between 0 and 1. ! ! This results in a uniform distribution of values in the unit circle. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer N, the number of values to compute. ! ! Input/output, integer SEED, the "seed" value, ! which should NOT be 0. ! On output, SEED has been updated. ! ! Output, complex ( kind = ck ) C(N), the pseudorandom complex vector. ! implicit none integer, parameter :: ck = kind ( ( 1.0D+00, 1.0D+00 ) ) integer, parameter :: rk = kind ( 1.0D+00 ) integer n complex ( kind = ck ) c(n) integer i integer, parameter :: i4_huge = 2147483647 real ( kind = rk ) r integer k real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 integer seed real ( kind = rk ) theta if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'C8VEC_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r = sqrt ( real ( seed, kind = rk ) * 4.656612875D-10 ) k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if theta = 2.0D+00 * pi * ( real ( seed, kind = rk ) * 4.656612875D-10 ) c(i) = r * cmplx ( cos ( theta ), sin ( theta ), kind = ck ) end do return end function ch_uniform_ab ( a, b, seed ) !*****************************************************************************80 ! !! CH_UNIFORM_AB returns a scaled pseudorandom CH. ! ! Discussion: ! ! A CH is an alphabetic character value. ! ! The value is scaled to lie between characters A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 October 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character A, B, the minimum and maximum acceptable characters. ! ! Input/output, integer SEED, a seed for the random ! number generator. ! ! Output, character CH_UNIFORM_AB, the randomly chosen character. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character a character b character ch_uniform_ab real ( kind = rk ) r8_uniform_01 integer i integer ihi integer ilo integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CH_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if ilo = ichar ( a ) ihi = ichar ( b ) i = ilo + int ( r8_uniform_01 ( seed ) * real ( ihi + 1 - ilo ) ) i = max ( i, ilo ) i = min ( i, ihi ) ch_uniform_ab = char ( i ) return end subroutine congruence ( a, b, c, ierror, x ) !*****************************************************************************80 ! !! CONGRUENCE solves a congruence of the form ( A * X = C ) mod B. ! ! Discussion: ! ! A, B and C are given integers. The equation is solvable if and only ! if the greatest common divisor of A and B also divides C. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Weisstein, ! CRC Concise Encyclopedia of Mathematics, ! CRC Press, 2002, ! Second edition, ! ISBN: 1584883472, ! LC: QA5.W45. ! ! Parameters: ! ! Input, integer A, B, C, the coefficients of the ! Diophantine equation. ! ! Output, integer IERROR, error flag. ! 0, no error, X was computed. ! 1, A = B = 0, C is nonzero. ! 2, A = 0, B and C nonzero, but C is not a multiple of B. ! 3, A nonzero, B zero, C nonzero, but C is not a multiple of A. ! 4, A, B, C nonzero, but GCD of A and B does not divide C. ! 5, algorithm ran out of internal space. ! ! Output, integer X, the solution of the Diophantine equation. ! X will be between 0 and B-1. ! implicit none integer, parameter :: nmax = 100 integer a integer a_copy integer a_mag integer a_sign integer b integer b_copy integer b_mag integer b_sign integer c integer c_copy integer g integer i4_gcd integer ierror integer k integer n integer, parameter :: one = 1 integer q(nmax) logical swap integer x integer y integer z ! ! Defaults for output parameters. ! ierror = 0 x = 0 y = 0 ! ! Special cases. ! if ( a == 0 .and. b == 0 .and. c == 0 ) then x = 0 return else if ( a == 0 .and. b == 0 .and. c /= 0 ) then ierror = 1 x = 0 return else if ( a == 0 .and. b /= 0 .and. c == 0 ) then x = 0 return else if ( a == 0 .and. b /= 0 .and. c /= 0 ) then x = 0 if ( mod ( c, b ) /= 0 ) then ierror = 2 end if return else if ( a /= 0 .and. b == 0 .and. c == 0 ) then x = 0 return else if ( a /= 0 .and. b == 0 .and. c /= 0 ) then x = c / a if ( mod ( c, a ) /= 0 ) then ierror = 3 end if return else if ( a /= 0 .and. b /= 0 .and. c == 0 ) then ! g = i4_gcd ( a, b ) ! x = b / g x = 0 return end if ! ! Handle the "general" case: A, B and C are nonzero. ! ! Step 1: Compute the GCD of A and B, which must also divide C. ! g = i4_gcd ( a, b ) if ( mod ( c, g ) /= 0 ) then ierror = 4 return end if a_copy = a / g b_copy = b / g c_copy = c / g ! ! Step 2: Split A and B into sign and magnitude. ! a_mag = abs ( a_copy ) a_sign = sign ( one, a_copy ) b_mag = abs ( b_copy ) b_sign = sign ( one, b_copy ) ! ! Another special case, A_MAG = 1 or B_MAG = 1. ! if ( a_mag == 1 ) then x = a_sign * c_copy return else if ( b_mag == 1 ) then x = 0 return end if ! ! Step 3: Produce the Euclidean remainder sequence. ! if ( b_mag <= a_mag ) then swap = .false. q(1) = a_mag q(2) = b_mag else swap = .true. q(1) = b_mag q(2) = a_mag end if n = 3 do q(n) = mod ( q(n-2), q(n-1) ) if ( q(n) == 1 ) then exit end if n = n + 1 if ( nmax < n ) then ierror = 5 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CONGRUENCE - Fatal error!' write ( *, '(a)' ) ' Exceeded number of iterations.' stop 1 end if end do ! ! Step 4: Go backwards to solve X * A_MAG + Y * B_MAG = 1. ! y = 0 do k = n, 2, -1 x = y y = ( 1 - x * q(k-1) ) / q(k) end do ! ! Step 5: Undo the swapping. ! if ( swap ) then z = x x = y y = z end if ! ! Step 6: Apply signs to X and Y so that X * A + Y * B = 1. ! x = x * a_sign ! ! Step 7: Multiply by C, so that X * A + Y * B = C. ! x = x * c_copy ! ! Step 8: Force 0 <= X < B. ! x = mod ( x, b ) ! ! Step 9: Force positivity. ! if ( x < 0 ) then x = x + b end if return end subroutine get_seed ( seed ) !*****************************************************************************80 ! !! GET_SEED returns a seed for the random number generator. ! ! Discussion: ! ! The seed depends on the current time, and ought to be (slightly) ! different every millisecond. Once the seed is obtained, a random ! number generator should be called a few times to further process ! the seed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer SEED, a pseudorandom seed value. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: i4_huge = 2147483647 integer seed real ( kind = rk ) temp character ( len = 10 ) time character ( len = 8 ) today integer values(8) character ( len = 5 ) zone call date_and_time ( today, time, zone, values ) temp = 0.0D+00 temp = temp + real ( values(2) - 1, kind = rk ) / 11.0D+00 temp = temp + real ( values(3) - 1, kind = rk ) / 30.0D+00 temp = temp + real ( values(5), kind = rk ) / 23.0D+00 temp = temp + real ( values(6), kind = rk ) / 59.0D+00 temp = temp + real ( values(7), kind = rk ) / 59.0D+00 temp = temp + real ( values(8), kind = rk ) / 999.0D+00 temp = temp / 6.0D+00 ! ! Force 0 < TEMP <= 1. ! do while ( temp <= 0.0D+00 ) temp = temp + 1.0D+00 end do do while ( 1.0D+00 < temp ) temp = temp - 1.0D+00 end do seed = int ( real ( i4_huge, kind = rk ) * temp ) ! ! Never use a seed of 0 or maximum integer. ! if ( seed == 0 ) then seed = 1 end if if ( seed == i4_huge ) then seed = seed - 1 end if return end function i4_gcd ( i, j ) !*****************************************************************************80 ! !! I4_GCD finds the greatest common divisor of two I4's. ! ! Discussion: ! ! An I4 is an integer value. ! ! Only the absolute values of I and J are ! considered, so that the result is always nonnegative. ! ! If I or J is 0, I4_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). ! ! If I and J have no common factor, I4_GCD is returned as 1. ! ! Otherwise, using the Euclidean algorithm, I4_GCD is the ! largest common factor of I and J. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, two numbers whose greatest common ! divisor is desired. ! ! Output, integer I4_GCD, the greatest common divisor ! of I and J. ! implicit none integer i integer i4_gcd integer ip integer iq integer ir integer j i4_gcd = 1 ! ! Return immediately if either I or J is zero. ! if ( i == 0 ) then i4_gcd = max ( 1, abs ( j ) ) return else if ( j == 0 ) then i4_gcd = max ( 1, abs ( i ) ) return end if ! ! Set IP to the larger of I and J, IQ to the smaller. ! This way, we can alter IP and IQ as we go. ! ip = max ( abs ( i ), abs ( j ) ) iq = min ( abs ( i ), abs ( j ) ) ! ! Carry out the Euclidean algorithm. ! do ir = mod ( ip, iq ) if ( ir == 0 ) then exit end if ip = iq iq = ir end do i4_gcd = iq return end function i4_huge ( ) !*****************************************************************************80 ! !! I4_HUGE returns a "huge" I4. ! ! Discussion: ! ! An I4 is an integer value. ! ! On an IEEE 32 bit machine, I4_HUGE should be 2^31 - 1, and its ! bit pattern should be ! ! 01111111111111111111111111111111 ! ! In this case, its numerical value is 2147483647. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer I4_HUGE, a "huge" I4. ! implicit none integer i4_huge i4_huge = 2147483647 return end function i4_seed_advance ( seed ) !*****************************************************************************80 ! !! I4_SEED_ADVANCE "advances" the seed. ! ! Discussion: ! ! An I4 is an integer value. ! ! This routine implements one step of the recursion ! ! SEED = ( 16807 * SEED ) mod ( 2^31 - 1 ) ! ! This version of the routine does not check whether the input value of ! SEED is zero. If the input value is zero, the output value will be zero. ! ! If we repeatedly use the output of SEED_ADVANCE as the next input, ! and we start with SEED = 12345, then the first few iterates are: ! ! Input Output ! SEED SEED ! ! 12345 207482415 ! 207482415 1790989824 ! 1790989824 2035175616 ! 2035175616 77048696 ! 77048696 24794531 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 April 2013 ! ! 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, integer SEED, the seed value. ! ! Output, integer I4_SEED_ADVANCE, the "next" seed. ! implicit none integer, parameter :: i4_huge = 2147483647 integer i4_seed_advance integer k integer seed integer seed_new seed_new = seed if ( seed_new < 0 ) then seed_new = seed_new + i4_huge end if k = seed_new / 127773 seed_new = 16807 * ( seed_new - k * 127773 ) - k * 2836 if ( seed_new < 0 ) then seed_new = seed_new + i4_huge end if i4_seed_advance = seed_new return end function i4_uniform_0i ( seed ) !*****************************************************************************80 ! !! I4_UNIFORM_0I returns a pseudorandom I4. ! ! Discussion: ! ! An I4 is an integer value. ! ! SEED = ( SEED * 7^5 ) mod (2^31 - 1) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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/output, integer SEED, the integer "seed" used to ! generate the output value. SEED should not be 0. ! ! Output, integer I4_UNIFORM_0I, a uniform random value ! between 1 and 2^31-1. ! ! Local: ! ! IA = 7^5 ! IB = 2^15 ! IB16 = 2^16 ! IP = 2^31-1 ! implicit none integer, parameter :: i4_huge = 2147483647 integer i4_uniform_0i integer, parameter :: ia = 16807 integer, parameter :: ib15 = 32768 integer, parameter :: ib16 = 65536 integer, parameter :: ip = 2147483647 integer iprhi integer ixhi integer k integer leftlo integer loxa integer seed ! ! Don't let SEED be 0. ! if ( seed == 0 ) then seed = i4_huge end if ! ! Get the 15 high order bits of SEED. ! ixhi = seed / ib16 ! ! Get the 16 low bits of SEED and form the low product. ! loxa = ( seed - ixhi * ib16 ) * ia ! ! Get the 15 high order bits of the low product. ! leftlo = loxa / ib16 ! ! Form the 31 highest bits of the full product. ! iprhi = ixhi * ia + leftlo ! ! Get overflow past the 31st bit of full product. ! k = iprhi / ib15 ! ! Assemble all the parts and presubtract IP. The parentheses are ! essential. ! seed = ( ( ( loxa - leftlo * ib16 ) - ip ) & + ( iprhi - k * ib15 ) * ib16 ) + k ! ! Add IP back in if necessary. ! if ( seed < 0 ) then seed = seed + ip end if i4_uniform_0i = seed return end function i4_uniform_ab ( a, b ) !*****************************************************************************80 ! !! i4_uniform_ab() returns a scaled pseudorandom I4. ! ! Discussion: ! ! An I4 is an integer value. ! ! The pseudorandom number will be scaled to be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 August 2021 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer A, B, the limits of the interval. ! ! Output: ! ! integer I4_UNIFORM_AB, a number between A and B. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer a integer b integer i4_uniform_ab real ( kind = rk ) r call random_number ( harvest = r ) i4_uniform_ab = a + int ( ( b + 1 - a ) * r ) return end subroutine i4mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_PRINT prints an I4MAT. ! ! Discussion: ! ! An I4MAT is a rectangular array of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, integer A(M,N), the matrix to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer m integer n integer a(m,n) integer ihi integer ilo integer jhi integer jlo character ( len = * ) title ilo = 1 ihi = m jlo = 1 jhi = n call i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) return end subroutine i4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_PRINT_SOME prints some of an I4MAT. ! ! Discussion: ! ! An I4MAT is a rectangular array of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 8 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i8)' ) j end do write ( *, '('' Col '',10a8)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(i8)' ) a(i,j) end do write ( *, '(i5,a,10a8)' ) i, ':', ( ctemp(j), j = 1, inc ) end do end do return end subroutine i4mat_uniform_ab ( m, n, a, b, seed, x ) !*****************************************************************************80 ! !! I4MAT_UNIFORM_AB returns a scaled pseudorandom I4MAT. ! ! Discussion: ! ! An I4MAT is a matrix of I4's. ! ! The pseudorandom numbers will be scaled to be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer M, N, the row and column dimensions ! of the matrix. ! ! Input, integer A, B, the limits of the interval. ! ! Input/output, integer SEED, the "seed" value, ! which should NOT be 0. ! On output, SEED has been updated. ! ! Output, integer X(M,N), a matrix of values between A and B. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer a integer b integer i integer, parameter :: i4_huge = 2147483647 integer j integer k real ( kind = rk ) r integer seed integer value integer x(m,n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r = real ( seed ) * 4.656612875E-10 ! ! Scale R to lie between A-0.5 and B+0.5. ! r = ( 1.0E+00 - r ) * ( real ( min ( a, b ) ) - 0.5E+00 ) & + r * ( real ( max ( a, b ) ) + 0.5E+00 ) ! ! Use rounding to convert R to an integer between A and B. ! value = nint ( r ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) x(i,j) = value end do end do return end subroutine i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 May 2010 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, integer A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer n integer a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,2x,i12)' ) i, ':', a(i) end do return end subroutine i4vec_uniform_ab ( n, a, b, seed, x ) !*****************************************************************************80 ! !! I4VEC_UNIFORM_AB returns a scaled pseudorandom I4VEC. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! The pseudorandom numbers should be scaled to be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer N, the dimension of the vector. ! ! Input, integer A, B, the limits of the interval. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, integer X(N), a vector of numbers between A and B. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer a integer b integer i integer, parameter :: i4_huge = 2147483647 integer k real ( kind = rk ) r integer seed integer value integer x(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4VEC_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r = real ( seed ) * 4.656612875E-10 ! ! Scale R to lie between A-0.5 and B+0.5. ! r = ( 1.0E+00 - r ) * ( real ( min ( a, b ) ) - 0.5E+00 ) & + r * ( real ( max ( a, b ) ) + 0.5E+00 ) ! ! Use rounding to convert R to an integer between A and B. ! value = nint ( r ) value = max ( value, min ( a, b ) ) value = min ( value, max ( a, b ) ) x(i) = value end do return end function l4_uniform ( seed ) !*****************************************************************************80 ! !! L4_UNIFORM returns a pseudorandom L4. ! ! Discussion: ! ! An L4 is a LOGICAL value. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 December 2007 ! ! 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/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, logical L4_UNIFORM, a pseudorandom logical value. ! implicit none integer, parameter :: i4_huge = 2147483647 integer, parameter :: i4_huge_half = 1073741823 integer k logical l4_uniform integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'L4_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if l4_uniform = ( i4_huge_half < seed ) return end subroutine l4mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! L4MAT_PRINT prints an L4MAT. ! ! Discussion: ! ! An L4MAT is an array of L4 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 November 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, logical A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer m integer n logical a(m,n) character ( len = * ) title call l4mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine l4mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! L4MAT_PRINT_SOME prints some of an L4MAT. ! ! Discussion: ! ! An L4MAT is an array of L4 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 November 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, logical A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: incx = 35 integer m integer n logical a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' if ( 100 <= j2hi ) then do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(1x,i1)' ) j / 100 end do write ( *, '('' '',35a2)' ) ctemp(1:inc) end if if ( 10 <= j2hi ) then do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(1x,i1)' ) mod ( j / 10, 10 ) end do write ( *, '('' '',35a2)' ) ctemp(1:inc) end if do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(1x,i1)' ) mod ( j, 10 ) end do write ( *, '('' Col '',35a2)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi write ( *, '(i5,a1,35(1x,l1))' ) i, ':', a(i,j2lo:j2hi) end do end do return end subroutine l4mat_uniform ( m, n, seed, lmat ) !*****************************************************************************80 ! !! L4MAT_UNIFORM returns a pseudorandom L4MAT. ! ! Discussion: ! ! An L4MAT is a two dimensional array of L4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 December 2007 ! ! 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, integer M, N, the order of the matrix. ! ! Input/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, logical LMAT(M,N), a pseudorandom logical matrix. ! implicit none integer m integer n integer, parameter :: i4_huge = 2147483647 integer, parameter :: i4_huge_half = 1073741823 integer i integer j integer k logical lmat(m,n) integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'L4MAT_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if lmat(i,j) = ( i4_huge_half < seed ) end do end do return end subroutine l4vec_print ( n, a, title ) !*****************************************************************************80 ! !! L4VEC_PRINT prints an L4VEC. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 May 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, logical A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer n logical a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,l1)' ) i, ':', a(i) end do return end subroutine l4vec_uniform ( n, seed, lvec ) !*****************************************************************************80 ! !! L4VEC_UNIFORM returns a pseudorandom L4VEC. ! ! Discussion: ! ! An L4VEC is a vector of L4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 December 2007 ! ! 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, integer N, the order of the vector. ! ! Input/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, logical LVEC(N), a pseudorandom logical vector. ! implicit none integer n integer, parameter :: i4_huge = 2147483647 integer, parameter :: i4_huge_half = 1073741823 integer i integer k logical lvec(n) integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'L4VEC_UNIFORM - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if lvec(i) = ( i4_huge_half < seed ) end do return end subroutine lcrg_anbn ( a, b, c, n, an, bn ) !*****************************************************************************80 ! !! LCRG_ANBN computes the "N-th power" of a linear congruential generator. ! ! Discussion: ! ! We are considering a linear congruential random number generator. ! The LCRG takes as input an integer value called SEED, and returns ! an updated value of SEED, ! ! SEED(out) = ( a * SEED(in) + b ) mod c. ! ! and an associated pseudorandom real value ! ! U = SEED(out) / c. ! ! In most cases, a user is content to call the LCRG repeatedly, with ! the updating of SEED being taken care of automatically. ! ! The purpose of this routine is to determine the values of AN and BN ! that describe the LCRG that is equivalent to N applications of the ! original LCRG. ! ! One use for such a facility would be to do random number computations ! in parallel. If each of N processors is to compute many random values, ! you can guarantee that they work with distinct random values ! by starting with a single value of SEED, using the original LCRG to ! generate the first N-1 "iterates" of SEED, so that you now have N "seed" ! values, and from now on, applying the N-th power of the LCRG to the seeds. ! ! If the K-th processor starts from the K-th seed, it will essentially ! be computing every N-th entry of the original random number sequence, ! offset by K. Thus the individual processors will be using a random ! number stream as good as the original one, and without repeating, and ! without having to communicate. ! ! To evaluate the N-th value of SEED directly, we start by ignoring ! the modular arithmetic, and working out the sequence of calculations ! as follows: ! ! SEED(0) = SEED. ! SEED(1) = a * SEED + b ! SEED(2) = a * SEED(1) + b = a^2 * SEED + a * b + b ! SEED(3) = a * SEED(2) + b = a^3 * SEED + a^2 * b + a * b + b ! ... ! SEED(N-1) = a * SEED(N-2) + b ! ! SEED(N) = a * SEED(N-1) + b = a^N * SEED ! + ( a^(n-1) + a^(n-2) + ... + a + 1 ) * b ! ! or, using the geometric series, ! ! SEED(N) = a^N * SEED + ( a^N - 1) / ( a - 1 ) * b ! = AN * SEED + BN ! ! Thus, from any SEED, we can determine the result of N applications of the ! original LCRG directly if we can solve ! ! ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic, ! ! and evaluate: ! ! AN = a^N ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 April 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Barry Wilkinson, Michael Allen, ! Parallel Programming: ! Techniques and Applications Using Networked Workstations ! and Parallel Computers, ! Prentice Hall, ! ISBN: 0-13-140563-2, ! LC: QA76.642.W54. ! ! Parameters: ! ! Input, integer A, the multiplier for the LCRG. ! ! Input, integer B, the added value for the LCRG. ! ! Input, integer C, the base for the modular arithmetic. ! For 32 bit arithmetic, this is often 2^31 - 1, or 2147483647. It is ! required that 0 < C. ! ! Input, integer N, the "index", or number of times that the ! LCRG is to be applied. It is required that 0 <= N. ! ! Output, integer AN, BN, the multiplier and added value for ! the LCRG that represent N applications of the original LCRG. ! implicit none integer a integer am1 integer an integer anm1tb integer b integer bn integer c integer ierror integer n if ( n < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_ANBN - Fatal error!' write ( *, '(a,i12)' ) ' Illegal input value of N = ', n stop 1 end if if ( c <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_ANBN - Fatal error!' write ( *, '(a,i12)' ) ' Illegal input value of C = ', c stop 1 end if if ( n == 0 ) then an = 1 bn = 0 else if ( n == 1 ) then an = a bn = b else ! ! Compute A^N. ! call power_mod ( a, n, c, an ) ! ! Solve ! ( a - 1 ) * BN = ( a^N - 1 ) mod B ! for BN. ! am1 = a - 1 anm1tb = ( an - 1 ) * b call congruence ( am1, c, anm1tb, ierror, bn ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_ANBN - Fatal error!' write ( *, '(a)' ) ' An error occurred in the CONGRUENCE routine.' write ( *, '(a,i8)' ) ' The error code was IERROR = ', ierror stop 1 end if end if return end subroutine lcrg_evaluate ( a, b, c, x, y ) !*****************************************************************************80 ! !! LCRG_EVALUATE evaluates an LCRG, y = ( A * x + B ) mod C. ! ! Discussion: ! ! This routine cannot be recommended for production use. Because we want ! to do modular arithmetic, but the base is not a power of 2, we need to ! use "double precision" integers to keep accuracy. ! ! If we knew the base C, we could try to avoid overflow while not changing ! precision. ! ! If the base C was a power of 2, we could rely on the usual properties of ! integer arithmetic on computers, in which overflow bits, which are always ! ignored, don't actually matter. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 April 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the multiplier for the LCRG. ! ! Input, integer B, the added value for the LCRG. ! ! Input, integer C, the base for the modular arithmetic. ! For 32 bit arithmetic, this is often 2^31 - 1, or 2147483647. It is ! required that 0 < C. ! ! Input, integer X, the value to be processed. ! ! Output, integer Y, the processed value. ! implicit none integer, parameter :: ik = selected_int_kind ( 16 ) integer a integer ( kind = ik ) a8 integer b integer ( kind = ik ) b8 integer c integer ( kind = ik ) c8 integer x integer ( kind = ik ) x8 integer y integer ( kind = ik ) y8 ! ! To avoid roundoff issues, we need to go to "double precision" integers. ! (Not available on all planets.) ! a8 = a b8 = b c8 = c x8 = x y8 = mod ( a8 * x8 + b8, c8 ) y = int ( y8 ) if ( y < 0 ) then y = y + c end if return end subroutine lcrg_seed ( a, b, c, n, seed, seed_lcrg ) !*****************************************************************************80 ! !! LCRG_SEED computes the N-th seed of a linear congruential generator. ! ! Discussion: ! ! We are considering a linear congruential random number generator. ! The LCRG takes as input an integer value called SEED, and returns ! an updated value of SEED, ! ! SEED(out) = ( a * SEED(in) + b ) mod c. ! ! and an associated pseudorandom real value ! ! U = SEED(out) / c. ! ! In most cases, a user is content to call the LCRG repeatedly, with ! the updating of SEED being taken care of automatically. ! ! The purpose of this routine is to determine the value of SEED that ! would be output after N successive applications of the LCRG. This ! allows the user to know, in advance, what the 1000-th value of ! SEED would be, for instance. Obviously, one way to do this is to ! apply the LCRG formula 1,000 times. However, it is possible to ! do this in a more direct and efficient way. ! ! One use for such a facility would be to do random number computations ! in parallel. If each processor is to compute 1,000 values, you can ! guarantee that they work with distinct random values by starting the ! first processor with SEED, the second with the value of SEED after ! 1,000 applications of the LCRG, and so on. ! ! To evaluate the N-th value of SEED directly, we start by ignoring ! the modular arithmetic, and working out the sequence of calculations ! as follows: ! ! SEED(0) = SEED. ! SEED(1) = a * SEED + b ! SEED(2) = a * SEED(1) + b = a^2 * SEED + a * b + b ! SEED(3) = a * SEED(2) + b = a^3 * SEED + a^2 * b + a * b + b ! ... ! SEED(N) = a * SEED(N-1) + b = a^N * SEED ! + ( a^(n-1) + a^(n-2) + ... + a + 1 ) * b ! ! or, using the geometric series, ! ! SEED(N) = a^N * SEED + ( a^N - 1) / ( a - 1 ) * b ! ! Therefore, we can determine SEED(N) directly if we can solve ! ! ( a - 1 ) * BN = ( a^N - 1 ) * b in modular arithmetic, ! ! and evaluated: ! ! AN = a^N ! ! Using the formula: ! ! SEED(N) = ( AN * SEED + BN ) mod c ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 April 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the multiplier for the LCRG. ! ! Input, integer B, the added value for the LCRG. ! ! Input, integer C, the base for the modular arithmetic. ! For 32 bit arithmetic, this is often 2^31 - 1, or 2147483647. It is ! required that 0 < C. ! ! Input, integer N, the "index", or number of times that the ! LCRG is to be applied. It is required that 0 <= N. ! ! Input, integer SEED, the starting value of SEED. It is ! customary that 0 < SEED. ! ! Output, integer SEED_LCRG, the value of SEED that would ! be output if the LCRG were applied to the starting value N times. ! implicit none integer, parameter :: ik = selected_int_kind ( 16 ) integer a integer am1 integer an integer anm1tb integer b integer bn integer c integer ierror integer n integer seed integer seed_lcrg integer ( kind = ik ) value2 if ( n < 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_SEED - Fatal error!' write ( *, '(a,i12)' ) ' Illegal input value of N = ', n stop 1 end if if ( c <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_SEED - Fatal error!' write ( *, '(a,i12)' ) ' Illegal input value of C = ', c stop 1 end if if ( n == 0 ) then seed_lcrg = mod ( seed, c ) if ( seed_lcrg < 0 ) then seed_lcrg = seed_lcrg + c end if return end if ! ! Get A^N. ! call power_mod ( a, n, c, an ) ! ! Solve ( a - 1 ) * BN = ( a^N - 1 ) for BN. ! ! The LCRG I have been investigating uses B = 0, so this code ! has not been properly tested yet. ! am1 = a - 1 anm1tb = ( an - 1 ) * b call congruence ( am1, c, anm1tb, ierror, bn ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LCRG_SEED - Fatal error!' write ( *, '(a)' ) ' An error occurred in the CONGRUENCE routine.' write ( *, '(a,i8)' ) ' The error code was IERROR = ', ierror stop 1 end if ! ! Set the new SEED. ! value2 = int ( an, kind = ik ) * int ( seed, kind = ik ) & + int ( bn, kind = ik ) value2 = mod ( value2, int ( c, kind = ik ) ) ! ! Guarantee that the value is positive. ! if ( value2 < 0 ) then value2 = value2 + int ( c, kind = ik ) end if seed_lcrg = int ( value2 ) return end subroutine power_mod ( a, n, m, x ) !*****************************************************************************80 ! !! POWER_MOD computes ( A^N ) mod M. ! ! Discussion: ! ! Some programming tricks are used to speed up the computation, and to ! allow computations in which the value A^N is much too large to ! store in an integer word. ! ! First, for efficiency, the power A^N is computed by determining ! the binary expansion of N, then computing A, A^2, A^4, and so on ! by repeated squaring, and multiplying only those factors that ! contribute to A^N. ! ! Secondly, the intermediate products are immediately "mod'ed", which ! keeps them small. ! ! For instance, to compute ( A^13 ) mod 11, we essentially compute ! ! 13 = 1 + 4 + 8 ! ! A^13 = A * A^4 * A^8 ! ! A^13 ( mod 11 ) = A ( mod 11 ) * A^4 ( mod 11 ) * A^8 ( mod 11 ). ! ! Fermat's little theorem says that if P is prime, and A is not divisible ! by P, then ( A^(P-1) - 1 ) is divisible by P. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 May 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the base of the expression to be tested. ! 0 <= A is required. ! ! Input, integer N, the power to which the base is raised. ! 0 <= N is required. ! ! Input, integer M, the divisor against which the expression ! is tested. 0 < M is required. ! ! Output, integer X, the remainder when A^N is divided by M. ! If any input quantity is unacceptable, then the nonsensical value ! X = -1 is returned. ! implicit none integer, parameter :: ik = selected_int_kind ( 16 ) integer a integer ( kind = ik ) a_square2 integer d integer m integer ( kind = ik ) m2 integer n integer ncopy integer x integer ( kind = ik ) x2 if ( a < 0 ) then x = -1 return end if if ( m <= 0 ) then x = -1 return end if if ( n < 0 ) then x = -1 return end if ! ! A_SQUARE contains the successive squares of A. ! a_square2 = int ( a, kind = ik ) x2 = int ( 1, kind = ik ) m2 = int ( m, kind = ik ) ncopy = n do while ( 0 < ncopy ) d = mod ( ncopy, 2 ) if ( d == 1 ) then x2 = mod ( x2 * a_square2, m2 ) end if a_square2 = mod ( a_square2 * a_square2, m2 ) ncopy = ( ncopy - d ) / 2 end do ! ! Fix up X so that it is nonnegative. ! do while ( x2 < 0 ) x2 = x2 + m2 end do x = int ( x2 ) return end function r8_uniform_01 ( seed ) !*****************************************************************************80 ! !! R8_UNIFORM_01 returns a unit pseudorandom R8. ! ! Discussion: ! ! An R8 is a real ( kind = rk ) value. ! ! For now, the input quantity SEED is an integer variable. ! ! 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 MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R8_UNIFORM_01, a new pseudorandom variate, ! strictly between 0 and 1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: i4_huge = 2147483647 integer k real ( kind = rk ) r8_uniform_01 integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r8_uniform_01 = real ( seed, kind = rk ) * 4.656612875D-10 return end function r8_uniform_ab ( a, b, seed ) !*****************************************************************************80 ! !! R8_UNIFORM_AB returns a scaled pseudorandom R8. ! ! Discussion: ! ! An R8 is a real ( kind = rk ) value. ! ! For now, the input quantity SEED is an integer variable. ! ! The pseudorandom number should be uniformly distributed ! between A and B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, real ( kind = rk ) A, B, the limits of the interval. ! ! Input/output, integer SEED, the "seed" value, which should ! NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R8_UNIFORM_AB, a number strictly between A and B. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a real ( kind = rk ) b integer, parameter :: i4_huge = 2147483647 integer k real ( kind = rk ) r8_uniform_ab integer seed if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r8_uniform_ab = a + ( b - a ) * real ( seed, kind = rk ) * 4.656612875D-10 return end subroutine r8col_uniform_ab ( m, n, a, b, seed, r ) !*****************************************************************************80 ! !! R8COL_UNIFORM_AB fills an R8COL with scaled pseudorandom numbers. ! ! Discussion: ! ! An R8COL is an array of R8 values, regarded as a set of column vectors. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 December 2014 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! 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. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns in ! the array. ! ! Input, real ( kind = rk ) A, B, the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a real ( kind = rk ) b integer i integer, parameter :: i4_huge = 2147483647 integer j integer k integer seed real ( kind = rk ) r(m,n) do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = a & + ( b - a ) * real ( seed, kind = rk ) * 4.656612875D-10 end do end do return end subroutine r8col_uniform_abvec ( m, n, a, b, seed, r ) !*****************************************************************************80 ! !! R8COL_UNIFORM_ABVEC fills an R8COL with scaled pseudorandom numbers. ! ! Discussion: ! ! An R8COL is an array of R8 values, regarded as a set of column vectors. ! ! The user specifies a minimum and maximum value for each row. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 December 2014 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! 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. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns in ! the array. ! ! Input, real ( kind = rk ) A(M), B(M), the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m) real ( kind = rk ) b(m) integer i integer, parameter :: i4_huge = 2147483647 integer j integer k integer seed real ( kind = rk ) r(m,n) do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = a(i) & + ( b(i) - a(i) ) * real ( seed, kind = rk ) * 4.656612875D-10 end do end do return end subroutine r8mat_print ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_PRINT prints an R8MAT. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows in A. ! ! Input, integer N, the number of columns in A. ! ! Input, real ( kind = rk ) A(M,N), the matrix. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8mat_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_PRINT_SOME prints some of an R8MAT. ! ! Discussion: ! ! An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M]. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) if ( m <= 0 .or. n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (None)' return end if do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i8,6x)' ) j end do write ( *, '('' Col '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' i2lo = max ( ilo, 1 ) i2hi = min ( ihi, m ) do i = i2lo, i2hi do j2 = 1, inc j = j2lo - 1 + j2 write ( ctemp(j2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc ) end do end do return end subroutine r8mat_uniform_01 ( m, n, seed, r ) !*****************************************************************************80 ! !! R8MAT_UNIFORM_01 returns a unit pseudorandom R8MAT. ! ! Discussion: ! ! An R8MAT is an array of R8's. ! ! For now, the input quantity SEED is an integer variable. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer M, N, the number of rows and columns ! in the array. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer i integer, parameter :: i4_huge = 2147483647 integer j integer k integer seed real ( kind = rk ) r(m,n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = real ( seed, kind = rk ) * 4.656612875D-10 end do end do return end subroutine r8mat_uniform_ab ( m, n, a, b, seed, r ) !*****************************************************************************80 ! !! R8MAT_UNIFORM_AB returns a scaled pseudorandom R8MAT. ! ! Discussion: ! ! A <= R(I,J) <= B. ! ! An R8MAT is an array of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer M, N, the number of rows and columns ! in the array. ! ! Input, real ( kind = rk ) A, B, the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a real ( kind = rk ) b integer i integer, parameter :: i4_huge = 2147483647 integer j integer k integer seed real ( kind = rk ) r(m,n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do j = 1, n do i = 1, m k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = a + ( b - a ) * real ( seed, kind = rk ) * 4.656612875D-10 end do end do return end subroutine r8row_uniform_abvec ( m, n, a, b, seed, r ) !*****************************************************************************80 ! !! R8ROW_UNIFORM_AB fills an R8ROW with scaled pseudorandom numbers. ! ! Discussion: ! ! An R8ROW is an array of R8 values, regarded as a set of row vectors. ! ! The user specifies a minimum and maximum value for each column. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 January 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Paul Bratley, Bennett Fox, Linus Schrage, ! A Guide to Simulation, ! Springer Verlag, pages 201-202, 1983. ! ! 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. ! ! Peter Lewis, Allen Goodman, James Miller, ! A Pseudo-Random Number Generator for the System/360, ! IBM Systems Journal, ! Volume 8, pages 136-143, 1969. ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns in ! the array. ! ! Input, real ( kind = rk ) A(N), B(N), the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(M,N), the array of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(n) real ( kind = rk ) b(n) integer i integer, parameter :: i4_huge = 2147483647 integer j integer k integer seed real ( kind = rk ) r(m,n) do i = 1, m do j = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i,j) = a(j) & + ( b(j) - a(j) ) * real ( seed, kind = rk ) * 4.656612875D-10 end do end do return end subroutine r8vec_normal_01 ( n, seed, x ) !*****************************************************************************80 ! !! R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! The standard normal probability distribution function (PDF) has ! mean 0 and standard deviation 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of values desired. ! ! Input/output, integer SEED, a seed for the random ! number generator. ! ! Output, real ( kind = rk ) X(N), a sample of the standard normal PDF. ! ! Local parameters: ! ! Local, real ( kind = rk ) R(N+1), is used to store some uniform ! random values. Its dimension is N+1, but really it is only needed ! to be the smallest even number greater than or equal to N. ! ! Local, integer X_LO_INDEX, X_HI_INDEX, records the range ! of entries of X that we need to compute ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer m real ( kind = rk ) r(n+1) real ( kind = rk ), parameter :: r8_pi = 3.141592653589793D+00 real ( kind = rk ) r8_uniform_01 integer seed real ( kind = rk ) x(n) integer x_hi_index integer x_lo_index ! ! Record the range of X we need to fill in. ! x_lo_index = 1 x_hi_index = n ! ! If we need just one new value, do that here to avoid null arrays. ! if ( x_hi_index - x_lo_index + 1 == 1 ) then r(1) = r8_uniform_01 ( seed ) if ( r(1) == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_NORMAL_01 - Fatal error!' write ( *, '(a)' ) ' R8_UNIFORM_01 returned a value of 0.' stop 1 end if r(2) = r8_uniform_01 ( seed ) x(x_hi_index) = & sqrt ( - 2.0D+00 * log ( r(1) ) ) * cos ( 2.0D+00 * r8_pi * r(2) ) ! ! If we require an even number of values, that's easy. ! else if ( mod ( x_hi_index - x_lo_index, 2 ) == 1 ) then m = ( x_hi_index - x_lo_index + 1 ) / 2 call r8vec_uniform_01 ( 2*m, seed, r ) x(x_lo_index:x_hi_index-1:2) = & sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & * cos ( 2.0D+00 * r8_pi * r(2:2*m:2) ) x(x_lo_index+1:x_hi_index:2) = & sqrt ( - 2.0D+00 * log ( r(1:2*m-1:2) ) ) & * sin ( 2.0D+00 * r8_pi * r(2:2*m:2) ) ! ! If we require an odd number of values, we generate an even number, ! and handle the last pair specially, storing one in X(N), and ! saving the other for later. ! else x_hi_index = x_hi_index - 1 m = ( x_hi_index - x_lo_index + 1 ) / 2 + 1 call r8vec_uniform_01 ( 2*m, seed, r ) x(x_lo_index:x_hi_index-1:2) = & sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & * cos ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) x(x_lo_index+1:x_hi_index:2) = & sqrt ( - 2.0D+00 * log ( r(1:2*m-3:2) ) ) & * sin ( 2.0D+00 * r8_pi * r(2:2*m-2:2) ) x(n) = sqrt ( - 2.0D+00 * log ( r(2*m-1) ) ) & * cos ( 2.0D+00 * r8_pi * r(2*m) ) end if return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i) end do return end subroutine r8vec_uniform_01 ( n, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer N, the number of entries in the vector. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(N), the vector of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer, parameter :: i4_huge = 2147483647 integer k integer seed real ( kind = rk ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_01 - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i) = real ( seed, kind = rk ) * 4.656612875D-10 end do return end subroutine r8vec_uniform_ab ( n, a, b, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_AB returns a scaled pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Each dimension ranges from A to B. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2007 ! ! 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, integer N, the number of entries in the vector. ! ! Input, real ( kind = rk ) A, B, the lower and upper limits. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(N), the vector of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a real ( kind = rk ) b integer i integer, parameter :: i4_huge = 2147483647 integer k integer seed real ( kind = rk ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i) = a + ( b - a ) * real ( seed, kind = rk ) * 4.656612875D-10 end do return end subroutine r8vec_uniform_abvec ( n, a, b, seed, r ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_ABVEC returns a scaled pseudorandom R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Dimension I ranges from A(I) to B(I). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 October 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, integer N, the number of entries in the vector. ! ! Input, real ( kind = rk ) A(N), B(N), the lower and upper limits ! for each dimension. ! ! Input/output, integer SEED, the "seed" value, which ! should NOT be 0. On output, SEED has been updated. ! ! Output, real ( kind = rk ) R(N), the vector of pseudorandom values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) real ( kind = rk ) b(n) integer i integer, parameter :: i4_huge = 2147483647 integer k integer seed real ( kind = rk ) r(n) if ( seed == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_ABVEC - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed < 0 ) then seed = seed + i4_huge end if r(i) = a(i) + ( b(i) - a(i) ) * real ( seed, kind = rk ) * 4.656612875D-10 end do return end subroutine r8vec_uniform_unit ( m, seed, w ) !*****************************************************************************80 ! !! R8VEC_UNIFORM_UNIT generates a uniformly random unit vector. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 October 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the spatial dimension. ! ! Input/output, integer SEED, a seed for the random number ! generator. ! ! Output, real ( kind = rk ) W(M), a random direction vector, ! with unit norm. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m real ( kind = rk ) norm integer seed real ( kind = rk ) w(m) ! ! Get M values from a standard normal distribution. ! call r8vec_normal_01 ( m, seed, w ) ! ! Compute the length of the vector. ! norm = sqrt ( sum ( w(1:m)**2 ) ) ! ! Normalize the vector. ! w(1:m) = w(1:m) / norm return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end