quadcc.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2010-2012 Pedro Gonnet
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include <stdlib.h>
00028 #include "lo-math.h"
00029 #include "lo-ieee.h"
00030 #include "oct.h"
00031 #include "parse.h"
00032 #include "ov-fcn-handle.h"
00033 
00034 /* Define the size of the interval heap. */
00035 #define cquad_heapsize                  200
00036 
00037 
00038 /* Data of a single interval */
00039 typedef struct
00040 {
00041   double a, b;
00042   double c[64];
00043   double fx[33];
00044   double igral, err;
00045   int depth, rdepth, ndiv;
00046 } cquad_ival;
00047 
00048 /* Some constants and matrices that we'll need.  */
00049 
00050 static const double xi[33] = {
00051   -1., -0.99518472667219688624, -0.98078528040323044912,
00052   -0.95694033573220886493, -0.92387953251128675612,
00053   -0.88192126434835502970, -0.83146961230254523708,
00054   -0.77301045336273696082, -0.70710678118654752440,
00055   -0.63439328416364549822, -0.55557023301960222475,
00056   -0.47139673682599764857, -0.38268343236508977173,
00057   -0.29028467725446236764, -0.19509032201612826785,
00058   -0.098017140329560601995, 0., 0.098017140329560601995,
00059   0.19509032201612826785, 0.29028467725446236764, 0.38268343236508977173,
00060   0.47139673682599764857, 0.55557023301960222475, 0.63439328416364549822,
00061   0.70710678118654752440, 0.77301045336273696082, 0.83146961230254523708,
00062   0.88192126434835502970, 0.92387953251128675612, 0.95694033573220886493,
00063   0.98078528040323044912, 0.99518472667219688624, 1.
00064 };
00065 
00066 static const double bee[68] = {
00067   0.00000000000000e+00, 2.28868854108532e-01, 0.00000000000000e+00,
00068   -8.15740215243451e-01, 0.00000000000000e+00, 5.31212715259731e-01,
00069   0.00000000000000e+00, 1.38538036812454e-02, 0.00000000000000e+00,
00070   3.74405228908818e-02, 0.00000000000000e+00, 2.12224115039342e-01,
00071   0.00000000000000e+00, -8.16362644507898e-01, 0.00000000000000e+00,
00072   5.35648426691481e-01, 0.00000000000000e+00, 1.52417902753662e-03,
00073   0.00000000000000e+00, 2.63058840550873e-03, 0.00000000000000e+00,
00074   4.15292106318904e-03, 0.00000000000000e+00, 6.97106011119775e-03,
00075   0.00000000000000e+00, 1.35535708431058e-02, 0.00000000000000e+00,
00076   3.52132898424856e-02, 0.00000000000000e+00, 2.06946714741884e-01,
00077   0.00000000000000e+00, -8.15674251283876e-01, 0.00000000000000e+00,
00078   5.38841175520580e-01, 0.00000000000000e+00, 1.84909689577590e-04,
00079   0.00000000000000e+00, 2.90936325007499e-04, 0.00000000000000e+00,
00080   3.84877750950089e-04, 0.00000000000000e+00, 4.86436656735046e-04,
00081   0.00000000000000e+00, 6.08688640346879e-04, 0.00000000000000e+00,
00082   7.66732830740331e-04, 0.00000000000000e+00, 9.82753336104205e-04,
00083   0.00000000000000e+00, 1.29359957505615e-03, 0.00000000000000e+00,
00084   1.76616363801885e-03, 0.00000000000000e+00, 2.53323433039089e-03,
00085   0.00000000000000e+00, 3.88872172121956e-03, 0.00000000000000e+00,
00086   6.58635106468291e-03, 0.00000000000000e+00, 1.30326736343254e-02,
00087   0.00000000000000e+00, 3.44353850696714e-02, 0.00000000000000e+00,
00088   2.05025409531915e-01, 0.00000000000000e+00, -8.14985893995401e-01,
00089   0.00000000000000e+00, 5.40679930965238e-01
00090 };
00091 
00092 static const double Lalpha[33] = {
00093   5.77350269189626e-01, 5.16397779494322e-01, 5.07092552837110e-01,
00094   5.03952630678970e-01, 5.02518907629606e-01, 5.01745206004255e-01,
00095   5.01280411827603e-01, 5.00979432868120e-01, 5.00773395667191e-01,
00096   5.00626174321759e-01, 5.00517330712619e-01, 5.00434593736979e-01,
00097   5.00370233297676e-01, 5.00319182924304e-01, 5.00278009473803e-01,
00098   5.00244319584578e-01, 5.00216403386025e-01, 5.00193012939056e-01,
00099   5.00173220168024e-01, 5.00156323280355e-01, 5.00141783641018e-01,
00100   5.00129182278347e-01, 5.00118189340972e-01, 5.00108542278496e-01,
00101   5.00100030010004e-01, 5.00092481273333e-01, 5.00085755939229e-01,
00102   5.00079738458365e-01, 5.00074332862969e-01, 5.00069458915387e-01,
00103   5.00065049112355e-01, 5.00061046334395e-01, 5.00057401986298e-01
00104 };
00105 
00106 static const double Lgamma[33] = {
00107   0.0, 0.0, 5.16397779494322e-01, 5.07092552837110e-01, 5.03952630678970e-01,
00108   5.02518907629606e-01, 5.01745206004255e-01, 5.01280411827603e-01,
00109   5.00979432868120e-01, 5.00773395667191e-01, 5.00626174321759e-01,
00110   5.00517330712619e-01, 5.00434593736979e-01, 5.00370233297676e-01,
00111   5.00319182924304e-01, 5.00278009473803e-01, 5.00244319584578e-01,
00112   5.00216403386025e-01, 5.00193012939056e-01, 5.00173220168024e-01,
00113   5.00156323280355e-01, 5.00141783641018e-01, 5.00129182278347e-01,
00114   5.00118189340972e-01, 5.00108542278496e-01, 5.00100030010003e-01,
00115   5.00092481273333e-01, 5.00085755939229e-01, 5.00079738458365e-01,
00116   5.00074332862969e-01, 5.00069458915387e-01, 5.00065049112355e-01,
00117   5.00061046334395e-01
00118 };
00119 
00120 static const double V1inv[5 * 5] = {
00121   .47140452079103168293e-1, .37712361663282534635, .56568542494923801952,
00122   .37712361663282534635, .47140452079103168293e-1,
00123   -.81649658092772603273e-1, -.46188021535170061160, 0,
00124   .46188021535170061160, .81649658092772603273e-1, .15058465048420853962,
00125   .12046772038736683169, -.54210474174315074262, .12046772038736683169,
00126   .15058465048420853962, -.21380899352993950775, .30237157840738178177, -0.,
00127   -.30237157840738178177, .21380899352993950775, .10774960475223581324,
00128   -.21549920950447162648, .21549920950447162648, -.21549920950447162648,
00129   .10774960475223581324
00130 };
00131 
00132 static const double V2inv[9 * 9] = {
00133   .11223917161691230546e-1, .10339219839658349826, .19754094204576565761,
00134   .25577315077753587922, .27835314560994251755, .25577315077753587922,
00135   .19754094204576565761, .10339219839658349826, .11223917161691230546e-1,
00136   -.19440394783993476970e-1, -.16544884625069155470, -.24193725566041460608,
00137   -.16953338808305493604, 0.0, .16953338808305493604, .24193725566041460608,
00138   .16544884625069155470, .19440394783993476970e-1, .26466393115406349388e-1,
00139   .17766815796285469394, .11316664642449611462, -.16306601003711325980,
00140   -.30847037493128779631, -.16306601003711325980, .11316664642449611462,
00141   .17766815796285469394, .26466393115406349388e-1,
00142   -.32395302049990834508e-1, -.15521142532414866547,
00143   .88573492664788602740e-1, .29570405784974857322, 0.0,
00144   -.29570405784974857322, -.88573492664788602740e-1, .15521142532414866547,
00145   .32395302049990834508e-1, .41442155673936851246e-1,
00146   .98186757907405608245e-1, -.23056908429499411784,
00147   -.68047008326360625520e-1, .31797435808002456774,
00148   -.68047008326360625520e-1, -.23056908429499411784,
00149   .98186757907405608245e-1, .41442155673936851246e-1,
00150   -.49981120317798783134e-1, -.24861810572835756217e-1,
00151   .23561326072010832539, -.24472785656448415351, 0.0, .24472785656448415351,
00152   -.23561326072010832539, .24861810572835756217e-1,
00153   .49981120317798783134e-1, .79691635865674781228e-1,
00154   -.95725617891693941833e-1, -.57957553356854386344e-1,
00155   .21164072460540271452, -.27529837844505833514, .21164072460540271452,
00156   -.57957553356854386344e-1, -.95725617891693941833e-1,
00157   .79691635865674781228e-1,
00158   -.10894869830716590913, .20131094491947531782, -.15407672674888869038,
00159   .83385723639789791384e-1, 0.0, -.83385723639789791384e-1,
00160   .15407672674888869038, -.20131094491947531782, .10894869830716590913,
00161   .54581057089643838221e-1, -.10916211417928767644, .10916211417928767644,
00162   -.10916211417928767644, .10916211417928767644, -.10916211417928767644,
00163   .10916211417928767644, -.10916211417928767644, .54581057089643838221e-1
00164 };
00165 
00166 static const double V3inv[17 * 17] = {
00167   .27729677693590098996e-2, .26423663180333065153e-1,
00168   .53374068493933898312e-1, .77007854739523195947e-1,
00169   .98257061072911596869e-1, .11538049741786835604, .12832134344120884559,
00170   .13612785914022865001, .13888293186236181317, .13612785914022865001,
00171   .12832134344120884559, .11538049741786835604, .98257061072911596869e-1,
00172   .77007854739523195947e-1, .53374068493933898312e-1,
00173   .26423663180333065153e-1, .27729677693590098996e-2,
00174   -.48029210642807413690e-2, -.44887724635478800254e-1,
00175   -.85409520147301089416e-1, -.11090267822061423050, -.12033983162705862441,
00176   -.11102786862182788886, -.85054870109799336515e-1,
00177   -.45998467987742225160e-1, 0.0, .45998467987742225160e-1,
00178   .85054870109799336515e-1, .11102786862182788886, .12033983162705862441,
00179   .11090267822061423050, .85409520147301089416e-1, .44887724635478800254e-1,
00180   .48029210642807413690e-2, .62758546879582030087e-2,
00181   .55561297093529155869e-1,
00182   .93281491021051539742e-1, .92320151237493695139e-1,
00183   .55077987469605684531e-1,
00184   -.96998141716497488255e-2, -.80285961895427405567e-1,
00185   -.13496839655913850224,
00186   -.15512521776684524331, -.13496839655913850224, -.80285961895427405567e-1,
00187   -.96998141716497488255e-2, .55077987469605684531e-1,
00188   .92320151237493695139e-1, .93281491021051539742e-1,
00189   .55561297093529155869e-1, .62758546879582030087e-2,
00190   -.74850969394858555939e-2, -.61751608943839234096e-1,
00191   -.82974150437304275958e-1, -.38437763431942633378e-1,
00192   .45745502025779701366e-1, .12369235652734542162, .14720439712852868239,
00193   .98768034347019704401e-1, 0.0,
00194   -.98768034347019704401e-1, -.14720439712852868239, -.12369235652734542162,
00195   -.45745502025779701366e-1, .38437763431942633378e-1,
00196   .82974150437304275958e-1, .61751608943839234096e-1,
00197   .74850969394858555939e-2, .86710099994384056338e-2,
00198   .64006230103659573344e-1, .58517426396091675690e-1,
00199   -.29743410528985802680e-1,
00200   -.11934127779157114754, -.12686773515361299409, -.30729137153877447035e-1,
00201   .97307836256600731568e-1, .15635811574451401023, .97307836256600731568e-1,
00202   -.30729137153877447035e-1, -.12686773515361299409, -.11934127779157114754,
00203   -.29743410528985802680e-1, .58517426396091675690e-1,
00204   .64006230103659573344e-1, .86710099994384056338e-2,
00205   -.97486395666294840165e-2, -.62995604908060224672e-1,
00206   -.24373234450275529219e-1, .87760984413626872730e-1,
00207   .12205204576993351394,
00208   .16216004196864002088e-1, -.12422320942156845775, -.13682714580929614678,
00209   0.0, .13682714580929614678, .12422320942156845775,
00210   -.16216004196864002088e-1, -.12205204576993351394,
00211   -.87760984413626872730e-1, .24373234450275529219e-1,
00212   .62995604908060224672e-1, .97486395666294840165e-2,
00213   .10956271233750488468e-1, .58613204255294358939e-1,
00214   -.13306063940736618859e-1, -.11606666444978454399,
00215   -.52059598001115805639e-1, .10868540217796151849, .12594452879014618005,
00216   -.44678658254872910434e-1, -.15617684362128533405,
00217   -.44678658254872910434e-1, .12594452879014618005, .10868540217796151849,
00218   -.52059598001115805639e-1, -.11606666444978454399,
00219   -.13306063940736618859e-1, .58613204255294358939e-1,
00220   .10956271233750488468e-1, -.12098893000863087230e-1,
00221   -.51626244709126208453e-1, .48919433304746979330e-1,
00222   .10467644465949427090,
00223   -.48729879523084673782e-1, -.13668732103524749234, .28190838706814496438e-1,
00224   .15434223333238741600, 0.0, -.15434223333238741600,
00225   -.28190838706814496438e-1, .13668732103524749234,
00226   .48729879523084673782e-1, -.10467644465949427090,
00227   -.48919433304746979330e-1, .51626244709126208453e-1,
00228   .12098893000863087230e-1, .13542668300437944822e-1,
00229   .41712033418258689308e-1,
00230   -.76190463272803434388e-1, -.58303943170068132010e-1, .12158068748245606853,
00231   .42121099930651007882e-1, -.14684425840766337756,
00232   -.16108203535058647043e-1, .15698075850757976092,
00233   -.16108203535058647043e-1, -.14684425840766337756,
00234   .42121099930651007882e-1, .12158068748245606853,
00235   -.58303943170068132010e-1, -.76190463272803434388e-1,
00236   .41712033418258689308e-1, .13542668300437944822e-1,
00237   -.14939634995117694417e-1, -.30047246373341564039e-1,
00238   .91624635082546425678e-1, -.79133374319110026377e-2,
00239   -.12292558212072233355, .90013382617762643524e-1,
00240   .84013717196539593395e-1, -.14813033309980695856, 0.0,
00241   .14813033309980695856, -.84013717196539593395e-1,
00242   -.90013382617762643524e-1,
00243   .12292558212072233355, .79133374319110026377e-2, -.91624635082546425678e-1,
00244   .30047246373341564039e-1, .14939634995117694417e-1,
00245   .16986031342807474208e-1,
00246   .15760203882617033601e-1, -.91494054040950941996e-1,
00247   .70082459207876130806e-1,
00248   .53390713710144539104e-1, -.14340746778352039430, .84048122493418898508e-1,
00249   .72456667788091316868e-1, -.15564535320096811360,
00250   .72456667788091316868e-1, .84048122493418898508e-1,
00251   -.14340746778352039430, .53390713710144539104e-1,
00252   .70082459207876130806e-1, -.91494054040950941996e-1,
00253   .15760203882617033601e-1,
00254   .16986031342807474208e-1, -.18994065631858742028e-1,
00255   -.82901821370405592927e-3, .77239669773015192888e-1,
00256   -.10850735431039424680, .47524484622086496464e-1,
00257   .69148184871588737021e-1, -.14829314646228194928, .11992057742398672066,
00258   0.0, -.11992057742398672066, .14829314646228194928,
00259   -.69148184871588737021e-1, -.47524484622086496464e-1,
00260   .10850735431039424680, -.77239669773015192888e-1,
00261   .82901821370405592927e-3, .18994065631858742028e-1,
00262   .22761703826371535132e-1, -.17728848711449643358e-1,
00263   -.47496371572480503788e-1, .10659958402328690063, -.11696013966166296514,
00264   .63073750910894244526e-1, .32928881123602721303e-1,
00265   -.12280950532497593683, .15926189077282729505, -.12280950532497593683,
00266   .32928881123602721303e-1, .63073750910894244526e-1,
00267   -.11696013966166296514, .10659958402328690063, -.47496371572480503788e-1,
00268   -.17728848711449643358e-1, .22761703826371535132e-1,
00269   -.26493215276042203434e-1, .35579780856128386192e-1,
00270   .10447309718398935122e-1, -.68616154085314996709e-1,
00271   .11775363082763954214, -.13918901977011837274, .12312819418827395690,
00272   -.72053565748259077905e-1, 0.0, .72053565748259077905e-1,
00273   -.12312819418827395690, .13918901977011837274, -.11775363082763954214,
00274   .68616154085314996709e-1, -.10447309718398935122e-1,
00275   -.35579780856128386192e-1,
00276   .26493215276042203434e-1, .40742523354399706918e-1,
00277   -.73124912999529117195e-1, .49317266444153837821e-1,
00278   -.13686605413876015320e-1, -.28342624942191100464e-1,
00279   .70371855298258216249e-1, -.10600251632853603875, .12981016288391131812,
00280   -.13817029659318161476, .12981016288391131812, -.10600251632853603875,
00281   .70371855298258216249e-1, -.28342624942191100464e-1,
00282   -.13686605413876015320e-1,
00283   .49317266444153837821e-1, -.73124912999529117195e-1,
00284   .40742523354399706918e-1, -.54944368958699908688e-1,
00285   .10777725663147408190, -.10152395581538265428, .91369146312596428468e-1,
00286   -.77703071757424700773e-1, .61050911730999815031e-1,
00287   -.42052599404498348871e-1, .21438229266251454773e-1, 0.0,
00288   -.21438229266251454773e-1, .42052599404498348871e-1,
00289   -.61050911730999815031e-1, .77703071757424700773e-1,
00290   -.91369146312596428468e-1,
00291   .10152395581538265428, -.10777725663147408190, .54944368958699908688e-1,
00292   .27485608464748840573e-1, -.54971216929497681146e-1,
00293   .54971216929497681146e-1,
00294   -.54971216929497681146e-1, .54971216929497681146e-1,
00295   -.54971216929497681146e-1, .54971216929497681146e-1,
00296   -.54971216929497681146e-1, .54971216929497681146e-1,
00297   -.54971216929497681146e-1, .54971216929497681146e-1,
00298   -.54971216929497681146e-1, .54971216929497681146e-1,
00299   -.54971216929497681146e-1, .54971216929497681146e-1,
00300   -.54971216929497681146e-1, .27485608464748840573e-1
00301 };
00302 
00303 static const double V4inv[33 * 33] = {
00304   .69120897476690862600e-3, .66419939766331555194e-2,
00305   .13600665164323186111e-1, .20122785860913684493e-1,
00306   .26583214101668429944e-1, .32712713318999268739e-1,
00307   .38576221976287138036e-1, .44033030938268925133e-1,
00308   .49092709529622799673e-1, .53657949874312515646e-1,
00309   .57724533144734311859e-1, .61219564530655179096e-1,
00310   .64138907503837875026e-1, .66427905189318792009e-1,
00311   .68088956652280022887e-1, .69083051391555695878e-1,
00312   .69422738116739271449e-1, .69083051391555695878e-1,
00313   .68088956652280022887e-1, .66427905189318792009e-1,
00314   .64138907503837875026e-1, .61219564530655179096e-1,
00315   .57724533144734311859e-1, .53657949874312515646e-1,
00316   .49092709529622799673e-1, .44033030938268925133e-1,
00317   .38576221976287138036e-1, .32712713318999268739e-1,
00318   .26583214101668429944e-1, .20122785860913684493e-1,
00319   .13600665164323186111e-1, .66419939766331555194e-2,
00320   .69120897476690862600e-3, -.11972090629438798134e-2,
00321   -.11448874821643225573e-1, -.23104401104002905904e-1,
00322   -.33352899418646530133e-1, -.42538626424075425908e-1,
00323   -.49969730733911825941e-1, -.55555454015360728353e-1,
00324   -.58955533624852604918e-1, -.60126044219122513907e-1,
00325   -.58959430451175833624e-1, -.55546925396227130606e-1,
00326   -.49984739749347973762e-1, -.42513009141170294365e-1,
00327   -.33399140950669746346e-1, -.23007690803851790829e-1,
00328   -.11728275717520066169e-1, 0.0, .11728275717520066169e-1,
00329   .23007690803851790829e-1, .33399140950669746346e-1,
00330   .42513009141170294365e-1, .49984739749347973762e-1,
00331   .55546925396227130606e-1, .58959430451175833624e-1,
00332   .60126044219122513907e-1, .58955533624852604918e-1,
00333   .55555454015360728353e-1, .49969730733911825941e-1,
00334   .42538626424075425908e-1, .33352899418646530133e-1,
00335   .23104401104002905904e-1, .11448874821643225573e-1,
00336   .11972090629438798134e-2, .15501585012936019146e-2,
00337   .14628781502199620482e-1, .28684915921474815271e-1,
00338   .39299396074628048026e-1, .46393418975496284204e-1,
00339   .48756902531094699526e-1, .46331333488337494692e-1,
00340   .39012645376980228775e-1, .27452795421085791153e-1,
00341   .12430953621169863781e-1, -.47682978056024928800e-2,
00342   -.22825828045428973853e-1,
00343   -.40195512090720278312e-1, -.55503004262826221955e-1,
00344   -.67424537752827046308e-1, -.75020199300113606452e-1,
00345   -.77607844312483656131e-1, -.75020199300113606452e-1,
00346   -.67424537752827046308e-1, -.55503004262826221955e-1,
00347   -.40195512090720278312e-1, -.22825828045428973853e-1,
00348   -.47682978056024928800e-2, .12430953621169863781e-1,
00349   .27452795421085791153e-1, .39012645376980228775e-1,
00350   .46331333488337494692e-1, .48756902531094699526e-1,
00351   .46393418975496284204e-1, .39299396074628048026e-1,
00352   .28684915921474815271e-1, .14628781502199620482e-1,
00353   .15501585012936019146e-2, -.18377757558949194214e-2,
00354   -.17050470050949761565e-1, -.31952119564923250836e-1,
00355   -.40197423449026348155e-1,
00356   -.41205649520281371624e-1, -.33909965817492272248e-1,
00357   -.19393664422115332144e-1, .56661049630886784692e-3,
00358   .22948272173686561721e-1, .44489719570904738207e-1,
00359   .61790363672287920596e-1, .72121014727028013894e-1,
00360   .73627151185287858579e-1, .65784665375961398923e-1,
00361   .49369676372333667559e-1, .26444326317059715065e-1, 0.0,
00362   -.26444326317059715065e-1, -.49369676372333667559e-1,
00363   -.65784665375961398923e-1, -.73627151185287858579e-1,
00364   -.72121014727028013894e-1, -.61790363672287920596e-1,
00365   -.44489719570904738207e-1, -.22948272173686561721e-1,
00366   -.56661049630886784692e-3, .19393664422115332144e-1,
00367   .33909965817492272248e-1, .41205649520281371624e-1,
00368   .40197423449026348155e-1, .31952119564923250836e-1,
00369   .17050470050949761565e-1, .18377757558949194214e-2,
00370   .20942714740729767769e-2, .18935902405146518232e-1,
00371   .33335840852491735126e-1, .36770680999102286065e-1,
00372   .28873194534132768509e-1, .10267303017729535513e-1,
00373   -.14607738306201572890e-1, -.40139568545572305818e-1,
00374   -.59808326733858291561e-1, -.68528358823372627506e-1,
00375   -.63306535387619244879e-1, -.44508601817574921056e-1,
00376   -.15449116105605395357e-1, .17941083795006546367e-1,
00377   .48747356011657242123e-1, .70329553984201665523e-1,
00378   .78106117292526169663e-1, .70329553984201665523e-1,
00379   .48747356011657242123e-1, .17941083795006546367e-1,
00380   -.15449116105605395357e-1, -.44508601817574921056e-1,
00381   -.63306535387619244879e-1, -.68528358823372627506e-1,
00382   -.59808326733858291561e-1,
00383   -.40139568545572305818e-1, -.14607738306201572890e-1,
00384   .10267303017729535513e-1, .28873194534132768509e-1,
00385   .36770680999102286065e-1, .33335840852491735126e-1,
00386   .18935902405146518232e-1, .20942714740729767769e-2,
00387   -.23245285491878278419e-2, -.20401404737639389919e-1,
00388   -.33019548231022514097e-1, -.29709828426463720091e-1,
00389   -.11760070922697422156e-1, .15987584743850393793e-1,
00390   .43619012891472813485e-1, .61177322409671487721e-1,
00391   .61144030218486655594e-1,
00392   .41895377620089086167e-1, .80232011820644308033e-2,
00393   -.30574701186675900915e-1,
00394   -.62072243008844865848e-1, -.76336186183574765586e-1,
00395   -.68435466095345537115e-1, -.40237669208466966207e-1, 0.0,
00396   .40237669208466966207e-1, .68435466095345537115e-1,
00397   .76336186183574765586e-1, .62072243008844865848e-1,
00398   .30574701186675900915e-1, -.80232011820644308033e-2,
00399   -.41895377620089086167e-1, -.61144030218486655594e-1,
00400   -.61177322409671487721e-1, -.43619012891472813485e-1,
00401   -.15987584743850393793e-1, .11760070922697422156e-1,
00402   .29709828426463720091e-1, .33019548231022514097e-1,
00403   .20401404737639389919e-1, .23245285491878278419e-2,
00404   .25451717261579269307e-2, .21480418595666878775e-1,
00405   .31177212469293007998e-1, .19816333607013379373e-1,
00406   -.72439496274458793681e-2, -.38404203906598342397e-1,
00407   -.57633632255322221046e-1, -.54070547403585392952e-1,
00408   -.26249823354368866005e-1, .15643058212336881516e-1,
00409   .54539832735118677194e-1, .73283028002473989724e-1,
00410   .62835303524135936213e-1, .26175977027801048141e-1,
00411   -.22193636309998606610e-1, -.62597049956093311234e-1,
00412   -.78206986173170212505e-1, -.62597049956093311234e-1,
00413   -.22193636309998606610e-1, .26175977027801048141e-1,
00414   .62835303524135936213e-1,
00415   .73283028002473989724e-1, .54539832735118677194e-1,
00416   .15643058212336881516e-1,
00417   -.26249823354368866005e-1, -.54070547403585392952e-1,
00418   -.57633632255322221046e-1, -.38404203906598342397e-1,
00419   -.72439496274458793681e-2, .19816333607013379373e-1,
00420   .31177212469293007998e-1, .21480418595666878775e-1,
00421   .25451717261579269307e-2, -.27506573922483820005e-2,
00422   -.22224442095099251870e-1, -.27949927254215773020e-1,
00423   -.80918481053370034987e-2, .25121859354449306916e-1,
00424   .51563535009373061074e-1, .51936965107145960512e-1,
00425   .22146626648171527753e-1,
00426   -.24172689882103382748e-1, -.61731229104853568296e-1,
00427   -.68477262429344201201e-1, -.38311232728303704742e-1,
00428   .14160578713659552679e-1, .61248813427564184033e-1,
00429   .77136328841293031805e-1, .52514801765183697988e-1, 0.0,
00430   -.52514801765183697988e-1, -.77136328841293031805e-1,
00431   -.61248813427564184033e-1, -.14160578713659552679e-1,
00432   .38311232728303704742e-1,
00433   .68477262429344201201e-1, .61731229104853568296e-1,
00434   .24172689882103382748e-1,
00435   -.22146626648171527753e-1, -.51936965107145960512e-1,
00436   -.51563535009373061074e-1, -.25121859354449306916e-1,
00437   .80918481053370034987e-2, .27949927254215773020e-1,
00438   .22224442095099251870e-1, .27506573922483820005e-2,
00439   .29562461131654311467e-2, .22630271480554450613e-1,
00440   .23547399831373800971e-1, -.43964593440902476642e-2,
00441   -.39055315767504970597e-1, -.52369643937940066804e-1,
00442   -.28506131614971613422e-1, .19906048093338832322e-1,
00443   .60408880866392420279e-1, .62493397473656883090e-1,
00444   .21391278377641297859e-1, -.37302864786623254746e-1,
00445   -.73665127933539496872e-1, -.61706142476854010202e-1,
00446   -.78065168882546327888e-2, .52335307373945544428e-1,
00447   .78278746279419264777e-1, .52335307373945544428e-1,
00448   -.78065168882546327888e-2, -.61706142476854010202e-1,
00449   -.73665127933539496872e-1, -.37302864786623254746e-1,
00450   .21391278377641297859e-1, .62493397473656883090e-1,
00451   .60408880866392420279e-1, .19906048093338832322e-1,
00452   -.28506131614971613422e-1, -.52369643937940066804e-1,
00453   -.39055315767504970597e-1, -.43964593440902476642e-2,
00454   .23547399831373800971e-1, .22630271480554450613e-1,
00455   .29562461131654311467e-2, -.31515718415504761303e-2,
00456   -.22739451096655080673e-1, -.18157123602272119779e-1,
00457   .16496480897167303621e-1, .46921166788569301124e-1,
00458   .40644395739978416354e-1, -.46275803430732216900e-2,
00459   -.52883375891308909486e-1, -.61116483226324111734e-1,
00460   -.17411698764545629853e-1, .44773430013166822765e-1,
00461   .73441577962383869198e-1, .42127368371995472815e-1,
00462   -.25504645957196772465e-1, -.74126818045972742488e-1,
00463   -.62780077864719287317e-1, 0.0, .62780077864719287317e-1,
00464   .74126818045972742488e-1, .25504645957196772465e-1,
00465   -.42127368371995472815e-1, -.73441577962383869198e-1,
00466   -.44773430013166822765e-1, .17411698764545629853e-1,
00467   .61116483226324111734e-1, .52883375891308909486e-1,
00468   .46275803430732216900e-2, -.40644395739978416354e-1,
00469   -.46921166788569301124e-1, -.16496480897167303621e-1,
00470   .18157123602272119779e-1, .22739451096655080673e-1,
00471   .31515718415504761303e-2, .33536559294882188208e-2,
00472   .22535348942792006185e-1,
00473   .12048629300953560767e-1, -.27166076791299493403e-1,
00474   -.47492745604230978367e-1, -.19246623430993153174e-1,
00475   .36231297307556299322e-1, .61713617181636122004e-1,
00476   .25928029734266134490e-1, -.40478700752883602818e-1,
00477   -.71053889866326412049e-1, -.31870824482961751482e-1,
00478   .41515251100219081281e-1, .76481960760098381651e-1,
00479   .36726509155999912440e-1, -.40090067032627055969e-1,
00480   -.78270742903374539397e-1, -.40090067032627055969e-1,
00481   .36726509155999912440e-1, .76481960760098381651e-1,
00482   .41515251100219081281e-1, -.31870824482961751482e-1,
00483   -.71053889866326412049e-1, -.40478700752883602818e-1,
00484   .25928029734266134490e-1, .61713617181636122004e-1,
00485   .36231297307556299322e-1, -.19246623430993153174e-1,
00486   -.47492745604230978367e-1, -.27166076791299493403e-1,
00487   .12048629300953560767e-1, .22535348942792006185e-1,
00488   .33536559294882188208e-2,
00489   -.35481220456925318865e-2, -.22062913693073191150e-1,
00490   -.54487362861834144999e-2, .35438821865804087489e-1,
00491   .40733077820527411302e-1, -.67403098138950720914e-2,
00492   -.55559584405239171054e-1, -.42417050790865158745e-1,
00493   .24499901971884704925e-1, .68721232891705409302e-1,
00494   .34086082787461126592e-1, -.43441000373118474002e-1,
00495   -.73878085292669148950e-1, -.18846995664706657127e-1,
00496   .59827776178286834498e-1, .70644634584085901794e-1, 0.0,
00497   -.70644634584085901794e-1, -.59827776178286834498e-1,
00498   .18846995664706657127e-1, .73878085292669148950e-1,
00499   .43441000373118474002e-1, -.34086082787461126592e-1,
00500   -.68721232891705409302e-1, -.24499901971884704925e-1,
00501   .42417050790865158745e-1, .55559584405239171054e-1,
00502   .67403098138950720914e-2, -.40733077820527411302e-1,
00503   -.35438821865804087489e-1, .54487362861834144999e-2,
00504   .22062913693073191150e-1, .35481220456925318865e-2,
00505   .37554176816665075631e-2, .21297045781589919482e-1,
00506   -.13327293083183431816e-2,
00507   -.40635299172764596484e-1, -.27659860508374175359e-1,
00508   .31089232744083445986e-1, .56113781541334176109e-1,
00509   .37577840643257763400e-2, -.60511227350664590865e-1,
00510   -.46670556446129053853e-1, .33263195878575888247e-1,
00511   .72757324720645228775e-1, .15011712351692283635e-1,
00512   -.65601212994924119078e-1, -.60016855838843789772e-1,
00513   .26220858553188665966e-1, .78322776605833552980e-1,
00514   .26220858553188665966e-1, -.60016855838843789772e-1,
00515   -.65601212994924119078e-1,
00516   .15011712351692283635e-1, .72757324720645228775e-1,
00517   .33263195878575888247e-1,
00518   -.46670556446129053853e-1, -.60511227350664590865e-1,
00519   .37577840643257763400e-2, .56113781541334176109e-1,
00520   .31089232744083445986e-1, -.27659860508374175359e-1,
00521   -.40635299172764596484e-1, -.13327293083183431816e-2,
00522   .21297045781589919482e-1, .37554176816665075631e-2,
00523   -.39566995305720591229e-2, -.20291873414438919995e-1,
00524   .80617453830770930551e-2, .42270189157016547906e-1,
00525   .10332624526759093004e-1, -.48054759547616142024e-1,
00526   -.37678032941171643972e-1,
00527   .36617192625732482394e-1, .61009425973424865714e-1,
00528   -.95589113168026591466e-2,
00529   -.71023202645076922361e-1, -.25097788086808784456e-1,
00530   .62406621963267050244e-1, .56907293171100693511e-1,
00531   -.36435383083882206257e-1, -.75790105119208756348e-1, 0.0,
00532   .75790105119208756348e-1, .36435383083882206257e-1,
00533   -.56907293171100693511e-1, -.62406621963267050244e-1,
00534   .25097788086808784456e-1, .71023202645076922361e-1,
00535   .95589113168026591466e-2,
00536   -.61009425973424865714e-1, -.36617192625732482394e-1,
00537   .37678032941171643972e-1, .48054759547616142024e-1,
00538   -.10332624526759093004e-1, -.42270189157016547906e-1,
00539   -.80617453830770930551e-2, .20291873414438919995e-1,
00540   .39566995305720591229e-2, .41776092289182138591e-2,
00541   .19013221163904414395e-1, -.14420609729849899876e-1,
00542   -.40259160586844441220e-1, .86327811113710831649e-2,
00543   .53564430703021034399e-1, .65469185402150431933e-2,
00544   -.60383116311280629856e-1,
00545   -.25657793784058876939e-1, .58745680576829226900e-1,
00546   .45649937869034420296e-1,
00547   -.49167932056844167772e-1, -.62696614328552187977e-1,
00548   .32540234556426699997e-1, .74280410383464269758e-1,
00549   -.11425672633410999870e-1, -.78280649404686404903e-1,
00550   -.11425672633410999870e-1, .74280410383464269758e-1,
00551   .32540234556426699997e-1, -.62696614328552187977e-1,
00552   -.49167932056844167772e-1, .45649937869034420296e-1,
00553   .58745680576829226900e-1, -.25657793784058876939e-1,
00554   -.60383116311280629856e-1, .65469185402150431933e-2,
00555   .53564430703021034399e-1,
00556   .86327811113710831649e-2, -.40259160586844441220e-1,
00557   -.14420609729849899876e-1, .19013221163904414395e-1,
00558   .41776092289182138591e-2, -.43935502082478059199e-2,
00559   -.17528761237509401631e-1, .20208915249153872535e-1,
00560   .34734743119040669109e-1, -.26275910172353637955e-1,
00561   -.46368003346018878786e-1,
00562   .26800056330709381025e-1, .56681476464606609921e-1,
00563   -.24749011438127255898e-1,
00564   -.64934612189056658992e-1, .20333742247679279535e-1,
00565   .71429299070059318651e-1,
00566   -.14452513210428671266e-1, -.75793341281736586582e-1,
00567   .74717094137184935270e-2, .78034921554757317374e-1, 0.0,
00568   -.78034921554757317374e-1, -.74717094137184935270e-2,
00569   .75793341281736586582e-1, .14452513210428671266e-1,
00570   -.71429299070059318651e-1, -.20333742247679279535e-1,
00571   .64934612189056658992e-1, .24749011438127255898e-1,
00572   -.56681476464606609921e-1,
00573   -.26800056330709381025e-1, .46368003346018878786e-1,
00574   .26275910172353637955e-1,
00575   -.34734743119040669109e-1, -.20208915249153872535e-1,
00576   .17528761237509401631e-1, .43935502082478059199e-2,
00577   .46379089482818671473e-2, .15791188144791287229e-1,
00578   -.25134290048737455284e-1, -.26249795071946841205e-1,
00579   .39960457575789924651e-1, .28111892450146525404e-1,
00580   -.51026476400767918226e-1,
00581   -.27266747278681831364e-1, .60708796647861610865e-1,
00582   .23532306960642115854e-1,
00583   -.68169639871532441111e-1, -.18204924701958312032e-1,
00584   .73822890510656128485e-1, .11373392486424717019e-1,
00585   -.77133324017644609416e-1, -.39295877480342619961e-2,
00586   .78351902829418987960e-1, -.39295877480342619961e-2,
00587   -.77133324017644609416e-1, .11373392486424717019e-1,
00588   .73822890510656128485e-1, -.18204924701958312032e-1,
00589   -.68169639871532441111e-1, .23532306960642115854e-1,
00590   .60708796647861610865e-1, -.27266747278681831364e-1,
00591   -.51026476400767918226e-1, .28111892450146525404e-1,
00592   .39960457575789924651e-1, -.26249795071946841205e-1,
00593   -.25134290048737455284e-1, .15791188144791287229e-1,
00594   .46379089482818671473e-2, -.48780095920069827068e-2,
00595   -.13886961667516983541e-1, .29071311049368895844e-1,
00596   .15480559452075811600e-1, -.47527977686242313065e-1,
00597   -.31929089844361042178e-2, .58015667638415922967e-1,
00598   -.14547915466597622925e-1, -.61067668299848923244e-1,
00599   .35093678009090186851e-1, .55378399159800654657e-1,
00600   -.54277226474891610385e-1, -.42023830782434076509e-1,
00601   .69197384645944912066e-1, .22610783557709586445e-1,
00602   -.77269275900637030185e-1, 0.0, .77269275900637030185e-1,
00603   -.22610783557709586445e-1,
00604   -.69197384645944912066e-1, .42023830782434076509e-1,
00605   .54277226474891610385e-1,
00606   -.55378399159800654657e-1, -.35093678009090186851e-1,
00607   .61067668299848923244e-1, .14547915466597622925e-1,
00608   -.58015667638415922967e-1, .31929089844361042178e-2,
00609   .47527977686242313065e-1, -.15480559452075811600e-1,
00610   -.29071311049368895844e-1, .13886961667516983541e-1,
00611   .48780095920069827068e-2, .51591759101720291381e-2,
00612   .11747497650231330965e-1, -.31777863364694653331e-1,
00613   -.34555825499804605557e-2, .47914131921157015198e-1,
00614   -.22573685920142225247e-1, -.45320344390022666738e-1,
00615   .49660630547172186418e-1, .25707858143963615736e-1,
00616   -.68132707341917233933e-1, .67534860185243140399e-2,
00617   .69268150370037450063e-1, -.41585011920451477177e-1,
00618   -.51622397460510041271e-1, .68408139576363036148e-1,
00619   .18981259024768933323e-1, -.78265472429342305554e-1,
00620   .18981259024768933323e-1, .68408139576363036148e-1,
00621   -.51622397460510041271e-1,
00622   -.41585011920451477177e-1, .69268150370037450063e-1,
00623   .67534860185243140399e-2,
00624   -.68132707341917233933e-1, .25707858143963615736e-1,
00625   .49660630547172186418e-1,
00626   -.45320344390022666738e-1, -.22573685920142225247e-1,
00627   .47914131921157015198e-1, -.34555825499804605557e-2,
00628   -.31777863364694653331e-1, .11747497650231330965e-1,
00629   .51591759101720291381e-2, -.54365757412741340377e-2,
00630   -.94862516619529080191e-2, .33240472093448190877e-1,
00631   -.88698898099681552229e-2,
00632   -.40973252097216337576e-1, .42995673349795657065e-1,
00633   .17320914507876958783e-1,
00634   -.62201292691914856803e-1, .24726274174637346693e-1,
00635   .51320859246515407288e-1,
00636   -.62882063373810501763e-1, -.11003569131725622672e-1,
00637   .73842261324108943465e-1, -.39240120294802923208e-1,
00638   -.49293966443941122807e-1, .73552644778818223475e-1, 0.0,
00639   -.73552644778818223475e-1, .49293966443941122807e-1,
00640   .39240120294802923208e-1, -.73842261324108943465e-1,
00641   .11003569131725622672e-1, .62882063373810501763e-1,
00642   -.51320859246515407288e-1,
00643   -.24726274174637346693e-1, .62201292691914856803e-1,
00644   -.17320914507876958783e-1, -.42995673349795657065e-1,
00645   .40973252097216337576e-1, .88698898099681552229e-2,
00646   -.33240472093448190877e-1, .94862516619529080191e-2,
00647   .54365757412741340377e-2, .57750194549356126240e-2,
00648   .69981166020044116791e-2, -.33274982140403110792e-1,
00649   .20297071020698356116e-1, .27898517839646066582e-1,
00650   -.53368678853282030262e-1, .16656482990394548343e-1,
00651   .46342901447260614255e-1,
00652   -.60536796508149003365e-1, .29109107483842596340e-2,
00653   .63224486124385124504e-1,
00654   -.59028872851312033411e-1, -.14783105962696191734e-1,
00655   .74269399241069253865e-1, -.49053677339382384625e-1,
00656   -.33525466624811186739e-1, .78397349622515386647e-1,
00657   -.33525466624811186739e-1, -.49053677339382384625e-1,
00658   .74269399241069253865e-1, -.14783105962696191734e-1,
00659   -.59028872851312033411e-1,
00660   .63224486124385124504e-1, .29109107483842596340e-2,
00661   -.60536796508149003365e-1,
00662   .46342901447260614255e-1, .16656482990394548343e-1,
00663   -.53368678853282030262e-1,
00664   .27898517839646066582e-1, .20297071020698356116e-1,
00665   -.33274982140403110792e-1,
00666   .69981166020044116791e-2, .57750194549356126240e-2,
00667   -.61100308370519200637e-2, -.44383614355738148616e-2,
00668   .32011283412619094811e-1, -.29965011866372897633e-1,
00669   -.10560682331349193348e-1, .51110336443392506342e-1,
00670   -.45012284729681775492e-1, -.94236825555873320102e-2,
00671   .60860695783141264746e-1,
00672   -.55014628647083368926e-1, -.73474782382499482121e-2,
00673   .66640148475243034781e-1, -.62533116045749887988e-1,
00674   -.38650525912400102585e-2, .68429769005837003777e-1,
00675   -.66984505412544901945e-1, 0.0, .66984505412544901945e-1,
00676   -.68429769005837003777e-1, .38650525912400102585e-2,
00677   .62533116045749887988e-1, -.66640148475243034781e-1,
00678   .73474782382499482121e-2,
00679   .55014628647083368926e-1, -.60860695783141264746e-1,
00680   .94236825555873320102e-2,
00681   .45012284729681775492e-1, -.51110336443392506342e-1,
00682   .10560682331349193348e-1,
00683   .29965011866372897633e-1, -.32011283412619094811e-1,
00684   .44383614355738148616e-2,
00685   .61100308370519200637e-2, .65409373892036191538e-2,
00686   .16350101107071157065e-2, -.29301957285983144319e-1,
00687   .36838667173388832579e-1, -.81922703976491586393e-2,
00688   -.36955670021050133434e-1, .58374851095540469865e-1,
00689   -.31977016246946181856e-1, -.25311073698658094646e-1,
00690   .66674413950106952577e-1,
00691   -.54865713324521039571e-1, -.39797027891537985440e-2,
00692   .62830285264808449064e-1, -.72226313251296100676e-1,
00693   .22560232697133353980e-1, .46455784709904033738e-1,
00694   -.78200930751070349956e-1, .46455784709904033738e-1,
00695   .22560232697133353980e-1, -.72226313251296100676e-1,
00696   .62830285264808449064e-1, -.39797027891537985440e-2,
00697   -.54865713324521039571e-1, .66674413950106952577e-1,
00698   -.25311073698658094646e-1, -.31977016246946181856e-1,
00699   .58374851095540469865e-1, -.36955670021050133434e-1,
00700   -.81922703976491586393e-2, .36838667173388832579e-1,
00701   -.29301957285983144319e-1, .16350101107071157065e-2,
00702   .65409373892036191538e-2, -.69686180931868703196e-2,
00703   .11849538727632789870e-2, .25452286414610537766e-1,
00704   -.40522480651713943230e-1, .25694679053362813183e-1,
00705   .14057118113748390637e-1, -.52037614725803488893e-1,
00706   .58849342223684035589e-1,
00707   -.25075229077361409271e-1, -.29559771094034181083e-1,
00708   .68296746944165720199e-1, -.62890462146423984955e-1,
00709   .14457636466274596445e-1, .45787612031322361496e-1,
00710   -.77231759014655809742e-1, .57881203613910543657e-1, 0.0,
00711   -.57881203613910543657e-1, .77231759014655809742e-1,
00712   -.45787612031322361496e-1, -.14457636466274596445e-1,
00713   .62890462146423984955e-1,
00714   -.68296746944165720199e-1, .29559771094034181083e-1,
00715   .25075229077361409271e-1,
00716   -.58849342223684035589e-1, .52037614725803488893e-1,
00717   -.14057118113748390637e-1, -.25694679053362813183e-1,
00718   .40522480651713943230e-1, -.25452286414610537766e-1,
00719   -.11849538727632789870e-2, .69686180931868703196e-2,
00720   .75611653617520254845e-2, -.43290610418608409141e-2,
00721   -.20277062025115566914e-1,
00722   .40362947027704828926e-1, -.38938808024132120254e-1,
00723   .11831186195916702262e-1,
00724   .28476667401744525357e-1, -.59320969056617684621e-1,
00725   .61101629747436200186e-1,
00726   -.29514834848355389223e-1, -.20668001885001084821e-1,
00727   .62923592802445122793e-1, -.73558456263588833115e-1,
00728   .45314556330160999776e-1, .79031645918426015574e-2,
00729   -.58136953576334689357e-1, .78538474524006405758e-1,
00730   -.58136953576334689357e-1, .79031645918426015574e-2,
00731   .45314556330160999776e-1, -.73558456263588833115e-1,
00732   .62923592802445122793e-1, -.20668001885001084821e-1,
00733   -.29514834848355389223e-1, .61101629747436200186e-1,
00734   -.59320969056617684621e-1, .28476667401744525357e-1,
00735   .11831186195916702262e-1, -.38938808024132120254e-1,
00736   .40362947027704828926e-1, -.20277062025115566914e-1,
00737   -.43290610418608409141e-2, .75611653617520254845e-2,
00738   -.81505692478987769484e-2, .74297333588288568430e-2,
00739   .14314212513540223314e-1, -.36711242251332751607e-1,
00740   .46240027755503814626e-1, -.34921532671769023773e-1,
00741   .46930051972353714773e-2,
00742   .32842770336385381562e-1, -.61317813706529588466e-1,
00743   .67000809902468893103e-1,
00744   -.45337449655535622885e-1, .35794459576271920867e-2,
00745   .41830061526027213385e-1,
00746   -.72091371931944711708e-1, .74150028530317793195e-1,
00747   -.46487632538609942002e-1, 0.0, .46487632538609942002e-1,
00748   -.74150028530317793195e-1, .72091371931944711708e-1,
00749   -.41830061526027213385e-1, -.35794459576271920867e-2,
00750   .45337449655535622885e-1, -.67000809902468893103e-1,
00751   .61317813706529588466e-1, -.32842770336385381562e-1,
00752   -.46930051972353714773e-2, .34921532671769023773e-1,
00753   -.46240027755503814626e-1, .36711242251332751607e-1,
00754   -.14314212513540223314e-1, -.74297333588288568430e-2,
00755   .81505692478987769484e-2, .90693182942442189743e-2,
00756   -.11121000903959576737e-1, -.71308296141317458546e-2,
00757   .29219439765986671645e-1, -.45820286629778129593e-1,
00758   .49088381175879124421e-1, -.35614888785023038938e-1,
00759   .78906970900092777895e-2,
00760   .26262843038404929480e-1, -.56143674270125757857e-1,
00761   .71700220472378350694e-1,
00762   -.66963544500697307945e-1, .42215091779892228883e-1,
00763   -.41338867413966866997e-2, -.36164891772995367321e-1,
00764   .66584367783847858225e-1, -.77874712365070098328e-1,
00765   .66584367783847858225e-1, -.36164891772995367321e-1,
00766   -.41338867413966866997e-2, .42215091779892228883e-1,
00767   -.66963544500697307945e-1,
00768   .71700220472378350694e-1, -.56143674270125757857e-1,
00769   .26262843038404929480e-1,
00770   .78906970900092777895e-2, -.35614888785023038938e-1,
00771   .49088381175879124421e-1,
00772   -.45820286629778129593e-1, .29219439765986671645e-1,
00773   -.71308296141317458546e-2, -.11121000903959576737e-1,
00774   .90693182942442189743e-2, -.99848472706332791043e-2,
00775   .14701271465939718856e-1, -.32917820356048383366e-3,
00776   -.19201195309873585230e-1, .38409681836626963278e-1,
00777   -.51647324405878909521e-1, .54522171113149311354e-1,
00778   -.45040302741689006270e-1, .24183738595685990149e-1,
00779   .42204134165479735097e-2, -.34317295181348742251e-1,
00780   .59542472465494579941e-1, -.74135115907618101263e-1,
00781   .74491937840566532596e-1, -.60042604725161994304e-1,
00782   .33437677409000083169e-1, 0.0,
00783   -.33437677409000083169e-1, .60042604725161994304e-1,
00784   -.74491937840566532596e-1, .74135115907618101263e-1,
00785   -.59542472465494579941e-1, .34317295181348742251e-1,
00786   -.42204134165479735097e-2, -.24183738595685990149e-1,
00787   .45040302741689006270e-1, -.54522171113149311354e-1,
00788   .51647324405878909521e-1, -.38409681836626963278e-1,
00789   .19201195309873585230e-1, .32917820356048383366e-3,
00790   -.14701271465939718856e-1, .99848472706332791043e-2,
00791   .11775579274769383373e-1, -.19892153937316935880e-1,
00792   .95335114477449041055e-2, .57661528440359081617e-2,
00793   -.23382690532380910781e-1, .40237257037170725321e-1,
00794   -.53280289903551636474e-1, .59974361806023689068e-1,
00795   -.58701684061992853224e-1, .49033407111597129616e-1,
00796   -.31818835267847249219e-1, .90800541261162098886e-2,
00797   .16272906819312603838e-1, -.40863896581186229487e-1,
00798   .61346046297517367703e-1,
00799   -.74896047554167268919e-1, .79632642148310325817e-1,
00800   -.74896047554167268919e-1, .61346046297517367703e-1,
00801   -.40863896581186229487e-1, .16272906819312603838e-1,
00802   .90800541261162098886e-2, -.31818835267847249219e-1,
00803   .49033407111597129616e-1, -.58701684061992853224e-1,
00804   .59974361806023689068e-1, -.53280289903551636474e-1,
00805   .40237257037170725321e-1, -.23382690532380910781e-1,
00806   .57661528440359081617e-2, .95335114477449041055e-2,
00807   -.19892153937316935880e-1,
00808   .11775579274769383373e-1, -.13562702617218467450e-1,
00809   .24885419969649845849e-1, -.18368693901908875583e-1,
00810   .81673147806084084638e-2, .47890591326129587131e-2,
00811   -.19313752945227974024e-1, .34065953398362954708e-1,
00812   -.47667045133463415672e-1, .58820377816690514309e-1,
00813   -.66424139824618415970e-1,
00814   .69667606260856092515e-1, -.68102459384364543253e-1,
00815   .61683024923302547971e-1,
00816   -.50771943476441639136e-1, .36110771847327189215e-1,
00817   -.18758028464284563358e-1, 0.0, .18758028464284563358e-1,
00818   -.36110771847327189215e-1, .50771943476441639136e-1,
00819   -.61683024923302547971e-1, .68102459384364543253e-1,
00820   -.69667606260856092515e-1, .66424139824618415970e-1,
00821   -.58820377816690514309e-1, .47667045133463415672e-1,
00822   -.34065953398362954708e-1, .19313752945227974024e-1,
00823   -.47890591326129587131e-2, -.81673147806084084638e-2,
00824   .18368693901908875583e-1, -.24885419969649845849e-1,
00825   .13562702617218467450e-1, .20576545037980523979e-1,
00826   -.40093155172981004337e-1, .36954083167944054826e-1,
00827   -.31856506837591907746e-1, .24996323181546255126e-1,
00828   -.16637165210473614136e-1, .71002706773325085237e-2,
00829   .32478629093205201133e-2,
00830   -.14009562579050569518e-1, .24771262248780618922e-1,
00831   -.35119395835433647559e-1, .44656290368574753171e-1,
00832   -.53015448339647394161e-1, .59875631995693046782e-1,
00833   -.64973208326045193862e-1, .68112280331082143373e-1,
00834   -.69172215234062186994e-1, .68112280331082143373e-1,
00835   -.64973208326045193862e-1, .59875631995693046782e-1,
00836   -.53015448339647394161e-1, .44656290368574753171e-1,
00837   -.35119395835433647559e-1, .24771262248780618922e-1,
00838   -.14009562579050569518e-1, .32478629093205201133e-2,
00839   .71002706773325085237e-2, -.16637165210473614136e-1,
00840   .24996323181546255126e-1, -.31856506837591907746e-1,
00841   .36954083167944054826e-1, -.40093155172981004337e-1,
00842   .20576545037980523979e-1, -.27584914609096156163e-1,
00843   .54904171411058497973e-1, -.54109756419563083153e-1,
00844   .52794234894345577483e-1, -.50970276026831042415e-1,
00845   .48655445537990983379e-1,
00846   -.45872036510847994332e-1, .42646854695899611372e-1,
00847   -.39010960357087507670e-1, .34999369144476467749e-1,
00848   -.30650714874402762189e-1, .26006877464703437057e-1,
00849   -.21112579608213651273e-1, .16014956068786763273e-1,
00850   -.10763099747751940252e-1, .54075888924374485533e-2, 0.0,
00851   -.54075888924374485533e-2, .10763099747751940252e-1,
00852   -.16014956068786763273e-1,
00853   .21112579608213651273e-1, -.26006877464703437057e-1,
00854   .30650714874402762189e-1,
00855   -.34999369144476467749e-1, .39010960357087507670e-1,
00856   -.42646854695899611372e-1, .45872036510847994332e-1,
00857   -.48655445537990983379e-1, .50970276026831042415e-1,
00858   -.52794234894345577483e-1, .54109756419563083153e-1,
00859   -.54904171411058497973e-1, .27584914609096156163e-1,
00860   .13794141262469565740e-1, -.27588282524939131481e-1,
00861   .27588282524939131481e-1, -.27588282524939131481e-1,
00862   .27588282524939131481e-1, -.27588282524939131481e-1,
00863   .27588282524939131481e-1,
00864   -.27588282524939131481e-1, .27588282524939131481e-1,
00865   -.27588282524939131481e-1, .27588282524939131481e-1,
00866   -.27588282524939131481e-1, .27588282524939131481e-1,
00867   -.27588282524939131481e-1, .27588282524939131481e-1,
00868   -.27588282524939131481e-1, .27588282524939131481e-1,
00869   -.27588282524939131481e-1, .27588282524939131481e-1,
00870   -.27588282524939131481e-1, .27588282524939131481e-1,
00871   -.27588282524939131481e-1, .27588282524939131481e-1,
00872   -.27588282524939131481e-1, .27588282524939131481e-1,
00873   -.27588282524939131481e-1, .27588282524939131481e-1,
00874   -.27588282524939131481e-1, .27588282524939131481e-1,
00875   -.27588282524939131481e-1, .27588282524939131481e-1,
00876   -.27588282524939131481e-1, .13794141262469565740e-1
00877 };
00878 
00879 static const double Tleft[33 * 33] = {
00880   1., -.86602540378443864678, 0., .33071891388307382381, 0.,
00881   -.20728904939721249057, 0., .15128841196122722208, 0.,
00882   -.11918864298744029244, 0., .98352013661686631224e-1, 0.,
00883   -.83727065404940845733e-1, 0., .72893399403505841203e-1, 0.,
00884   -.64544632643375022436e-1, 0., .57913170372415565639e-1, 0.,
00885   -.52518242575729562263e-1, 0., .48043311993977520457e-1, 0.,
00886   -.44271433659733990243e-1, 0., .41048928022856771981e-1, 0.,
00887   -.38263878662008271459e-1, 0., .35832844026365304501e-1, 0., 0.,
00888   .50000000000000000000, -.96824583655185422130, .57282196186948000082,
00889   .21650635094610966169, -.35903516540862679125, -.97578093724974971969e-1,
00890   .26203921611325660506, .55792409597991015609e-1, -.20644078533943456204,
00891   -.36172381205961199479e-1, .17035068468874958194,
00892   .25371838001497225980e-1, -.14501953125000000000,
00893   -.18786835250972344757e-1, .12625507130328301066,
00894   .14473795929590520582e-1, -.11179458309419422675,
00895   -.11494434254897626155e-1, .10030855351241635862,
00896   .93498556820544479096e-2, -.90964264465390582629e-1,
00897   -.77546391824364392762e-2, .83213457337452292745e-1,
00898   .65358085945588638605e-2, -.76680372422574234569e-1,
00899   -.55835321940047427169e-2, .71098828931825789428e-1,
00900   .48253327982967591019e-2, -.66274981937248958553e-1,
00901   -.42118078245337801387e-2, .62064306433355646267e-1,
00902   .37083386598903548973e-2, 0., 0., .25000000000000000000,
00903   -.73950997288745200531, .83852549156242113615, -.23175620272173946716,
00904   -.37791833195149451496, .25710129174850522325, .21608307321780204633,
00905   -.22844049245646009157, -.14009503000335388415, .19897685605518413847,
00906   .98264706042471226893e-1, -.17445445004279014046,
00907   -.72761100054958328401e-1, .15463589893742108388,
00908   .56056770591708784481e-1, -.13855313872640495158,
00909   -.44517752443294564781e-1, .12534277657695128850,
00910   .36211835346039665762e-1, -.11434398255136139683,
00911   -.30033588409423828125e-1, .10506705408753910481,
00912   .25313077840725783008e-1, -.97149327637744872155e-1,
00913   -.21624927200393328444e-1, .90319582367202122625e-1,
00914   .18688433567711780666e-1, -.84372291635345108584e-1,
00915   -.16312261561845420752e-1, .79149526894804751586e-1,
00916   .14362333871852474757e-1, 0., 0., 0., .12500000000000000000,
00917   -.49607837082461073572, .82265291131801144317, -.59621200088559103072,
00918   -.80054302859059362371e-1, .42612156697795759420,
00919   -.90098145270865592887e-1, -.29769623255090078484, .13630307904779758221,
00920   .21638835185708931831, -.14600247270306082052, -.16348801804014290453,
00921   .14340708728599057249, .12755243353979286190, -.13661523715071346961,
00922   -.10215585947881057394, .12864248070157166547, .83592528025348693602e-1,
00923   -.12066728689302565222, -.69633728678718053052e-1, .11314245177331919532,
00924   .58882939251410088028e-1, -.10621835858758221487,
00925   -.50432266865187597572e-1, .99916834723527771581e-1,
00926   .43672094283057258509e-1, -.94206380251950852413e-1,
00927   -.38181356812697746418e-1, .89035739656537771225e-1,
00928   .33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1,
00929   -.31093357409581873586, .67604086414949799246, -.75644205980613611039,
00930   .28990586430124175741, .30648508196770360914, -.35801372616842500052,
00931   -.91326869828709014708e-1, .31127929687500000000,
00932   -.90915752838698393094e-2, -.25637381283965534330,
00933   .57601077850322797594e-1, .21019685709225757945,
00934   -.81244992138514014256e-1, -.17375078516720988858,
00935   .92289437277967051125e-1, .14527351914265391374,
00936   -.96675340792832019889e-1, -.12289485697108543415,
00937   .97448175340011084006e-1, .10511755943298339844,
00938   -.96242247086378239657e-1, -.90822942272780513537e-1,
00939   .93966350452322132384e-1, .79189411876493712558e-1,
00940   -.91139307067989309325e-1, -.69613039934383197265e-1,
00941   .88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0.,
00942   .31250000000000000000e-1, -.18684782411095934408, .50176689760410660236,
00943   -.74784031498626095398, .56472001151566251186, .14842464993721351203e-1,
00944   -.41162920273003120936, .20243071230196532282, .23772054897172750436,
00945   -.24963810923972235950, -.12116179938394678936, .24330535483519110663,
00946   .47903849781124471359e-1, -.22133299683101224293,
00947   -.20542915138527200983e-2, .19653465717678146728,
00948   -.26818172626509178444e-1, -.17319122357631210944,
00949   .45065391411065545445e-1, .15253391395444065941,
00950   -.56543897711725408302e-1, -.13469154928743585367,
00951   .63632471400208840155e-1, .11941684923913523817,
00952   -.67828850207933293098e-1, -.10636309084510652670,
00953   .70095786922999181504e-1, .95187373095150709082e-1, 0., 0., 0., 0., 0.,
00954   0., .15625000000000000000e-1, -.10909562534194485289,
00955   .34842348626527747318, -.64461114561628111443, .69382480527334683659,
00956   -.29551102358528827763, -.25527584713978439819, .38878771718544715394,
00957   -.82956185835347407489e-2, -.31183177761966943912, .12831420840372374767,
00958   .22067618205599434368, -.17569196937129496961, -.14598057000132284135,
00959   .18864406621763419484, .89921002550386645767e-1, -.18571835020187122114,
00960   -.48967672227195481777e-1, .17584685670380332798,
00961   .19267984545067426324e-1, -.16335437520503462738,
00962   .22598055455032407594e-2, .15032800884170631129,
00963   -.17883358353754640871e-1, -.13774837869432209951,
00964   .29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0.,
00965   0., .78125000000000000000e-2, -.62377810244809812496e-1,
00966   .23080781467370883845, -.50841310636012325368, .69834547012574056043,
00967   -.52572723156526459672, .11464215704954976471e-1, .38698869011491210342,
00968   -.26125646622255207507, -.16951698812361607510, .29773875898928782269,
00969   .20130501202570367491e-1, -.26332493149159310198,
00970   .67734613690401207009e-1, .21207315477103762715, -.11541543390889415193,
00971   -.16249634759782417533, .13885887405041735068, .11996491328010275427,
00972   -.14810432001630926895, -.85177658352556243411e-1, .14918860659904380587,
00973   .57317789510444151564e-1, -.14569827645586660151,
00974   -.35213090145965327390e-1, .13975998126844578198, 0., 0., 0., 0., 0., 0.,
00975   0., 0., .39062500000000000000e-2, -.35101954600803571207e-1,
00976   .14761284084133737720, -.37655033076080192966, .62410290231517322776,
00977   -.64335622317683389875, .28188168266139524244, .22488495672137010675,
00978   -.39393811089283576186, .75184777995770096714e-1, .28472023119398293003,
00979   -.20410910833705899572, -.15590046962908511750, .23814567544617953125,
00980   .54442805556829031204e-1, -.22855930338589720954,
00981   .16303223615756629897e-1, .20172722433875559213,
00982   -.62723406421217419404e-1, -.17012230831020922010,
00983   .91754642766136561612e-1, .13927644821381121197, -.10886600968068418181,
00984   -.11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0.,
00985   0., 0., .19531250000000000000e-2, -.19506820659607596598e-1,
00986   .91865676095362231937e-1, -.26604607809696493849, .51425874205091288223,
00987   -.66047561132505329292, .48660109511591303851, -.17575661168678285615e-1,
00988   -.36594333408055703366, .29088854695378694533, .11318677346656537927,
00989   -.31110645235730182168, .60733219161008787341e-1, .24333848233620420826,
00990   -.15254312332655419708, -.15995968483455388613, .19010344455215289289,
00991   .86040636766440260000e-1, -.19652589954665259945,
00992   -.27633388517205837713e-1, .18660848552712880387,
00993   -.15942583868416775867e-1, -.16902042462382064786,
00994   .47278526495327740646e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
00995   .97656250000000000000e-3, -.10731084460857378207e-1,
00996   .55939644713816406331e-1, -.18118487371914493668, .39914857299829864263,
00997   -.60812322949933902435, .60011887183061967583, -.26002695805835928795,
00998   -.20883922404786010096, .38988130966114638081, -.11797833550782589082,
00999   -.25231824756239520077, .24817859972953934712, .90516417677868996417e-1,
01000   -.26079073291293066798, .30259468817169480161e-1, .22178195264114178432,
01001   -.10569877864302048175, -.16679648389266977455, .14637718550245050850,
01002   .11219272032739559870, -.16359363640525750353, -.64358194509092101393e-1,
01003   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3,
01004   -.58542865274813470967e-2, .33461741635290096452e-1,
01005   -.11979993155896201271, .29580223766987206958, -.51874761979436016742,
01006   .62861483498014306968, -.44868895761051453296, .12567502628371529386e-1,
01007   .35040366183235474275, -.30466868455569500886, -.70903913601490112666e-1,
01008   .30822791893032512740, -.11969443264190207736, -.20764760317621313946,
01009   .20629838355452128532, .95269702915334718507e-1, -.22432624768705133300,
01010   -.33103381593477797101e-2, .20570036048155716333,
01011   -.62208282720094518964e-1, -.17095309330441436348, 0., 0., 0., 0., 0., 0.,
01012   0., 0., 0., 0., 0., 0., .24414062500000000000e-3,
01013   -.31714797501871532475e-2, .19721062526127334100e-1,
01014   -.77311181185536498246e-1, .21124871792841566575, -.41777980401893650886,
01015   .59401977834943551650, -.56132417807488349048, .23433675061367565951,
01016   .20222775295220942126, -.38280372496506190127, .14443804214023095767,
01017   .22268950939178466797, -.27211314150777981984, -.34184876506180717313e-1,
01018   .26006498895669734842, -.97650425186005090107e-1, -.19024527660129101293,
01019   .16789164198044635671, .10875811641651905252, -.19276785058805921298, 0.,
01020   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3,
01021   -.17078941137247586143e-2, .11477733754843910060e-1,
01022   -.48887017020924625462e-1, .14634927241421789683, -.32156282683019547854,
01023   .52165811920227223937, -.60001958466396926460, .41208501541480733755,
01024   -.11366945503190350975e-2, -.33968093962672089159, .30955190935923386766,
01025   .40657421856578262210e-1, -.29873400409871531764, .16094481791768257440,
01026   .16876122436206497694, -.23650217045022161255, -.33070260090574765012e-1,
01027   .22985258456375907796, -.68645651043827097771e-1, 0., 0., 0., 0., 0., 0.,
01028   0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4,
01029   -.91501857608428649078e-3, .66085179496951987952e-2,
01030   -.30383171695850355404e-1, .98840838845366876117e-1,
01031   -.23855447246420318989, .43322017468145613917, -.58049033744876107191,
01032   .52533893203742699346, -.20681056202371946180, -.20180000924562504384,
01033   .37503922291962681797, -.15988102869837429062, -.19823558102762374094,
01034   .28393023878803799622, -.11188133439357510403e-1, -.24730368377168229255,
01035   .14731529061377942839, .14878558042884266021, 0., 0., 0., 0., 0., 0., 0.,
01036   0., 0., 0., 0., 0., 0., 0., 0., .30517578125000000000e-4,
01037   -.48804277318479845551e-3, .37696080990601968396e-2,
01038   -.18603912108994738255e-1, .65325006755649582964e-1,
01039   -.17162960707938819795, .34411527956476971322, -.52289350347082497959,
01040   .57319653625674910592, -.37662253421045430413, -.14099055105384663902e-1,
01041   .33265570610216904208, -.30921265572647566661, -.19911390594166455281e-1,
01042   .28738590811031797718, -.18912130469738472647, -.13235936203215819193,
01043   .25076406142356675279, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01044   0., 0., 0., .15258789062500000000e-4, -.25928719280954633249e-3,
01045   .21327398937568540428e-2, -.11244626133630732010e-1,
01046   .42375605740664331966e-1, -.12031130345907846211, .26352562258934426830,
01047   -.44590628258512682078, .56682835613700749379, -.49116715128261660395,
01048   .17845943097110339078, .20541650677432497477, -.36739803642257458221,
01049   .16776034069210108273, .17920950989905112908, -.28867732805385066532,
01050   .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01051   0., 0., 0., 0., 0., .76293945312500000000e-5, -.13727610943181290891e-3,
01052   .11979683091449349286e-2, -.67195313034570709806e-2,
01053   .27044920779931968175e-1, -.82472196498517457862e-1,
01054   .19570475044896150093, -.36391620788543817693, .52241392782736588032,
01055   -.54727504974907879912, .34211551468813581183, .31580472732719957762e-1,
01056   -.32830006549176759667, .30563797665254420769, .64905014620683140120e-2,
01057   -.27642986248995073032, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01058   0., 0., 0., 0., 0., 0., .38146972656250000000e-5,
01059   -.72454147007837596854e-4, .66859847582761390285e-3,
01060   -.39751311980366118437e-2, .17015198650201528366e-1,
01061   -.55443621868993855715e-1, .14157060481641692131, -.28641242619559616836,
01062   .45610665490966615415, -.55262786406029265394, .45818352706035500108,
01063   -.14984403004611673047, -.21163807462970713245, .36007252928843413718,
01064   -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01065   0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5,
01066   -.38135049864067468562e-4, .37101393638555730015e-3,
01067   -.23305339886279723213e-2, .10569913448297127219e-1,
01068   -.36640175162216897547e-1, .10010476414320235508, -.21860074212675559892,
01069   .38124757096345313719, -.52020999209879669177, .52172632730659212045,
01070   -.30841620620308814614, -.50322546186721500184e-1, .32577618885114899053,
01071   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01072   0., 0., .95367431640625000000e-6, -.20021483206955925244e-4,
01073   .20481807322420625431e-3, -.13553476938058909882e-2,
01074   .64919676350791905019e-2, -.23848725425069251903e-1,
01075   .69384632678886421292e-1, -.16249711393618776934, .30736618106830314788,
01076   -.46399909601971539157, .53765031034002467225, -.42598991476520183929,
01077   .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01078   0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6,
01079   -.10487707828484902486e-4, .11254146162337528943e-3,
01080   -.78248929534271987118e-3, .39468337145306794566e-2,
01081   -.15313546659475671763e-1, .47249070825218564146e-1,
01082   -.11804374107101480543, .24031796927792491122, -.39629215049166341285,
01083   .51629108968402548545, -.49622372075429782915, 0., 0., 0., 0., 0., 0., 0.,
01084   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01085   .23841857910156250000e-6, -.54823314130625337326e-5,
01086   .61575377321535518154e-4, -.44877834366497538134e-3,
01087   .23774612048621955857e-2, -.97136347645161687796e-2,
01088   .31671599547606636717e-1, -.84028665767000747480e-1,
01089   .18298487576742964949, -.32647878537696945218, .46970971486488895077, 0.,
01090   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01091   0., 0., 0., 0., .11920928955078125000e-6, -.28604020001177375838e-5,
01092   .33559227978295551013e-4, -.25583821662860610560e-3,
01093   .14201552747787302339e-2, -.60938046986874414969e-2,
01094   .20930869247951926793e-1, -.58745021125678072911e-1,
01095   .13613725780285953720, -.26083988356030237586, 0., 0., 0., 0., 0., 0., 0.,
01096   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01097   .59604644775390625000e-7, -.14898180663526043291e-5,
01098   .18224991282807693921e-4, -.14504433444608833821e-3,
01099   .84184722720281809548e-3, -.37846965430000478789e-2,
01100   .13656355548211376864e-1, -.40409541997718853934e-1,
01101   .99226988101858325902e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01102   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01103   .29802322387695312500e-7, -.77471708843445529468e-6,
01104   .98649879372606876995e-5, -.81814934772838523887e-4,
01105   .49554483992403011328e-3, -.23290922072351413938e-2,
01106   .88068134250844034186e-2, -.27393666952485719070e-1, 0., 0., 0., 0., 0.,
01107   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01108   0., 0., 0., .14901161193847656250e-7, -.40226235946098233685e-6,
01109   .53236418690561306700e-5, -.45933829691164002269e-4,
01110   .28982005232838857913e-3, -.14212974043211018374e-2,
01111   .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01112   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01113   .74505805969238281250e-8, -.20858299254133430408e-6,
01114   .28648457300134381744e-5, -.25677535898258910850e-4,
01115   .16849420429491355445e-3, -.86062824010315834002e-3, 0., 0., 0., 0., 0.,
01116   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01117   0., 0., 0., 0., 0., .37252902984619140625e-8, -.10801736017613096861e-6,
01118   .15376606719887104015e-5, -.14296523739727437959e-4,
01119   .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01120   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01121   .18626451492309570312e-8, -.55871592916438890146e-7,
01122   .82331193828137454068e-6, -.79302250528382787666e-5, 0., 0., 0., 0., 0.,
01123   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01124   0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9,
01125   -.28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0.,
01126   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01127   0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9,
01128   -.14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01129   0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01130   0., 0., .23283064365386962891e-9
01131 };
01132 
01133 static const double Tright[33 * 33] = {
01134   1., .86602540378443864678, 0., -.33071891388307382381, 0.,
01135     .20728904939721249057, 0., -.15128841196122722208, 0.,
01136     .11918864298744029244, 0., -.98352013661686631224e-1, 0.,
01137     .83727065404940845733e-1, 0., -.72893399403505841203e-1, 0.,
01138     .64544632643375022436e-1, 0., -.57913170372415565639e-1, 0.,
01139     .52518242575729562263e-1, 0., -.48043311993977520457e-1, 0.,
01140     .44271433659733990243e-1, 0., -.41048928022856771981e-1, 0.,
01141     .38263878662008271459e-1, 0., -.35832844026365304501e-1, 0., 0.,
01142     .50000000000000000000, .96824583655185422130, .57282196186948000082,
01143     -.21650635094610966169, -.35903516540862679125, .97578093724974971969e-1,
01144     .26203921611325660506, -.55792409597991015609e-1, -.20644078533943456204,
01145     .36172381205961199479e-1, .17035068468874958194,
01146     -.25371838001497225980e-1, -.14501953125000000000,
01147     .18786835250972344757e-1, .12625507130328301066,
01148     -.14473795929590520582e-1, -.11179458309419422675,
01149     .11494434254897626155e-1, .10030855351241635862,
01150     -.93498556820544479096e-2, -.90964264465390582629e-1,
01151     .77546391824364392762e-2, .83213457337452292745e-1,
01152     -.65358085945588638605e-2, -.76680372422574234569e-1,
01153     .55835321940047427169e-2, .71098828931825789428e-1,
01154     -.48253327982967591019e-2, -.66274981937248958553e-1,
01155     .42118078245337801387e-2, .62064306433355646267e-1,
01156     -.37083386598903548973e-2, 0., 0., .25000000000000000000,
01157     .73950997288745200531, .83852549156242113615, .23175620272173946716,
01158     -.37791833195149451496, -.25710129174850522325, .21608307321780204633,
01159     .22844049245646009157, -.14009503000335388415, -.19897685605518413847,
01160     .98264706042471226893e-1, .17445445004279014046,
01161     -.72761100054958328401e-1, -.15463589893742108388,
01162     .56056770591708784481e-1, .13855313872640495158,
01163     -.44517752443294564781e-1, -.12534277657695128850,
01164     .36211835346039665762e-1, .11434398255136139683,
01165     -.30033588409423828125e-1, -.10506705408753910481,
01166     .25313077840725783008e-1, .97149327637744872155e-1,
01167     -.21624927200393328444e-1, -.90319582367202122625e-1,
01168     .18688433567711780666e-1, .84372291635345108584e-1,
01169     -.16312261561845420752e-1, -.79149526894804751586e-1,
01170     .14362333871852474757e-1, 0., 0., 0., .12500000000000000000,
01171     .49607837082461073572, .82265291131801144317, .59621200088559103072,
01172     -.80054302859059362371e-1, -.42612156697795759420,
01173     -.90098145270865592887e-1, .29769623255090078484, .13630307904779758221,
01174     -.21638835185708931831, -.14600247270306082052, .16348801804014290453,
01175     .14340708728599057249, -.12755243353979286190, -.13661523715071346961,
01176     .10215585947881057394, .12864248070157166547, -.83592528025348693602e-1,
01177     -.12066728689302565222, .69633728678718053052e-1, .11314245177331919532,
01178     -.58882939251410088028e-1, -.10621835858758221487,
01179     .50432266865187597572e-1, .99916834723527771581e-1,
01180     -.43672094283057258509e-1, -.94206380251950852413e-1,
01181     .38181356812697746418e-1, .89035739656537771225e-1,
01182     -.33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1,
01183     .31093357409581873586, .67604086414949799246, .75644205980613611039,
01184     .28990586430124175741, -.30648508196770360914, -.35801372616842500052,
01185     .91326869828709014708e-1, .31127929687500000000, .90915752838698393094e-2,
01186     -.25637381283965534330, -.57601077850322797594e-1, .21019685709225757945,
01187     .81244992138514014256e-1, -.17375078516720988858,
01188     -.92289437277967051125e-1, .14527351914265391374,
01189     .96675340792832019889e-1, -.12289485697108543415,
01190     -.97448175340011084006e-1, .10511755943298339844,
01191     .96242247086378239657e-1, -.90822942272780513537e-1,
01192     -.93966350452322132384e-1, .79189411876493712558e-1,
01193     .91139307067989309325e-1, -.69613039934383197265e-1,
01194     -.88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0.,
01195     .31250000000000000000e-1, .18684782411095934408, .50176689760410660236,
01196     .74784031498626095398, .56472001151566251186, -.14842464993721351203e-1,
01197     -.41162920273003120936, -.20243071230196532282, .23772054897172750436,
01198     .24963810923972235950, -.12116179938394678936, -.24330535483519110663,
01199     .47903849781124471359e-1, .22133299683101224293,
01200     -.20542915138527200983e-2, -.19653465717678146728,
01201     -.26818172626509178444e-1, .17319122357631210944,
01202     .45065391411065545445e-1, -.15253391395444065941,
01203     -.56543897711725408302e-1, .13469154928743585367,
01204     .63632471400208840155e-1, -.11941684923913523817,
01205     -.67828850207933293098e-1, .10636309084510652670,
01206     .70095786922999181504e-1, -.95187373095150709082e-1, 0., 0., 0., 0., 0.,
01207     0., .15625000000000000000e-1, .10909562534194485289,
01208     .34842348626527747318, .64461114561628111443, .69382480527334683659,
01209     .29551102358528827763, -.25527584713978439819, -.38878771718544715394,
01210     -.82956185835347407489e-2, .31183177761966943912, .12831420840372374767,
01211     -.22067618205599434368, -.17569196937129496961, .14598057000132284135,
01212     .18864406621763419484, -.89921002550386645767e-1, -.18571835020187122114,
01213     .48967672227195481777e-1, .17584685670380332798,
01214     -.19267984545067426324e-1, -.16335437520503462738,
01215     -.22598055455032407594e-2, .15032800884170631129,
01216     .17883358353754640871e-1, -.13774837869432209951,
01217     -.29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0.,
01218     0., .78125000000000000000e-2, .62377810244809812496e-1,
01219     .23080781467370883845, .50841310636012325368, .69834547012574056043,
01220     .52572723156526459672, .11464215704954976471e-1, -.38698869011491210342,
01221     -.26125646622255207507, .16951698812361607510, .29773875898928782269,
01222     -.20130501202570367491e-1, -.26332493149159310198,
01223     -.67734613690401207009e-1, .21207315477103762715, .11541543390889415193,
01224     -.16249634759782417533, -.13885887405041735068, .11996491328010275427,
01225     .14810432001630926895, -.85177658352556243411e-1, -.14918860659904380587,
01226     .57317789510444151564e-1, .14569827645586660151,
01227     -.35213090145965327390e-1, -.13975998126844578198, 0., 0., 0., 0., 0., 0.,
01228     0., 0., .39062500000000000000e-2, .35101954600803571207e-1,
01229     .14761284084133737720, .37655033076080192966, .62410290231517322776,
01230     .64335622317683389875, .28188168266139524244, -.22488495672137010675,
01231     -.39393811089283576186, -.75184777995770096714e-1, .28472023119398293003,
01232     .20410910833705899572, -.15590046962908511750, -.23814567544617953125,
01233     .54442805556829031204e-1, .22855930338589720954, .16303223615756629897e-1,
01234     -.20172722433875559213, -.62723406421217419404e-1, .17012230831020922010,
01235     .91754642766136561612e-1, -.13927644821381121197, -.10886600968068418181,
01236     .11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0.,
01237     0., 0., .19531250000000000000e-2, .19506820659607596598e-1,
01238     .91865676095362231937e-1, .26604607809696493849, .51425874205091288223,
01239     .66047561132505329292, .48660109511591303851, .17575661168678285615e-1,
01240     -.36594333408055703366, -.29088854695378694533, .11318677346656537927,
01241     .31110645235730182168, .60733219161008787341e-1, -.24333848233620420826,
01242     -.15254312332655419708, .15995968483455388613, .19010344455215289289,
01243     -.86040636766440260000e-1, -.19652589954665259945,
01244     .27633388517205837713e-1, .18660848552712880387, .15942583868416775867e-1,
01245     -.16902042462382064786, -.47278526495327740646e-1, 0., 0., 0., 0., 0., 0.,
01246     0., 0., 0., 0., .97656250000000000000e-3, .10731084460857378207e-1,
01247     .55939644713816406331e-1, .18118487371914493668, .39914857299829864263,
01248     .60812322949933902435, .60011887183061967583, .26002695805835928795,
01249     -.20883922404786010096, -.38988130966114638081, -.11797833550782589082,
01250     .25231824756239520077, .24817859972953934712, -.90516417677868996417e-1,
01251     -.26079073291293066798, -.30259468817169480161e-1, .22178195264114178432,
01252     .10569877864302048175, -.16679648389266977455, -.14637718550245050850,
01253     .11219272032739559870, .16359363640525750353, -.64358194509092101393e-1,
01254     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3,
01255     .58542865274813470967e-2, .33461741635290096452e-1, .11979993155896201271,
01256     .29580223766987206958, .51874761979436016742, .62861483498014306968,
01257     .44868895761051453296, .12567502628371529386e-1, -.35040366183235474275,
01258     -.30466868455569500886, .70903913601490112666e-1, .30822791893032512740,
01259     .11969443264190207736, -.20764760317621313946, -.20629838355452128532,
01260     .95269702915334718507e-1, .22432624768705133300,
01261     -.33103381593477797101e-2, -.20570036048155716333,
01262     -.62208282720094518964e-1, .17095309330441436348, 0., 0., 0., 0., 0., 0.,
01263     0., 0., 0., 0., 0., 0., .24414062500000000000e-3,
01264     .31714797501871532475e-2, .19721062526127334100e-1,
01265     .77311181185536498246e-1, .21124871792841566575, .41777980401893650886,
01266     .59401977834943551650, .56132417807488349048, .23433675061367565951,
01267     -.20222775295220942126, -.38280372496506190127, -.14443804214023095767,
01268     .22268950939178466797, .27211314150777981984, -.34184876506180717313e-1,
01269     -.26006498895669734842, -.97650425186005090107e-1, .19024527660129101293,
01270     .16789164198044635671, -.10875811641651905252, -.19276785058805921298, 0.,
01271     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3,
01272     .17078941137247586143e-2, .11477733754843910060e-1,
01273     .48887017020924625462e-1, .14634927241421789683, .32156282683019547854,
01274     .52165811920227223937, .60001958466396926460, .41208501541480733755,
01275     .11366945503190350975e-2, -.33968093962672089159, -.30955190935923386766,
01276     .40657421856578262210e-1, .29873400409871531764, .16094481791768257440,
01277     -.16876122436206497694, -.23650217045022161255, .33070260090574765012e-1,
01278     .22985258456375907796, .68645651043827097771e-1, 0., 0., 0., 0., 0., 0.,
01279     0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4,
01280     .91501857608428649078e-3, .66085179496951987952e-2,
01281     .30383171695850355404e-1, .98840838845366876117e-1, .23855447246420318989,
01282     .43322017468145613917, .58049033744876107191, .52533893203742699346,
01283     .20681056202371946180, -.20180000924562504384, -.37503922291962681797,
01284     -.15988102869837429062, .19823558102762374094, .28393023878803799622,
01285     .11188133439357510403e-1, -.24730368377168229255, -.14731529061377942839,
01286     .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01287     0., 0., .30517578125000000000e-4, .48804277318479845551e-3,
01288     .37696080990601968396e-2, .18603912108994738255e-1,
01289     .65325006755649582964e-1, .17162960707938819795, .34411527956476971322,
01290     .52289350347082497959, .57319653625674910592, .37662253421045430413,
01291     -.14099055105384663902e-1, -.33265570610216904208, -.30921265572647566661,
01292     .19911390594166455281e-1, .28738590811031797718, .18912130469738472647,
01293     -.13235936203215819193, -.25076406142356675279, 0., 0., 0., 0., 0., 0.,
01294     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .15258789062500000000e-4,
01295     .25928719280954633249e-3, .21327398937568540428e-2,
01296     .11244626133630732010e-1, .42375605740664331966e-1, .12031130345907846211,
01297     .26352562258934426830, .44590628258512682078, .56682835613700749379,
01298     .49116715128261660395, .17845943097110339078, -.20541650677432497477,
01299     -.36739803642257458221, -.16776034069210108273, .17920950989905112908,
01300     .28867732805385066532, .46473465543376206337e-1, 0., 0., 0., 0., 0., 0.,
01301     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .76293945312500000000e-5,
01302     .13727610943181290891e-3, .11979683091449349286e-2,
01303     .67195313034570709806e-2, .27044920779931968175e-1,
01304     .82472196498517457862e-1, .19570475044896150093, .36391620788543817693,
01305     .52241392782736588032, .54727504974907879912, .34211551468813581183,
01306     -.31580472732719957762e-1, -.32830006549176759667, -.30563797665254420769,
01307     .64905014620683140120e-2, .27642986248995073032, 0., 0., 0., 0., 0., 0.,
01308     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .38146972656250000000e-5,
01309     .72454147007837596854e-4, .66859847582761390285e-3,
01310     .39751311980366118437e-2, .17015198650201528366e-1,
01311     .55443621868993855715e-1, .14157060481641692131, .28641242619559616836,
01312     .45610665490966615415, .55262786406029265394, .45818352706035500108,
01313     .14984403004611673047, -.21163807462970713245, -.36007252928843413718,
01314     -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01315     0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5,
01316     .38135049864067468562e-4, .37101393638555730015e-3,
01317     .23305339886279723213e-2, .10569913448297127219e-1,
01318     .36640175162216897547e-1, .10010476414320235508, .21860074212675559892,
01319     .38124757096345313719, .52020999209879669177, .52172632730659212045,
01320     .30841620620308814614, -.50322546186721500184e-1, -.32577618885114899053,
01321     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01322     0., 0., .95367431640625000000e-6, .20021483206955925244e-4,
01323     .20481807322420625431e-3, .13553476938058909882e-2,
01324     .64919676350791905019e-2, .23848725425069251903e-1,
01325     .69384632678886421292e-1, .16249711393618776934, .30736618106830314788,
01326     .46399909601971539157, .53765031034002467225, .42598991476520183929,
01327     .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01328     0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6,
01329     .10487707828484902486e-4, .11254146162337528943e-3,
01330     .78248929534271987118e-3, .39468337145306794566e-2,
01331     .15313546659475671763e-1, .47249070825218564146e-1, .11804374107101480543,
01332     .24031796927792491122, .39629215049166341285, .51629108968402548545,
01333     .49622372075429782915, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01334     0., 0., 0., 0., 0., 0., 0., 0., 0., .23841857910156250000e-6,
01335     .54823314130625337326e-5, .61575377321535518154e-4,
01336     .44877834366497538134e-3, .23774612048621955857e-2,
01337     .97136347645161687796e-2, .31671599547606636717e-1,
01338     .84028665767000747480e-1, .18298487576742964949, .32647878537696945218,
01339     .46970971486488895077, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01340     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .11920928955078125000e-6,
01341     .28604020001177375838e-5, .33559227978295551013e-4,
01342     .25583821662860610560e-3, .14201552747787302339e-2,
01343     .60938046986874414969e-2, .20930869247951926793e-1,
01344     .58745021125678072911e-1, .13613725780285953720, .26083988356030237586,
01345     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01346     0., 0., 0., 0., 0., 0., .59604644775390625000e-7,
01347     .14898180663526043291e-5, .18224991282807693921e-4,
01348     .14504433444608833821e-3, .84184722720281809548e-3,
01349     .37846965430000478789e-2, .13656355548211376864e-1,
01350     .40409541997718853934e-1, .99226988101858325902e-1, 0., 0., 0., 0., 0.,
01351     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01352     0., 0., .29802322387695312500e-7, .77471708843445529468e-6,
01353     .98649879372606876995e-5, .81814934772838523887e-4,
01354     .49554483992403011328e-3, .23290922072351413938e-2,
01355     .88068134250844034186e-2, .27393666952485719070e-1, 0., 0., 0., 0., 0.,
01356     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01357     0., 0., 0., .14901161193847656250e-7, .40226235946098233685e-6,
01358     .53236418690561306700e-5, .45933829691164002269e-4,
01359     .28982005232838857913e-3, .14212974043211018374e-2,
01360     .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01361     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01362     .74505805969238281250e-8, .20858299254133430408e-6,
01363     .28648457300134381744e-5, .25677535898258910850e-4,
01364     .16849420429491355445e-3, .86062824010315834002e-3, 0., 0., 0., 0., 0.,
01365     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01366     0., 0., 0., 0., 0., .37252902984619140625e-8, .10801736017613096861e-6,
01367     .15376606719887104015e-5, .14296523739727437959e-4,
01368     .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01369     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01370     .18626451492309570312e-8, .55871592916438890146e-7,
01371     .82331193828137454068e-6, .79302250528382787666e-5, 0., 0., 0., 0., 0.,
01372     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01373     0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9,
01374     .28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0.,
01375     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01376     0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9,
01377     .14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01378     0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
01379     0., 0., .23283064365386962891e-9
01380 };
01381 
01382 /* Allocates a workspace for the given maximum number of intervals.
01383     Note that if the workspace gets filled, the intervals with the
01384     lowest error estimates are dropped. The maximum number of
01385     intervals is therefore not the maximum number of intervals
01386     that will be computed, but merely the size of the buffer.
01387     */
01388 
01389 /* Compute the product of the fx with one of the inverse
01390     Vandermonde-like matrices. */
01391 
01392 void
01393 Vinvfx (const double *fx, double *c, const int d)
01394 {
01395 
01396   int i, j;
01397 
01398   switch (d)
01399     {
01400     case 0:
01401       for (i = 0; i <= 4; i++)
01402         {
01403           c[i] = 0.0;
01404           for (j = 0; j <= 4; j++)
01405             c[i] += V1inv[i * 5 + j] * fx[j * 8];
01406         }
01407       break;
01408     case 1:
01409       for (i = 0; i <= 8; i++)
01410         {
01411           c[i] = 0.0;
01412           for (j = 0; j <= 8; j++)
01413             c[i] += V2inv[i * 9 + j] * fx[j * 4];
01414         }
01415       break;
01416     case 2:
01417       for (i = 0; i <= 16; i++)
01418         {
01419           c[i] = 0.0;
01420           for (j = 0; j <= 16; j++)
01421             c[i] += V3inv[i * 17 + j] * fx[j * 2];
01422         }
01423       break;
01424     case 3:
01425       for (i = 0; i <= 32; i++)
01426         {
01427           c[i] = 0.0;
01428           for (j = 0; j <= 32; j++)
01429             c[i] += V4inv[i * 33 + j] * fx[j];
01430         }
01431       break;
01432     }
01433 
01434 }
01435 
01436 
01437 /* Downdate the interpolation given by the n coefficients c
01438     by removing the nodes with indices in nans. */
01439 
01440 void
01441 downdate (double *c, int n, int d, int *nans, int nnans)
01442 {
01443 
01444   static const int bidx[4] = { 0, 6, 16, 34 };
01445   double b_new[34], alpha;
01446   int i, j;
01447 
01448   for (i = 0; i <= n + 1; i++)
01449     b_new[i] = bee[bidx[d] + i];
01450   for (i = 0; i < nnans; i++)
01451     {
01452       b_new[n + 1] = b_new[n + 1] / Lalpha[n];
01453       b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1];
01454       for (j = n - 1; j > 0; j--)
01455         b_new[j] =
01456           (b_new[j] + xi[nans[i]] * b_new[j + 1] -
01457            Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1];
01458       for (j = 0; j <= n; j++)
01459         b_new[j] = b_new[j + 1];
01460       alpha = c[n] / b_new[n];
01461       for (j = 0; j < n; j++)
01462         c[j] -= alpha * b_new[j];
01463       c[n] = 0;
01464       n--;
01465     }
01466 
01467 }
01468 
01469 
01470 /* The actual integration routine.  */
01471 
01472 DEFUN_DLD (quadcc, args, nargout,
01473 "-*- texinfo -*-\n\
01474 @deftypefn  {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b})\n\
01475 @deftypefnx {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})\n\
01476 @deftypefnx {Function File} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\
01477 @deftypefnx {Function File} {[@var{q}, @var{err}, @var{nr_points}] =} quadcc (@dots{})\n\
01478 Numerically evaluate the integral of @var{f} from @var{a} to @var{b}\n\
01479 using the doubly-adaptive Clenshaw-Curtis quadrature described by P. Gonnet\n\
01480 in @cite{Increasing the Reliability of Adaptive Quadrature Using Explicit\n\
01481 Interpolants}.\n\
01482 @var{f} is a function handle, inline function, or string\n\
01483 containing the name of the function to evaluate.\n\
01484 The function @var{f} must be vectorized and must return a vector of output\n\
01485 values if given a vector of input values.  For example,\n\
01486 \n\
01487 @example\n\
01488    f = @@(x) x .* sin (1./x) .* sqrt (abs (1 - x));\n\
01489 @end example\n\
01490 \n\
01491 @noindent\n\
01492 which uses the element-by-element `dot' form for all operators.\n\
01493 \n\
01494 @var{a} and @var{b} are the lower and upper limits of integration.  Either\n\
01495 or both limits may be infinite.  @code{quadcc} handles an inifinite limit\n\
01496 by substituting the variable of integration with @code{x=tan(pi/2*u)}.\n\
01497 \n\
01498 The optional argument @var{tol} defines the relative tolerance used to stop\n\
01499 the integration procedure.  The default value is @math{1e^{-6}}.\n\
01500 \n\
01501 The optional argument @var{sing} contains a list of points where the\n\
01502 integrand has known singularities, or discontinuities\n\
01503 in any of its derivatives, inside the integration interval.\n\
01504 For the example above, which has a discontinuity at x=1, the call to\n\
01505 @code{quadcc} would be as follows\n\
01506 \n\
01507 @example\n\
01508    int = quadcc (f, a, b, 1.0e-6, [ 1 ]);\n\
01509 @end example\n\
01510 \n\
01511 The result of the integration is returned in @var{q}.\n\
01512 @var{err} is an estimate of the absolute integration error and\n\
01513 @var{nr_points} is the number of points at which the integrand was evaluated.\n\
01514 If the adaptive integration did not converge, the value of\n\
01515 @var{err} will be larger than the requested tolerance.  Therefore, it is\n\
01516 recommended to verify this value for difficult integrands.\n\
01517 \n\
01518 @code{quadcc} is capable of dealing with non-numeric\n\
01519 values of the integrand such as @code{NaN} or @code{Inf}.\n\
01520 If the integral diverges, and @code{quadcc} detects this,\n\
01521 then a warning is issued and @code{Inf} or @code{-Inf} is returned.\n\
01522 \n\
01523 Note: @code{quadcc} is a general purpose quadrature algorithm\n\
01524 and, as such, may be less efficient for a smooth or otherwise\n\
01525 well-behaved integrand than other methods such as @code{quadgk}.\n\
01526 \n\
01527 The algorithm uses Clenshaw-Curtis quadrature rules of increasing\n\
01528 degree in each interval and bisects the interval if either the\n\
01529 function does not appear to be smooth or a rule of maximum\n\
01530 degree has been reached.  The error estimate is computed from the\n\
01531 L2-norm of the difference between two successive interpolations\n\
01532 of the integrand over the nodes of the respective quadrature rules.\n\
01533 \n\
01534 Reference: P. Gonnet, @cite{Increasing the Reliability of Adaptive\n\
01535 Quadrature Using Explicit Interpolants}, ACM Transactions on\n\
01536 Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010.\n\
01537 @seealso{quad, quadv, quadl, quadgk, trapz, dblquad, triplequad}\n\
01538 @end deftypefn")
01539 {
01540   octave_value_list retval;
01541 
01542   /* Some constants that we will need. */
01543   static const int n[4] = { 4, 8, 16, 32 };
01544   static const int skip[4] = { 8, 4, 2, 1 };
01545   static const int idx[4] = { 0, 5, 14, 31 };
01546   static const double w = M_SQRT2 / 2;
01547   static const int ndiv_max = 20;
01548 
01549   /* The interval heap. */
01550   cquad_ival ivals[cquad_heapsize];
01551   int heap[cquad_heapsize];
01552 
01553   /* Arguments left and right */
01554   int nargin = args.length ();
01555   octave_function *fcn;
01556   double a, b, tol, iivals[cquad_heapsize], *sing;
01557 
01558   /* Variables needed for transforming the integrand. */
01559   bool wrap = false;
01560   double xw;
01561 
01562   /* Stuff we will need to call the integrand. */
01563   octave_value_list fargs, fvals;
01564 
01565   /* Actual variables (as opposed to constants above). */
01566   double m, h, ml, hl, mr, hr, temp;
01567   double igral, err, igral_final, err_final;
01568   int nivals, neval = 0;
01569   int i, j, d, split, t;
01570   int nnans, nans[33];
01571   cquad_ival *iv, *ivl, *ivr;
01572   double nc, ncdiff;
01573 
01574 
01575   /* Parse the input arguments. */
01576   if (nargin < 3)
01577     {
01578       print_usage ();
01579       return retval;
01580     }
01581 
01582   if (args(0).is_function_handle () || args(0).is_inline_function ())
01583     fcn = args(0).function_value ();
01584   else
01585     {
01586        std::string fcn_name = unique_symbol_name ("__quadcc_fcn_");
01587        std::string fname = "function y = ";
01588        fname.append (fcn_name);
01589        fname.append ("(x) y = ");
01590        fcn = extract_function (args(0), "quadcc", fcn_name, fname,
01591                                "; endfunction");
01592     }
01593 
01594   if (!args(1).is_real_scalar ())
01595     {
01596       error ("quadcc: lower limit of integration (A) must be a single real scalar");
01597       return retval;
01598     }
01599   else
01600     a = args(1).double_value ();
01601 
01602   if (!args(2).is_real_scalar ())
01603     {
01604       error ("quadcc: upper limit of integration (B) must be a single real scalar");
01605       return retval;
01606     }
01607   else
01608     b = args(2).double_value ();
01609 
01610   if (nargin < 4 || args(3).is_empty ())
01611     tol = 1.0e-6;
01612   else if (!args(3).is_real_scalar () || args(3).double_value () <= 0)
01613     {
01614       error ("quadcc: tolerance (TOL) must be a single real scalar > 0");
01615       return retval;
01616     }
01617   else
01618     tol = args(3).double_value ();
01619 
01620   if (nargin < 5)
01621     {
01622       nivals = 1;
01623       iivals[0] = a;
01624       iivals[1] = b;
01625     }
01626   else if (!(args(4).is_real_scalar () || args(4).is_real_matrix ()))
01627     {
01628       error ("quadcc: list of singularities (SING) must be a vector of real values");
01629       return retval;
01630     }
01631   else
01632     {
01633       nivals = 1 + args(4).length ();
01634       if (nivals > cquad_heapsize)
01635         {
01636           error ("quadcc: maximum number of singular points is limited to %i",
01637                  cquad_heapsize-1);
01638           return retval;
01639         }
01640       sing = args(4).array_value ().fortran_vec ();
01641       iivals[0] = a;
01642       for (i = 0; i < nivals - 2; i++)
01643         iivals[i + 1] = sing[i];
01644       iivals[nivals] = b;
01645     }
01646 
01647   /* If a or b are +/-Inf, transform the integral. */
01648   if (xisinf (a) || xisinf (b))
01649     {
01650       wrap = true;
01651       for (i = 0; i <= nivals; i++)
01652         if (xisinf (iivals[i]))
01653           iivals[i] = gnulib::copysign (1.0, iivals[i]);
01654         else
01655           iivals[i] = 2.0 * atan (iivals[i]) / M_PI;
01656     }
01657 
01658 
01659   /* Initialize the heaps. */
01660   for (i = 0; i < cquad_heapsize; i++)
01661     heap[i] = i;
01662 
01663 
01664   /* Create the first interval(s). */
01665   igral = 0.0;
01666   err = 0.0;
01667   for (j = 0; j < nivals; j++)
01668     {
01669 
01670       /* Initialize the interval. */
01671       iv = &(ivals[heap[j]]);
01672       m = (iivals[j] + iivals[j + 1]) / 2;
01673       h = (iivals[j + 1] - iivals[j]) / 2;
01674       nnans = 0;
01675       ColumnVector ex (33);
01676       if (wrap)
01677         {
01678           for (i = 0; i <= n[3]; i++)
01679             ex (i) = tan (M_PI / 2 * (m + xi[i] * h));
01680         }
01681       else
01682         {
01683           for (i = 0; i <= n[3]; i++)
01684             ex (i) = m + xi[i] * h;
01685         }
01686       fargs(0) = ex;
01687       fvals = feval (fcn, fargs, 1);
01688       if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
01689         {
01690           error ("quadcc: integrand F must return a single, real-valued vector");
01691           return retval;
01692         }
01693       Matrix effex = fvals(0).matrix_value ();
01694       if (effex.length () != ex.length ())
01695         {
01696           error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
01697           return retval;
01698         }
01699       for (i = 0; i <= n[3]; i++)
01700         {
01701           iv->fx[i] = effex (i);
01702           if (wrap)
01703             {
01704               xw = ex(i);
01705               iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2;
01706             }
01707           neval++;
01708           if (!xfinite (iv->fx[i]))
01709             {
01710               nans[nnans++] = i;
01711               iv->fx[i] = 0.0;
01712             }
01713         }
01714       Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
01715       Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
01716       Vinvfx (iv->fx, &(iv->c[0]), 0);
01717       for (i = 0; i < nnans; i++)
01718         iv->fx[nans[i]] = octave_NaN;
01719       iv->a = iivals[j];
01720       iv->b = iivals[j + 1];
01721       iv->depth = 3;
01722       iv->rdepth = 1;
01723       iv->ndiv = 0;
01724       iv->igral = 2 * h * iv->c[idx[3]] * w;
01725       nc = 0.0;
01726       for (i = n[2] + 1; i <= n[3]; i++)
01727         {
01728           temp = iv->c[idx[3] + i];
01729           nc += temp * temp;
01730         }
01731       ncdiff = nc;
01732       for (i = 0; i <= n[2]; i++)
01733         {
01734           temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
01735           ncdiff += temp * temp;
01736           nc += iv->c[idx[3] + i] * iv->c[idx[3] + i];
01737         }
01738       ncdiff = sqrt (ncdiff);
01739       nc = sqrt (nc);
01740       iv->err = ncdiff * 2 * h;
01741       if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
01742         iv->err = 2 * h * nc;
01743 
01744       /* Tabulate this interval's data. */
01745       igral += iv->igral;
01746       err += iv->err;
01747 
01748       /* Sift it up the heap. */
01749       i = j;
01750       while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err)
01751         {
01752           temp = heap[i];
01753           heap[i] = heap[i / 2];
01754           heap[i / 2] = temp;
01755           i /= 2;
01756         }
01757 
01758     }
01759 
01760 
01761   /* Initialize some global values. */
01762   igral_final = 0.0;
01763   err_final = 0.0;
01764 
01765 
01766   /* Main loop. */
01767   while (nivals > 0 && err > 0.0 && err > fabs (igral) * tol
01768          && !(err_final > fabs (igral) * tol
01769               && err - err_final < fabs (igral) * tol))
01770     {
01771 
01772       /* Allow the user to interrupt. */
01773       OCTAVE_QUIT;
01774 
01775       /* Put our finger on the interval with the largest error. */
01776       iv = &(ivals[heap[0]]);
01777       m = (iv->a + iv->b) / 2;
01778       h = (iv->b - iv->a) / 2;
01779 
01780 /*      printf
01781         ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
01782          heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
01783 */
01784       /* Should we try to increase the degree? */
01785       if (iv->depth < 3)
01786         {
01787 
01788           /* Keep tabs on some variables. */
01789           d = ++iv->depth;
01790 
01791           /* Get the new (missing) function values */
01792           {
01793             ColumnVector ex (n[d] / 2);
01794             if (wrap)
01795               {
01796                 for (i = 0; i < n[d] / 2; i++)
01797                   ex (i) =
01798                     tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h));
01799               }
01800             else
01801               {
01802                 for (i = 0; i < n[d] / 2; i++)
01803                   ex (i) = m + xi[(2 * i + 1) * skip[d]] * h;
01804               }
01805             fargs(0) = ex;
01806             fvals = feval (fcn, fargs, 1);
01807             if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
01808               {
01809                 error ("quadcc: integrand F must return a single, real-valued vector");
01810                 return retval;
01811               }
01812             Matrix effex = fvals(0).matrix_value ();
01813             if (effex.length () != ex.length ())
01814               {
01815                 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
01816                 return retval;
01817               }
01818             neval += effex.length ();
01819             for (i = 0; i < n[d] / 2; i++)
01820               {
01821                 j = (2 * i + 1) * skip[d];
01822                 iv->fx[j] = effex (i);
01823                 if (wrap)
01824                   {
01825                     xw = ex(i);
01826                     iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
01827                   }
01828               }
01829           }
01830           nnans = 0;
01831           for (i = 0; i <= 32; i += skip[d])
01832             {
01833               if (!xfinite (iv->fx[i]))
01834                 {
01835                   nans[nnans++] = i;
01836                   iv->fx[i] = 0.0;
01837                 }
01838             }
01839 
01840           /* Compute the new coefficients. */
01841           Vinvfx (iv->fx, &(iv->c[idx[d]]), d);
01842           /* Downdate any NaNs. */
01843           if (nnans > 0)
01844             {
01845               downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
01846               for (i = 0; i < nnans; i++)
01847                 iv->fx[nans[i]] = octave_NaN;
01848             }
01849 
01850           /* Compute the error estimate. */
01851           nc = 0.0;
01852           for (i = n[d - 1] + 1; i <= n[d]; i++)
01853             {
01854               temp = iv->c[idx[d] + i];
01855               nc += temp * temp;
01856             }
01857           ncdiff = nc;
01858           for (i = 0; i <= n[d - 1]; i++)
01859             {
01860               temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
01861               ncdiff += temp * temp;
01862               nc += iv->c[idx[d] + i] * iv->c[idx[d] + i];
01863             }
01864           ncdiff = sqrt (ncdiff);
01865           nc = sqrt (nc);
01866           iv->err = ncdiff * 2 * h;
01867           /* Compute the local integral. */
01868           iv->igral = 2 * h * w * iv->c[idx[d]];
01869           /* Split the interval prematurely? */
01870           split = (nc > 0 && ncdiff / nc > 0.1);
01871         }
01872 
01873       /* Maximum degree reached, just split. */
01874       else
01875         {
01876           split = 1;
01877         }
01878 
01879 
01880       /* Should we drop this interval? */
01881       if ((m + h * xi[0]) >= (m + h * xi[1])
01882           || (m + h * xi[31]) >= (m + h * xi[32])
01883           || iv->err < fabs (iv->igral) * DBL_EPSILON * 10)
01884         {
01885 
01886 /*          printf
01887             ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
01888              heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
01889              iv->depth);
01890 */
01891           /* Keep this interval's contribution */
01892           err_final += iv->err;
01893           igral_final += iv->igral;
01894           /* Swap with the last element on the heap */
01895           t = heap[nivals - 1];
01896           heap[nivals - 1] = heap[0];
01897           heap[0] = t;
01898           nivals--;
01899           /* Fix up the heap */
01900           i = 0;
01901           while (2 * i + 1 < nivals)
01902             {
01903 
01904               /* Get the kids */
01905               j = 2 * i + 1;
01906               /* If the j+1st entry exists and is larger than the jth,
01907                  use it instead. */
01908               if (j + 1 < nivals
01909                   && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
01910                 j++;
01911               /* Do we need to move the ith entry up? */
01912               if (ivals[heap[j]].err <= ivals[heap[i]].err)
01913                 break;
01914               else
01915                 {
01916                   t = heap[j];
01917                   heap[j] = heap[i];
01918                   heap[i] = t;
01919                   i = j;
01920                 }
01921             }
01922 
01923         }
01924 
01925       /* Do we need to split this interval? */
01926       else if (split)
01927         {
01928 
01929           /* Some values we will need often... */
01930           d = iv->depth;
01931           /* Generate the interval on the left */
01932           ivl = &(ivals[heap[nivals++]]);
01933           ivl->a = iv->a;
01934           ivl->b = m;
01935           ml = (ivl->a + ivl->b) / 2;
01936           hl = h / 2;
01937           ivl->depth = 0;
01938           ivl->rdepth = iv->rdepth + 1;
01939           ivl->fx[0] = iv->fx[0];
01940           ivl->fx[32] = iv->fx[16];
01941           {
01942             ColumnVector ex (n[0] - 1);
01943             if (wrap)
01944               {
01945                 for (i = 0; i < n[0] - 1; i++)
01946                   ex (i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl));
01947               }
01948             else
01949               {
01950                 for (i = 0; i < n[0] - 1; i++)
01951                   ex (i) = ml + xi[(i + 1) * skip[0]] * hl;
01952               }
01953             fargs(0) = ex;
01954             fvals = feval (fcn, fargs, 1);
01955             if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
01956               {
01957                 error ("quadcc: integrand F must return a single, real-valued vector");
01958                 return retval;
01959               }
01960             Matrix effex = fvals(0).matrix_value ();
01961             if (effex.length () != ex.length ())
01962               {
01963                 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
01964                 return retval;
01965               }
01966             neval += effex.length ();
01967             for (i = 0; i < n[0] - 1; i++)
01968               {
01969                 j = (i + 1) * skip[0];
01970                 ivl->fx[j] = effex (i);
01971                 if (wrap)
01972                   {
01973                     xw = ex(i);
01974                     ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
01975                   }
01976               }
01977           }
01978           nnans = 0;
01979           for (i = 0; i <= 32; i += skip[0])
01980             {
01981               if (!xfinite (ivl->fx[i]))
01982                 {
01983                   nans[nnans++] = i;
01984                   ivl->fx[i] = 0.0;
01985                 }
01986             }
01987           Vinvfx (ivl->fx, ivl->c, 0);
01988           if (nnans > 0)
01989             {
01990               downdate (ivl->c, n[0], 0, nans, nnans);
01991               for (i = 0; i < nnans; i++)
01992                 ivl->fx[nans[i]] = octave_NaN;
01993             }
01994           for (i = 0; i <= n[d]; i++)
01995             {
01996               ivl->c[idx[d] + i] = 0.0;
01997               for (j = i; j <= n[d]; j++)
01998                 ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j];
01999             }
02000           ncdiff = 0.0;
02001           for (i = 0; i <= n[0]; i++)
02002             {
02003               temp = ivl->c[i] - ivl->c[idx[d] + i];
02004               ncdiff += temp * temp;
02005             }
02006           for (i = n[0] + 1; i <= n[d]; i++)
02007             {
02008               temp = ivl->c[idx[d] + i];
02009               ncdiff += temp * temp;
02010             }
02011           ncdiff = sqrt (ncdiff);
02012           ivl->err = ncdiff * h;
02013           /* Check for divergence. */
02014           ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
02015                                   && ivl->c[0] / iv->c[0] > 2);
02016           if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth)
02017             {
02018               igral = gnulib::copysign (octave_Inf, igral);
02019               warning ("quadcc: divergent integral detected");
02020               break;
02021             }
02022 
02023           /* Compute the local integral. */
02024           ivl->igral = h * w * ivl->c[0];
02025 
02026 
02027           /* Generate the interval on the right */
02028           ivr = &(ivals[heap[nivals++]]);
02029           ivr->a = m;
02030           ivr->b = iv->b;
02031           mr = (ivr->a + ivr->b) / 2;
02032           hr = h / 2;
02033           ivr->depth = 0;
02034           ivr->rdepth = iv->rdepth + 1;
02035           ivr->fx[0] = iv->fx[16];
02036           ivr->fx[32] = iv->fx[32];
02037           {
02038             ColumnVector ex (n[0] - 1);
02039             if (wrap)
02040               {
02041                 for (i = 0; i < n[0] - 1; i++)
02042                   ex (i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr));
02043               }
02044             else
02045               {
02046                 for (i = 0; i < n[0] - 1; i++)
02047                   ex (i) = mr + xi[(i + 1) * skip[0]] * hr;
02048               }
02049             fargs(0) = ex;
02050             fvals = feval (fcn, fargs, 1);
02051             if (fvals.length () != 1 || !fvals(0).is_real_matrix ())
02052               {
02053                 error ("quadcc: integrand F must return a single, real-valued vector");
02054                 return retval;
02055               }
02056             Matrix effex = fvals(0).matrix_value ();
02057             if (effex.length () != ex.length ())
02058               {
02059                 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
02060                 return retval;
02061               }
02062             neval += effex.length ();
02063             for (i = 0; i < n[0] - 1; i++)
02064               {
02065                 j = (i + 1) * skip[0];
02066                 ivr->fx[j] = effex (i);
02067                 if (wrap)
02068                   {
02069                     xw = ex(i);
02070                     ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
02071                   }
02072               }
02073           }
02074           nnans = 0;
02075           for (i = 0; i <= 32; i += skip[0])
02076             {
02077               if (!xfinite (ivr->fx[i]))
02078                 {
02079                   nans[nnans++] = i;
02080                   ivr->fx[i] = 0.0;
02081                 }
02082             }
02083           Vinvfx (ivr->fx, ivr->c, 0);
02084           if (nnans > 0)
02085             {
02086               downdate (ivr->c, n[0], 0, nans, nnans);
02087               for (i = 0; i < nnans; i++)
02088                 ivr->fx[nans[i]] = octave_NaN;
02089             }
02090           for (i = 0; i <= n[d]; i++)
02091             {
02092               ivr->c[idx[d] + i] = 0.0;
02093               for (j = i; j <= n[d]; j++)
02094                 ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j];
02095             }
02096           ncdiff = 0.0;
02097           for (i = 0; i <= n[0]; i++)
02098             {
02099               temp = ivr->c[i] - ivr->c[idx[d] + i];
02100               ncdiff += temp * temp;
02101             }
02102           for (i = n[0] + 1; i <= n[d]; i++)
02103             {
02104               temp = ivr->c[idx[d] + i];
02105               ncdiff += temp * temp;
02106             }
02107           ncdiff = sqrt (ncdiff);
02108           ivr->err = ncdiff * h;
02109           /* Check for divergence. */
02110           ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
02111                                   && ivr->c[0] / iv->c[0] > 2);
02112           if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth)
02113             {
02114               igral = gnulib::copysign (octave_Inf, igral);
02115               warning ("quadcc: divergent integral detected");
02116               break;
02117             }
02118 
02119           /* Compute the local integral. */
02120           ivr->igral = h * w * ivr->c[0];
02121 
02122 
02123           /* Fix-up the heap: we now have one interval on top
02124              that we don't need any more and two new, unsorted
02125              ones at the bottom. */
02126           /* Flip the last interval to the top of the heap and
02127              sift down. */
02128           t = heap[nivals - 1];
02129           heap[nivals - 1] = heap[0];
02130           heap[0] = t;
02131           nivals--;
02132           /* Sift this interval back down the heap. */
02133           i = 0;
02134           while (2 * i + 1 < nivals - 1)
02135             {
02136               j = 2 * i + 1;
02137               if (j + 1 < nivals - 1
02138                   && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
02139                 j++;
02140               if (ivals[heap[j]].err <= ivals[heap[i]].err)
02141                 break;
02142               else
02143                 {
02144                   t = heap[j];
02145                   heap[j] = heap[i];
02146                   heap[i] = t;
02147                   i = j;
02148                 }
02149             }
02150 
02151           /* Now grab the last interval and sift it up the heap. */
02152           i = nivals - 1;
02153           while (i > 0)
02154             {
02155               j = (i - 1) / 2;
02156               if (ivals[heap[j]].err < ivals[heap[i]].err)
02157                 {
02158                   t = heap[j];
02159                   heap[j] = heap[i];
02160                   heap[i] = t;
02161                   i = j;
02162                 }
02163               else
02164                 break;
02165             }
02166 
02167 
02168         }
02169 
02170       /* Otherwise, just fix-up the heap. */
02171       else
02172         {
02173           i = 0;
02174           while (2 * i + 1 < nivals)
02175             {
02176               j = 2 * i + 1;
02177               if (j + 1 < nivals
02178                   && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
02179                 j++;
02180               if (ivals[heap[j]].err <= ivals[heap[i]].err)
02181                 break;
02182               else
02183                 {
02184                   t = heap[j];
02185                   heap[j] = heap[i];
02186                   heap[i] = t;
02187                   i = j;
02188                 }
02189             }
02190 
02191         }
02192 
02193       /* If the heap is about to overflow, remove the last two
02194          intervals. */
02195       while (nivals > cquad_heapsize - 2)
02196         {
02197           iv = &(ivals[heap[nivals - 1]]);
02198 /*          printf
02199             ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
02200              heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
02201              iv->depth);
02202 */
02203           err_final += iv->err;
02204           igral_final += iv->igral;
02205           nivals--;
02206         }
02207 
02208       /* Collect the value of the integral and error. */
02209       igral = igral_final;
02210       err = err_final;
02211       for (i = 0; i < nivals; i++)
02212         {
02213           igral += ivals[heap[i]].igral;
02214           err += ivals[heap[i]].err;
02215         }
02216 
02217     }
02218 
02219   /* Dump the contents of the heap. */
02220 /*  for (i = 0; i < nivals; i++)
02221     {
02222       iv = &(ivals[heap[i]]);
02223       printf
02224         ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n",
02225          i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
02226          iv->rdepth, iv->ndiv);
02227     }
02228 */
02229   /* Clean up and present the results. */
02230   if (nargout > 2)
02231     retval(2) = neval;
02232   if (nargout > 1)
02233     retval(1) = err;
02234   retval(0) = igral;
02235   /* All is well that ends well. */
02236   return retval;
02237 }
02238 
02239 
02240 /*
02241 
02242 %!assert (quadcc(@sin,-pi,pi), 0, 1e-6)
02243 %!assert (quadcc(inline('sin'),-pi,pi), 0, 1e-6)
02244 %!assert (quadcc('sin',-pi,pi), 0, 1e-6)
02245 
02246 %!assert (quadcc(@sin,-pi,0), -2, 1e-6)
02247 %!assert (quadcc(@sin,0,pi), 2, 1e-6)
02248 %!assert (quadcc(@(x) 1./sqrt(x), 0, 1), 2, 1e-6)
02249 %!assert (quadcc(@(x) 1./(sqrt(x).*(x+1)), 0, Inf), pi, 1e-6)
02250 
02251 %!assert (quadcc (@(x) exp(-x .^ 2), -Inf, Inf), sqrt(pi), 1e-6)
02252 %!assert (quadcc (@(x) exp(-x .^ 2), -Inf, 0), sqrt(pi)/2, 1e-6)
02253 
02254 ## Test function with NaNs in interval 
02255 %!function y = __nansin (x)
02256 %!  nan_locs = [-3*pi/4, -pi/4, 0, pi/3, pi/2, pi];
02257 %!  y = sin (x);
02258 %!  idx = min (abs (bsxfun (@minus, x(:), nan_locs)), [], 2); 
02259 %!  y(idx < 1e-10) = NaN;
02260 %!endfunction 
02261 
02262 %!test
02263 %! [q, err, npoints] = quadcc ("__nansin", -pi, pi); 
02264 %! assert (q, 0, eps);
02265 %! assert (err, 0, 15*eps);
02266 
02267 %% Test input validation
02268 %!error (quadcc ())
02269 %!error (quadcc (@sin))
02270 %!error (quadcc (@sin, 0))
02271 %!error (quadcc (@sin, ones(2), pi))
02272 %!error (quadcc (@sin, -i, pi))
02273 %!error (quadcc (@sin, 0, ones(2)))
02274 %!error (quadcc (@sin, 0, i))
02275 %!error (quadcc (@sin, 0, pi, 0))
02276 %!error (quadcc (@sin, 0, pi, 1e-6, [ i ]))
02277 
02278 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Thu Mar 7 2013 16:03:23 for GNU Octave by  doxygen 1.7.1