subroutine legendre_set ( n, x, w ) !*****************************************************************************80 ! !! LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature. ! ! Discussion: ! ! The integral: ! ! integral ( -1 <= x <= 1 ) f(x) dx ! ! The quadrature rule: ! ! sum ( 1 <= i <= n ) w(i) * f ( x(i) ) ! ! The quadrature rule is exact for polynomials through degree 2*N-1. ! ! The abscissas are the zeroes of the Legendre polynomial P(N)(X). ! ! Mathematica can compute the abscissas and weights of a Gauss-Legendre ! rule of order N for the interval [A,B] with P digits of precision ! by the commands: ! ! Needs["NumericalDifferentialEquationAnalysis`"] ! GaussianQuadratureWeights [ n, a, b, p ] ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2011 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Vladimir Krylov, ! Approximate Calculation of Integrals, ! Dover, 2006, ! ISBN: 0486445798, ! LC: QA311.K713. ! ! Arthur Stroud, Don Secrest, ! Gaussian Quadrature Formulas, ! Prentice Hall, 1966, ! LC: QA299.4G3S7. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Cambridge University Press, 1999, ! ISBN: 0-521-64314-7, ! LC: QA76.95.W65. ! ! Daniel Zwillinger, editor, ! CRC Standard Mathematical Tables and Formulae, ! 30th Edition, ! CRC Press, 1996, ! ISBN: 0-8493-2479-3, ! LC: QA47.M315. ! ! Parameters: ! ! Input, integer N, the order. ! N must be 1, 3, 7, 15, 31, 63, 127, or 255. ! ! Output, real ( kind = rk ) X(N), the abscissas. ! ! Output, real ( kind = rk ) W(N), the weights. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) w(n) real ( kind = rk ) x(n) if ( n == 1 ) then x(1) = 0.000000000000000000000000000000D+00 w(1) = 2.000000000000000000000000000000D+00 else if ( n == 3 ) then x(1) = -0.774596669241483377035853079956D+00 x(2) = 0.000000000000000000000000000000D+00 x(3) = 0.774596669241483377035853079956D+00 w(1) = 0.555555555555555555555555555556D+00 w(2) = 0.888888888888888888888888888889D+00 w(3) = 0.555555555555555555555555555556D+00 else if ( n == 7 ) then x(1) = -0.949107912342758524526189684048D+00 x(2) = -0.741531185599394439863864773281D+00 x(3) = -0.405845151377397166906606412077D+00 x(4) = 0.000000000000000000000000000000D+00 x(5) = 0.405845151377397166906606412077D+00 x(6) = 0.741531185599394439863864773281D+00 x(7) = 0.949107912342758524526189684048D+00 w(1) = 0.129484966168869693270611432679D+00 w(2) = 0.279705391489276667901467771424D+00 w(3) = 0.381830050505118944950369775489D+00 w(4) = 0.417959183673469387755102040816D+00 w(5) = 0.381830050505118944950369775489D+00 w(6) = 0.279705391489276667901467771424D+00 w(7) = 0.129484966168869693270611432679D+00 else if ( n == 15 ) then x(1) = -0.987992518020485428489565719D+00 x(2) = -0.937273392400705904307758948D+00 x(3) = -0.848206583410427216200648321D+00 x(4) = -0.724417731360170047416186055D+00 x(5) = -0.570972172608538847537226737D+00 x(6) = -0.394151347077563369897207371D+00 x(7) = -0.201194093997434522300628303D+00 x(8) = 0.00000000000000000000000000D+00 x(9) = 0.20119409399743452230062830D+00 x(10) = 0.39415134707756336989720737D+00 x(11) = 0.57097217260853884753722674D+00 x(12) = 0.72441773136017004741618605D+00 x(13) = 0.84820658341042721620064832D+00 x(14) = 0.93727339240070590430775895D+00 x(15) = 0.98799251802048542848956572D+00 w(1) = 0.030753241996117268354628393577D+00 w(2) = 0.070366047488108124709267416451D+00 w(3) = 0.107159220467171935011869546686D+00 w(4) = 0.13957067792615431444780479451D+00 w(5) = 0.16626920581699393355320086048D+00 w(6) = 0.18616100001556221102680056187D+00 w(7) = 0.19843148532711157645611832644D+00 w(8) = 0.20257824192556127288062019997D+00 w(9) = 0.19843148532711157645611832644D+00 w(10) = 0.18616100001556221102680056187D+00 w(11) = 0.16626920581699393355320086048D+00 w(12) = 0.13957067792615431444780479451D+00 w(13) = 0.107159220467171935011869546686D+00 w(14) = 0.070366047488108124709267416451D+00 w(15) = 0.030753241996117268354628393577D+00 else if ( n == 31 ) then x(1) = -0.99708748181947707405562655D+00 x(2) = -0.98468590966515248400246517D+00 x(3) = -0.96250392509294966178905240D+00 x(4) = -0.93075699789664816495694576D+00 x(5) = -0.88976002994827104337419201D+00 x(6) = -0.83992032014626734008690454D+00 x(7) = -0.78173314841662494040636002D+00 x(8) = -0.71577678458685328390597087D+00 x(9) = -0.64270672292426034618441820D+00 x(10) = -0.56324916140714926272094492D+00 x(11) = -0.47819378204490248044059404D+00 x(12) = -0.38838590160823294306135146D+00 x(13) = -0.29471806998170161661790390D+00 x(14) = -0.19812119933557062877241300D+00 x(15) = -0.09955531215234152032517479D+00 x(16) = 0.00000000000000000000000000D+00 x(17) = 0.09955531215234152032517479D+00 x(18) = 0.19812119933557062877241300D+00 x(19) = 0.29471806998170161661790390D+00 x(20) = 0.38838590160823294306135146D+00 x(21) = 0.47819378204490248044059404D+00 x(22) = 0.56324916140714926272094492D+00 x(23) = 0.64270672292426034618441820D+00 x(24) = 0.71577678458685328390597087D+00 x(25) = 0.78173314841662494040636002D+00 x(26) = 0.83992032014626734008690454D+00 x(27) = 0.88976002994827104337419201D+00 x(28) = 0.93075699789664816495694576D+00 x(29) = 0.96250392509294966178905240D+00 x(30) = 0.98468590966515248400246517D+00 x(31) = 0.99708748181947707405562655D+00 w(1) = 0.0074708315792487758586968750322D+00 w(2) = 0.017318620790310582463157996087D+00 w(3) = 0.027009019184979421800608708092D+00 w(4) = 0.036432273912385464024392010468D+00 w(5) = 0.045493707527201102902315857895D+00 w(6) = 0.054103082424916853711666259087D+00 w(7) = 0.062174786561028426910343543687D+00 w(8) = 0.069628583235410366167756126255D+00 w(9) = 0.076390386598776616426357674901D+00 w(10) = 0.082392991761589263903823367432D+00 w(11) = 0.087576740608477876126198069695D+00 w(12) = 0.091890113893641478215362871607D+00 w(13) = 0.095290242912319512807204197488D+00 w(14) = 0.097743335386328725093474010979D+00 w(15) = 0.099225011226672307874875514429D+00 w(16) = 0.09972054479342645142753383373D+00 w(17) = 0.099225011226672307874875514429D+00 w(18) = 0.097743335386328725093474010979D+00 w(19) = 0.095290242912319512807204197488D+00 w(20) = 0.091890113893641478215362871607D+00 w(21) = 0.087576740608477876126198069695D+00 w(22) = 0.082392991761589263903823367432D+00 w(23) = 0.076390386598776616426357674901D+00 w(24) = 0.069628583235410366167756126255D+00 w(25) = 0.062174786561028426910343543687D+00 w(26) = 0.054103082424916853711666259087D+00 w(27) = 0.045493707527201102902315857895D+00 w(28) = 0.036432273912385464024392010468D+00 w(29) = 0.027009019184979421800608708092D+00 w(30) = 0.017318620790310582463157996087D+00 w(31) = 0.0074708315792487758586968750322D+00 else if ( n == 63 ) then x(1) = -0.99928298402912378037893614D+00 x(2) = -0.99622401277797010860219336D+00 x(3) = -0.99072854689218946681089467D+00 x(4) = -0.98280881059372723486251141D+00 x(5) = -0.97248403469757002280196068D+00 x(6) = -0.95977944975894192707035417D+00 x(7) = -0.94472613404100980296637532D+00 x(8) = -0.92736092062184320544703138D+00 x(9) = -0.90772630277853155803695313D+00 x(10) = -0.88587032850785342629029846D+00 x(11) = -0.86184648236412371953961184D+00 x(12) = -0.83571355431950284347180777D+00 x(13) = -0.80753549577345676005146599D+00 x(14) = -0.7773812629903723355633302D+00 x(15) = -0.7453246483178474178293217D+00 x(16) = -0.7114440995848458078514315D+00 x(17) = -0.6758225281149860901311033D+00 x(18) = -0.6385471058213653850003070D+00 x(19) = -0.5997090518776252357390089D+00 x(20) = -0.5594034094862850132676978D+00 x(21) = -0.5177288132900332481244776D+00 x(22) = -0.4747872479948043999222123D+00 x(23) = -0.4306837987951116006620889D+00 x(24) = -0.3855263942122478924776150D+00 x(25) = -0.3394255419745844024688344D+00 x(26) = -0.2924940585862514400361572D+00 x(27) = -0.2448467932459533627484046D+00 x(28) = -0.1966003467915066845576275D+00 x(29) = -0.1478727863578719685698391D+00 x(30) = -0.0987833564469452795297037D+00 x(31) = -0.0494521871161596272342338D+00 x(32) = 0.0000000000000000000000000D+00 x(33) = 0.0494521871161596272342338D+00 x(34) = 0.0987833564469452795297037D+00 x(35) = 0.1478727863578719685698391D+00 x(36) = 0.1966003467915066845576275D+00 x(37) = 0.2448467932459533627484046D+00 x(38) = 0.2924940585862514400361572D+00 x(39) = 0.3394255419745844024688344D+00 x(40) = 0.3855263942122478924776150D+00 x(41) = 0.4306837987951116006620889D+00 x(42) = 0.4747872479948043999222123D+00 x(43) = 0.5177288132900332481244776D+00 x(44) = 0.5594034094862850132676978D+00 x(45) = 0.5997090518776252357390089D+00 x(46) = 0.6385471058213653850003070D+00 x(47) = 0.6758225281149860901311033D+00 x(48) = 0.7114440995848458078514315D+00 x(49) = 0.7453246483178474178293217D+00 x(50) = 0.7773812629903723355633302D+00 x(51) = 0.8075354957734567600514660D+00 x(52) = 0.8357135543195028434718078D+00 x(53) = 0.8618464823641237195396118D+00 x(54) = 0.8858703285078534262902985D+00 x(55) = 0.9077263027785315580369531D+00 x(56) = 0.9273609206218432054470314D+00 x(57) = 0.9447261340410098029663753D+00 x(58) = 0.9597794497589419270703542D+00 x(59) = 0.9724840346975700228019607D+00 x(60) = 0.9828088105937272348625114D+00 x(61) = 0.9907285468921894668108947D+00 x(62) = 0.9962240127779701086021934D+00 x(63) = 0.9992829840291237803789361D+00 w(1) = 0.0018398745955770841170924455540D+00 w(2) = 0.0042785083468637618660784110826D+00 w(3) = 0.0067102917659601362519069307298D+00 w(4) = 0.0091259686763266563540586454218D+00 w(5) = 0.011519376076880041750750606149D+00 w(6) = 0.013884612616115610824866086368D+00 w(7) = 0.016215878410338338882283672975D+00 w(8) = 0.018507464460161270409260545805D+00 w(9) = 0.020753761258039090775341953421D+00 w(10) = 0.022949271004889933148942319562D+00 w(11) = 0.025088620553344986618630138068D+00 w(12) = 0.027166574359097933225189839439D+00 w(13) = 0.029178047208280526945551502154D+00 w(14) = 0.031118116622219817508215988557D+00 w(15) = 0.032982034883779341765683179672D+00 w(16) = 0.034765240645355877697180504643D+00 w(17) = 0.036463370085457289630452409788D+00 w(18) = 0.038072267584349556763638324928D+00 w(19) = 0.039587995891544093984807928149D+00 w(20) = 0.041006845759666398635110037009D+00 w(21) = 0.042325345020815822982505485403D+00 w(22) = 0.043540267083027590798964315704D+00 w(23) = 0.044648638825941395370332669517D+00 w(24) = 0.045647747876292608685885992609D+00 w(25) = 0.046535149245383696510395418747D+00 w(26) = 0.047308671312268919080604988339D+00 w(27) = 0.047966421137995131411052756195D+00 w(28) = 0.048506789097883847864090099146D+00 w(29) = 0.048928452820511989944709361549D+00 w(30) = 0.049230380423747560785043116988D+00 w(31) = 0.049411833039918178967039646117D+00 w(32) = 0.04947236662393102088866936042D+00 w(33) = 0.049411833039918178967039646117D+00 w(34) = 0.049230380423747560785043116988D+00 w(35) = 0.048928452820511989944709361549D+00 w(36) = 0.048506789097883847864090099146D+00 w(37) = 0.047966421137995131411052756195D+00 w(38) = 0.047308671312268919080604988339D+00 w(39) = 0.046535149245383696510395418747D+00 w(40) = 0.045647747876292608685885992609D+00 w(41) = 0.044648638825941395370332669517D+00 w(42) = 0.043540267083027590798964315704D+00 w(43) = 0.042325345020815822982505485403D+00 w(44) = 0.041006845759666398635110037009D+00 w(45) = 0.039587995891544093984807928149D+00 w(46) = 0.038072267584349556763638324928D+00 w(47) = 0.036463370085457289630452409788D+00 w(48) = 0.034765240645355877697180504643D+00 w(49) = 0.032982034883779341765683179672D+00 w(50) = 0.031118116622219817508215988557D+00 w(51) = 0.029178047208280526945551502154D+00 w(52) = 0.027166574359097933225189839439D+00 w(53) = 0.025088620553344986618630138068D+00 w(54) = 0.022949271004889933148942319562D+00 w(55) = 0.020753761258039090775341953421D+00 w(56) = 0.018507464460161270409260545805D+00 w(57) = 0.016215878410338338882283672975D+00 w(58) = 0.013884612616115610824866086368D+00 w(59) = 0.011519376076880041750750606149D+00 w(60) = 0.0091259686763266563540586454218D+00 w(61) = 0.0067102917659601362519069307298D+00 w(62) = 0.0042785083468637618660784110826D+00 w(63) = 0.0018398745955770841170924455540D+00 else if ( n == 127 ) then x(1) = -0.9998221304153061462673512D+00 x(2) = -0.9990629343553118951383159D+00 x(3) = -0.9976975661898046210744170D+00 x(4) = -0.9957265513520272266354334D+00 x(5) = -0.9931510492545171473611308D+00 x(6) = -0.9899726145914841576077867D+00 x(7) = -0.9861931740169316667104383D+00 x(8) = -0.9818150208038141100334631D+00 x(9) = -0.9768408123430703268174439D+00 x(10) = -0.9712735681615291922889469D+00 x(11) = -0.9651166679452921210908251D+00 x(12) = -0.9583738494252387711491029D+00 x(13) = -0.9510492060778803105479076D+00 x(14) = -0.9431471846248148273454496D+00 x(15) = -0.9346725823247379685736349D+00 x(16) = -0.9256305440562338491274647D+00 x(17) = -0.9160265591914658093130886D+00 x(18) = -0.9058664582618213828024613D+00 x(19) = -0.8951564094170837089690438D+00 x(20) = -0.8839029146800265699452579D+00 x(21) = -0.8721128059985607114196375D+00 x(22) = -0.8597932410977408098120313D+00 x(23) = -0.8469516991340975984533393D+00 x(24) = -0.8335959761548995143795572D+00 x(25) = -0.8197341803650786741551191D+00 x(26) = -0.8053747272046802146665608D+00 x(27) = -0.7905263342398137999454500D+00 x(28) = -0.7751980158702023824449628D+00 x(29) = -0.7593990778565366715566637D+00 x(30) = -0.7431391116709545129205669D+00 x(31) = -0.7264279886740726855356929D+00 x(32) = -0.7092758541221045609994446D+00 x(33) = -0.6916931210077006701564414D+00 x(34) = -0.6736904637382504853466825D+00 x(35) = -0.6552788116554826302767651D+00 x(36) = -0.6364693424002972413476082D+00 x(37) = -0.6172734751268582838576392D+00 x(38) = -0.5977028635700652293844120D+00 x(39) = -0.5777693889706125800032517D+00 x(40) = -0.5574851528619322329218619D+00 x(41) = -0.5368624697233975674581664D+00 x(42) = -0.5159138595042493572772773D+00 x(43) = -0.4946520400227821173949402D+00 x(44) = -0.4730899192454052416450999D+00 x(45) = -0.4512405874502662273318986D+00 x(46) = -0.4291173092801933762625441D+00 x(47) = -0.4067335156897825634086729D+00 x(48) = -0.3841027957915169357790778D+00 x(49) = -0.3612388886058697060709248D+00 x(50) = -0.3381556747203985013760003D+00 x(51) = -0.3148671678628949814860148D+00 x(52) = -0.2913875063937056207945188D+00 x(53) = -0.2677309447223886208883435D+00 x(54) = -0.2439118446539178579707132D+00 x(55) = -0.2199446666696875424545234D+00 x(56) = -0.1958439611486108515042816D+00 x(57) = -0.1716243595336421650083449D+00 x(58) = -0.1473005654490856693893293D+00 x(59) = -0.1228873457740829717260337D+00 x(60) = -0.0983995216776989707510918D+00 x(61) = -0.0738519596210485452734404D+00 x(62) = -0.0492595623319266303153793D+00 x(63) = -0.0246372597574209446148971D+00 x(64) = 0.0000000000000000000000000D+00 x(65) = 0.0246372597574209446148971D+00 x(66) = 0.0492595623319266303153793D+00 x(67) = 0.0738519596210485452734404D+00 x(68) = 0.0983995216776989707510918D+00 x(69) = 0.1228873457740829717260337D+00 x(70) = 0.1473005654490856693893293D+00 x(71) = 0.1716243595336421650083449D+00 x(72) = 0.1958439611486108515042816D+00 x(73) = 0.2199446666696875424545234D+00 x(74) = 0.2439118446539178579707132D+00 x(75) = 0.2677309447223886208883435D+00 x(76) = 0.2913875063937056207945188D+00 x(77) = 0.3148671678628949814860148D+00 x(78) = 0.3381556747203985013760003D+00 x(79) = 0.3612388886058697060709248D+00 x(80) = 0.3841027957915169357790778D+00 x(81) = 0.4067335156897825634086729D+00 x(82) = 0.4291173092801933762625441D+00 x(83) = 0.4512405874502662273318986D+00 x(84) = 0.4730899192454052416450999D+00 x(85) = 0.4946520400227821173949402D+00 x(86) = 0.5159138595042493572772773D+00 x(87) = 0.5368624697233975674581664D+00 x(88) = 0.5574851528619322329218619D+00 x(89) = 0.5777693889706125800032517D+00 x(90) = 0.5977028635700652293844120D+00 x(91) = 0.6172734751268582838576392D+00 x(92) = 0.6364693424002972413476082D+00 x(93) = 0.6552788116554826302767651D+00 x(94) = 0.6736904637382504853466825D+00 x(95) = 0.6916931210077006701564414D+00 x(96) = 0.7092758541221045609994446D+00 x(97) = 0.7264279886740726855356929D+00 x(98) = 0.7431391116709545129205669D+00 x(99) = 0.7593990778565366715566637D+00 x(100) = 0.7751980158702023824449628D+00 x(101) = 0.7905263342398137999454500D+00 x(102) = 0.8053747272046802146665608D+00 x(103) = 0.8197341803650786741551191D+00 x(104) = 0.8335959761548995143795572D+00 x(105) = 0.8469516991340975984533393D+00 x(106) = 0.8597932410977408098120313D+00 x(107) = 0.8721128059985607114196375D+00 x(108) = 0.8839029146800265699452579D+00 x(109) = 0.8951564094170837089690438D+00 x(110) = 0.9058664582618213828024613D+00 x(111) = 0.9160265591914658093130886D+00 x(112) = 0.9256305440562338491274647D+00 x(113) = 0.9346725823247379685736349D+00 x(114) = 0.9431471846248148273454496D+00 x(115) = 0.9510492060778803105479076D+00 x(116) = 0.9583738494252387711491029D+00 x(117) = 0.965116667945292121090825D+00 x(118) = 0.971273568161529192288947D+00 x(119) = 0.976840812343070326817444D+00 x(120) = 0.981815020803814110033463D+00 x(121) = 0.986193174016931666710438D+00 x(122) = 0.989972614591484157607787D+00 x(123) = 0.993151049254517147361131D+00 x(124) = 0.995726551352027226635433D+00 x(125) = 0.997697566189804621074417D+00 x(126) = 0.999062934355311895138316D+00 x(127) = 0.999822130415306146267351D+00 w(1) = 0.00045645726109586662791936519265D+00 w(2) = 0.00106227668695384869596523598532D+00 w(3) = 0.0016683488125171936761028862915D+00 w(4) = 0.0022734860707492547802810840776D+00 w(5) = 0.0028772587656289004082883197514D+00 w(6) = 0.0034792893810051465908910894100D+00 w(7) = 0.0040792095178254605327114733457D+00 w(8) = 0.0046766539777779034772638165663D+00 w(9) = 0.0052712596565634400891303815906D+00 w(10) = 0.0058626653903523901033648343751D+00 w(11) = 0.0064505120486899171845442463869D+00 w(12) = 0.0070344427036681608755685893033D+00 w(13) = 0.0076141028256526859356393930849D+00 w(14) = 0.0081891404887415730817235884719D+00 w(15) = 0.0087592065795403145773316804234D+00 w(16) = 0.0093239550065309714787536985834D+00 w(17) = 0.0098830429087554914716648010900D+00 w(18) = 0.0104361308631410052256731719977D+00 w(19) = 0.0109828830900689757887996573761D+00 w(20) = 0.011522967656921087154811609735D+00 w(21) = 0.012056056679400848183529562145D+00 w(22) = 0.012581826520465013101514365424D+00 w(23) = 0.013099957986718627426172681913D+00 w(24) = 0.013610136522139249906034237534D+00 w(25) = 0.014112052399003395774044161634D+00 w(26) = 0.014605400905893418351737288079D+00 w(27) = 0.015089882532666922992635733981D+00 w(28) = 0.015565203152273955098532590263D+00 w(29) = 0.016031074199309941802254151843D+00 w(30) = 0.016487212845194879399346060358D+00 w(31) = 0.016933342169871654545878815295D+00 w(32) = 0.017369191329918731922164721250D+00 w(33) = 0.017794495722974774231027912900D+00 w(34) = 0.018208997148375106468721469154D+00 w(35) = 0.018612443963902310429440419899D+00 w(36) = 0.019004591238555646611148901045D+00 w(37) = 0.019385200901246454628112623489D+00 w(38) = 0.019754041885329183081815217323D+00 w(39) = 0.020110890268880247225644623956D+00 w(40) = 0.020455529410639508279497065713D+00 w(41) = 0.020787750081531811812652137291D+00 w(42) = 0.021107350591688713643523847922D+00 w(43) = 0.021414136912893259295449693234D+00 w(44) = 0.021707922796373466052301324695D+00 w(45) = 0.021988529885872983756478409759D+00 w(46) = 0.022255787825930280235631416460D+00 w(47) = 0.022509534365300608085694429903D+00 w(48) = 0.022749615455457959852242553241D+00 w(49) = 0.022975885344117206754377437839D+00 w(50) = 0.023188206663719640249922582982D+00 w(51) = 0.023386450514828194170722043497D+00 w(52) = 0.023570496544381716050033676844D+00 w(53) = 0.023740233018760777777714726703D+00 w(54) = 0.023895556891620665983864481754D+00 w(55) = 0.024036373866450369675132086026D+00 w(56) = 0.024162598453819584716522917711D+00 w(57) = 0.024274154023278979833195063937D+00 w(58) = 0.024370972849882214952813561907D+00 w(59) = 0.024452996155301467956140198472D+00 w(60) = 0.024520174143511508275183033290D+00 w(61) = 0.024572466031020653286354137335D+00 w(62) = 0.024609840071630254092545634003D+00 w(63) = 0.024632273575707679066033370218D+00 w(64) = 0.02463975292396109441957941748D+00 w(65) = 0.024632273575707679066033370218D+00 w(66) = 0.024609840071630254092545634003D+00 w(67) = 0.024572466031020653286354137335D+00 w(68) = 0.024520174143511508275183033290D+00 w(69) = 0.024452996155301467956140198472D+00 w(70) = 0.024370972849882214952813561907D+00 w(71) = 0.024274154023278979833195063937D+00 w(72) = 0.024162598453819584716522917711D+00 w(73) = 0.024036373866450369675132086026D+00 w(74) = 0.023895556891620665983864481754D+00 w(75) = 0.023740233018760777777714726703D+00 w(76) = 0.023570496544381716050033676844D+00 w(77) = 0.023386450514828194170722043497D+00 w(78) = 0.023188206663719640249922582982D+00 w(79) = 0.022975885344117206754377437839D+00 w(80) = 0.022749615455457959852242553241D+00 w(81) = 0.022509534365300608085694429903D+00 w(82) = 0.022255787825930280235631416460D+00 w(83) = 0.021988529885872983756478409759D+00 w(84) = 0.021707922796373466052301324695D+00 w(85) = 0.021414136912893259295449693234D+00 w(86) = 0.021107350591688713643523847922D+00 w(87) = 0.020787750081531811812652137291D+00 w(88) = 0.020455529410639508279497065713D+00 w(89) = 0.020110890268880247225644623956D+00 w(90) = 0.019754041885329183081815217323D+00 w(91) = 0.019385200901246454628112623489D+00 w(92) = 0.019004591238555646611148901045D+00 w(93) = 0.018612443963902310429440419899D+00 w(94) = 0.018208997148375106468721469154D+00 w(95) = 0.017794495722974774231027912900D+00 w(96) = 0.017369191329918731922164721250D+00 w(97) = 0.016933342169871654545878815295D+00 w(98) = 0.016487212845194879399346060358D+00 w(99) = 0.016031074199309941802254151843D+00 w(100) = 0.015565203152273955098532590263D+00 w(101) = 0.015089882532666922992635733981D+00 w(102) = 0.014605400905893418351737288079D+00 w(103) = 0.014112052399003395774044161634D+00 w(104) = 0.013610136522139249906034237534D+00 w(105) = 0.013099957986718627426172681913D+00 w(106) = 0.012581826520465013101514365424D+00 w(107) = 0.012056056679400848183529562145D+00 w(108) = 0.011522967656921087154811609735D+00 w(109) = 0.0109828830900689757887996573761D+00 w(110) = 0.0104361308631410052256731719977D+00 w(111) = 0.0098830429087554914716648010900D+00 w(112) = 0.0093239550065309714787536985834D+00 w(113) = 0.0087592065795403145773316804234D+00 w(114) = 0.0081891404887415730817235884719D+00 w(115) = 0.0076141028256526859356393930849D+00 w(116) = 0.0070344427036681608755685893033D+00 w(117) = 0.0064505120486899171845442463869D+00 w(118) = 0.0058626653903523901033648343751D+00 w(119) = 0.0052712596565634400891303815906D+00 w(120) = 0.0046766539777779034772638165663D+00 w(121) = 0.0040792095178254605327114733457D+00 w(122) = 0.0034792893810051465908910894100D+00 w(123) = 0.0028772587656289004082883197514D+00 w(124) = 0.0022734860707492547802810840776D+00 w(125) = 0.0016683488125171936761028862915D+00 w(126) = 0.00106227668695384869596523598532D+00 w(127) = 0.00045645726109586662791936519265D+00 else if ( n == 255 ) then x(1) = -0.999955705317563751730191D+00 x(2) = -0.999766621312000569367063D+00 x(3) = -0.999426474680169959344386D+00 x(4) = -0.998935241284654635142155D+00 x(5) = -0.998292986136967889228248D+00 x(6) = -0.997499804126615814044844D+00 x(7) = -0.996555814435198617028738D+00 x(8) = -0.995461159480026294089975D+00 x(9) = -0.994216004616630164799381D+00 x(10) = -0.992820538021989138984811D+00 x(11) = -0.991274970630385567164523D+00 x(12) = -0.989579536085920123498574D+00 x(13) = -0.987734490699732356281248D+00 x(14) = -0.985740113407419277752900D+00 x(15) = -0.983596705724776358640192D+00 x(16) = -0.981304591701017185126565D+00 x(17) = -0.978864117869068155239121D+00 x(18) = -0.976275653192735980815246D+00 x(19) = -0.973539589010643617645393D+00 x(20) = -0.970656338976880365477697D+00 x(21) = -0.967626338998338798105523D+00 x(22) = -0.964450047168726298761719D+00 x(23) = -0.961127943699247839572910D+00 x(24) = -0.957660530845962076295490D+00 x(25) = -0.954048332833816317950921D+00 x(26) = -0.950291895777368285733522D+00 x(27) = -0.946391787598204251752103D+00 x(28) = -0.942348597939064408301480D+00 x(29) = -0.938162938074687317626793D+00 x(30) = -0.933835440819386124349338D+00 x(31) = -0.929366760431369935739045D+00 x(32) = -0.924757572513824425220425D+00 x(33) = -0.920008573912766315142721D+00 x(34) = -0.915120482611686961035103D+00 x(35) = -0.910094037623000801254172D+00 x(36) = -0.904929998876314959753358D+00 x(37) = -0.899629147103536800144342D+00 x(38) = -0.894192283720836729335637D+00 x(39) = -0.888620230707484040924981D+00 x(40) = -0.882913830481574073645470D+00 x(41) = -0.877073945772665439532627D+00 x(42) = -0.871101459491346550796200D+00 x(43) = -0.864997274595751144137121D+00 x(44) = -0.858762313955042966785823D+00 x(45) = -0.852397520209890250084237D+00 x(46) = -0.845903855629951054143931D+00 x(47) = -0.839282301968391021084600D+00 x(48) = -0.832533860313455524647230D+00 x(49) = -0.825659550937118650611534D+00 x(50) = -0.818660413140831885432406D+00 x(51) = -0.811537505098395829833580D+00 x(52) = -0.804291903695978689734633D+00 x(53) = -0.796924704369305728807154D+00 x(54) = -0.789437020938044295117764D+00 x(55) = -0.781829985437409458675147D+00 x(56) = -0.774104747947015717207115D+00 x(57) = -0.766262476417000644100858D+00 x(58) = -0.758304356491446765092016D+00 x(59) = -0.750231591329128358931528D+00 x(60) = -0.742045401421610281838045D+00 x(61) = -0.733747024408726316001889D+00 x(62) = -0.725337714891464938687812D+00 x(63) = -0.716818744242290800531501D+00 x(64) = -0.708191400412930589382399D+00 x(65) = -0.699456987739652339456557D+00 x(66) = -0.690616826746067624571761D+00 x(67) = -0.681672253943486448787259D+00 x(68) = -0.672624621628855017806731D+00 x(69) = -0.663475297680306939970658D+00 x(70) = -0.654225665350358766508700D+00 x(71) = -0.644877123056781136890077D+00 x(72) = -0.635431084171177146547142D+00 x(73) = -0.625888976805299900901619D+00 x(74) = -0.616252243595141561442344D+00 x(75) = -0.606522341482826526536576D+00 x(76) = -0.596700741496341721653202D+00 x(77) = -0.586788928527137300685706D+00 x(78) = -0.576788401105631382036211D+00 x(79) = -0.566700671174652760010815D+00 x(80) = -0.556527263860855843833077D+00 x(81) = -0.546269717244142383159817D+00 x(82) = -0.535929582125124840335150D+00 x(83) = -0.525508421790666565699453D+00 x(84) = -0.515007811777534223035005D+00 x(85) = -0.504429339634198197635551D+00 x(86) = -0.493774604680816999489812D+00 x(87) = -0.483045217767441948626854D+00 x(88) = -0.472242801030478698742627D+00 x(89) = -0.461368987647442418771401D+00 x(90) = -0.450425421590043710043279D+00 x(91) = -0.439413757375642589040685D+00 x(92) = -0.428335659817108112494341D+00 x(93) = -0.417192803771121462605751D+00 x(94) = -0.405986873884960545511889D+00 x(95) = -0.394719564341804385683361D+00 x(96) = -0.383392578604595822734854D+00 x(97) = -0.372007629158501235092510D+00 x(98) = -0.360566437252006227074021D+00 x(99) = -0.349070732636686422161576D+00 x(100) = -0.337522253305692705554261D+00 x(101) = -0.325922745230990453444769D+00 x(102) = -0.314273962099392474845918D+00 x(103) = -0.302577665047425574167140D+00 x(104) = -0.290835622395070819082047D+00 x(105) = -0.279049609378417768508970D+00 x(106) = -0.267221407881273079721012D+00 x(107) = -0.255352806165764071686080D+00 x(108) = -0.243445598601977973686482D+00 x(109) = -0.231501585396677734059116D+00 x(110) = -0.219522572321135403508985D+00 x(111) = -0.207510370438124240859625D+00 x(112) = -0.195466795828110816293869D+00 x(113) = -0.183393669314688508087976D+00 x(114) = -0.171292816189293903533225D+00 x(115) = -0.159166065935247723154292D+00 x(116) = -0.147015251951161989456661D+00 x(117) = -0.134842211273755257250625D+00 x(118) = -0.122648784300117812092492D+00 x(119) = -0.110436814509468826540991D+00 x(120) = -0.098208148184447540736015D+00 x(121) = -0.085964634131980604256000D+00 x(122) = -0.073708123403767780288977D+00 x(123) = -0.061440469016428270850728D+00 x(124) = -0.049163525671349973093019D+00 x(125) = -0.036879149474284021657652D+00 x(126) = -0.024589197654727010541405D+00 x(127) = -0.012295528285133320036860D+00 x(128) = 0.000000000000000000000000D+00 x(129) = 0.012295528285133320036860D+00 x(130) = 0.024589197654727010541405D+00 x(131) = 0.036879149474284021657652D+00 x(132) = 0.049163525671349973093019D+00 x(133) = 0.061440469016428270850728D+00 x(134) = 0.073708123403767780288977D+00 x(135) = 0.085964634131980604256000D+00 x(136) = 0.098208148184447540736015D+00 x(137) = 0.110436814509468826540991D+00 x(138) = 0.122648784300117812092492D+00 x(139) = 0.134842211273755257250625D+00 x(140) = 0.147015251951161989456661D+00 x(141) = 0.159166065935247723154292D+00 x(142) = 0.171292816189293903533225D+00 x(143) = 0.183393669314688508087976D+00 x(144) = 0.195466795828110816293869D+00 x(145) = 0.207510370438124240859625D+00 x(146) = 0.219522572321135403508985D+00 x(147) = 0.231501585396677734059116D+00 x(148) = 0.243445598601977973686482D+00 x(149) = 0.255352806165764071686080D+00 x(150) = 0.267221407881273079721012D+00 x(151) = 0.279049609378417768508970D+00 x(152) = 0.290835622395070819082047D+00 x(153) = 0.302577665047425574167140D+00 x(154) = 0.314273962099392474845918D+00 x(155) = 0.325922745230990453444769D+00 x(156) = 0.337522253305692705554261D+00 x(157) = 0.349070732636686422161576D+00 x(158) = 0.360566437252006227074021D+00 x(159) = 0.372007629158501235092510D+00 x(160) = 0.383392578604595822734854D+00 x(161) = 0.394719564341804385683361D+00 x(162) = 0.405986873884960545511889D+00 x(163) = 0.417192803771121462605751D+00 x(164) = 0.428335659817108112494341D+00 x(165) = 0.439413757375642589040685D+00 x(166) = 0.450425421590043710043279D+00 x(167) = 0.461368987647442418771401D+00 x(168) = 0.472242801030478698742627D+00 x(169) = 0.483045217767441948626854D+00 x(170) = 0.493774604680816999489812D+00 x(171) = 0.504429339634198197635551D+00 x(172) = 0.515007811777534223035005D+00 x(173) = 0.525508421790666565699453D+00 x(174) = 0.535929582125124840335150D+00 x(175) = 0.546269717244142383159817D+00 x(176) = 0.556527263860855843833077D+00 x(177) = 0.566700671174652760010815D+00 x(178) = 0.576788401105631382036211D+00 x(179) = 0.586788928527137300685706D+00 x(180) = 0.596700741496341721653202D+00 x(181) = 0.606522341482826526536576D+00 x(182) = 0.616252243595141561442344D+00 x(183) = 0.625888976805299900901619D+00 x(184) = 0.635431084171177146547142D+00 x(185) = 0.644877123056781136890077D+00 x(186) = 0.654225665350358766508700D+00 x(187) = 0.663475297680306939970658D+00 x(188) = 0.672624621628855017806731D+00 x(189) = 0.681672253943486448787259D+00 x(190) = 0.690616826746067624571761D+00 x(191) = 0.699456987739652339456557D+00 x(192) = 0.708191400412930589382399D+00 x(193) = 0.716818744242290800531501D+00 x(194) = 0.725337714891464938687812D+00 x(195) = 0.733747024408726316001889D+00 x(196) = 0.742045401421610281838045D+00 x(197) = 0.750231591329128358931528D+00 x(198) = 0.758304356491446765092016D+00 x(199) = 0.766262476417000644100858D+00 x(200) = 0.774104747947015717207115D+00 x(201) = 0.781829985437409458675147D+00 x(202) = 0.789437020938044295117764D+00 x(203) = 0.796924704369305728807154D+00 x(204) = 0.804291903695978689734633D+00 x(205) = 0.811537505098395829833580D+00 x(206) = 0.818660413140831885432406D+00 x(207) = 0.825659550937118650611534D+00 x(208) = 0.832533860313455524647230D+00 x(209) = 0.839282301968391021084600D+00 x(210) = 0.845903855629951054143931D+00 x(211) = 0.852397520209890250084237D+00 x(212) = 0.858762313955042966785823D+00 x(213) = 0.864997274595751144137121D+00 x(214) = 0.871101459491346550796200D+00 x(215) = 0.877073945772665439532627D+00 x(216) = 0.882913830481574073645470D+00 x(217) = 0.888620230707484040924981D+00 x(218) = 0.894192283720836729335637D+00 x(219) = 0.899629147103536800144342D+00 x(220) = 0.904929998876314959753358D+00 x(221) = 0.910094037623000801254172D+00 x(222) = 0.915120482611686961035103D+00 x(223) = 0.920008573912766315142721D+00 x(224) = 0.924757572513824425220425D+00 x(225) = 0.929366760431369935739045D+00 x(226) = 0.933835440819386124349338D+00 x(227) = 0.938162938074687317626793D+00 x(228) = 0.942348597939064408301480D+00 x(229) = 0.946391787598204251752103D+00 x(230) = 0.950291895777368285733522D+00 x(231) = 0.954048332833816317950921D+00 x(232) = 0.957660530845962076295490D+00 x(233) = 0.961127943699247839572910D+00 x(234) = 0.964450047168726298761719D+00 x(235) = 0.967626338998338798105523D+00 x(236) = 0.970656338976880365477697D+00 x(237) = 0.973539589010643617645393D+00 x(238) = 0.976275653192735980815246D+00 x(239) = 0.978864117869068155239121D+00 x(240) = 0.981304591701017185126565D+00 x(241) = 0.983596705724776358640192D+00 x(242) = 0.985740113407419277752900D+00 x(243) = 0.987734490699732356281248D+00 x(244) = 0.989579536085920123498574D+00 x(245) = 0.991274970630385567164523D+00 x(246) = 0.992820538021989138984811D+00 x(247) = 0.994216004616630164799381D+00 x(248) = 0.995461159480026294089975D+00 x(249) = 0.996555814435198617028738D+00 x(250) = 0.997499804126615814044844D+00 x(251) = 0.998292986136967889228248D+00 x(252) = 0.998935241284654635142155D+00 x(253) = 0.999426474680169959344386D+00 x(254) = 0.999766621312000569367063D+00 x(255) = 0.999955705317563751730191D+00 w(1) = 0.00011367361999142272115645954414D+00 w(2) = 0.00026459387119083065532790838855D+00 w(3) = 0.00041569762526823913616284210066D+00 w(4) = 0.00056675794564824918946626058353D+00 w(5) = 0.00071773647800611087798371518325D+00 w(6) = 0.00086860766611945667949717690640D+00 w(7) = 0.00101934797642732530281229369360D+00 w(8) = 0.0011699343729388079886897709773D+00 w(9) = 0.0013203439900221692090523602144D+00 w(10) = 0.0014705540427783843160097204304D+00 w(11) = 0.0016205417990415653896921100325D+00 w(12) = 0.0017702845706603213070421243905D+00 w(13) = 0.0019197597117132050055085980675D+00 w(14) = 0.0020689446195015801533643667413D+00 w(15) = 0.0022178167367540171700373764020D+00 w(16) = 0.0023663535543962867157201855305D+00 w(17) = 0.0025145326145997073931298921370D+00 w(18) = 0.0026623315139717112732749157331D+00 w(19) = 0.0028097279068204407457332299361D+00 w(20) = 0.0029566995084575002760043344138D+00 w(21) = 0.0031032240985191112621977893133D+00 w(22) = 0.0032492795242943133198690930777D+00 w(23) = 0.0033948437040533928255056951665D+00 w(24) = 0.0035398946303722552150296713510D+00 w(25) = 0.0036844103734499176530742235517D+00 w(26) = 0.0038283690844171626400743524999D+00 w(27) = 0.0039717489986349171988699773906D+00 w(28) = 0.0041145284389812475901826468094D+00 w(29) = 0.0042566858191260658425395494472D+00 w(30) = 0.0043981996467927779838546384780D+00 w(31) = 0.0045390485270061921259394035112D+00 w(32) = 0.0046792111653260640506279893190D+00 w(33) = 0.0048186663710656988918572043815D+00 w(34) = 0.0049573930604950563104281084148D+00 w(35) = 0.0050953702600278273039420404117D+00 w(36) = 0.0052325771093919661294970523234D+00 w(37) = 0.0053689928647831724787741258653D+00 w(38) = 0.0055045969020008281904902120813D+00 w(39) = 0.0056393687195659001929970994675D+00 w(40) = 0.0057732879418203275712033691864D+00 w(41) = 0.0059063343220074160130475409466D+00 w(42) = 0.0060384877453327676663371666884D+00 w(43) = 0.0061697282320052788060812561217D+00 w(44) = 0.0063000359402577418025981070425D+00 w(45) = 0.0064293911693465917826140832500D+00 w(46) = 0.0065577743625303421548456356354D+00 w(47) = 0.0066851661100262568757892743568D+00 w(48) = 0.0068115471519448109954345674817D+00 w(49) = 0.0069368983812014946719507501243D+00 w(50) = 0.0070612008464055194979848418291D+00 w(51) = 0.0071844357547249896530757997058D+00 w(52) = 0.0073065844747281040972736443146D+00 w(53) = 0.0074276285391999597581348419714D+00 w(54) = 0.0075475496479345294426435656724D+00 w(55) = 0.0076663296705013920315933272426D+00 w(56) = 0.0077839506489867963897419914623D+00 w(57) = 0.0079003948007086443529587296692D+00 w(58) = 0.0080156445209049821352946484008D+00 w(59) = 0.0081296823853955935356080649925D+00 w(60) = 0.0082424911532162924158504385939D+00 w(61) = 0.0083540537692255160718568405530D+00 w(62) = 0.0084643533666828253227353760036D+00 w(63) = 0.0085733732697989214067758505840D+00 w(64) = 0.0086810969962567940901133439612D+00 w(65) = 0.0087875082597036197689825483144D+00 w(66) = 0.0088925909722130327769834298578D+00 w(67) = 0.0089963292467173975949700110383D+00 w(68) = 0.0090987073994097142025303711406D+00 w(69) = 0.0091997099521147934060534414075D+00 w(70) = 0.0092993216346293436285393234867D+00 w(71) = 0.0093975273870306153500305317074D+00 w(72) = 0.0094943123619532541442165010292D+00 w(73) = 0.0095896619268340180657610209655D+00 w(74) = 0.0096835616661240200035669970076D+00 w(75) = 0.0097759973834681605268499842249D+00 w(76) = 0.0098669551038514217128483481814D+00 w(77) = 0.0099564210757116974565448593910D+00 w(78) = 0.0100443817730188408231888789497D+00 w(79) = 0.0101308238973196141129538950955D+00 w(80) = 0.0102157343797482324629939488415D+00 w(81) = 0.0102991003830021970147153502911D+00 w(82) = 0.0103809093032831189224876935085D+00 w(83) = 0.0104611487722022407735015844669D+00 w(84) = 0.0105398066586503673262517188088D+00 w(85) = 0.0106168710706319228563864391054D+00 w(86) = 0.0106923303570628578226139809571D+00 w(87) = 0.0107661731095321330311788312990D+00 w(88) = 0.0108383881640265149842990798832D+00 w(89) = 0.0109089646026184216450603134401D+00 w(90) = 0.0109778917551165634377595759712D+00 w(91) = 0.0110451592006791299277436662993D+00 w(92) = 0.0111107567693892782875426356195D+00 w(93) = 0.0111746745437926853557086684962D+00 w(94) = 0.0112369028603969308303734810332D+00 w(95) = 0.0112974323111324849102690558722D+00 w(96) = 0.0113562537447750795009464486204D+00 w(97) = 0.011413358268329247942299599697D+00 w(98) = 0.011468737248372824084374355981D+00 w(99) = 0.011522382312362197440930930031D+00 w(100) = 0.011574285349898127083439539046D+00 w(101) = 0.011624438513951922901227922331D+00 w(102) = 0.011672834222051808845465154244D+00 w(103) = 0.011719465157429288794653489478D+00 w(104) = 0.011764324270125341726399410909D+00 w(105) = 0.011807404778056278953532930501D+00 w(106) = 0.011848700168039102281222824051D+00 w(107) = 0.011888204196776208064673282076D+00 w(108) = 0.011925910891799288293359117699D+00 w(109) = 0.011961814552372285996633285380D+00 w(110) = 0.011995909750353268455989686823D+00 w(111) = 0.012028191331015087920350431142D+00 w(112) = 0.012058654413824705751531083631D+00 w(113) = 0.012087294393181062176578184854D+00 w(114) = 0.012114106939111380091025793650D+00 w(115) = 0.012139087997925797641334635250D+00 w(116) = 0.012162233792830230614908682534D+00 w(117) = 0.012183540824497371981177306326D+00 w(118) = 0.012203005871595742256331865516D+00 w(119) = 0.012220625991276710706457005806D+00 w(120) = 0.012236398519619413758040249691D+00 w(121) = 0.012250321072033503350218104906D+00 w(122) = 0.012262391543619664338660618398D+00 w(123) = 0.012272608109487846445745237751D+00 w(124) = 0.012280969225033162644659793962D+00 w(125) = 0.012287473626169412265336919908D+00 w(126) = 0.012292120329520193516690694701D+00 w(127) = 0.012294908632567576531532225710D+00 w(128) = 0.01229583811375831445681490730D+00 w(129) = 0.012294908632567576531532225710D+00 w(130) = 0.012292120329520193516690694701D+00 w(131) = 0.012287473626169412265336919908D+00 w(132) = 0.012280969225033162644659793962D+00 w(133) = 0.012272608109487846445745237751D+00 w(134) = 0.012262391543619664338660618398D+00 w(135) = 0.012250321072033503350218104906D+00 w(136) = 0.012236398519619413758040249691D+00 w(137) = 0.012220625991276710706457005806D+00 w(138) = 0.012203005871595742256331865516D+00 w(139) = 0.012183540824497371981177306326D+00 w(140) = 0.012162233792830230614908682534D+00 w(141) = 0.012139087997925797641334635250D+00 w(142) = 0.012114106939111380091025793650D+00 w(143) = 0.012087294393181062176578184854D+00 w(144) = 0.012058654413824705751531083631D+00 w(145) = 0.012028191331015087920350431142D+00 w(146) = 0.011995909750353268455989686823D+00 w(147) = 0.011961814552372285996633285380D+00 w(148) = 0.011925910891799288293359117699D+00 w(149) = 0.011888204196776208064673282076D+00 w(150) = 0.011848700168039102281222824051D+00 w(151) = 0.011807404778056278953532930501D+00 w(152) = 0.011764324270125341726399410909D+00 w(153) = 0.011719465157429288794653489478D+00 w(154) = 0.011672834222051808845465154244D+00 w(155) = 0.011624438513951922901227922331D+00 w(156) = 0.011574285349898127083439539046D+00 w(157) = 0.011522382312362197440930930031D+00 w(158) = 0.011468737248372824084374355981D+00 w(159) = 0.011413358268329247942299599697D+00 w(160) = 0.0113562537447750795009464486204D+00 w(161) = 0.0112974323111324849102690558722D+00 w(162) = 0.0112369028603969308303734810332D+00 w(163) = 0.0111746745437926853557086684962D+00 w(164) = 0.0111107567693892782875426356195D+00 w(165) = 0.0110451592006791299277436662993D+00 w(166) = 0.0109778917551165634377595759712D+00 w(167) = 0.0109089646026184216450603134401D+00 w(168) = 0.0108383881640265149842990798832D+00 w(169) = 0.0107661731095321330311788312990D+00 w(170) = 0.0106923303570628578226139809571D+00 w(171) = 0.0106168710706319228563864391054D+00 w(172) = 0.0105398066586503673262517188088D+00 w(173) = 0.0104611487722022407735015844669D+00 w(174) = 0.0103809093032831189224876935085D+00 w(175) = 0.0102991003830021970147153502911D+00 w(176) = 0.0102157343797482324629939488415D+00 w(177) = 0.0101308238973196141129538950955D+00 w(178) = 0.0100443817730188408231888789497D+00 w(179) = 0.0099564210757116974565448593910D+00 w(180) = 0.0098669551038514217128483481814D+00 w(181) = 0.0097759973834681605268499842249D+00 w(182) = 0.0096835616661240200035669970076D+00 w(183) = 0.0095896619268340180657610209655D+00 w(184) = 0.0094943123619532541442165010292D+00 w(185) = 0.0093975273870306153500305317074D+00 w(186) = 0.0092993216346293436285393234867D+00 w(187) = 0.0091997099521147934060534414075D+00 w(188) = 0.0090987073994097142025303711406D+00 w(189) = 0.0089963292467173975949700110383D+00 w(190) = 0.0088925909722130327769834298578D+00 w(191) = 0.0087875082597036197689825483144D+00 w(192) = 0.0086810969962567940901133439612D+00 w(193) = 0.0085733732697989214067758505840D+00 w(194) = 0.0084643533666828253227353760036D+00 w(195) = 0.0083540537692255160718568405530D+00 w(196) = 0.0082424911532162924158504385939D+00 w(197) = 0.0081296823853955935356080649925D+00 w(198) = 0.0080156445209049821352946484008D+00 w(199) = 0.0079003948007086443529587296692D+00 w(200) = 0.0077839506489867963897419914623D+00 w(201) = 0.0076663296705013920315933272426D+00 w(202) = 0.0075475496479345294426435656724D+00 w(203) = 0.0074276285391999597581348419714D+00 w(204) = 0.0073065844747281040972736443146D+00 w(205) = 0.0071844357547249896530757997058D+00 w(206) = 0.0070612008464055194979848418291D+00 w(207) = 0.0069368983812014946719507501243D+00 w(208) = 0.0068115471519448109954345674817D+00 w(209) = 0.0066851661100262568757892743568D+00 w(210) = 0.0065577743625303421548456356354D+00 w(211) = 0.0064293911693465917826140832500D+00 w(212) = 0.0063000359402577418025981070425D+00 w(213) = 0.0061697282320052788060812561217D+00 w(214) = 0.0060384877453327676663371666884D+00 w(215) = 0.0059063343220074160130475409466D+00 w(216) = 0.0057732879418203275712033691864D+00 w(217) = 0.0056393687195659001929970994675D+00 w(218) = 0.0055045969020008281904902120813D+00 w(219) = 0.0053689928647831724787741258653D+00 w(220) = 0.0052325771093919661294970523234D+00 w(221) = 0.0050953702600278273039420404117D+00 w(222) = 0.0049573930604950563104281084148D+00 w(223) = 0.0048186663710656988918572043815D+00 w(224) = 0.0046792111653260640506279893190D+00 w(225) = 0.0045390485270061921259394035112D+00 w(226) = 0.0043981996467927779838546384780D+00 w(227) = 0.0042566858191260658425395494472D+00 w(228) = 0.0041145284389812475901826468094D+00 w(229) = 0.0039717489986349171988699773906D+00 w(230) = 0.0038283690844171626400743524999D+00 w(231) = 0.0036844103734499176530742235517D+00 w(232) = 0.0035398946303722552150296713510D+00 w(233) = 0.0033948437040533928255056951665D+00 w(234) = 0.0032492795242943133198690930777D+00 w(235) = 0.0031032240985191112621977893133D+00 w(236) = 0.0029566995084575002760043344138D+00 w(237) = 0.0028097279068204407457332299361D+00 w(238) = 0.0026623315139717112732749157331D+00 w(239) = 0.0025145326145997073931298921370D+00 w(240) = 0.0023663535543962867157201855305D+00 w(241) = 0.0022178167367540171700373764020D+00 w(242) = 0.0020689446195015801533643667413D+00 w(243) = 0.0019197597117132050055085980675D+00 w(244) = 0.0017702845706603213070421243905D+00 w(245) = 0.0016205417990415653896921100325D+00 w(246) = 0.0014705540427783843160097204304D+00 w(247) = 0.0013203439900221692090523602144D+00 w(248) = 0.0011699343729388079886897709773D+00 w(249) = 0.00101934797642732530281229369360D+00 w(250) = 0.00086860766611945667949717690640D+00 w(251) = 0.00071773647800611087798371518325D+00 w(252) = 0.00056675794564824918946626058353D+00 w(253) = 0.00041569762526823913616284210066D+00 w(254) = 0.00026459387119083065532790838855D+00 w(255) = 0.00011367361999142272115645954414D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEGENDRE_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1, 3, 7, 15, 31, 63, 127, and 255.' stop end if return end subroutine p00_exact ( problem, exact ) !*****************************************************************************80 ! !! P00_EXACT returns the exact integral for any problem. ! ! Discussion: ! ! This routine provides a "generic" interface to the exact integral ! routines for the various problems, and allows a problem to be called ! by index rather than by name. ! ! In some cases, the "exact" value of the integral is in fact ! merely a respectable approximation. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer PROBLEM, the problem index. ! ! Output, real ( kind = rk ) EXACT, the exact value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact integer problem if ( problem == 1 ) then call p01_exact ( exact ) else if ( problem == 2 ) then call p02_exact ( exact ) else if ( problem == 3 ) then call p03_exact ( exact ) else if ( problem == 4 ) then call p04_exact ( exact ) else if ( problem == 5 ) then call p05_exact ( exact ) else if ( problem == 6 ) then call p06_exact ( exact ) else if ( problem == 7 ) then call p07_exact ( exact ) else if ( problem == 8 ) then call p08_exact ( exact ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_EXACT - Fatal error!' write ( *, '(a,i8)' ) ' Illegal problem index = ', problem stop end if return end subroutine p00_fun ( problem, n, x, fx ) !*****************************************************************************80 ! !! P00_FUN evaluates the integrand for any problem. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer PROBLEM, the problem index. ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) integer problem real ( kind = rk ) x(2,n) if ( problem == 1 ) then call p01_fun ( n, x, fx ) else if ( problem == 2 ) then call p02_fun ( n, x, fx ) else if ( problem == 3 ) then call p03_fun ( n, x, fx ) else if ( problem == 4 ) then call p04_fun ( n, x, fx ) else if ( problem == 5 ) then call p05_fun ( n, x, fx ) else if ( problem == 6 ) then call p06_fun ( n, x, fx ) else if ( problem == 7 ) then call p07_fun ( n, x, fx ) else if ( problem == 8 ) then call p08_fun ( n, x, fx ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_FUN - Fatal error!' write ( *, '(a,i8)' ) ' Illegal problem index = ', problem stop end if return end subroutine p00_lim ( problem, a, b ) !*****************************************************************************80 ! !! P00_LIM returns the integration limits for any problem. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer PROBLEM, the problem index. ! ! Output, real ( kind = rk ) A(2), B(2), the lower and upper limits ! of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) integer problem if ( problem == 1 ) then call p01_lim ( a, b ) else if ( problem == 2 ) then call p02_lim ( a, b ) else if ( problem == 3 ) then call p03_lim ( a, b ) else if ( problem == 4 ) then call p04_lim ( a, b ) else if ( problem == 5 ) then call p05_lim ( a, b ) else if ( problem == 6 ) then call p06_lim ( a, b ) else if ( problem == 7 ) then call p07_lim ( a, b ) else if ( problem == 8 ) then call p08_lim ( a, b ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'P00_LIM - Fatal error!' write ( *, '(a,i8)' ) ' Illegal problem index = ', problem stop end if return end subroutine p00_problem_num ( problem_num ) !*****************************************************************************80 ! !! P00_PROBLEM_NUM returns the number of test integration problems. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer PROBLEM_NUM, the number of problems. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer problem_num problem_num = 8 return end subroutine p01_exact ( exact ) !*****************************************************************************80 ! !! P01_EXACT returns the exact integral for problem 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 exact = pi * pi / 6.0D+00 return end subroutine p01_fun ( n, x, fx ) !*****************************************************************************80 ! !! P01_FUN evaluates the integrand for problem 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = 1.0D+00 / ( 1.0D+00 - x(1,1:n) * x(2,1:n) ) return end subroutine p01_lim ( a, b ) !*****************************************************************************80 ! !! P01_LIM returns the integration limits for problem 1. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = 0.0D+00 b(1:2) = 1.0D+00 return end subroutine p02_exact ( exact ) !*****************************************************************************80 ! !! P02_EXACT returns the exact integral for problem 2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 exact = 2.0D+00 * pi * log ( 2.0D+00 ) return end subroutine p02_fun ( n, x, fx ) !*****************************************************************************80 ! !! P02_FUN evaluates the integrand for problem 2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = 1.0D+00 / sqrt ( 1.0D+00 - x(1,1:n)**2 * x(2,1:n)**2 ) return end subroutine p02_lim ( a, b ) !*****************************************************************************80 ! !! P02_LIM returns the integration limits for problem 2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = -1.0D+00 b(1:2) = 1.0D+00 return end subroutine p03_exact ( exact ) !*****************************************************************************80 ! !! P03_EXACT returns the exact integral for problem 3. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact exact = ( 16.0D+00 / 3.0D+00 ) * ( 2.0D+00 - sqrt ( 2.0D+00 ) ) return end subroutine p03_fun ( n, x, fx ) !*****************************************************************************80 ! !! P03_FUN evaluates the integrand for problem 3. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = 1.0D+00 / sqrt ( 2.0D+00 - x(1,1:n) - x(2,1:n) ) return end subroutine p03_lim ( a, b ) !*****************************************************************************80 ! !! P03_LIM returns the integration limits for problem 3. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = -1.0D+00 b(1:2) = 1.0D+00 return end subroutine p04_exact ( exact ) !*****************************************************************************80 ! !! P04_EXACT returns the exact integral for problem 4. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact exact = ( sqrt ( 32.0D+00 ) / 3.0D+00 ) & * ( sqrt ( 27.0D+00 ) - sqrt ( 8.0D+00 ) - 1.0D+00 ) return end subroutine p04_fun ( n, x, fx ) !*****************************************************************************80 ! !! P04_FUN evaluates the integrand for problem 4. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = 1.0D+00 / sqrt ( 3.0D+00 - x(1,1:n) - 2.0D+00 * x(2,1:n) ) return end subroutine p04_lim ( a, b ) !*****************************************************************************80 ! !! P04_LIM returns the integration limits for problem 4. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = -1.0D+00 b(1:2) = 1.0D+00 return end subroutine p05_exact ( exact ) !*****************************************************************************80 ! !! P05_EXACT returns the exact integral for problem 5. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact exact = 4.0D+00 / 9.0D+00 return end subroutine p05_fun ( n, x, fx ) !*****************************************************************************80 ! !! P05_FUN evaluates the integrand for problem 5. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = sqrt ( x(1,1:n) * x(2,1:n) ) return end subroutine p05_lim ( a, b ) !*****************************************************************************80 ! !! P05_LIM returns the integration limits for problem 5. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = 0.0D+00 b(1:2) = 1.0D+00 return end subroutine p06_exact ( exact ) !*****************************************************************************80 ! !! P06_EXACT returns the exact integral for problem 6. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 exact = ( 5.0D+00 / 3.0D+00 ) + ( pi / 16.0D+00 ) return end subroutine p06_fun ( n, x, fx ) !*****************************************************************************80 ! !! P06_FUN evaluates the integrand for problem 6. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = abs ( x(1,1:n)**2 + x(2,1:n)**2 - 0.25D+00 ) return end subroutine p06_lim ( a, b ) !*****************************************************************************80 ! !! P06_LIM returns the integration limits for problem 6. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = -1.0D+00 b(1:2) = 1.0D+00 return end subroutine p07_exact ( exact ) !*****************************************************************************80 ! !! P07_EXACT returns the exact integral for problem 7. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact exact = 8.0D+00 / 15.0D+00 return end subroutine p07_fun ( n, x, fx ) !*****************************************************************************80 ! !! P07_FUN evaluates the integrand for problem 7. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Gwynne Evans, ! Practical Numerical Integration, ! Wiley, 1993, ! ISBN: 047193898X, ! LC: QA299.3E93. ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = sqrt ( abs ( x(1,1:n) - x(2,1:n) ) ) return end subroutine p07_lim ( a, b ) !*****************************************************************************80 ! !! P07_LIM returns the integration limits for problem 7. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 January 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = 0.0D+00 b(1:2) = 1.0D+00 return end subroutine p08_exact ( exact ) !*****************************************************************************80 ! !! P08_EXACT returns the exact integral for problem 8. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) EXACT, the value of the integral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) exact real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) r8_erf exact = 0.25D+00 * pi * ( r8_erf ( 1.0D+00 ) + r8_erf ( 4.0D+00 ) )**2 return end subroutine p08_fun ( n, x, fx ) !*****************************************************************************80 ! !! P08_FUN evaluates the integrand for problem 8. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of evaluation points. ! ! Input, real ( kind = rk ) X(2,N), the evaluation points. ! ! Output, real ( kind = rk ) FX(N), the integrand values. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) fx(n) real ( kind = rk ) x(2,n) fx(1:n) = exp ( - ( ( x(1,1:n) - 4.0D+00 )**2 + ( x(2,1:n) - 1.0D+00 )**2 ) ) return end subroutine p08_lim ( a, b ) !*****************************************************************************80 ! !! P08_LIM returns the integration limits for problem 8. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 16 September 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = rk ) A(2), B(2), the limits of integration. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) a(2) real ( kind = rk ) b(2) a(1:2) = 0.0D+00 b(1:2) = 5.0D+00 return end function r8_csevl ( x, a, n ) !*****************************************************************************80 ! !! R8_CSEVL evaluates a Chebyshev series. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2011 ! ! Author: ! ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Roger Broucke, ! Algorithm 446: ! Ten Subroutines for the Manipulation of Chebyshev Series, ! Communications of the ACM, ! Volume 16, Number 4, April 1973, pages 254-256. ! ! Parameters: ! ! Input, real ( kind = rk ) X, the evaluation point. ! ! Input, real ( kind = rk ) CS(N), the Chebyshev coefficients. ! ! Input, integer N, the number of Chebyshev coefficients. ! ! Output, real ( kind = rk ) R8_CSEVL, the Chebyshev series evaluated at X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) real ( kind = rk ) b0 real ( kind = rk ) b1 real ( kind = rk ) b2 integer i real ( kind = rk ) r8_csevl real ( kind = rk ) twox real ( kind = rk ) x if ( n < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_CSEVL - Fatal error!' write ( *, '(a)' ) ' Number of terms <= 0.' stop end if if ( 1000 < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_CSEVL - Fatal error!' write ( *, '(a)' ) ' Number of terms > 1000.' stop end if if ( x < -1.1D+00 .or. 1.1D+00 < x ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_CSEVL - Fatal error!' write ( *, '(a)' ) ' X outside (-1,+1)' write ( *, '(a,g14.6)' ) ' X = ', x stop end if twox = 2.0D+00 * x b1 = 0.0D+00 b0 = 0.0D+00 do i = n, 1, -1 b2 = b1 b1 = b0 b0 = twox * b1 - b2 + a(i) end do r8_csevl = 0.5D+00 * ( b0 - b2 ) return end function r8_erf ( x ) !*****************************************************************************80 ! !! R8_ERF evaluates the error function of an R8 argument. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2011 ! ! Author: ! ! Original FORTRAN77 version by Wayne Fullerton. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Wayne Fullerton, ! Portable Special Function Routines, ! in Portability of Numerical Software, ! edited by Wayne Cowell, ! Lecture Notes in Computer Science, Volume 57, ! Springer 1977, ! ISBN: 978-3-540-08446-4, ! LC: QA297.W65. ! ! Parameters: ! ! Input, real ( kind = rk ) X, the argument. ! ! Output, real ( kind = rk ) R8_ERF, the error function of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) erfcs(21) integer nterf real ( kind = rk ) r8_csevl real ( kind = rk ) r8_erf real ( kind = rk ) r8_erfc integer r8_inits real ( kind = rk ) r8_mach real ( kind = rk ), save :: sqeps = 0.0D+00 real ( kind = rk ) sqrtpi real ( kind = rk ) value real ( kind = rk ) x real ( kind = rk ) xbig real ( kind = rk ) y save erfcs save nterf save sqrtpi save xbig data erfcs( 1) / -0.49046121234691808039984544033376D-01 / data erfcs( 2) / -0.14226120510371364237824741899631D+00 / data erfcs( 3) / +0.10035582187599795575754676712933D-01 / data erfcs( 4) / -0.57687646997674847650827025509167D-03 / data erfcs( 5) / +0.27419931252196061034422160791471D-04 / data erfcs( 6) / -0.11043175507344507604135381295905D-05 / data erfcs( 7) / +0.38488755420345036949961311498174D-07 / data erfcs( 8) / -0.11808582533875466969631751801581D-08 / data erfcs( 9) / +0.32334215826050909646402930953354D-10 / data erfcs( 10) / -0.79910159470045487581607374708595D-12 / data erfcs( 11) / +0.17990725113961455611967245486634D-13 / data erfcs( 12) / -0.37186354878186926382316828209493D-15 / data erfcs( 13) / +0.71035990037142529711689908394666D-17 / data erfcs( 14) / -0.12612455119155225832495424853333D-18 / data erfcs( 15) / +0.20916406941769294369170500266666D-20 / data erfcs( 16) / -0.32539731029314072982364160000000D-22 / data erfcs( 17) / +0.47668672097976748332373333333333D-24 / data erfcs( 18) / -0.65980120782851343155199999999999D-26 / data erfcs( 19) / +0.86550114699637626197333333333333D-28 / data erfcs( 20) / -0.10788925177498064213333333333333D-29 / data erfcs( 21) / +0.12811883993017002666666666666666D-31 / data nterf / 0 / data sqrtpi / 1.77245385090551602729816748334115D+00 / data xbig / 0.0D+00 / if ( nterf == 0 ) then nterf = r8_inits ( erfcs, 21, 0.1D+00 * r8_mach ( 3 ) ) xbig = sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) ) sqeps = sqrt ( 2.0D+00 * r8_mach ( 3 ) ) end if y = abs ( x ) if ( y <= sqeps ) then value = 2.0D+00 * x / sqrtpi else if ( y <= 1.0D+00 ) then value = x * ( 1.0D+00 & + r8_csevl ( 2.0D+00 * x * x - 1.0D+00, erfcs, nterf ) ) else if ( y <= xbig ) then value = 1.0D+00 - r8_erfc ( y ) if ( x < 0.0D+00 ) then value = - value end if else value = 1.0D+00 if ( x < 0.0D+00 ) then value = - value end if end if r8_erf = value return end function r8_erfc ( x ) !*****************************************************************************80 ! !! R8_ERFC evaluates the co-error function of an R8 argument. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2011 ! ! Author: ! ! Original FORTRAN77 version by Wayne Fullerton. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Wayne Fullerton, ! Portable Special Function Routines, ! in Portability of Numerical Software, ! edited by Wayne Cowell, ! Lecture Notes in Computer Science, Volume 57, ! Springer 1977, ! ISBN: 978-3-540-08446-4, ! LC: QA297.W65. ! ! Parameters: ! ! Input, real ( kind = rk ) X, the argument. ! ! Output, real ( kind = rk ) R8_ERFC, the co-error function of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) erc2cs(49) real ( kind = rk ) erfccs(59) real ( kind = rk ) erfcs(21) real ( kind = rk ) eta integer nterc2 integer nterf integer nterfc real ( kind = rk ) r8_csevl real ( kind = rk ) r8_erfc integer r8_inits real ( kind = rk ) r8_mach real ( kind = rk ), save :: sqeps = 0.0D+00 real ( kind = rk ) sqrtpi real ( kind = rk ) x real ( kind = rk ) xmax real ( kind = rk ) xsml real ( kind = rk ) y save erfccs save erfcs save erc2cs save nterc2 save nterf save nterfc save sqrtpi save xmax save xsml data erfcs( 1) / -0.49046121234691808039984544033376D-01 / data erfcs( 2) / -0.14226120510371364237824741899631D+00 / data erfcs( 3) / +0.10035582187599795575754676712933D-01 / data erfcs( 4) / -0.57687646997674847650827025509167D-03 / data erfcs( 5) / +0.27419931252196061034422160791471D-04 / data erfcs( 6) / -0.11043175507344507604135381295905D-05 / data erfcs( 7) / +0.38488755420345036949961311498174D-07 / data erfcs( 8) / -0.11808582533875466969631751801581D-08 / data erfcs( 9) / +0.32334215826050909646402930953354D-10 / data erfcs( 10) / -0.79910159470045487581607374708595D-12 / data erfcs( 11) / +0.17990725113961455611967245486634D-13 / data erfcs( 12) / -0.37186354878186926382316828209493D-15 / data erfcs( 13) / +0.71035990037142529711689908394666D-17 / data erfcs( 14) / -0.12612455119155225832495424853333D-18 / data erfcs( 15) / +0.20916406941769294369170500266666D-20 / data erfcs( 16) / -0.32539731029314072982364160000000D-22 / data erfcs( 17) / +0.47668672097976748332373333333333D-24 / data erfcs( 18) / -0.65980120782851343155199999999999D-26 / data erfcs( 19) / +0.86550114699637626197333333333333D-28 / data erfcs( 20) / -0.10788925177498064213333333333333D-29 / data erfcs( 21) / +0.12811883993017002666666666666666D-31 / data erc2cs( 1) / -0.6960134660230950112739150826197D-01 / data erc2cs( 2) / -0.4110133936262089348982212084666D-01 / data erc2cs( 3) / +0.3914495866689626881561143705244D-02 / data erc2cs( 4) / -0.4906395650548979161280935450774D-03 / data erc2cs( 5) / +0.7157479001377036380760894141825D-04 / data erc2cs( 6) / -0.1153071634131232833808232847912D-04 / data erc2cs( 7) / +0.1994670590201997635052314867709D-05 / data erc2cs( 8) / -0.3642666471599222873936118430711D-06 / data erc2cs( 9) / +0.6944372610005012589931277214633D-07 / data erc2cs( 10) / -0.1371220902104366019534605141210D-07 / data erc2cs( 11) / +0.2788389661007137131963860348087D-08 / data erc2cs( 12) / -0.5814164724331161551864791050316D-09 / data erc2cs( 13) / +0.1238920491752753181180168817950D-09 / data erc2cs( 14) / -0.2690639145306743432390424937889D-10 / data erc2cs( 15) / +0.5942614350847910982444709683840D-11 / data erc2cs( 16) / -0.1332386735758119579287754420570D-11 / data erc2cs( 17) / +0.3028046806177132017173697243304D-12 / data erc2cs( 18) / -0.6966648814941032588795867588954D-13 / data erc2cs( 19) / +0.1620854541053922969812893227628D-13 / data erc2cs( 20) / -0.3809934465250491999876913057729D-14 / data erc2cs( 21) / +0.9040487815978831149368971012975D-15 / data erc2cs( 22) / -0.2164006195089607347809812047003D-15 / data erc2cs( 23) / +0.5222102233995854984607980244172D-16 / data erc2cs( 24) / -0.1269729602364555336372415527780D-16 / data erc2cs( 25) / +0.3109145504276197583836227412951D-17 / data erc2cs( 26) / -0.7663762920320385524009566714811D-18 / data erc2cs( 27) / +0.1900819251362745202536929733290D-18 / data erc2cs( 28) / -0.4742207279069039545225655999965D-19 / data erc2cs( 29) / +0.1189649200076528382880683078451D-19 / data erc2cs( 30) / -0.3000035590325780256845271313066D-20 / data erc2cs( 31) / +0.7602993453043246173019385277098D-21 / data erc2cs( 32) / -0.1935909447606872881569811049130D-21 / data erc2cs( 33) / +0.4951399124773337881000042386773D-22 / data erc2cs( 34) / -0.1271807481336371879608621989888D-22 / data erc2cs( 35) / +0.3280049600469513043315841652053D-23 / data erc2cs( 36) / -0.8492320176822896568924792422399D-24 / data erc2cs( 37) / +0.2206917892807560223519879987199D-24 / data erc2cs( 38) / -0.5755617245696528498312819507199D-25 / data erc2cs( 39) / +0.1506191533639234250354144051199D-25 / data erc2cs( 40) / -0.3954502959018796953104285695999D-26 / data erc2cs( 41) / +0.1041529704151500979984645051733D-26 / data erc2cs( 42) / -0.2751487795278765079450178901333D-27 / data erc2cs( 43) / +0.7290058205497557408997703680000D-28 / data erc2cs( 44) / -0.1936939645915947804077501098666D-28 / data erc2cs( 45) / +0.5160357112051487298370054826666D-29 / data erc2cs( 46) / -0.1378419322193094099389644800000D-29 / data erc2cs( 47) / +0.3691326793107069042251093333333D-30 / data erc2cs( 48) / -0.9909389590624365420653226666666D-31 / data erc2cs( 49) / +0.2666491705195388413323946666666D-31 / data erfccs( 1) / +0.715179310202924774503697709496D-01 / data erfccs( 2) / -0.265324343376067157558893386681D-01 / data erfccs( 3) / +0.171115397792085588332699194606D-02 / data erfccs( 4) / -0.163751663458517884163746404749D-03 / data erfccs( 5) / +0.198712935005520364995974806758D-04 / data erfccs( 6) / -0.284371241276655508750175183152D-05 / data erfccs( 7) / +0.460616130896313036969379968464D-06 / data erfccs( 8) / -0.822775302587920842057766536366D-07 / data erfccs( 9) / +0.159214187277090112989358340826D-07 / data erfccs( 10) / -0.329507136225284321486631665072D-08 / data erfccs( 11) / +0.722343976040055546581261153890D-09 / data erfccs( 12) / -0.166485581339872959344695966886D-09 / data erfccs( 13) / +0.401039258823766482077671768814D-10 / data erfccs( 14) / -0.100481621442573113272170176283D-10 / data erfccs( 15) / +0.260827591330033380859341009439D-11 / data erfccs( 16) / -0.699111056040402486557697812476D-12 / data erfccs( 17) / +0.192949233326170708624205749803D-12 / data erfccs( 18) / -0.547013118875433106490125085271D-13 / data erfccs( 19) / +0.158966330976269744839084032762D-13 / data erfccs( 20) / -0.472689398019755483920369584290D-14 / data erfccs( 21) / +0.143587337678498478672873997840D-14 / data erfccs( 22) / -0.444951056181735839417250062829D-15 / data erfccs( 23) / +0.140481088476823343737305537466D-15 / data erfccs( 24) / -0.451381838776421089625963281623D-16 / data erfccs( 25) / +0.147452154104513307787018713262D-16 / data erfccs( 26) / -0.489262140694577615436841552532D-17 / data erfccs( 27) / +0.164761214141064673895301522827D-17 / data erfccs( 28) / -0.562681717632940809299928521323D-18 / data erfccs( 29) / +0.194744338223207851429197867821D-18 / data erfccs( 30) / -0.682630564294842072956664144723D-19 / data erfccs( 31) / +0.242198888729864924018301125438D-19 / data erfccs( 32) / -0.869341413350307042563800861857D-20 / data erfccs( 33) / +0.315518034622808557122363401262D-20 / data erfccs( 34) / -0.115737232404960874261239486742D-20 / data erfccs( 35) / +0.428894716160565394623737097442D-21 / data erfccs( 36) / -0.160503074205761685005737770964D-21 / data erfccs( 37) / +0.606329875745380264495069923027D-22 / data erfccs( 38) / -0.231140425169795849098840801367D-22 / data erfccs( 39) / +0.888877854066188552554702955697D-23 / data erfccs( 40) / -0.344726057665137652230718495566D-23 / data erfccs( 41) / +0.134786546020696506827582774181D-23 / data erfccs( 42) / -0.531179407112502173645873201807D-24 / data erfccs( 43) / +0.210934105861978316828954734537D-24 / data erfccs( 44) / -0.843836558792378911598133256738D-25 / data erfccs( 45) / +0.339998252494520890627359576337D-25 / data erfccs( 46) / -0.137945238807324209002238377110D-25 / data erfccs( 47) / +0.563449031183325261513392634811D-26 / data erfccs( 48) / -0.231649043447706544823427752700D-26 / data erfccs( 49) / +0.958446284460181015263158381226D-27 / data erfccs( 50) / -0.399072288033010972624224850193D-27 / data erfccs( 51) / +0.167212922594447736017228709669D-27 / data erfccs( 52) / -0.704599152276601385638803782587D-28 / data erfccs( 53) / +0.297976840286420635412357989444D-28 / data erfccs( 54) / -0.126252246646061929722422632994D-28 / data erfccs( 55) / +0.539543870454248793985299653154D-29 / data erfccs( 56) / -0.238099288253145918675346190062D-29 / data erfccs( 57) / +0.109905283010276157359726683750D-29 / data erfccs( 58) / -0.486771374164496572732518677435D-30 / data erfccs( 59) / +0.152587726411035756763200828211D-30 / data nterc2 / 0 / data nterf / 0 / data nterfc / 0 / data sqrtpi / 1.77245385090551602729816748334115D+00 / data xmax / 0.0D+00 / data xsml / 0.0D+00 / if ( nterf == 0 ) then eta = 0.1D+00 * r8_mach ( 3 ) nterf = r8_inits ( erfcs, 21, eta ) nterfc = r8_inits ( erfccs, 59, eta ) nterc2 = r8_inits ( erc2cs, 49, eta ) xsml = - sqrt ( - log ( sqrtpi * r8_mach ( 3 ) ) ) xmax = sqrt (- log ( sqrtpi * r8_mach ( 1 ) ) ) xmax = xmax - 0.5D+00 * log ( xmax ) / xmax - 0.01D+00 sqeps = sqrt ( 2.0D+00 * r8_mach ( 3 ) ) end if if ( x <= xsml ) then r8_erfc = 2.0D+00 return end if if ( xmax < x ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_ERFC - Warning!' write ( *, '(a)' ) ' X so big that ERFC underflows.' r8_erfc = 0.0D+00 return end if y = abs ( x ) if ( y < sqeps ) then r8_erfc = 1.0D+00 - 2.0D+00 * x / sqrtpi return else if ( y <= 1.0D+00 ) then r8_erfc = 1.0D+00 - x * ( 1.0D+00 & + r8_csevl ( 2.0D+00 * x * x - 1.0D+00, erfcs, nterf ) ) return end if y = y * y if ( y <= 4.0D+00 ) then r8_erfc = exp ( - y ) / abs ( x ) * ( 0.5D+00 & + r8_csevl ( ( 8.0D+00 / y - 5.0D+00 ) / 3.0D+00, erc2cs, nterc2 ) ) else r8_erfc = exp ( - y ) / abs ( x ) * ( 0.5D+00 & + r8_csevl ( 8.0D+00 / y - 1.0D+00, erfccs, nterfc ) ) end if if ( x < 0.0D+00 ) then r8_erfc = 2.0D+00 - r8_erfc end if return end function r8_inits ( dos, nos, eta ) !*****************************************************************************80 ! !! R8_INITS initializes a Chebyshev series. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 September 2011 ! ! Author: ! ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Roger Broucke, ! Algorithm 446: ! Ten Subroutines for the Manipulation of Chebyshev Series, ! Communications of the ACM, ! Volume 16, Number 4, April 1973, pages 254-256. ! ! Parameters: ! ! Input, real ( kind = rk ) DOS(NOS), the Chebyshev coefficients. ! ! Input, integer NOS, the number of coefficients. ! ! Input, real ( kind = rk ) ETA, the desired accuracy. ! ! Output, integer R8_INITS, the number of terms of the series ! needed to ensure the requested accuracy. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nos real ( kind = rk ) dos(nos) real ( kind = rk ) err real ( kind = rk ) eta integer i integer r8_inits if ( nos < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_INITS - Fatal error!' write ( *, '(a)' ) ' Number of coefficients < 1.' stop end if err = 0.0D+00 do i = nos, 1, -1 err = err + abs ( dos(i) ) if ( eta < err ) then r8_inits = i return end if end do r8_inits = nos write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_INITS - Warning!' write ( *, '(a)' ) ' ETA may be too small.' return end function r8_mach ( i ) !*****************************************************************************80 ! !! R8_MACH returns real ( kind = rk ) real machine-dependent constants. ! ! Discussion: ! ! R8_MACH can be used to obtain machine-dependent parameters ! for the local machine environment. It is a function ! with one input argument, and can be called as follows: ! ! D = R8_MACH ( I ) ! ! where I=1,...,5. The output value of D above is ! determined by the input value of I:. ! ! R8_MACH ( 1) = B^(EMIN-1), the smallest positive magnitude. ! R8_MACH ( 2) = B^EMAX*(1 - B**(-T)), the largest magnitude. ! R8_MACH ( 3) = B^(-T), the smallest relative spacing. ! R8_MACH ( 4) = B^(1-T), the largest relative spacing. ! R8_MACH ( 5) = LOG10(B) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 April 2007 ! ! Author: ! ! Original FORTRAN77 version by Phyllis Fox, Andrew Hall, Norman Schryer. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Phyllis Fox, Andrew Hall, Norman Schryer, ! Algorithm 528: ! Framework for a Portable Library, ! ACM Transactions on Mathematical Software, ! Volume 4, Number 2, June 1978, page 176-188. ! ! Parameters: ! ! Input, integer I, the index of the desired constant. ! ! Output, real ( kind = rk ) R8_MACH, the value of the constant. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) r8_mach integer i if ( i < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i r8_mach = 0.0D+00 stop else if ( i == 1 ) then r8_mach = 4.450147717014403D-308 else if ( i == 2 ) then r8_mach = 8.988465674311579D+307 else if ( i == 3 ) then r8_mach = 1.110223024625157D-016 else if ( i == 4 ) then r8_mach = 2.220446049250313D-016 else if ( i == 5 ) then r8_mach = 0.301029995663981D+000 else if ( 5 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8_MACH - Fatal error!' write ( *, '(a)' ) ' The input argument I is out of bounds.' write ( *, '(a)' ) ' Legal values satisfy 1 <= I <= 5.' write ( *, '(a,i12)' ) ' I = ', i r8_mach = 0.0D+00 stop end if 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 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 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 integer, parameter :: rk = kind ( 1.0D+00 ) 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.2,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