Mercurial > octave
changeset 11562:1811f4f8b3b5
Rename cquad to quadcc and add to documentation.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Mon, 17 Jan 2011 20:18:07 -0800 |
parents | 1140bd9f9b21 |
children | 3c6e8aaa9555 |
files | ChangeLog NEWS doc/ChangeLog doc/interpreter/quad.txi scripts/ChangeLog scripts/general/dblquad.m scripts/general/quadgk.m scripts/general/quadl.m scripts/general/quadv.m scripts/general/triplequad.m src/ChangeLog src/DLD-FUNCTIONS/cquad.cc src/DLD-FUNCTIONS/module-files src/DLD-FUNCTIONS/quad.cc src/DLD-FUNCTIONS/quadcc.cc |
diffstat | 15 files changed, 2297 insertions(+), 2273 deletions(-) [+] |
line wrap: on
line diff
--- a/ChangeLog Mon Jan 17 17:25:42 2011 -0500 +++ b/ChangeLog Mon Jan 17 20:18:07 2011 -0800 @@ -1,3 +1,7 @@ +2010-01-17 Rik <octave@nomad.inbox5.com> + + * NEWS: Add quadcc to list of new functions added. + 2011-01-17 John W. Eaton <jwe@octave.org> * mkoctfile.cc.in (main): Add + missing from previous change.
--- a/NEWS Mon Jan 17 17:25:42 2011 -0500 +++ b/NEWS Mon Jan 17 20:18:07 2011 -0800 @@ -62,7 +62,7 @@ `issymmetric' has been changed for better consistency. * The `ismatrix' function now returns true for all numeric, - logical and character 2d or Nd matrices. Previously, `ismatrix' + logical and character 2-D or N-D matrices. Previously, `ismatrix' returned false if the first or second dimension was zero. Hence, `ismatrix ([])' was false, while `ismatrix (zeros (1,2,0))' was true. @@ -76,7 +76,7 @@ now yields true. * The `issymmetric' function now checks for symmetry instead of - hermitianness. For the latter, ishermitian was created. Also, + Hermitianness. For the latter, ishermitian was created. Also, logical scalar is returned rather than the dimension, so `issymmetric ([])' is now true. @@ -313,15 +313,15 @@ ** The following functions are new in Octave 3.4: - accumdim erfcx nfields pqpnonneg uigetfile - bitpack fileread nth_element randi uiputfile - bitunpack fminbnd onCleanup repelems uimenu - blkmm fskipl pbaspect reset whitebg - cbrt ifelse pie3 rsf2csf - curl ishermitian powerset saveas - chop isindex ppder strread - daspect luupdate ppint textread - divergence merge ppjumps uigetdir + accumdim erfcx nfields pqpnonneg uigetdir + bitpack fileread nth_element quadcc uigetfile + bitunpack fminbnd onCleanup randi uiputfile + blkmm fskipl pbaspect repelems uimenu + cbrt ifelse pie3 reset whitebg + curl ishermitian powerset rsf2csf + chop isindex ppder saveas + daspect luupdate ppint strread + divergence merge ppjumps textread ** Using the image function to view images with external programs such as display, xv, and xloadimage is no longer supported. The
--- a/doc/ChangeLog Mon Jan 17 17:25:42 2011 -0500 +++ b/doc/ChangeLog Mon Jan 17 20:18:07 2011 -0800 @@ -1,3 +1,7 @@ +2011-01-17 Rik <octave@nomad.inbox5.com> + + * interpreter/quad.txi: Add quadcc to documentation. + 2011-01-15 Rik <octave@nomad.inbox5.com> * interpreter/expr.txi: Add merge/ifelse function to documentation.
--- a/doc/interpreter/quad.txi Mon Jan 17 17:25:42 2011 -0500 +++ b/doc/interpreter/quad.txi Mon Jan 17 20:18:07 2011 -0800 @@ -54,6 +54,9 @@ @item quadv Numerical integration using an adaptive vectorized Simpson's rule. +@item quadcc +Numerical integration using adaptive Clenshaw-Curtis rules. + @item trapz Numerical integration using the trapezoidal method. @end table @@ -124,6 +127,8 @@ @DOCSTRING(quadv) +@DOCSTRING(quadcc) + @DOCSTRING(trapz) @DOCSTRING(cumtrapz)
--- a/scripts/ChangeLog Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/ChangeLog Mon Jan 17 20:18:07 2011 -0800 @@ -1,3 +1,9 @@ +2011-01-17 Rik <octave@nomad.inbox5.com> + + * general/dblquad.m, general/quadgk.m, general/quadl.m, + general/quadv.m, general/triplequad.m: Improve docstring with seealso + links to quadcc. + 2011-01-17 John W. Eaton <jwe@octave.org> * miscellaneous/isdeployed.m: New function.
--- a/scripts/general/dblquad.m Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/general/dblquad.m Mon Jan 17 20:18:07 2011 -0800 @@ -29,7 +29,7 @@ ## ## Additional arguments, are passed directly to @var{f}. To use the default ## value for @var{tol} one may pass an empty matrix. -## @seealso{triplequad, quad, quadv, quadl, quadgk, trapz} +## @seealso{triplequad,quad,quadv,quadl,quadgk,quadcc,trapz} ## @end deftypefn function q = dblquad(f, xa, xb, ya, yb, tol, quadf, varargin)
--- a/scripts/general/quadgk.m Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/general/quadgk.m Mon Jan 17 20:18:07 2011 -0800 @@ -112,7 +112,7 @@ ## approximate bounds on the error in the integral @code{abs (@var{q} - ## @var{i})}, where @var{i} is the exact value of the integral. ## -## @seealso{triplequad, dblquad, quad, quadl, quadv, trapz} +## @seealso{quad,quadv,quadl,quadcc,trapz,dblquad,triplequad} ## @end deftypefn function [q, err] = quadgk (f, a, b, varargin)
--- a/scripts/general/quadl.m Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/general/quadl.m Mon Jan 17 20:18:07 2011 -0800 @@ -41,7 +41,7 @@ ## Reference: W. Gander and W. Gautschi, @cite{Adaptive Quadrature - ## Revisited}, BIT Vol. 40, No. 1, March 2000, pp. 84--101. ## @url{http://www.inf.ethz.ch/personal/gander/} -## +## @seealso{quad,quadv,quadgk,quadcc,trapz,dblquad,triplequad} ## @end deftypefn ## Author: Walter Gautschi
--- a/scripts/general/quadv.m Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/general/quadv.m Mon Jan 17 20:18:07 2011 -0800 @@ -42,7 +42,7 @@ ## To use default values for @var{tol} and @var{trace}, one may pass ## empty matrices. ## -## @seealso{triplequad, dblquad, quad, quadl, quadgk, trapz} +##@seealso{quad,quadl,quadgk,quadcc,trapz,dblquad,triplequad} ## @end deftypefn function [q, fcnt] = quadv (f, a, b, tol, trace, varargin)
--- a/scripts/general/triplequad.m Mon Jan 17 17:25:42 2011 -0500 +++ b/scripts/general/triplequad.m Mon Jan 17 20:18:07 2011 -0800 @@ -30,7 +30,7 @@ ## ## Additional arguments, are passed directly to @var{f}. To use the default ## value for @var{tol} one may pass an empty matrix. -## @seealso{dblquad, quad, quadv, quadl, quadgk, trapz} +## @seealso{dblquad,quad,quadv,quadl,quadgk,quadcc,trapz} ## @end deftypefn function Q = triplequad(f, xa, xb, ya, yb, za, zb, tol, quadf, varargin)
--- a/src/ChangeLog Mon Jan 17 17:25:42 2011 -0500 +++ b/src/ChangeLog Mon Jan 17 20:18:07 2011 -0800 @@ -1,3 +1,8 @@ +2011-01-17 Rik <octave@nomad.inbox5.com> + + * DLD-FUNCTIONS/module-files: Add quadcc.cc to list of files. + * DLD-FUNCTIONS/quad.cc: Add Seealso links to quadcc. + 2011-01-17 Jaroslav Hajek <highegg@gmail.com> * DLD-FUNCTION/lookup.cc (Flookup): Validate option string.
--- a/src/DLD-FUNCTIONS/cquad.cc Mon Jan 17 17:25:42 2011 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,2255 +0,0 @@ -/* - -Copyright (C) 2010-2011 Pedro Gonnet - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdlib.h> -#include "lo-math.h" -#include "lo-ieee.h" -#include "oct.h" -#include "parse.h" -#include "ov-fcn-handle.h" - -#if ! defined (HAVE_COPYSIGN) && defined (HAVE__COPYSIGN) -#define copysign _copysign -#define HAVE_COPYSIGN 1 -#endif - -/* Define the size of the interval heap. */ -#define cquad_heapsize 200 - - -/* Data of a single interval */ -typedef struct -{ - double a, b; - double c[64]; - double fx[33]; - double igral, err; - int depth, rdepth, ndiv; -} cquad_ival; - -/* Some constants and matrices that we'll need. - */ - -static const double xi[33] = { - -1., -0.99518472667219688624, -0.98078528040323044912, - -0.95694033573220886493, -0.92387953251128675612, - -0.88192126434835502970, -0.83146961230254523708, - -0.77301045336273696082, -0.70710678118654752440, - -0.63439328416364549822, -0.55557023301960222475, - -0.47139673682599764857, -0.38268343236508977173, - -0.29028467725446236764, -0.19509032201612826785, - -0.098017140329560601995, 0., 0.098017140329560601995, - 0.19509032201612826785, 0.29028467725446236764, 0.38268343236508977173, - 0.47139673682599764857, 0.55557023301960222475, 0.63439328416364549822, - 0.70710678118654752440, 0.77301045336273696082, 0.83146961230254523708, - 0.88192126434835502970, 0.92387953251128675612, 0.95694033573220886493, - 0.98078528040323044912, 0.99518472667219688624, 1. -}; - -static const double bee[68] = { - 0.00000000000000e+00, 2.28868854108532e-01, 0.00000000000000e+00, - -8.15740215243451e-01, 0.00000000000000e+00, 5.31212715259731e-01, - 0.00000000000000e+00, 1.38538036812454e-02, 0.00000000000000e+00, - 3.74405228908818e-02, 0.00000000000000e+00, 2.12224115039342e-01, - 0.00000000000000e+00, -8.16362644507898e-01, 0.00000000000000e+00, - 5.35648426691481e-01, 0.00000000000000e+00, 1.52417902753662e-03, - 0.00000000000000e+00, 2.63058840550873e-03, 0.00000000000000e+00, - 4.15292106318904e-03, 0.00000000000000e+00, 6.97106011119775e-03, - 0.00000000000000e+00, 1.35535708431058e-02, 0.00000000000000e+00, - 3.52132898424856e-02, 0.00000000000000e+00, 2.06946714741884e-01, - 0.00000000000000e+00, -8.15674251283876e-01, 0.00000000000000e+00, - 5.38841175520580e-01, 0.00000000000000e+00, 1.84909689577590e-04, - 0.00000000000000e+00, 2.90936325007499e-04, 0.00000000000000e+00, - 3.84877750950089e-04, 0.00000000000000e+00, 4.86436656735046e-04, - 0.00000000000000e+00, 6.08688640346879e-04, 0.00000000000000e+00, - 7.66732830740331e-04, 0.00000000000000e+00, 9.82753336104205e-04, - 0.00000000000000e+00, 1.29359957505615e-03, 0.00000000000000e+00, - 1.76616363801885e-03, 0.00000000000000e+00, 2.53323433039089e-03, - 0.00000000000000e+00, 3.88872172121956e-03, 0.00000000000000e+00, - 6.58635106468291e-03, 0.00000000000000e+00, 1.30326736343254e-02, - 0.00000000000000e+00, 3.44353850696714e-02, 0.00000000000000e+00, - 2.05025409531915e-01, 0.00000000000000e+00, -8.14985893995401e-01, - 0.00000000000000e+00, 5.40679930965238e-01 -}; - -static const double Lalpha[33] = { - 5.77350269189626e-01, 5.16397779494322e-01, 5.07092552837110e-01, - 5.03952630678970e-01, 5.02518907629606e-01, 5.01745206004255e-01, - 5.01280411827603e-01, 5.00979432868120e-01, 5.00773395667191e-01, - 5.00626174321759e-01, 5.00517330712619e-01, 5.00434593736979e-01, - 5.00370233297676e-01, 5.00319182924304e-01, 5.00278009473803e-01, - 5.00244319584578e-01, 5.00216403386025e-01, 5.00193012939056e-01, - 5.00173220168024e-01, 5.00156323280355e-01, 5.00141783641018e-01, - 5.00129182278347e-01, 5.00118189340972e-01, 5.00108542278496e-01, - 5.00100030010004e-01, 5.00092481273333e-01, 5.00085755939229e-01, - 5.00079738458365e-01, 5.00074332862969e-01, 5.00069458915387e-01, - 5.00065049112355e-01, 5.00061046334395e-01, 5.00057401986298e-01 -}; - -static const double Lgamma[33] = { - 0.0, 0.0, 5.16397779494322e-01, 5.07092552837110e-01, 5.03952630678970e-01, - 5.02518907629606e-01, 5.01745206004255e-01, 5.01280411827603e-01, - 5.00979432868120e-01, 5.00773395667191e-01, 5.00626174321759e-01, - 5.00517330712619e-01, 5.00434593736979e-01, 5.00370233297676e-01, - 5.00319182924304e-01, 5.00278009473803e-01, 5.00244319584578e-01, - 5.00216403386025e-01, 5.00193012939056e-01, 5.00173220168024e-01, - 5.00156323280355e-01, 5.00141783641018e-01, 5.00129182278347e-01, - 5.00118189340972e-01, 5.00108542278496e-01, 5.00100030010003e-01, - 5.00092481273333e-01, 5.00085755939229e-01, 5.00079738458365e-01, - 5.00074332862969e-01, 5.00069458915387e-01, 5.00065049112355e-01, - 5.00061046334395e-01 -}; - -static const double V1inv[5 * 5] = { - .47140452079103168293e-1, .37712361663282534635, .56568542494923801952, - .37712361663282534635, .47140452079103168293e-1, - -.81649658092772603273e-1, -.46188021535170061160, 0, - .46188021535170061160, .81649658092772603273e-1, .15058465048420853962, - .12046772038736683169, -.54210474174315074262, .12046772038736683169, - .15058465048420853962, -.21380899352993950775, .30237157840738178177, -0., - -.30237157840738178177, .21380899352993950775, .10774960475223581324, - -.21549920950447162648, .21549920950447162648, -.21549920950447162648, - .10774960475223581324 -}; - -static const double V2inv[9 * 9] = { - .11223917161691230546e-1, .10339219839658349826, .19754094204576565761, - .25577315077753587922, .27835314560994251755, .25577315077753587922, - .19754094204576565761, .10339219839658349826, .11223917161691230546e-1, - -.19440394783993476970e-1, -.16544884625069155470, -.24193725566041460608, - -.16953338808305493604, 0.0, .16953338808305493604, .24193725566041460608, - .16544884625069155470, .19440394783993476970e-1, .26466393115406349388e-1, - .17766815796285469394, .11316664642449611462, -.16306601003711325980, - -.30847037493128779631, -.16306601003711325980, .11316664642449611462, - .17766815796285469394, .26466393115406349388e-1, - -.32395302049990834508e-1, -.15521142532414866547, - .88573492664788602740e-1, .29570405784974857322, 0.0, - -.29570405784974857322, -.88573492664788602740e-1, .15521142532414866547, - .32395302049990834508e-1, .41442155673936851246e-1, - .98186757907405608245e-1, -.23056908429499411784, - -.68047008326360625520e-1, .31797435808002456774, - -.68047008326360625520e-1, -.23056908429499411784, - .98186757907405608245e-1, .41442155673936851246e-1, - -.49981120317798783134e-1, -.24861810572835756217e-1, - .23561326072010832539, -.24472785656448415351, 0.0, .24472785656448415351, - -.23561326072010832539, .24861810572835756217e-1, - .49981120317798783134e-1, .79691635865674781228e-1, - -.95725617891693941833e-1, -.57957553356854386344e-1, - .21164072460540271452, -.27529837844505833514, .21164072460540271452, - -.57957553356854386344e-1, -.95725617891693941833e-1, - .79691635865674781228e-1, - -.10894869830716590913, .20131094491947531782, -.15407672674888869038, - .83385723639789791384e-1, 0.0, -.83385723639789791384e-1, - .15407672674888869038, -.20131094491947531782, .10894869830716590913, - .54581057089643838221e-1, -.10916211417928767644, .10916211417928767644, - -.10916211417928767644, .10916211417928767644, -.10916211417928767644, - .10916211417928767644, -.10916211417928767644, .54581057089643838221e-1 -}; - -static const double V3inv[17 * 17] = { - .27729677693590098996e-2, .26423663180333065153e-1, - .53374068493933898312e-1, .77007854739523195947e-1, - .98257061072911596869e-1, .11538049741786835604, .12832134344120884559, - .13612785914022865001, .13888293186236181317, .13612785914022865001, - .12832134344120884559, .11538049741786835604, .98257061072911596869e-1, - .77007854739523195947e-1, .53374068493933898312e-1, - .26423663180333065153e-1, .27729677693590098996e-2, - -.48029210642807413690e-2, -.44887724635478800254e-1, - -.85409520147301089416e-1, -.11090267822061423050, -.12033983162705862441, - -.11102786862182788886, -.85054870109799336515e-1, - -.45998467987742225160e-1, 0.0, .45998467987742225160e-1, - .85054870109799336515e-1, .11102786862182788886, .12033983162705862441, - .11090267822061423050, .85409520147301089416e-1, .44887724635478800254e-1, - .48029210642807413690e-2, .62758546879582030087e-2, - .55561297093529155869e-1, - .93281491021051539742e-1, .92320151237493695139e-1, - .55077987469605684531e-1, - -.96998141716497488255e-2, -.80285961895427405567e-1, - -.13496839655913850224, - -.15512521776684524331, -.13496839655913850224, -.80285961895427405567e-1, - -.96998141716497488255e-2, .55077987469605684531e-1, - .92320151237493695139e-1, .93281491021051539742e-1, - .55561297093529155869e-1, .62758546879582030087e-2, - -.74850969394858555939e-2, -.61751608943839234096e-1, - -.82974150437304275958e-1, -.38437763431942633378e-1, - .45745502025779701366e-1, .12369235652734542162, .14720439712852868239, - .98768034347019704401e-1, 0.0, - -.98768034347019704401e-1, -.14720439712852868239, -.12369235652734542162, - -.45745502025779701366e-1, .38437763431942633378e-1, - .82974150437304275958e-1, .61751608943839234096e-1, - .74850969394858555939e-2, .86710099994384056338e-2, - .64006230103659573344e-1, .58517426396091675690e-1, - -.29743410528985802680e-1, - -.11934127779157114754, -.12686773515361299409, -.30729137153877447035e-1, - .97307836256600731568e-1, .15635811574451401023, .97307836256600731568e-1, - -.30729137153877447035e-1, -.12686773515361299409, -.11934127779157114754, - -.29743410528985802680e-1, .58517426396091675690e-1, - .64006230103659573344e-1, .86710099994384056338e-2, - -.97486395666294840165e-2, -.62995604908060224672e-1, - -.24373234450275529219e-1, .87760984413626872730e-1, - .12205204576993351394, - .16216004196864002088e-1, -.12422320942156845775, -.13682714580929614678, - 0.0, .13682714580929614678, .12422320942156845775, - -.16216004196864002088e-1, -.12205204576993351394, - -.87760984413626872730e-1, .24373234450275529219e-1, - .62995604908060224672e-1, .97486395666294840165e-2, - .10956271233750488468e-1, .58613204255294358939e-1, - -.13306063940736618859e-1, -.11606666444978454399, - -.52059598001115805639e-1, .10868540217796151849, .12594452879014618005, - -.44678658254872910434e-1, -.15617684362128533405, - -.44678658254872910434e-1, .12594452879014618005, .10868540217796151849, - -.52059598001115805639e-1, -.11606666444978454399, - -.13306063940736618859e-1, .58613204255294358939e-1, - .10956271233750488468e-1, -.12098893000863087230e-1, - -.51626244709126208453e-1, .48919433304746979330e-1, - .10467644465949427090, - -.48729879523084673782e-1, -.13668732103524749234, .28190838706814496438e-1, - .15434223333238741600, 0.0, -.15434223333238741600, - -.28190838706814496438e-1, .13668732103524749234, - .48729879523084673782e-1, -.10467644465949427090, - -.48919433304746979330e-1, .51626244709126208453e-1, - .12098893000863087230e-1, .13542668300437944822e-1, - .41712033418258689308e-1, - -.76190463272803434388e-1, -.58303943170068132010e-1, .12158068748245606853, - .42121099930651007882e-1, -.14684425840766337756, - -.16108203535058647043e-1, .15698075850757976092, - -.16108203535058647043e-1, -.14684425840766337756, - .42121099930651007882e-1, .12158068748245606853, - -.58303943170068132010e-1, -.76190463272803434388e-1, - .41712033418258689308e-1, .13542668300437944822e-1, - -.14939634995117694417e-1, -.30047246373341564039e-1, - .91624635082546425678e-1, -.79133374319110026377e-2, - -.12292558212072233355, .90013382617762643524e-1, - .84013717196539593395e-1, -.14813033309980695856, 0.0, - .14813033309980695856, -.84013717196539593395e-1, - -.90013382617762643524e-1, - .12292558212072233355, .79133374319110026377e-2, -.91624635082546425678e-1, - .30047246373341564039e-1, .14939634995117694417e-1, - .16986031342807474208e-1, - .15760203882617033601e-1, -.91494054040950941996e-1, - .70082459207876130806e-1, - .53390713710144539104e-1, -.14340746778352039430, .84048122493418898508e-1, - .72456667788091316868e-1, -.15564535320096811360, - .72456667788091316868e-1, .84048122493418898508e-1, - -.14340746778352039430, .53390713710144539104e-1, - .70082459207876130806e-1, -.91494054040950941996e-1, - .15760203882617033601e-1, - .16986031342807474208e-1, -.18994065631858742028e-1, - -.82901821370405592927e-3, .77239669773015192888e-1, - -.10850735431039424680, .47524484622086496464e-1, - .69148184871588737021e-1, -.14829314646228194928, .11992057742398672066, - 0.0, -.11992057742398672066, .14829314646228194928, - -.69148184871588737021e-1, -.47524484622086496464e-1, - .10850735431039424680, -.77239669773015192888e-1, - .82901821370405592927e-3, .18994065631858742028e-1, - .22761703826371535132e-1, -.17728848711449643358e-1, - -.47496371572480503788e-1, .10659958402328690063, -.11696013966166296514, - .63073750910894244526e-1, .32928881123602721303e-1, - -.12280950532497593683, .15926189077282729505, -.12280950532497593683, - .32928881123602721303e-1, .63073750910894244526e-1, - -.11696013966166296514, .10659958402328690063, -.47496371572480503788e-1, - -.17728848711449643358e-1, .22761703826371535132e-1, - -.26493215276042203434e-1, .35579780856128386192e-1, - .10447309718398935122e-1, -.68616154085314996709e-1, - .11775363082763954214, -.13918901977011837274, .12312819418827395690, - -.72053565748259077905e-1, 0.0, .72053565748259077905e-1, - -.12312819418827395690, .13918901977011837274, -.11775363082763954214, - .68616154085314996709e-1, -.10447309718398935122e-1, - -.35579780856128386192e-1, - .26493215276042203434e-1, .40742523354399706918e-1, - -.73124912999529117195e-1, .49317266444153837821e-1, - -.13686605413876015320e-1, -.28342624942191100464e-1, - .70371855298258216249e-1, -.10600251632853603875, .12981016288391131812, - -.13817029659318161476, .12981016288391131812, -.10600251632853603875, - .70371855298258216249e-1, -.28342624942191100464e-1, - -.13686605413876015320e-1, - .49317266444153837821e-1, -.73124912999529117195e-1, - .40742523354399706918e-1, -.54944368958699908688e-1, - .10777725663147408190, -.10152395581538265428, .91369146312596428468e-1, - -.77703071757424700773e-1, .61050911730999815031e-1, - -.42052599404498348871e-1, .21438229266251454773e-1, 0.0, - -.21438229266251454773e-1, .42052599404498348871e-1, - -.61050911730999815031e-1, .77703071757424700773e-1, - -.91369146312596428468e-1, - .10152395581538265428, -.10777725663147408190, .54944368958699908688e-1, - .27485608464748840573e-1, -.54971216929497681146e-1, - .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .54971216929497681146e-1, - -.54971216929497681146e-1, .27485608464748840573e-1 -}; - -static const double V4inv[33 * 33] = { - .69120897476690862600e-3, .66419939766331555194e-2, - .13600665164323186111e-1, .20122785860913684493e-1, - .26583214101668429944e-1, .32712713318999268739e-1, - .38576221976287138036e-1, .44033030938268925133e-1, - .49092709529622799673e-1, .53657949874312515646e-1, - .57724533144734311859e-1, .61219564530655179096e-1, - .64138907503837875026e-1, .66427905189318792009e-1, - .68088956652280022887e-1, .69083051391555695878e-1, - .69422738116739271449e-1, .69083051391555695878e-1, - .68088956652280022887e-1, .66427905189318792009e-1, - .64138907503837875026e-1, .61219564530655179096e-1, - .57724533144734311859e-1, .53657949874312515646e-1, - .49092709529622799673e-1, .44033030938268925133e-1, - .38576221976287138036e-1, .32712713318999268739e-1, - .26583214101668429944e-1, .20122785860913684493e-1, - .13600665164323186111e-1, .66419939766331555194e-2, - .69120897476690862600e-3, -.11972090629438798134e-2, - -.11448874821643225573e-1, -.23104401104002905904e-1, - -.33352899418646530133e-1, -.42538626424075425908e-1, - -.49969730733911825941e-1, -.55555454015360728353e-1, - -.58955533624852604918e-1, -.60126044219122513907e-1, - -.58959430451175833624e-1, -.55546925396227130606e-1, - -.49984739749347973762e-1, -.42513009141170294365e-1, - -.33399140950669746346e-1, -.23007690803851790829e-1, - -.11728275717520066169e-1, 0.0, .11728275717520066169e-1, - .23007690803851790829e-1, .33399140950669746346e-1, - .42513009141170294365e-1, .49984739749347973762e-1, - .55546925396227130606e-1, .58959430451175833624e-1, - .60126044219122513907e-1, .58955533624852604918e-1, - .55555454015360728353e-1, .49969730733911825941e-1, - .42538626424075425908e-1, .33352899418646530133e-1, - .23104401104002905904e-1, .11448874821643225573e-1, - .11972090629438798134e-2, .15501585012936019146e-2, - .14628781502199620482e-1, .28684915921474815271e-1, - .39299396074628048026e-1, .46393418975496284204e-1, - .48756902531094699526e-1, .46331333488337494692e-1, - .39012645376980228775e-1, .27452795421085791153e-1, - .12430953621169863781e-1, -.47682978056024928800e-2, - -.22825828045428973853e-1, - -.40195512090720278312e-1, -.55503004262826221955e-1, - -.67424537752827046308e-1, -.75020199300113606452e-1, - -.77607844312483656131e-1, -.75020199300113606452e-1, - -.67424537752827046308e-1, -.55503004262826221955e-1, - -.40195512090720278312e-1, -.22825828045428973853e-1, - -.47682978056024928800e-2, .12430953621169863781e-1, - .27452795421085791153e-1, .39012645376980228775e-1, - .46331333488337494692e-1, .48756902531094699526e-1, - .46393418975496284204e-1, .39299396074628048026e-1, - .28684915921474815271e-1, .14628781502199620482e-1, - .15501585012936019146e-2, -.18377757558949194214e-2, - -.17050470050949761565e-1, -.31952119564923250836e-1, - -.40197423449026348155e-1, - -.41205649520281371624e-1, -.33909965817492272248e-1, - -.19393664422115332144e-1, .56661049630886784692e-3, - .22948272173686561721e-1, .44489719570904738207e-1, - .61790363672287920596e-1, .72121014727028013894e-1, - .73627151185287858579e-1, .65784665375961398923e-1, - .49369676372333667559e-1, .26444326317059715065e-1, 0.0, - -.26444326317059715065e-1, -.49369676372333667559e-1, - -.65784665375961398923e-1, -.73627151185287858579e-1, - -.72121014727028013894e-1, -.61790363672287920596e-1, - -.44489719570904738207e-1, -.22948272173686561721e-1, - -.56661049630886784692e-3, .19393664422115332144e-1, - .33909965817492272248e-1, .41205649520281371624e-1, - .40197423449026348155e-1, .31952119564923250836e-1, - .17050470050949761565e-1, .18377757558949194214e-2, - .20942714740729767769e-2, .18935902405146518232e-1, - .33335840852491735126e-1, .36770680999102286065e-1, - .28873194534132768509e-1, .10267303017729535513e-1, - -.14607738306201572890e-1, -.40139568545572305818e-1, - -.59808326733858291561e-1, -.68528358823372627506e-1, - -.63306535387619244879e-1, -.44508601817574921056e-1, - -.15449116105605395357e-1, .17941083795006546367e-1, - .48747356011657242123e-1, .70329553984201665523e-1, - .78106117292526169663e-1, .70329553984201665523e-1, - .48747356011657242123e-1, .17941083795006546367e-1, - -.15449116105605395357e-1, -.44508601817574921056e-1, - -.63306535387619244879e-1, -.68528358823372627506e-1, - -.59808326733858291561e-1, - -.40139568545572305818e-1, -.14607738306201572890e-1, - .10267303017729535513e-1, .28873194534132768509e-1, - .36770680999102286065e-1, .33335840852491735126e-1, - .18935902405146518232e-1, .20942714740729767769e-2, - -.23245285491878278419e-2, -.20401404737639389919e-1, - -.33019548231022514097e-1, -.29709828426463720091e-1, - -.11760070922697422156e-1, .15987584743850393793e-1, - .43619012891472813485e-1, .61177322409671487721e-1, - .61144030218486655594e-1, - .41895377620089086167e-1, .80232011820644308033e-2, - -.30574701186675900915e-1, - -.62072243008844865848e-1, -.76336186183574765586e-1, - -.68435466095345537115e-1, -.40237669208466966207e-1, 0.0, - .40237669208466966207e-1, .68435466095345537115e-1, - .76336186183574765586e-1, .62072243008844865848e-1, - .30574701186675900915e-1, -.80232011820644308033e-2, - -.41895377620089086167e-1, -.61144030218486655594e-1, - -.61177322409671487721e-1, -.43619012891472813485e-1, - -.15987584743850393793e-1, .11760070922697422156e-1, - .29709828426463720091e-1, .33019548231022514097e-1, - .20401404737639389919e-1, .23245285491878278419e-2, - .25451717261579269307e-2, .21480418595666878775e-1, - .31177212469293007998e-1, .19816333607013379373e-1, - -.72439496274458793681e-2, -.38404203906598342397e-1, - -.57633632255322221046e-1, -.54070547403585392952e-1, - -.26249823354368866005e-1, .15643058212336881516e-1, - .54539832735118677194e-1, .73283028002473989724e-1, - .62835303524135936213e-1, .26175977027801048141e-1, - -.22193636309998606610e-1, -.62597049956093311234e-1, - -.78206986173170212505e-1, -.62597049956093311234e-1, - -.22193636309998606610e-1, .26175977027801048141e-1, - .62835303524135936213e-1, - .73283028002473989724e-1, .54539832735118677194e-1, - .15643058212336881516e-1, - -.26249823354368866005e-1, -.54070547403585392952e-1, - -.57633632255322221046e-1, -.38404203906598342397e-1, - -.72439496274458793681e-2, .19816333607013379373e-1, - .31177212469293007998e-1, .21480418595666878775e-1, - .25451717261579269307e-2, -.27506573922483820005e-2, - -.22224442095099251870e-1, -.27949927254215773020e-1, - -.80918481053370034987e-2, .25121859354449306916e-1, - .51563535009373061074e-1, .51936965107145960512e-1, - .22146626648171527753e-1, - -.24172689882103382748e-1, -.61731229104853568296e-1, - -.68477262429344201201e-1, -.38311232728303704742e-1, - .14160578713659552679e-1, .61248813427564184033e-1, - .77136328841293031805e-1, .52514801765183697988e-1, 0.0, - -.52514801765183697988e-1, -.77136328841293031805e-1, - -.61248813427564184033e-1, -.14160578713659552679e-1, - .38311232728303704742e-1, - .68477262429344201201e-1, .61731229104853568296e-1, - .24172689882103382748e-1, - -.22146626648171527753e-1, -.51936965107145960512e-1, - -.51563535009373061074e-1, -.25121859354449306916e-1, - .80918481053370034987e-2, .27949927254215773020e-1, - .22224442095099251870e-1, .27506573922483820005e-2, - .29562461131654311467e-2, .22630271480554450613e-1, - .23547399831373800971e-1, -.43964593440902476642e-2, - -.39055315767504970597e-1, -.52369643937940066804e-1, - -.28506131614971613422e-1, .19906048093338832322e-1, - .60408880866392420279e-1, .62493397473656883090e-1, - .21391278377641297859e-1, -.37302864786623254746e-1, - -.73665127933539496872e-1, -.61706142476854010202e-1, - -.78065168882546327888e-2, .52335307373945544428e-1, - .78278746279419264777e-1, .52335307373945544428e-1, - -.78065168882546327888e-2, -.61706142476854010202e-1, - -.73665127933539496872e-1, -.37302864786623254746e-1, - .21391278377641297859e-1, .62493397473656883090e-1, - .60408880866392420279e-1, .19906048093338832322e-1, - -.28506131614971613422e-1, -.52369643937940066804e-1, - -.39055315767504970597e-1, -.43964593440902476642e-2, - .23547399831373800971e-1, .22630271480554450613e-1, - .29562461131654311467e-2, -.31515718415504761303e-2, - -.22739451096655080673e-1, -.18157123602272119779e-1, - .16496480897167303621e-1, .46921166788569301124e-1, - .40644395739978416354e-1, -.46275803430732216900e-2, - -.52883375891308909486e-1, -.61116483226324111734e-1, - -.17411698764545629853e-1, .44773430013166822765e-1, - .73441577962383869198e-1, .42127368371995472815e-1, - -.25504645957196772465e-1, -.74126818045972742488e-1, - -.62780077864719287317e-1, 0.0, .62780077864719287317e-1, - .74126818045972742488e-1, .25504645957196772465e-1, - -.42127368371995472815e-1, -.73441577962383869198e-1, - -.44773430013166822765e-1, .17411698764545629853e-1, - .61116483226324111734e-1, .52883375891308909486e-1, - .46275803430732216900e-2, -.40644395739978416354e-1, - -.46921166788569301124e-1, -.16496480897167303621e-1, - .18157123602272119779e-1, .22739451096655080673e-1, - .31515718415504761303e-2, .33536559294882188208e-2, - .22535348942792006185e-1, - .12048629300953560767e-1, -.27166076791299493403e-1, - -.47492745604230978367e-1, -.19246623430993153174e-1, - .36231297307556299322e-1, .61713617181636122004e-1, - .25928029734266134490e-1, -.40478700752883602818e-1, - -.71053889866326412049e-1, -.31870824482961751482e-1, - .41515251100219081281e-1, .76481960760098381651e-1, - .36726509155999912440e-1, -.40090067032627055969e-1, - -.78270742903374539397e-1, -.40090067032627055969e-1, - .36726509155999912440e-1, .76481960760098381651e-1, - .41515251100219081281e-1, -.31870824482961751482e-1, - -.71053889866326412049e-1, -.40478700752883602818e-1, - .25928029734266134490e-1, .61713617181636122004e-1, - .36231297307556299322e-1, -.19246623430993153174e-1, - -.47492745604230978367e-1, -.27166076791299493403e-1, - .12048629300953560767e-1, .22535348942792006185e-1, - .33536559294882188208e-2, - -.35481220456925318865e-2, -.22062913693073191150e-1, - -.54487362861834144999e-2, .35438821865804087489e-1, - .40733077820527411302e-1, -.67403098138950720914e-2, - -.55559584405239171054e-1, -.42417050790865158745e-1, - .24499901971884704925e-1, .68721232891705409302e-1, - .34086082787461126592e-1, -.43441000373118474002e-1, - -.73878085292669148950e-1, -.18846995664706657127e-1, - .59827776178286834498e-1, .70644634584085901794e-1, 0.0, - -.70644634584085901794e-1, -.59827776178286834498e-1, - .18846995664706657127e-1, .73878085292669148950e-1, - .43441000373118474002e-1, -.34086082787461126592e-1, - -.68721232891705409302e-1, -.24499901971884704925e-1, - .42417050790865158745e-1, .55559584405239171054e-1, - .67403098138950720914e-2, -.40733077820527411302e-1, - -.35438821865804087489e-1, .54487362861834144999e-2, - .22062913693073191150e-1, .35481220456925318865e-2, - .37554176816665075631e-2, .21297045781589919482e-1, - -.13327293083183431816e-2, - -.40635299172764596484e-1, -.27659860508374175359e-1, - .31089232744083445986e-1, .56113781541334176109e-1, - .37577840643257763400e-2, -.60511227350664590865e-1, - -.46670556446129053853e-1, .33263195878575888247e-1, - .72757324720645228775e-1, .15011712351692283635e-1, - -.65601212994924119078e-1, -.60016855838843789772e-1, - .26220858553188665966e-1, .78322776605833552980e-1, - .26220858553188665966e-1, -.60016855838843789772e-1, - -.65601212994924119078e-1, - .15011712351692283635e-1, .72757324720645228775e-1, - .33263195878575888247e-1, - -.46670556446129053853e-1, -.60511227350664590865e-1, - .37577840643257763400e-2, .56113781541334176109e-1, - .31089232744083445986e-1, -.27659860508374175359e-1, - -.40635299172764596484e-1, -.13327293083183431816e-2, - .21297045781589919482e-1, .37554176816665075631e-2, - -.39566995305720591229e-2, -.20291873414438919995e-1, - .80617453830770930551e-2, .42270189157016547906e-1, - .10332624526759093004e-1, -.48054759547616142024e-1, - -.37678032941171643972e-1, - .36617192625732482394e-1, .61009425973424865714e-1, - -.95589113168026591466e-2, - -.71023202645076922361e-1, -.25097788086808784456e-1, - .62406621963267050244e-1, .56907293171100693511e-1, - -.36435383083882206257e-1, -.75790105119208756348e-1, 0.0, - .75790105119208756348e-1, .36435383083882206257e-1, - -.56907293171100693511e-1, -.62406621963267050244e-1, - .25097788086808784456e-1, .71023202645076922361e-1, - .95589113168026591466e-2, - -.61009425973424865714e-1, -.36617192625732482394e-1, - .37678032941171643972e-1, .48054759547616142024e-1, - -.10332624526759093004e-1, -.42270189157016547906e-1, - -.80617453830770930551e-2, .20291873414438919995e-1, - .39566995305720591229e-2, .41776092289182138591e-2, - .19013221163904414395e-1, -.14420609729849899876e-1, - -.40259160586844441220e-1, .86327811113710831649e-2, - .53564430703021034399e-1, .65469185402150431933e-2, - -.60383116311280629856e-1, - -.25657793784058876939e-1, .58745680576829226900e-1, - .45649937869034420296e-1, - -.49167932056844167772e-1, -.62696614328552187977e-1, - .32540234556426699997e-1, .74280410383464269758e-1, - -.11425672633410999870e-1, -.78280649404686404903e-1, - -.11425672633410999870e-1, .74280410383464269758e-1, - .32540234556426699997e-1, -.62696614328552187977e-1, - -.49167932056844167772e-1, .45649937869034420296e-1, - .58745680576829226900e-1, -.25657793784058876939e-1, - -.60383116311280629856e-1, .65469185402150431933e-2, - .53564430703021034399e-1, - .86327811113710831649e-2, -.40259160586844441220e-1, - -.14420609729849899876e-1, .19013221163904414395e-1, - .41776092289182138591e-2, -.43935502082478059199e-2, - -.17528761237509401631e-1, .20208915249153872535e-1, - .34734743119040669109e-1, -.26275910172353637955e-1, - -.46368003346018878786e-1, - .26800056330709381025e-1, .56681476464606609921e-1, - -.24749011438127255898e-1, - -.64934612189056658992e-1, .20333742247679279535e-1, - .71429299070059318651e-1, - -.14452513210428671266e-1, -.75793341281736586582e-1, - .74717094137184935270e-2, .78034921554757317374e-1, 0.0, - -.78034921554757317374e-1, -.74717094137184935270e-2, - .75793341281736586582e-1, .14452513210428671266e-1, - -.71429299070059318651e-1, -.20333742247679279535e-1, - .64934612189056658992e-1, .24749011438127255898e-1, - -.56681476464606609921e-1, - -.26800056330709381025e-1, .46368003346018878786e-1, - .26275910172353637955e-1, - -.34734743119040669109e-1, -.20208915249153872535e-1, - .17528761237509401631e-1, .43935502082478059199e-2, - .46379089482818671473e-2, .15791188144791287229e-1, - -.25134290048737455284e-1, -.26249795071946841205e-1, - .39960457575789924651e-1, .28111892450146525404e-1, - -.51026476400767918226e-1, - -.27266747278681831364e-1, .60708796647861610865e-1, - .23532306960642115854e-1, - -.68169639871532441111e-1, -.18204924701958312032e-1, - .73822890510656128485e-1, .11373392486424717019e-1, - -.77133324017644609416e-1, -.39295877480342619961e-2, - .78351902829418987960e-1, -.39295877480342619961e-2, - -.77133324017644609416e-1, .11373392486424717019e-1, - .73822890510656128485e-1, -.18204924701958312032e-1, - -.68169639871532441111e-1, .23532306960642115854e-1, - .60708796647861610865e-1, -.27266747278681831364e-1, - -.51026476400767918226e-1, .28111892450146525404e-1, - .39960457575789924651e-1, -.26249795071946841205e-1, - -.25134290048737455284e-1, .15791188144791287229e-1, - .46379089482818671473e-2, -.48780095920069827068e-2, - -.13886961667516983541e-1, .29071311049368895844e-1, - .15480559452075811600e-1, -.47527977686242313065e-1, - -.31929089844361042178e-2, .58015667638415922967e-1, - -.14547915466597622925e-1, -.61067668299848923244e-1, - .35093678009090186851e-1, .55378399159800654657e-1, - -.54277226474891610385e-1, -.42023830782434076509e-1, - .69197384645944912066e-1, .22610783557709586445e-1, - -.77269275900637030185e-1, 0.0, .77269275900637030185e-1, - -.22610783557709586445e-1, - -.69197384645944912066e-1, .42023830782434076509e-1, - .54277226474891610385e-1, - -.55378399159800654657e-1, -.35093678009090186851e-1, - .61067668299848923244e-1, .14547915466597622925e-1, - -.58015667638415922967e-1, .31929089844361042178e-2, - .47527977686242313065e-1, -.15480559452075811600e-1, - -.29071311049368895844e-1, .13886961667516983541e-1, - .48780095920069827068e-2, .51591759101720291381e-2, - .11747497650231330965e-1, -.31777863364694653331e-1, - -.34555825499804605557e-2, .47914131921157015198e-1, - -.22573685920142225247e-1, -.45320344390022666738e-1, - .49660630547172186418e-1, .25707858143963615736e-1, - -.68132707341917233933e-1, .67534860185243140399e-2, - .69268150370037450063e-1, -.41585011920451477177e-1, - -.51622397460510041271e-1, .68408139576363036148e-1, - .18981259024768933323e-1, -.78265472429342305554e-1, - .18981259024768933323e-1, .68408139576363036148e-1, - -.51622397460510041271e-1, - -.41585011920451477177e-1, .69268150370037450063e-1, - .67534860185243140399e-2, - -.68132707341917233933e-1, .25707858143963615736e-1, - .49660630547172186418e-1, - -.45320344390022666738e-1, -.22573685920142225247e-1, - .47914131921157015198e-1, -.34555825499804605557e-2, - -.31777863364694653331e-1, .11747497650231330965e-1, - .51591759101720291381e-2, -.54365757412741340377e-2, - -.94862516619529080191e-2, .33240472093448190877e-1, - -.88698898099681552229e-2, - -.40973252097216337576e-1, .42995673349795657065e-1, - .17320914507876958783e-1, - -.62201292691914856803e-1, .24726274174637346693e-1, - .51320859246515407288e-1, - -.62882063373810501763e-1, -.11003569131725622672e-1, - .73842261324108943465e-1, -.39240120294802923208e-1, - -.49293966443941122807e-1, .73552644778818223475e-1, 0.0, - -.73552644778818223475e-1, .49293966443941122807e-1, - .39240120294802923208e-1, -.73842261324108943465e-1, - .11003569131725622672e-1, .62882063373810501763e-1, - -.51320859246515407288e-1, - -.24726274174637346693e-1, .62201292691914856803e-1, - -.17320914507876958783e-1, -.42995673349795657065e-1, - .40973252097216337576e-1, .88698898099681552229e-2, - -.33240472093448190877e-1, .94862516619529080191e-2, - .54365757412741340377e-2, .57750194549356126240e-2, - .69981166020044116791e-2, -.33274982140403110792e-1, - .20297071020698356116e-1, .27898517839646066582e-1, - -.53368678853282030262e-1, .16656482990394548343e-1, - .46342901447260614255e-1, - -.60536796508149003365e-1, .29109107483842596340e-2, - .63224486124385124504e-1, - -.59028872851312033411e-1, -.14783105962696191734e-1, - .74269399241069253865e-1, -.49053677339382384625e-1, - -.33525466624811186739e-1, .78397349622515386647e-1, - -.33525466624811186739e-1, -.49053677339382384625e-1, - .74269399241069253865e-1, -.14783105962696191734e-1, - -.59028872851312033411e-1, - .63224486124385124504e-1, .29109107483842596340e-2, - -.60536796508149003365e-1, - .46342901447260614255e-1, .16656482990394548343e-1, - -.53368678853282030262e-1, - .27898517839646066582e-1, .20297071020698356116e-1, - -.33274982140403110792e-1, - .69981166020044116791e-2, .57750194549356126240e-2, - -.61100308370519200637e-2, -.44383614355738148616e-2, - .32011283412619094811e-1, -.29965011866372897633e-1, - -.10560682331349193348e-1, .51110336443392506342e-1, - -.45012284729681775492e-1, -.94236825555873320102e-2, - .60860695783141264746e-1, - -.55014628647083368926e-1, -.73474782382499482121e-2, - .66640148475243034781e-1, -.62533116045749887988e-1, - -.38650525912400102585e-2, .68429769005837003777e-1, - -.66984505412544901945e-1, 0.0, .66984505412544901945e-1, - -.68429769005837003777e-1, .38650525912400102585e-2, - .62533116045749887988e-1, -.66640148475243034781e-1, - .73474782382499482121e-2, - .55014628647083368926e-1, -.60860695783141264746e-1, - .94236825555873320102e-2, - .45012284729681775492e-1, -.51110336443392506342e-1, - .10560682331349193348e-1, - .29965011866372897633e-1, -.32011283412619094811e-1, - .44383614355738148616e-2, - .61100308370519200637e-2, .65409373892036191538e-2, - .16350101107071157065e-2, -.29301957285983144319e-1, - .36838667173388832579e-1, -.81922703976491586393e-2, - -.36955670021050133434e-1, .58374851095540469865e-1, - -.31977016246946181856e-1, -.25311073698658094646e-1, - .66674413950106952577e-1, - -.54865713324521039571e-1, -.39797027891537985440e-2, - .62830285264808449064e-1, -.72226313251296100676e-1, - .22560232697133353980e-1, .46455784709904033738e-1, - -.78200930751070349956e-1, .46455784709904033738e-1, - .22560232697133353980e-1, -.72226313251296100676e-1, - .62830285264808449064e-1, -.39797027891537985440e-2, - -.54865713324521039571e-1, .66674413950106952577e-1, - -.25311073698658094646e-1, -.31977016246946181856e-1, - .58374851095540469865e-1, -.36955670021050133434e-1, - -.81922703976491586393e-2, .36838667173388832579e-1, - -.29301957285983144319e-1, .16350101107071157065e-2, - .65409373892036191538e-2, -.69686180931868703196e-2, - .11849538727632789870e-2, .25452286414610537766e-1, - -.40522480651713943230e-1, .25694679053362813183e-1, - .14057118113748390637e-1, -.52037614725803488893e-1, - .58849342223684035589e-1, - -.25075229077361409271e-1, -.29559771094034181083e-1, - .68296746944165720199e-1, -.62890462146423984955e-1, - .14457636466274596445e-1, .45787612031322361496e-1, - -.77231759014655809742e-1, .57881203613910543657e-1, 0.0, - -.57881203613910543657e-1, .77231759014655809742e-1, - -.45787612031322361496e-1, -.14457636466274596445e-1, - .62890462146423984955e-1, - -.68296746944165720199e-1, .29559771094034181083e-1, - .25075229077361409271e-1, - -.58849342223684035589e-1, .52037614725803488893e-1, - -.14057118113748390637e-1, -.25694679053362813183e-1, - .40522480651713943230e-1, -.25452286414610537766e-1, - -.11849538727632789870e-2, .69686180931868703196e-2, - .75611653617520254845e-2, -.43290610418608409141e-2, - -.20277062025115566914e-1, - .40362947027704828926e-1, -.38938808024132120254e-1, - .11831186195916702262e-1, - .28476667401744525357e-1, -.59320969056617684621e-1, - .61101629747436200186e-1, - -.29514834848355389223e-1, -.20668001885001084821e-1, - .62923592802445122793e-1, -.73558456263588833115e-1, - .45314556330160999776e-1, .79031645918426015574e-2, - -.58136953576334689357e-1, .78538474524006405758e-1, - -.58136953576334689357e-1, .79031645918426015574e-2, - .45314556330160999776e-1, -.73558456263588833115e-1, - .62923592802445122793e-1, -.20668001885001084821e-1, - -.29514834848355389223e-1, .61101629747436200186e-1, - -.59320969056617684621e-1, .28476667401744525357e-1, - .11831186195916702262e-1, -.38938808024132120254e-1, - .40362947027704828926e-1, -.20277062025115566914e-1, - -.43290610418608409141e-2, .75611653617520254845e-2, - -.81505692478987769484e-2, .74297333588288568430e-2, - .14314212513540223314e-1, -.36711242251332751607e-1, - .46240027755503814626e-1, -.34921532671769023773e-1, - .46930051972353714773e-2, - .32842770336385381562e-1, -.61317813706529588466e-1, - .67000809902468893103e-1, - -.45337449655535622885e-1, .35794459576271920867e-2, - .41830061526027213385e-1, - -.72091371931944711708e-1, .74150028530317793195e-1, - -.46487632538609942002e-1, 0.0, .46487632538609942002e-1, - -.74150028530317793195e-1, .72091371931944711708e-1, - -.41830061526027213385e-1, -.35794459576271920867e-2, - .45337449655535622885e-1, -.67000809902468893103e-1, - .61317813706529588466e-1, -.32842770336385381562e-1, - -.46930051972353714773e-2, .34921532671769023773e-1, - -.46240027755503814626e-1, .36711242251332751607e-1, - -.14314212513540223314e-1, -.74297333588288568430e-2, - .81505692478987769484e-2, .90693182942442189743e-2, - -.11121000903959576737e-1, -.71308296141317458546e-2, - .29219439765986671645e-1, -.45820286629778129593e-1, - .49088381175879124421e-1, -.35614888785023038938e-1, - .78906970900092777895e-2, - .26262843038404929480e-1, -.56143674270125757857e-1, - .71700220472378350694e-1, - -.66963544500697307945e-1, .42215091779892228883e-1, - -.41338867413966866997e-2, -.36164891772995367321e-1, - .66584367783847858225e-1, -.77874712365070098328e-1, - .66584367783847858225e-1, -.36164891772995367321e-1, - -.41338867413966866997e-2, .42215091779892228883e-1, - -.66963544500697307945e-1, - .71700220472378350694e-1, -.56143674270125757857e-1, - .26262843038404929480e-1, - .78906970900092777895e-2, -.35614888785023038938e-1, - .49088381175879124421e-1, - -.45820286629778129593e-1, .29219439765986671645e-1, - -.71308296141317458546e-2, -.11121000903959576737e-1, - .90693182942442189743e-2, -.99848472706332791043e-2, - .14701271465939718856e-1, -.32917820356048383366e-3, - -.19201195309873585230e-1, .38409681836626963278e-1, - -.51647324405878909521e-1, .54522171113149311354e-1, - -.45040302741689006270e-1, .24183738595685990149e-1, - .42204134165479735097e-2, -.34317295181348742251e-1, - .59542472465494579941e-1, -.74135115907618101263e-1, - .74491937840566532596e-1, -.60042604725161994304e-1, - .33437677409000083169e-1, 0.0, - -.33437677409000083169e-1, .60042604725161994304e-1, - -.74491937840566532596e-1, .74135115907618101263e-1, - -.59542472465494579941e-1, .34317295181348742251e-1, - -.42204134165479735097e-2, -.24183738595685990149e-1, - .45040302741689006270e-1, -.54522171113149311354e-1, - .51647324405878909521e-1, -.38409681836626963278e-1, - .19201195309873585230e-1, .32917820356048383366e-3, - -.14701271465939718856e-1, .99848472706332791043e-2, - .11775579274769383373e-1, -.19892153937316935880e-1, - .95335114477449041055e-2, .57661528440359081617e-2, - -.23382690532380910781e-1, .40237257037170725321e-1, - -.53280289903551636474e-1, .59974361806023689068e-1, - -.58701684061992853224e-1, .49033407111597129616e-1, - -.31818835267847249219e-1, .90800541261162098886e-2, - .16272906819312603838e-1, -.40863896581186229487e-1, - .61346046297517367703e-1, - -.74896047554167268919e-1, .79632642148310325817e-1, - -.74896047554167268919e-1, .61346046297517367703e-1, - -.40863896581186229487e-1, .16272906819312603838e-1, - .90800541261162098886e-2, -.31818835267847249219e-1, - .49033407111597129616e-1, -.58701684061992853224e-1, - .59974361806023689068e-1, -.53280289903551636474e-1, - .40237257037170725321e-1, -.23382690532380910781e-1, - .57661528440359081617e-2, .95335114477449041055e-2, - -.19892153937316935880e-1, - .11775579274769383373e-1, -.13562702617218467450e-1, - .24885419969649845849e-1, -.18368693901908875583e-1, - .81673147806084084638e-2, .47890591326129587131e-2, - -.19313752945227974024e-1, .34065953398362954708e-1, - -.47667045133463415672e-1, .58820377816690514309e-1, - -.66424139824618415970e-1, - .69667606260856092515e-1, -.68102459384364543253e-1, - .61683024923302547971e-1, - -.50771943476441639136e-1, .36110771847327189215e-1, - -.18758028464284563358e-1, 0.0, .18758028464284563358e-1, - -.36110771847327189215e-1, .50771943476441639136e-1, - -.61683024923302547971e-1, .68102459384364543253e-1, - -.69667606260856092515e-1, .66424139824618415970e-1, - -.58820377816690514309e-1, .47667045133463415672e-1, - -.34065953398362954708e-1, .19313752945227974024e-1, - -.47890591326129587131e-2, -.81673147806084084638e-2, - .18368693901908875583e-1, -.24885419969649845849e-1, - .13562702617218467450e-1, .20576545037980523979e-1, - -.40093155172981004337e-1, .36954083167944054826e-1, - -.31856506837591907746e-1, .24996323181546255126e-1, - -.16637165210473614136e-1, .71002706773325085237e-2, - .32478629093205201133e-2, - -.14009562579050569518e-1, .24771262248780618922e-1, - -.35119395835433647559e-1, .44656290368574753171e-1, - -.53015448339647394161e-1, .59875631995693046782e-1, - -.64973208326045193862e-1, .68112280331082143373e-1, - -.69172215234062186994e-1, .68112280331082143373e-1, - -.64973208326045193862e-1, .59875631995693046782e-1, - -.53015448339647394161e-1, .44656290368574753171e-1, - -.35119395835433647559e-1, .24771262248780618922e-1, - -.14009562579050569518e-1, .32478629093205201133e-2, - .71002706773325085237e-2, -.16637165210473614136e-1, - .24996323181546255126e-1, -.31856506837591907746e-1, - .36954083167944054826e-1, -.40093155172981004337e-1, - .20576545037980523979e-1, -.27584914609096156163e-1, - .54904171411058497973e-1, -.54109756419563083153e-1, - .52794234894345577483e-1, -.50970276026831042415e-1, - .48655445537990983379e-1, - -.45872036510847994332e-1, .42646854695899611372e-1, - -.39010960357087507670e-1, .34999369144476467749e-1, - -.30650714874402762189e-1, .26006877464703437057e-1, - -.21112579608213651273e-1, .16014956068786763273e-1, - -.10763099747751940252e-1, .54075888924374485533e-2, 0.0, - -.54075888924374485533e-2, .10763099747751940252e-1, - -.16014956068786763273e-1, - .21112579608213651273e-1, -.26006877464703437057e-1, - .30650714874402762189e-1, - -.34999369144476467749e-1, .39010960357087507670e-1, - -.42646854695899611372e-1, .45872036510847994332e-1, - -.48655445537990983379e-1, .50970276026831042415e-1, - -.52794234894345577483e-1, .54109756419563083153e-1, - -.54904171411058497973e-1, .27584914609096156163e-1, - .13794141262469565740e-1, -.27588282524939131481e-1, - .27588282524939131481e-1, -.27588282524939131481e-1, - .27588282524939131481e-1, -.27588282524939131481e-1, - .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .27588282524939131481e-1, - -.27588282524939131481e-1, .13794141262469565740e-1 -}; - -static const double Tleft[33 * 33] = { - 1., -.86602540378443864678, 0., .33071891388307382381, 0., - -.20728904939721249057, 0., .15128841196122722208, 0., - -.11918864298744029244, 0., .98352013661686631224e-1, 0., - -.83727065404940845733e-1, 0., .72893399403505841203e-1, 0., - -.64544632643375022436e-1, 0., .57913170372415565639e-1, 0., - -.52518242575729562263e-1, 0., .48043311993977520457e-1, 0., - -.44271433659733990243e-1, 0., .41048928022856771981e-1, 0., - -.38263878662008271459e-1, 0., .35832844026365304501e-1, 0., 0., - .50000000000000000000, -.96824583655185422130, .57282196186948000082, - .21650635094610966169, -.35903516540862679125, -.97578093724974971969e-1, - .26203921611325660506, .55792409597991015609e-1, -.20644078533943456204, - -.36172381205961199479e-1, .17035068468874958194, - .25371838001497225980e-1, -.14501953125000000000, - -.18786835250972344757e-1, .12625507130328301066, - .14473795929590520582e-1, -.11179458309419422675, - -.11494434254897626155e-1, .10030855351241635862, - .93498556820544479096e-2, -.90964264465390582629e-1, - -.77546391824364392762e-2, .83213457337452292745e-1, - .65358085945588638605e-2, -.76680372422574234569e-1, - -.55835321940047427169e-2, .71098828931825789428e-1, - .48253327982967591019e-2, -.66274981937248958553e-1, - -.42118078245337801387e-2, .62064306433355646267e-1, - .37083386598903548973e-2, 0., 0., .25000000000000000000, - -.73950997288745200531, .83852549156242113615, -.23175620272173946716, - -.37791833195149451496, .25710129174850522325, .21608307321780204633, - -.22844049245646009157, -.14009503000335388415, .19897685605518413847, - .98264706042471226893e-1, -.17445445004279014046, - -.72761100054958328401e-1, .15463589893742108388, - .56056770591708784481e-1, -.13855313872640495158, - -.44517752443294564781e-1, .12534277657695128850, - .36211835346039665762e-1, -.11434398255136139683, - -.30033588409423828125e-1, .10506705408753910481, - .25313077840725783008e-1, -.97149327637744872155e-1, - -.21624927200393328444e-1, .90319582367202122625e-1, - .18688433567711780666e-1, -.84372291635345108584e-1, - -.16312261561845420752e-1, .79149526894804751586e-1, - .14362333871852474757e-1, 0., 0., 0., .12500000000000000000, - -.49607837082461073572, .82265291131801144317, -.59621200088559103072, - -.80054302859059362371e-1, .42612156697795759420, - -.90098145270865592887e-1, -.29769623255090078484, .13630307904779758221, - .21638835185708931831, -.14600247270306082052, -.16348801804014290453, - .14340708728599057249, .12755243353979286190, -.13661523715071346961, - -.10215585947881057394, .12864248070157166547, .83592528025348693602e-1, - -.12066728689302565222, -.69633728678718053052e-1, .11314245177331919532, - .58882939251410088028e-1, -.10621835858758221487, - -.50432266865187597572e-1, .99916834723527771581e-1, - .43672094283057258509e-1, -.94206380251950852413e-1, - -.38181356812697746418e-1, .89035739656537771225e-1, - .33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1, - -.31093357409581873586, .67604086414949799246, -.75644205980613611039, - .28990586430124175741, .30648508196770360914, -.35801372616842500052, - -.91326869828709014708e-1, .31127929687500000000, - -.90915752838698393094e-2, -.25637381283965534330, - .57601077850322797594e-1, .21019685709225757945, - -.81244992138514014256e-1, -.17375078516720988858, - .92289437277967051125e-1, .14527351914265391374, - -.96675340792832019889e-1, -.12289485697108543415, - .97448175340011084006e-1, .10511755943298339844, - -.96242247086378239657e-1, -.90822942272780513537e-1, - .93966350452322132384e-1, .79189411876493712558e-1, - -.91139307067989309325e-1, -.69613039934383197265e-1, - .88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0., - .31250000000000000000e-1, -.18684782411095934408, .50176689760410660236, - -.74784031498626095398, .56472001151566251186, .14842464993721351203e-1, - -.41162920273003120936, .20243071230196532282, .23772054897172750436, - -.24963810923972235950, -.12116179938394678936, .24330535483519110663, - .47903849781124471359e-1, -.22133299683101224293, - -.20542915138527200983e-2, .19653465717678146728, - -.26818172626509178444e-1, -.17319122357631210944, - .45065391411065545445e-1, .15253391395444065941, - -.56543897711725408302e-1, -.13469154928743585367, - .63632471400208840155e-1, .11941684923913523817, - -.67828850207933293098e-1, -.10636309084510652670, - .70095786922999181504e-1, .95187373095150709082e-1, 0., 0., 0., 0., 0., - 0., .15625000000000000000e-1, -.10909562534194485289, - .34842348626527747318, -.64461114561628111443, .69382480527334683659, - -.29551102358528827763, -.25527584713978439819, .38878771718544715394, - -.82956185835347407489e-2, -.31183177761966943912, .12831420840372374767, - .22067618205599434368, -.17569196937129496961, -.14598057000132284135, - .18864406621763419484, .89921002550386645767e-1, -.18571835020187122114, - -.48967672227195481777e-1, .17584685670380332798, - .19267984545067426324e-1, -.16335437520503462738, - .22598055455032407594e-2, .15032800884170631129, - -.17883358353754640871e-1, -.13774837869432209951, - .29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0., - 0., .78125000000000000000e-2, -.62377810244809812496e-1, - .23080781467370883845, -.50841310636012325368, .69834547012574056043, - -.52572723156526459672, .11464215704954976471e-1, .38698869011491210342, - -.26125646622255207507, -.16951698812361607510, .29773875898928782269, - .20130501202570367491e-1, -.26332493149159310198, - .67734613690401207009e-1, .21207315477103762715, -.11541543390889415193, - -.16249634759782417533, .13885887405041735068, .11996491328010275427, - -.14810432001630926895, -.85177658352556243411e-1, .14918860659904380587, - .57317789510444151564e-1, -.14569827645586660151, - -.35213090145965327390e-1, .13975998126844578198, 0., 0., 0., 0., 0., 0., - 0., 0., .39062500000000000000e-2, -.35101954600803571207e-1, - .14761284084133737720, -.37655033076080192966, .62410290231517322776, - -.64335622317683389875, .28188168266139524244, .22488495672137010675, - -.39393811089283576186, .75184777995770096714e-1, .28472023119398293003, - -.20410910833705899572, -.15590046962908511750, .23814567544617953125, - .54442805556829031204e-1, -.22855930338589720954, - .16303223615756629897e-1, .20172722433875559213, - -.62723406421217419404e-1, -.17012230831020922010, - .91754642766136561612e-1, .13927644821381121197, -.10886600968068418181, - -.11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0., - 0., 0., .19531250000000000000e-2, -.19506820659607596598e-1, - .91865676095362231937e-1, -.26604607809696493849, .51425874205091288223, - -.66047561132505329292, .48660109511591303851, -.17575661168678285615e-1, - -.36594333408055703366, .29088854695378694533, .11318677346656537927, - -.31110645235730182168, .60733219161008787341e-1, .24333848233620420826, - -.15254312332655419708, -.15995968483455388613, .19010344455215289289, - .86040636766440260000e-1, -.19652589954665259945, - -.27633388517205837713e-1, .18660848552712880387, - -.15942583868416775867e-1, -.16902042462382064786, - .47278526495327740646e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .97656250000000000000e-3, -.10731084460857378207e-1, - .55939644713816406331e-1, -.18118487371914493668, .39914857299829864263, - -.60812322949933902435, .60011887183061967583, -.26002695805835928795, - -.20883922404786010096, .38988130966114638081, -.11797833550782589082, - -.25231824756239520077, .24817859972953934712, .90516417677868996417e-1, - -.26079073291293066798, .30259468817169480161e-1, .22178195264114178432, - -.10569877864302048175, -.16679648389266977455, .14637718550245050850, - .11219272032739559870, -.16359363640525750353, -.64358194509092101393e-1, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3, - -.58542865274813470967e-2, .33461741635290096452e-1, - -.11979993155896201271, .29580223766987206958, -.51874761979436016742, - .62861483498014306968, -.44868895761051453296, .12567502628371529386e-1, - .35040366183235474275, -.30466868455569500886, -.70903913601490112666e-1, - .30822791893032512740, -.11969443264190207736, -.20764760317621313946, - .20629838355452128532, .95269702915334718507e-1, -.22432624768705133300, - -.33103381593477797101e-2, .20570036048155716333, - -.62208282720094518964e-1, -.17095309330441436348, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., .24414062500000000000e-3, - -.31714797501871532475e-2, .19721062526127334100e-1, - -.77311181185536498246e-1, .21124871792841566575, -.41777980401893650886, - .59401977834943551650, -.56132417807488349048, .23433675061367565951, - .20222775295220942126, -.38280372496506190127, .14443804214023095767, - .22268950939178466797, -.27211314150777981984, -.34184876506180717313e-1, - .26006498895669734842, -.97650425186005090107e-1, -.19024527660129101293, - .16789164198044635671, .10875811641651905252, -.19276785058805921298, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3, - -.17078941137247586143e-2, .11477733754843910060e-1, - -.48887017020924625462e-1, .14634927241421789683, -.32156282683019547854, - .52165811920227223937, -.60001958466396926460, .41208501541480733755, - -.11366945503190350975e-2, -.33968093962672089159, .30955190935923386766, - .40657421856578262210e-1, -.29873400409871531764, .16094481791768257440, - .16876122436206497694, -.23650217045022161255, -.33070260090574765012e-1, - .22985258456375907796, -.68645651043827097771e-1, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4, - -.91501857608428649078e-3, .66085179496951987952e-2, - -.30383171695850355404e-1, .98840838845366876117e-1, - -.23855447246420318989, .43322017468145613917, -.58049033744876107191, - .52533893203742699346, -.20681056202371946180, -.20180000924562504384, - .37503922291962681797, -.15988102869837429062, -.19823558102762374094, - .28393023878803799622, -.11188133439357510403e-1, -.24730368377168229255, - .14731529061377942839, .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .30517578125000000000e-4, - -.48804277318479845551e-3, .37696080990601968396e-2, - -.18603912108994738255e-1, .65325006755649582964e-1, - -.17162960707938819795, .34411527956476971322, -.52289350347082497959, - .57319653625674910592, -.37662253421045430413, -.14099055105384663902e-1, - .33265570610216904208, -.30921265572647566661, -.19911390594166455281e-1, - .28738590811031797718, -.18912130469738472647, -.13235936203215819193, - .25076406142356675279, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., .15258789062500000000e-4, -.25928719280954633249e-3, - .21327398937568540428e-2, -.11244626133630732010e-1, - .42375605740664331966e-1, -.12031130345907846211, .26352562258934426830, - -.44590628258512682078, .56682835613700749379, -.49116715128261660395, - .17845943097110339078, .20541650677432497477, -.36739803642257458221, - .16776034069210108273, .17920950989905112908, -.28867732805385066532, - .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., .76293945312500000000e-5, -.13727610943181290891e-3, - .11979683091449349286e-2, -.67195313034570709806e-2, - .27044920779931968175e-1, -.82472196498517457862e-1, - .19570475044896150093, -.36391620788543817693, .52241392782736588032, - -.54727504974907879912, .34211551468813581183, .31580472732719957762e-1, - -.32830006549176759667, .30563797665254420769, .64905014620683140120e-2, - -.27642986248995073032, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., .38146972656250000000e-5, - -.72454147007837596854e-4, .66859847582761390285e-3, - -.39751311980366118437e-2, .17015198650201528366e-1, - -.55443621868993855715e-1, .14157060481641692131, -.28641242619559616836, - .45610665490966615415, -.55262786406029265394, .45818352706035500108, - -.14984403004611673047, -.21163807462970713245, .36007252928843413718, - -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5, - -.38135049864067468562e-4, .37101393638555730015e-3, - -.23305339886279723213e-2, .10569913448297127219e-1, - -.36640175162216897547e-1, .10010476414320235508, -.21860074212675559892, - .38124757096345313719, -.52020999209879669177, .52172632730659212045, - -.30841620620308814614, -.50322546186721500184e-1, .32577618885114899053, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .95367431640625000000e-6, -.20021483206955925244e-4, - .20481807322420625431e-3, -.13553476938058909882e-2, - .64919676350791905019e-2, -.23848725425069251903e-1, - .69384632678886421292e-1, -.16249711393618776934, .30736618106830314788, - -.46399909601971539157, .53765031034002467225, -.42598991476520183929, - .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6, - -.10487707828484902486e-4, .11254146162337528943e-3, - -.78248929534271987118e-3, .39468337145306794566e-2, - -.15313546659475671763e-1, .47249070825218564146e-1, - -.11804374107101480543, .24031796927792491122, -.39629215049166341285, - .51629108968402548545, -.49622372075429782915, 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .23841857910156250000e-6, -.54823314130625337326e-5, - .61575377321535518154e-4, -.44877834366497538134e-3, - .23774612048621955857e-2, -.97136347645161687796e-2, - .31671599547606636717e-1, -.84028665767000747480e-1, - .18298487576742964949, -.32647878537696945218, .46970971486488895077, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., .11920928955078125000e-6, -.28604020001177375838e-5, - .33559227978295551013e-4, -.25583821662860610560e-3, - .14201552747787302339e-2, -.60938046986874414969e-2, - .20930869247951926793e-1, -.58745021125678072911e-1, - .13613725780285953720, -.26083988356030237586, 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .59604644775390625000e-7, -.14898180663526043291e-5, - .18224991282807693921e-4, -.14504433444608833821e-3, - .84184722720281809548e-3, -.37846965430000478789e-2, - .13656355548211376864e-1, -.40409541997718853934e-1, - .99226988101858325902e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .29802322387695312500e-7, -.77471708843445529468e-6, - .98649879372606876995e-5, -.81814934772838523887e-4, - .49554483992403011328e-3, -.23290922072351413938e-2, - .88068134250844034186e-2, -.27393666952485719070e-1, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., .14901161193847656250e-7, -.40226235946098233685e-6, - .53236418690561306700e-5, -.45933829691164002269e-4, - .28982005232838857913e-3, -.14212974043211018374e-2, - .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .74505805969238281250e-8, -.20858299254133430408e-6, - .28648457300134381744e-5, -.25677535898258910850e-4, - .16849420429491355445e-3, -.86062824010315834002e-3, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., .37252902984619140625e-8, -.10801736017613096861e-6, - .15376606719887104015e-5, -.14296523739727437959e-4, - .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .18626451492309570312e-8, -.55871592916438890146e-7, - .82331193828137454068e-6, -.79302250528382787666e-5, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9, - -.28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9, - -.14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .23283064365386962891e-9 -}; - -static const double Tright[33 * 33] = { - 1., .86602540378443864678, 0., -.33071891388307382381, 0., - .20728904939721249057, 0., -.15128841196122722208, 0., - .11918864298744029244, 0., -.98352013661686631224e-1, 0., - .83727065404940845733e-1, 0., -.72893399403505841203e-1, 0., - .64544632643375022436e-1, 0., -.57913170372415565639e-1, 0., - .52518242575729562263e-1, 0., -.48043311993977520457e-1, 0., - .44271433659733990243e-1, 0., -.41048928022856771981e-1, 0., - .38263878662008271459e-1, 0., -.35832844026365304501e-1, 0., 0., - .50000000000000000000, .96824583655185422130, .57282196186948000082, - -.21650635094610966169, -.35903516540862679125, .97578093724974971969e-1, - .26203921611325660506, -.55792409597991015609e-1, -.20644078533943456204, - .36172381205961199479e-1, .17035068468874958194, - -.25371838001497225980e-1, -.14501953125000000000, - .18786835250972344757e-1, .12625507130328301066, - -.14473795929590520582e-1, -.11179458309419422675, - .11494434254897626155e-1, .10030855351241635862, - -.93498556820544479096e-2, -.90964264465390582629e-1, - .77546391824364392762e-2, .83213457337452292745e-1, - -.65358085945588638605e-2, -.76680372422574234569e-1, - .55835321940047427169e-2, .71098828931825789428e-1, - -.48253327982967591019e-2, -.66274981937248958553e-1, - .42118078245337801387e-2, .62064306433355646267e-1, - -.37083386598903548973e-2, 0., 0., .25000000000000000000, - .73950997288745200531, .83852549156242113615, .23175620272173946716, - -.37791833195149451496, -.25710129174850522325, .21608307321780204633, - .22844049245646009157, -.14009503000335388415, -.19897685605518413847, - .98264706042471226893e-1, .17445445004279014046, - -.72761100054958328401e-1, -.15463589893742108388, - .56056770591708784481e-1, .13855313872640495158, - -.44517752443294564781e-1, -.12534277657695128850, - .36211835346039665762e-1, .11434398255136139683, - -.30033588409423828125e-1, -.10506705408753910481, - .25313077840725783008e-1, .97149327637744872155e-1, - -.21624927200393328444e-1, -.90319582367202122625e-1, - .18688433567711780666e-1, .84372291635345108584e-1, - -.16312261561845420752e-1, -.79149526894804751586e-1, - .14362333871852474757e-1, 0., 0., 0., .12500000000000000000, - .49607837082461073572, .82265291131801144317, .59621200088559103072, - -.80054302859059362371e-1, -.42612156697795759420, - -.90098145270865592887e-1, .29769623255090078484, .13630307904779758221, - -.21638835185708931831, -.14600247270306082052, .16348801804014290453, - .14340708728599057249, -.12755243353979286190, -.13661523715071346961, - .10215585947881057394, .12864248070157166547, -.83592528025348693602e-1, - -.12066728689302565222, .69633728678718053052e-1, .11314245177331919532, - -.58882939251410088028e-1, -.10621835858758221487, - .50432266865187597572e-1, .99916834723527771581e-1, - -.43672094283057258509e-1, -.94206380251950852413e-1, - .38181356812697746418e-1, .89035739656537771225e-1, - -.33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1, - .31093357409581873586, .67604086414949799246, .75644205980613611039, - .28990586430124175741, -.30648508196770360914, -.35801372616842500052, - .91326869828709014708e-1, .31127929687500000000, .90915752838698393094e-2, - -.25637381283965534330, -.57601077850322797594e-1, .21019685709225757945, - .81244992138514014256e-1, -.17375078516720988858, - -.92289437277967051125e-1, .14527351914265391374, - .96675340792832019889e-1, -.12289485697108543415, - -.97448175340011084006e-1, .10511755943298339844, - .96242247086378239657e-1, -.90822942272780513537e-1, - -.93966350452322132384e-1, .79189411876493712558e-1, - .91139307067989309325e-1, -.69613039934383197265e-1, - -.88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0., - .31250000000000000000e-1, .18684782411095934408, .50176689760410660236, - .74784031498626095398, .56472001151566251186, -.14842464993721351203e-1, - -.41162920273003120936, -.20243071230196532282, .23772054897172750436, - .24963810923972235950, -.12116179938394678936, -.24330535483519110663, - .47903849781124471359e-1, .22133299683101224293, - -.20542915138527200983e-2, -.19653465717678146728, - -.26818172626509178444e-1, .17319122357631210944, - .45065391411065545445e-1, -.15253391395444065941, - -.56543897711725408302e-1, .13469154928743585367, - .63632471400208840155e-1, -.11941684923913523817, - -.67828850207933293098e-1, .10636309084510652670, - .70095786922999181504e-1, -.95187373095150709082e-1, 0., 0., 0., 0., 0., - 0., .15625000000000000000e-1, .10909562534194485289, - .34842348626527747318, .64461114561628111443, .69382480527334683659, - .29551102358528827763, -.25527584713978439819, -.38878771718544715394, - -.82956185835347407489e-2, .31183177761966943912, .12831420840372374767, - -.22067618205599434368, -.17569196937129496961, .14598057000132284135, - .18864406621763419484, -.89921002550386645767e-1, -.18571835020187122114, - .48967672227195481777e-1, .17584685670380332798, - -.19267984545067426324e-1, -.16335437520503462738, - -.22598055455032407594e-2, .15032800884170631129, - .17883358353754640871e-1, -.13774837869432209951, - -.29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0., - 0., .78125000000000000000e-2, .62377810244809812496e-1, - .23080781467370883845, .50841310636012325368, .69834547012574056043, - .52572723156526459672, .11464215704954976471e-1, -.38698869011491210342, - -.26125646622255207507, .16951698812361607510, .29773875898928782269, - -.20130501202570367491e-1, -.26332493149159310198, - -.67734613690401207009e-1, .21207315477103762715, .11541543390889415193, - -.16249634759782417533, -.13885887405041735068, .11996491328010275427, - .14810432001630926895, -.85177658352556243411e-1, -.14918860659904380587, - .57317789510444151564e-1, .14569827645586660151, - -.35213090145965327390e-1, -.13975998126844578198, 0., 0., 0., 0., 0., 0., - 0., 0., .39062500000000000000e-2, .35101954600803571207e-1, - .14761284084133737720, .37655033076080192966, .62410290231517322776, - .64335622317683389875, .28188168266139524244, -.22488495672137010675, - -.39393811089283576186, -.75184777995770096714e-1, .28472023119398293003, - .20410910833705899572, -.15590046962908511750, -.23814567544617953125, - .54442805556829031204e-1, .22855930338589720954, .16303223615756629897e-1, - -.20172722433875559213, -.62723406421217419404e-1, .17012230831020922010, - .91754642766136561612e-1, -.13927644821381121197, -.10886600968068418181, - .11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0., - 0., 0., .19531250000000000000e-2, .19506820659607596598e-1, - .91865676095362231937e-1, .26604607809696493849, .51425874205091288223, - .66047561132505329292, .48660109511591303851, .17575661168678285615e-1, - -.36594333408055703366, -.29088854695378694533, .11318677346656537927, - .31110645235730182168, .60733219161008787341e-1, -.24333848233620420826, - -.15254312332655419708, .15995968483455388613, .19010344455215289289, - -.86040636766440260000e-1, -.19652589954665259945, - .27633388517205837713e-1, .18660848552712880387, .15942583868416775867e-1, - -.16902042462382064786, -.47278526495327740646e-1, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., .97656250000000000000e-3, .10731084460857378207e-1, - .55939644713816406331e-1, .18118487371914493668, .39914857299829864263, - .60812322949933902435, .60011887183061967583, .26002695805835928795, - -.20883922404786010096, -.38988130966114638081, -.11797833550782589082, - .25231824756239520077, .24817859972953934712, -.90516417677868996417e-1, - -.26079073291293066798, -.30259468817169480161e-1, .22178195264114178432, - .10569877864302048175, -.16679648389266977455, -.14637718550245050850, - .11219272032739559870, .16359363640525750353, -.64358194509092101393e-1, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3, - .58542865274813470967e-2, .33461741635290096452e-1, .11979993155896201271, - .29580223766987206958, .51874761979436016742, .62861483498014306968, - .44868895761051453296, .12567502628371529386e-1, -.35040366183235474275, - -.30466868455569500886, .70903913601490112666e-1, .30822791893032512740, - .11969443264190207736, -.20764760317621313946, -.20629838355452128532, - .95269702915334718507e-1, .22432624768705133300, - -.33103381593477797101e-2, -.20570036048155716333, - -.62208282720094518964e-1, .17095309330441436348, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., .24414062500000000000e-3, - .31714797501871532475e-2, .19721062526127334100e-1, - .77311181185536498246e-1, .21124871792841566575, .41777980401893650886, - .59401977834943551650, .56132417807488349048, .23433675061367565951, - -.20222775295220942126, -.38280372496506190127, -.14443804214023095767, - .22268950939178466797, .27211314150777981984, -.34184876506180717313e-1, - -.26006498895669734842, -.97650425186005090107e-1, .19024527660129101293, - .16789164198044635671, -.10875811641651905252, -.19276785058805921298, 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3, - .17078941137247586143e-2, .11477733754843910060e-1, - .48887017020924625462e-1, .14634927241421789683, .32156282683019547854, - .52165811920227223937, .60001958466396926460, .41208501541480733755, - .11366945503190350975e-2, -.33968093962672089159, -.30955190935923386766, - .40657421856578262210e-1, .29873400409871531764, .16094481791768257440, - -.16876122436206497694, -.23650217045022161255, .33070260090574765012e-1, - .22985258456375907796, .68645651043827097771e-1, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4, - .91501857608428649078e-3, .66085179496951987952e-2, - .30383171695850355404e-1, .98840838845366876117e-1, .23855447246420318989, - .43322017468145613917, .58049033744876107191, .52533893203742699346, - .20681056202371946180, -.20180000924562504384, -.37503922291962681797, - -.15988102869837429062, .19823558102762374094, .28393023878803799622, - .11188133439357510403e-1, -.24730368377168229255, -.14731529061377942839, - .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .30517578125000000000e-4, .48804277318479845551e-3, - .37696080990601968396e-2, .18603912108994738255e-1, - .65325006755649582964e-1, .17162960707938819795, .34411527956476971322, - .52289350347082497959, .57319653625674910592, .37662253421045430413, - -.14099055105384663902e-1, -.33265570610216904208, -.30921265572647566661, - .19911390594166455281e-1, .28738590811031797718, .18912130469738472647, - -.13235936203215819193, -.25076406142356675279, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .15258789062500000000e-4, - .25928719280954633249e-3, .21327398937568540428e-2, - .11244626133630732010e-1, .42375605740664331966e-1, .12031130345907846211, - .26352562258934426830, .44590628258512682078, .56682835613700749379, - .49116715128261660395, .17845943097110339078, -.20541650677432497477, - -.36739803642257458221, -.16776034069210108273, .17920950989905112908, - .28867732805385066532, .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .76293945312500000000e-5, - .13727610943181290891e-3, .11979683091449349286e-2, - .67195313034570709806e-2, .27044920779931968175e-1, - .82472196498517457862e-1, .19570475044896150093, .36391620788543817693, - .52241392782736588032, .54727504974907879912, .34211551468813581183, - -.31580472732719957762e-1, -.32830006549176759667, -.30563797665254420769, - .64905014620683140120e-2, .27642986248995073032, 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .38146972656250000000e-5, - .72454147007837596854e-4, .66859847582761390285e-3, - .39751311980366118437e-2, .17015198650201528366e-1, - .55443621868993855715e-1, .14157060481641692131, .28641242619559616836, - .45610665490966615415, .55262786406029265394, .45818352706035500108, - .14984403004611673047, -.21163807462970713245, -.36007252928843413718, - -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5, - .38135049864067468562e-4, .37101393638555730015e-3, - .23305339886279723213e-2, .10569913448297127219e-1, - .36640175162216897547e-1, .10010476414320235508, .21860074212675559892, - .38124757096345313719, .52020999209879669177, .52172632730659212045, - .30841620620308814614, -.50322546186721500184e-1, -.32577618885114899053, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .95367431640625000000e-6, .20021483206955925244e-4, - .20481807322420625431e-3, .13553476938058909882e-2, - .64919676350791905019e-2, .23848725425069251903e-1, - .69384632678886421292e-1, .16249711393618776934, .30736618106830314788, - .46399909601971539157, .53765031034002467225, .42598991476520183929, - .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6, - .10487707828484902486e-4, .11254146162337528943e-3, - .78248929534271987118e-3, .39468337145306794566e-2, - .15313546659475671763e-1, .47249070825218564146e-1, .11804374107101480543, - .24031796927792491122, .39629215049166341285, .51629108968402548545, - .49622372075429782915, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., .23841857910156250000e-6, - .54823314130625337326e-5, .61575377321535518154e-4, - .44877834366497538134e-3, .23774612048621955857e-2, - .97136347645161687796e-2, .31671599547606636717e-1, - .84028665767000747480e-1, .18298487576742964949, .32647878537696945218, - .46970971486488895077, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .11920928955078125000e-6, - .28604020001177375838e-5, .33559227978295551013e-4, - .25583821662860610560e-3, .14201552747787302339e-2, - .60938046986874414969e-2, .20930869247951926793e-1, - .58745021125678072911e-1, .13613725780285953720, .26083988356030237586, - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., .59604644775390625000e-7, - .14898180663526043291e-5, .18224991282807693921e-4, - .14504433444608833821e-3, .84184722720281809548e-3, - .37846965430000478789e-2, .13656355548211376864e-1, - .40409541997718853934e-1, .99226988101858325902e-1, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .29802322387695312500e-7, .77471708843445529468e-6, - .98649879372606876995e-5, .81814934772838523887e-4, - .49554483992403011328e-3, .23290922072351413938e-2, - .88068134250844034186e-2, .27393666952485719070e-1, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., .14901161193847656250e-7, .40226235946098233685e-6, - .53236418690561306700e-5, .45933829691164002269e-4, - .28982005232838857913e-3, .14212974043211018374e-2, - .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .74505805969238281250e-8, .20858299254133430408e-6, - .28648457300134381744e-5, .25677535898258910850e-4, - .16849420429491355445e-3, .86062824010315834002e-3, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., .37252902984619140625e-8, .10801736017613096861e-6, - .15376606719887104015e-5, .14296523739727437959e-4, - .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - .18626451492309570312e-8, .55871592916438890146e-7, - .82331193828137454068e-6, .79302250528382787666e-5, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9, - .28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9, - .14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., .23283064365386962891e-9 -}; - -/* Allocates a workspace for the given maximum number of intervals. - Note that if the workspace gets filled, the intervals with the - lowest error estimates are dropped. The maximum number of - intervals is therefore not the maximum number of intervals - that will be computed, but merely the size of the buffer. - */ - -/* Compute the product of the fx with one of the inverse - Vandermonde-like matrices. */ - -void -Vinvfx (const double *fx, double *c, const int d) -{ - - int i, j; - - switch (d) - { - case 0: - for (i = 0; i <= 4; i++) - { - c[i] = 0.0; - for (j = 0; j <= 4; j++) - c[i] += V1inv[i * 5 + j] * fx[j * 8]; - } - break; - case 1: - for (i = 0; i <= 8; i++) - { - c[i] = 0.0; - for (j = 0; j <= 8; j++) - c[i] += V2inv[i * 9 + j] * fx[j * 4]; - } - break; - case 2: - for (i = 0; i <= 16; i++) - { - c[i] = 0.0; - for (j = 0; j <= 16; j++) - c[i] += V3inv[i * 17 + j] * fx[j * 2]; - } - break; - case 3: - for (i = 0; i <= 32; i++) - { - c[i] = 0.0; - for (j = 0; j <= 32; j++) - c[i] += V4inv[i * 33 + j] * fx[j]; - } - break; - } - -} - - -/* Downdate the interpolation given by the n coefficients c - by removing the nodes with indices in nans. */ - -void -downdate (double *c, int n, int d, int *nans, int nnans) -{ - - static const int bidx[4] = { 0, 6, 16, 34 }; - double b_new[34], alpha; - int i, j; - - for (i = 0; i <= n + 1; i++) - b_new[i] = bee[bidx[d] + i]; - for (i = 0; i < nnans; i++) - { - b_new[n + 1] = b_new[n + 1] / Lalpha[n]; - b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1]; - for (j = n - 1; j > 0; j--) - b_new[j] = - (b_new[j] + xi[nans[i]] * b_new[j + 1] - - Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1]; - for (j = 0; j <= n; j++) - b_new[j] = b_new[j + 1]; - alpha = c[n] / b_new[n]; - for (j = 0; j < n; j++) - c[j] -= alpha * b_new[j]; - c[n] = 0; - n--; - } - -} - - - -/* The actual integration routine. - */ - -DEFUN_DLD (cquad, args, nargout, -"-*- texinfo -*-\n\ -@deftypefn {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} cquad (@var{f}, @var{a}, @var{b}, @var{tol})\n\ -@deftypefnx {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} cquad (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ -Numerically evaluates an integral using the doubly-adaptive\n\ -quadrature described by P. Gonnet in @cite{Increasing the\n\ -Reliability of Adaptive Quadrature Using Explicit Interpolants,\n\ -ACM Transactions on Mathematical Software, in Press, 2010.\n\ -The algorithm uses Clenshaw-Curtis quadrature rules of increasing\n\ -degree in each interval and bisects the interval if either the\n\ -function does not appear to be smooth or a rule of maximum\n\ -degree has been reached. The error estimate is computed from the\n\ -L2-norm of the difference between two successive interpolations\n\ -of the integrand over the nodes of the respective quadrature rules.\n\ -\n\ -For example,\n\ -\n\ -@example\n\ - int = cquad (f, a, b, 1.0e-6);\n\ -@end example\n\ -\n\ -@noindent\n\ -computes the integral of a function @var{f} in the interval\n\ -[@var{a}, @var{b}] to the relative precision of six\n\ -decimal digits.\n\ -The integrand @var{f} should accept a vector argument and return a vector\n\ -result containing the integrand evaluated at each element of the\n\ -argument, for example:\n\ -\n\ -@example\n\ - f = @@(x) x .* sin (1 ./ x) .* sqrt (abs (1 - x));\n\ -@end example\n\ -\n\ -If the integrand has known singularities or discontinuities\n\ -in any of its derivatives inside the interval,\n\ -as does the above example at x=1, these can be specified in\n\ -the additional argument @var{sing} as follows\n\ -\n\ -@example\n\ - int = cquad (f, a, b, 1.0e-6, [ 1 ]);\n\ -@end example\n\ -\n\ -The two additional output variables @var{err} and @var{nr_points}\n\ -return an estimate of the absolute integration error and\n\ -the number of points at which the integrand was evaluated\n\ -respectively.\n\ -If the adaptive integration did not converge, the value of\n\ -@var{err} will be larger than the requested tolerance. It is\n\ -therefore recommended to verify this value for difficult\n\ -integrands.\n\ -\n\ -If either @var{a} or @var{b} are @code{+/-Inf}, @code{cquad}\n\ -integrates @var{f} by substituting the variable of integration\n\ -with @code{x=tan(pi/2*u)}.\n\ -\n\ -@code{cquad} is capable of dealing with non-numeric\n\ -values of the integrand such as @code{NaN}, @code{Inf}\n\ -or @code{-Inf}, as in the above example at x=0.\n\ -If the integral diverges and @code{cquad} detects this, \n\ -a warning is issued and @code{Inf} or @code{-Inf} is returned.\n\ -\n\ -Note that @code{cquad} is a general purpose quadrature algorithm\n\ -and as such may be less efficient for smooth or otherwise\n\ -well-behaved integrand than other methods such as\n\ -@code{quadgk} or @code{trapz}.\n\ -\n\ -@seealso{quad,quadv,quadl,quadgk,trapz}\n\ -@end deftypefn") -{ - - /* Some constants that we will need. */ - static const int n[4] = { 4, 8, 16, 32 }; - static const int skip[4] = { 8, 4, 2, 1 }; - static const int idx[4] = { 0, 5, 14, 31 }; - static const double w = M_SQRT2 / 2; - static const int ndiv_max = 20; - - /* The interval heap. */ - cquad_ival ivals[cquad_heapsize]; - int heap[cquad_heapsize]; - - /* Arguments left and right */ - int nargin = args.length (); - octave_function *fcn; - double a, b, tol, iivals[cquad_heapsize], *sing; - - /* Variables needed for transforming the integrand. */ - int wrap = 0; - double xw; - - /* Stuff we will need to call the integrand. */ - octave_value_list fargs, retval; - - /* Actual variables (as opposed to constants above). */ - double m, h, ml, hl, mr, hr, temp; - double igral, err, igral_final, err_final, err_excess; - int nivals, neval = 0; - int i, j, d, split, t; - int nnans, nans[33]; - cquad_ival *iv, *ivl, *ivr; - double nc, ncdiff; - - - /* Parse the input arguments. */ - if (nargin < 1) - { - error - ("cquad: first argument (integrand) of type function handle required."); - return octave_value_list (); - } - else - { - if (args (0).is_function_handle () || args (0).is_inline_function ()) - fcn = args (0).function_value (); - else { - error - ("cquad: first argument (integrand) must be a function handle or an inline function."); - return octave_value_list(); - } - } - - if (nargin < 2 || !args (1).is_real_scalar ()) - { - error - ("cquad: second argument (left interval edge) must be a single real scalar."); - return octave_value_list (); - } - else - a = args (1).double_value (); - - if (nargin < 3 || !args (2).is_real_scalar ()) - { - error - ("cquad: third argument (right interval edge) must be a single real scalar."); - return octave_value_list (); - } - else - b = args (2).double_value (); - - if (nargin < 4) - tol = 1.0e-6; - else if (!args (3).is_real_scalar ()) - { - error - ("cquad: fourth argument (tolerance) must be a single real scalar."); - return octave_value_list (); - } - else - tol = args (3).double_value (); - - if (nargin < 5) - { - nivals = 1; - iivals[0] = a; - iivals[1] = b; - } - else if (!(args (4).is_real_scalar () || args (4).is_real_matrix ())) - { - error - ("cquad: fifth argument (singularities) must be a vector of real values."); - return octave_value_list (); - } - else - { - nivals = 1 + args (4).length (); - if ( nivals > cquad_heapsize ) { - error("cquad: maximum number of singular points is limited to %i.",cquad_heapsize-1); - return octave_value_list(); - } - sing = args (4).array_value ().fortran_vec (); - iivals[0] = a; - for (i = 0; i < nivals - 2; i++) - iivals[i + 1] = sing[i]; - iivals[nivals] = b; - } - - /* If a or b are +/-Inf, transform the integral. */ - if (xisinf (a) || xisinf (b)) - { - wrap = 1; - for (i = 0; i <= nivals; i++) - if (xisinf (iivals[i])) - iivals[i] = copysign (1.0, iivals[i]); - else - iivals[i] = 2.0 * atan (iivals[i]) / M_PI; - } - - - /* Initialize the heaps. */ - for (i = 0; i < cquad_heapsize; i++) - heap[i] = i; - - - /* Create the first interval(s). */ - igral = 0.0; - err = 0.0; - for (j = 0; j < nivals; j++) - { - - /* Initialize the interval. */ - iv = &(ivals[heap[j]]); - m = (iivals[j] + iivals[j + 1]) / 2; - h = (iivals[j + 1] - iivals[j]) / 2; - nnans = 0; - ColumnVector ex (33); - if (wrap) - { - for (i = 0; i <= n[3]; i++) - ex (i) = tan (M_PI / 2 * (m + xi[i] * h)); - } - else - { - for (i = 0; i <= n[3]; i++) - ex (i) = m + xi[i] * h; - } - fargs (0) = ex; - retval = feval (fcn, fargs, 1); - if (retval.length () != 1 || !retval (0).is_real_matrix ()) - { - error - ("cquad: integrand must return a single, real-valued vector."); - return octave_value_list (); - } - Matrix effex = retval (0).matrix_value (); - if (effex.length () != ex.length ()) - { - error - ("cquad: integrand must return a single, real-valued vector of the same size as the input"); - return octave_value_list (); - } - for (i = 0; i <= n[3]; i++) - { - iv->fx[i] = effex (i); - if (wrap) - { - xw = ex(i); - iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2; - } - neval++; - if (!xfinite (iv->fx[i])) - { - nans[nnans++] = i; - iv->fx[i] = 0.0; - } - } - Vinvfx (iv->fx, &(iv->c[idx[3]]), 3); - Vinvfx (iv->fx, &(iv->c[idx[2]]), 2); - Vinvfx (iv->fx, &(iv->c[0]), 0); - for (i = 0; i < nnans; i++) - iv->fx[i] = octave_NaN; - iv->a = iivals[j]; - iv->b = iivals[j + 1]; - iv->depth = 3; - iv->rdepth = 1; - iv->ndiv = 0; - iv->igral = 2 * h * iv->c[idx[3]] * w; - nc = 0.0; - for (i = n[2] + 1; i <= n[3]; i++) - { - temp = iv->c[idx[3] + i]; - nc += temp * temp; - } - ncdiff = nc; - for (i = 0; i <= n[2]; i++) - { - temp = iv->c[idx[2] + i] - iv->c[idx[3] + i]; - ncdiff += temp * temp; - nc += iv->c[idx[3] + i] * iv->c[idx[3] + i]; - } - ncdiff = sqrt (ncdiff); - nc = sqrt (nc); - iv->err = ncdiff * 2 * h; - if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc) - iv->err = 2 * h * nc; - - /* Tabulate this interval's data. */ - igral += iv->igral; - err += iv->err; - - /* Sift it up the heap. */ - i = j; - while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err) - { - temp = heap[i]; - heap[i] = heap[i / 2]; - heap[i / 2] = temp; - i /= 2; - } - - } - - - /* Initialize some global values. */ - igral_final = 0.0; - err_final = 0.0; - err_excess = 0.0; - - - /* Main loop. */ - while (nivals > 0 && err > 0.0 && err > fabs (igral) * tol - && !(err_final > fabs (igral) * tol - && err - err_final < fabs (igral) * tol)) - { - - /* Allow the user to interrupt. */ - OCTAVE_QUIT; - - /* Put our finger on the interval with the largest error. */ - iv = &(ivals[heap[0]]); - m = (iv->a + iv->b) / 2; - h = (iv->b - iv->a) / 2; - -/* printf - ("cquad: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", - heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth); -*/ - /* Should we try to increase the degree? */ - if (iv->depth < 3) - { - - /* Keep tabs on some variables. */ - d = ++iv->depth; - - /* Get the new (missing) function values */ - { - ColumnVector ex (n[d] / 2); - if (wrap) - { - for (i = 0; i < n[d] / 2; i++) - ex (i) = - tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h)); - } - else - { - for (i = 0; i < n[d] / 2; i++) - ex (i) = m + xi[(2 * i + 1) * skip[d]] * h; - } - fargs (0) = ex; - retval = feval (fcn, fargs, 1); - if (retval.length () != 1 || !retval (0).is_real_matrix ()) - { - error - ("cquad: integrand must return a single, real-valued vector."); - return octave_value_list (); - } - Matrix effex = retval (0).matrix_value (); - if (effex.length () != ex.length ()) - { - error - ("cquad: integrand must return a single, real-valued vector of the same size as the input."); - return octave_value_list (); - } - neval += effex.length (); - for (i = 0; i < n[d] / 2; i++) - { - j = (2 * i + 1) * skip[d]; - iv->fx[j] = effex (i); - if (wrap) - { - xw = ex(i); - iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2; - } - } - } - nnans = 0; - for (i = 0; i <= 32; i += skip[d]) - { - if (!xfinite (iv->fx[i])) - { - nans[nnans++] = i; - iv->fx[i] = 0.0; - } - } - - /* Compute the new coefficients. */ - Vinvfx (iv->fx, &(iv->c[idx[d]]), d); - /* Downdate any NaNs. */ - if (nnans > 0) - { - downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans); - for (i = 0; i < nnans; i++) - iv->fx[i] = octave_NaN; - } - - /* Compute the error estimate. */ - nc = 0.0; - for (i = n[d - 1] + 1; i <= n[d]; i++) - { - temp = iv->c[idx[d] + i]; - nc += temp * temp; - } - ncdiff = nc; - for (i = 0; i <= n[d - 1]; i++) - { - temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i]; - ncdiff += temp * temp; - nc += iv->c[idx[d] + i] * iv->c[idx[d] + i]; - } - ncdiff = sqrt (ncdiff); - nc = sqrt (nc); - iv->err = ncdiff * 2 * h; - /* Compute the local integral. */ - iv->igral = 2 * h * w * iv->c[idx[d]]; - /* Split the interval prematurely? */ - split = (nc > 0 && ncdiff / nc > 0.1); - } - - /* Maximum degree reached, just split. */ - else - { - split = 1; - } - - - /* Should we drop this interval? */ - if ((m + h * xi[0]) >= (m + h * xi[1]) - || (m + h * xi[31]) >= (m + h * xi[32]) - || iv->err < fabs (iv->igral) * DBL_EPSILON * 10) - { - -/* printf - ("cquad: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", - heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, - iv->depth); -*/ - /* Keep this interval's contribution */ - err_final += iv->err; - igral_final += iv->igral; - /* Swap with the last element on the heap */ - t = heap[nivals - 1]; - heap[nivals - 1] = heap[0]; - heap[0] = t; - nivals--; - /* Fix up the heap */ - i = 0; - while (2 * i + 1 < nivals) - { - - /* Get the kids */ - j = 2 * i + 1; - /* If the j+1st entry exists and is larger than the jth, - use it instead. */ - if (j + 1 < nivals - && ivals[heap[j + 1]].err >= ivals[heap[j]].err) - j++; - /* Do we need to move the ith entry up? */ - if (ivals[heap[j]].err <= ivals[heap[i]].err) - break; - else - { - t = heap[j]; - heap[j] = heap[i]; - heap[i] = t; - i = j; - } - } - - } - - /* Do we need to split this interval? */ - else if (split) - { - - /* Some values we will need often... */ - d = iv->depth; - /* Generate the interval on the left */ - ivl = &(ivals[heap[nivals++]]); - ivl->a = iv->a; - ivl->b = m; - ml = (ivl->a + ivl->b) / 2; - hl = h / 2; - ivl->depth = 0; - ivl->rdepth = iv->rdepth + 1; - ivl->fx[0] = iv->fx[0]; - ivl->fx[32] = iv->fx[16]; - { - ColumnVector ex (n[0] - 1); - if (wrap) - { - for (i = 0; i < n[0] - 1; i++) - ex (i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl)); - } - else - { - for (i = 0; i < n[0] - 1; i++) - ex (i) = ml + xi[(i + 1) * skip[0]] * hl; - } - fargs (0) = ex; - retval = feval (fcn, fargs, 1); - if (retval.length () != 1 || !retval (0).is_real_matrix ()) - { - error - ("cquad: integrand must return a single, real-valued vector."); - return octave_value_list (); - } - Matrix effex = retval (0).matrix_value (); - if (effex.length () != ex.length ()) - { - error - ("cquad: integrand must return a single, real-valued vector of the same size as the input."); - return octave_value_list (); - } - neval += effex.length (); - for (i = 0; i < n[0] - 1; i++) - { - j = (i + 1) * skip[0]; - ivl->fx[j] = effex (i); - if (wrap) - { - xw = ex(i); - ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2; - } - } - } - nnans = 0; - for (i = 0; i <= 32; i += skip[0]) - { - if (!xfinite (ivl->fx[i])) - { - nans[nnans++] = i; - ivl->fx[i] = 0.0; - } - } - Vinvfx (ivl->fx, ivl->c, 0); - if (nnans > 0) - { - downdate (ivl->c, n[0], 0, nans, nnans); - for (i = 0; i < nnans; i++) - ivl->fx[i] = octave_NaN; - } - for (i = 0; i <= n[d]; i++) - { - ivl->c[idx[d] + i] = 0.0; - for (j = i; j <= n[d]; j++) - ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j]; - } - ncdiff = 0.0; - for (i = 0; i <= n[0]; i++) - { - temp = ivl->c[i] - ivl->c[idx[d] + i]; - ncdiff += temp * temp; - } - for (i = n[0] + 1; i <= n[d]; i++) - { - temp = ivl->c[idx[d] + i]; - ncdiff += temp * temp; - } - ncdiff = sqrt (ncdiff); - ivl->err = ncdiff * h; - /* Check for divergence. */ - ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0 - && ivl->c[0] / iv->c[0] > 2); - if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth) - { - igral = copysign (octave_Inf, igral); - warning ("cquad: divergent integral detected."); - break; - } - - /* Compute the local integral. */ - ivl->igral = h * w * ivl->c[0]; - - - /* Generate the interval on the right */ - ivr = &(ivals[heap[nivals++]]); - ivr->a = m; - ivr->b = iv->b; - mr = (ivr->a + ivr->b) / 2; - hr = h / 2; - ivr->depth = 0; - ivr->rdepth = iv->rdepth + 1; - ivr->fx[0] = iv->fx[16]; - ivr->fx[32] = iv->fx[32]; - { - ColumnVector ex (n[0] - 1); - if (wrap) - { - for (i = 0; i < n[0] - 1; i++) - ex (i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr)); - } - else - { - for (i = 0; i < n[0] - 1; i++) - ex (i) = mr + xi[(i + 1) * skip[0]] * hr; - } - fargs (0) = ex; - retval = feval (fcn, fargs, 1); - if (retval.length () != 1 || !retval (0).is_real_matrix ()) - { - error - ("cquad: integrand must return a single, real-valued vector."); - return octave_value_list (); - } - Matrix effex = retval (0).matrix_value (); - if (effex.length () != ex.length ()) - { - error - ("cquad: integrand must return a single, real-valued vector of the same size as the input."); - return octave_value_list (); - } - neval += effex.length (); - for (i = 0; i < n[0] - 1; i++) - { - j = (i + 1) * skip[0]; - ivr->fx[j] = effex (i); - if (wrap) - { - xw = ex(i); - ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2; - } - } - } - nnans = 0; - for (i = 0; i <= 32; i += skip[0]) - { - if (!xfinite (ivr->fx[i])) - { - nans[nnans++] = i; - ivr->fx[i] = 0.0; - } - } - Vinvfx (ivr->fx, ivr->c, 0); - if (nnans > 0) - { - downdate (ivr->c, n[0], 0, nans, nnans); - for (i = 0; i < nnans; i++) - ivr->fx[i] = octave_NaN; - } - for (i = 0; i <= n[d]; i++) - { - ivr->c[idx[d] + i] = 0.0; - for (j = i; j <= n[d]; j++) - ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j]; - } - ncdiff = 0.0; - for (i = 0; i <= n[0]; i++) - { - temp = ivr->c[i] - ivr->c[idx[d] + i]; - ncdiff += temp * temp; - } - for (i = n[0] + 1; i <= n[d]; i++) - { - temp = ivr->c[idx[d] + i]; - ncdiff += temp * temp; - } - ncdiff = sqrt (ncdiff); - ivr->err = ncdiff * h; - /* Check for divergence. */ - ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0 - && ivr->c[0] / iv->c[0] > 2); - if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth) - { - igral = copysign (octave_Inf, igral); - warning ("cquad: divergent integral detected."); - break; - } - - /* Compute the local integral. */ - ivr->igral = h * w * ivr->c[0]; - - - /* Fix-up the heap: we now have one interval on top - that we don't need any more and two new, unsorted - ones at the bottom. */ - /* Flip the last interval to the top of the heap and - sift down. */ - t = heap[nivals - 1]; - heap[nivals - 1] = heap[0]; - heap[0] = t; - nivals--; - /* Sift this interval back down the heap. */ - i = 0; - while (2 * i + 1 < nivals - 1) - { - j = 2 * i + 1; - if (j + 1 < nivals - 1 - && ivals[heap[j + 1]].err >= ivals[heap[j]].err) - j++; - if (ivals[heap[j]].err <= ivals[heap[i]].err) - break; - else - { - t = heap[j]; - heap[j] = heap[i]; - heap[i] = t; - i = j; - } - } - - /* Now grab the last interval and sift it up the heap. */ - i = nivals - 1; - while (i > 0) - { - j = (i - 1) / 2; - if (ivals[heap[j]].err < ivals[heap[i]].err) - { - t = heap[j]; - heap[j] = heap[i]; - heap[i] = t; - i = j; - } - else - break; - } - - - } - - /* Otherwise, just fix-up the heap. */ - else - { - i = 0; - while (2 * i + 1 < nivals) - { - j = 2 * i + 1; - if (j + 1 < nivals - && ivals[heap[j + 1]].err >= ivals[heap[j]].err) - j++; - if (ivals[heap[j]].err <= ivals[heap[i]].err) - break; - else - { - t = heap[j]; - heap[j] = heap[i]; - heap[i] = t; - i = j; - } - } - - } - - /* If the heap is about to overflow, remove the last two - intervals. */ - while (nivals > cquad_heapsize - 2) - { - iv = &(ivals[heap[nivals - 1]]); -/* printf - ("cquad: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", - heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, - iv->depth); -*/ - err_final += iv->err; - igral_final += iv->igral; - nivals--; - } - - /* Collect the value of the integral and error. */ - igral = igral_final; - err = err_final; - for (i = 0; i < nivals; i++) - { - igral += ivals[heap[i]].igral; - err += ivals[heap[i]].err; - } - - } - - /* Dump the contents of the heap. */ -/* for (i = 0; i < nivals; i++) - { - iv = &(ivals[heap[i]]); - printf - ("cquad: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n", - i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth, - iv->rdepth, iv->ndiv); - } -*/ - /* Clean up and present the results. */ - retval (0) = igral; - if (nargout > 1) - retval (1) = err; - if (nargout > 2) - retval (2) = neval; - /* All is well that ends well. */ - return retval; -}
--- a/src/DLD-FUNCTIONS/module-files Mon Jan 17 17:25:42 2011 -0500 +++ b/src/DLD-FUNCTIONS/module-files Mon Jan 17 20:18:07 2011 -0800 @@ -20,7 +20,6 @@ colloc.cc conv2.cc convhulln.cc -cquad.cc daspk.cc dasrt.cc dassl.cc @@ -58,6 +57,7 @@ pinv.cc qr.cc quad.cc +quadcc.cc qz.cc rand.cc rcond.cc
--- a/src/DLD-FUNCTIONS/quad.cc Mon Jan 17 17:25:42 2011 -0500 +++ b/src/DLD-FUNCTIONS/quad.cc Mon Jan 17 20:18:07 2011 -0800 @@ -214,7 +214,7 @@ \n\ Note: because @code{quad} is written in Fortran it\n\ cannot be called recursively.\n\ -@seealso{quad_options,quadv,quadl,quadgk,trapz}\n\ +@seealso{quad_options,quadv,quadl,quadgk,quadcc,trapz,dblquad,triplequad}\n\ @end deftypefn") { octave_value_list retval;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/quadcc.cc Mon Jan 17 20:18:07 2011 -0800 @@ -0,0 +1,2255 @@ +/* + +Copyright (C) 2010-2011 Pedro Gonnet + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include "lo-math.h" +#include "lo-ieee.h" +#include "oct.h" +#include "parse.h" +#include "ov-fcn-handle.h" + +#if ! defined (HAVE_COPYSIGN) && defined (HAVE__COPYSIGN) +#define copysign _copysign +#define HAVE_COPYSIGN 1 +#endif + +/* Define the size of the interval heap. */ +#define cquad_heapsize 200 + + +/* Data of a single interval */ +typedef struct +{ + double a, b; + double c[64]; + double fx[33]; + double igral, err; + int depth, rdepth, ndiv; +} cquad_ival; + +/* Some constants and matrices that we'll need. + */ + +static const double xi[33] = { + -1., -0.99518472667219688624, -0.98078528040323044912, + -0.95694033573220886493, -0.92387953251128675612, + -0.88192126434835502970, -0.83146961230254523708, + -0.77301045336273696082, -0.70710678118654752440, + -0.63439328416364549822, -0.55557023301960222475, + -0.47139673682599764857, -0.38268343236508977173, + -0.29028467725446236764, -0.19509032201612826785, + -0.098017140329560601995, 0., 0.098017140329560601995, + 0.19509032201612826785, 0.29028467725446236764, 0.38268343236508977173, + 0.47139673682599764857, 0.55557023301960222475, 0.63439328416364549822, + 0.70710678118654752440, 0.77301045336273696082, 0.83146961230254523708, + 0.88192126434835502970, 0.92387953251128675612, 0.95694033573220886493, + 0.98078528040323044912, 0.99518472667219688624, 1. +}; + +static const double bee[68] = { + 0.00000000000000e+00, 2.28868854108532e-01, 0.00000000000000e+00, + -8.15740215243451e-01, 0.00000000000000e+00, 5.31212715259731e-01, + 0.00000000000000e+00, 1.38538036812454e-02, 0.00000000000000e+00, + 3.74405228908818e-02, 0.00000000000000e+00, 2.12224115039342e-01, + 0.00000000000000e+00, -8.16362644507898e-01, 0.00000000000000e+00, + 5.35648426691481e-01, 0.00000000000000e+00, 1.52417902753662e-03, + 0.00000000000000e+00, 2.63058840550873e-03, 0.00000000000000e+00, + 4.15292106318904e-03, 0.00000000000000e+00, 6.97106011119775e-03, + 0.00000000000000e+00, 1.35535708431058e-02, 0.00000000000000e+00, + 3.52132898424856e-02, 0.00000000000000e+00, 2.06946714741884e-01, + 0.00000000000000e+00, -8.15674251283876e-01, 0.00000000000000e+00, + 5.38841175520580e-01, 0.00000000000000e+00, 1.84909689577590e-04, + 0.00000000000000e+00, 2.90936325007499e-04, 0.00000000000000e+00, + 3.84877750950089e-04, 0.00000000000000e+00, 4.86436656735046e-04, + 0.00000000000000e+00, 6.08688640346879e-04, 0.00000000000000e+00, + 7.66732830740331e-04, 0.00000000000000e+00, 9.82753336104205e-04, + 0.00000000000000e+00, 1.29359957505615e-03, 0.00000000000000e+00, + 1.76616363801885e-03, 0.00000000000000e+00, 2.53323433039089e-03, + 0.00000000000000e+00, 3.88872172121956e-03, 0.00000000000000e+00, + 6.58635106468291e-03, 0.00000000000000e+00, 1.30326736343254e-02, + 0.00000000000000e+00, 3.44353850696714e-02, 0.00000000000000e+00, + 2.05025409531915e-01, 0.00000000000000e+00, -8.14985893995401e-01, + 0.00000000000000e+00, 5.40679930965238e-01 +}; + +static const double Lalpha[33] = { + 5.77350269189626e-01, 5.16397779494322e-01, 5.07092552837110e-01, + 5.03952630678970e-01, 5.02518907629606e-01, 5.01745206004255e-01, + 5.01280411827603e-01, 5.00979432868120e-01, 5.00773395667191e-01, + 5.00626174321759e-01, 5.00517330712619e-01, 5.00434593736979e-01, + 5.00370233297676e-01, 5.00319182924304e-01, 5.00278009473803e-01, + 5.00244319584578e-01, 5.00216403386025e-01, 5.00193012939056e-01, + 5.00173220168024e-01, 5.00156323280355e-01, 5.00141783641018e-01, + 5.00129182278347e-01, 5.00118189340972e-01, 5.00108542278496e-01, + 5.00100030010004e-01, 5.00092481273333e-01, 5.00085755939229e-01, + 5.00079738458365e-01, 5.00074332862969e-01, 5.00069458915387e-01, + 5.00065049112355e-01, 5.00061046334395e-01, 5.00057401986298e-01 +}; + +static const double Lgamma[33] = { + 0.0, 0.0, 5.16397779494322e-01, 5.07092552837110e-01, 5.03952630678970e-01, + 5.02518907629606e-01, 5.01745206004255e-01, 5.01280411827603e-01, + 5.00979432868120e-01, 5.00773395667191e-01, 5.00626174321759e-01, + 5.00517330712619e-01, 5.00434593736979e-01, 5.00370233297676e-01, + 5.00319182924304e-01, 5.00278009473803e-01, 5.00244319584578e-01, + 5.00216403386025e-01, 5.00193012939056e-01, 5.00173220168024e-01, + 5.00156323280355e-01, 5.00141783641018e-01, 5.00129182278347e-01, + 5.00118189340972e-01, 5.00108542278496e-01, 5.00100030010003e-01, + 5.00092481273333e-01, 5.00085755939229e-01, 5.00079738458365e-01, + 5.00074332862969e-01, 5.00069458915387e-01, 5.00065049112355e-01, + 5.00061046334395e-01 +}; + +static const double V1inv[5 * 5] = { + .47140452079103168293e-1, .37712361663282534635, .56568542494923801952, + .37712361663282534635, .47140452079103168293e-1, + -.81649658092772603273e-1, -.46188021535170061160, 0, + .46188021535170061160, .81649658092772603273e-1, .15058465048420853962, + .12046772038736683169, -.54210474174315074262, .12046772038736683169, + .15058465048420853962, -.21380899352993950775, .30237157840738178177, -0., + -.30237157840738178177, .21380899352993950775, .10774960475223581324, + -.21549920950447162648, .21549920950447162648, -.21549920950447162648, + .10774960475223581324 +}; + +static const double V2inv[9 * 9] = { + .11223917161691230546e-1, .10339219839658349826, .19754094204576565761, + .25577315077753587922, .27835314560994251755, .25577315077753587922, + .19754094204576565761, .10339219839658349826, .11223917161691230546e-1, + -.19440394783993476970e-1, -.16544884625069155470, -.24193725566041460608, + -.16953338808305493604, 0.0, .16953338808305493604, .24193725566041460608, + .16544884625069155470, .19440394783993476970e-1, .26466393115406349388e-1, + .17766815796285469394, .11316664642449611462, -.16306601003711325980, + -.30847037493128779631, -.16306601003711325980, .11316664642449611462, + .17766815796285469394, .26466393115406349388e-1, + -.32395302049990834508e-1, -.15521142532414866547, + .88573492664788602740e-1, .29570405784974857322, 0.0, + -.29570405784974857322, -.88573492664788602740e-1, .15521142532414866547, + .32395302049990834508e-1, .41442155673936851246e-1, + .98186757907405608245e-1, -.23056908429499411784, + -.68047008326360625520e-1, .31797435808002456774, + -.68047008326360625520e-1, -.23056908429499411784, + .98186757907405608245e-1, .41442155673936851246e-1, + -.49981120317798783134e-1, -.24861810572835756217e-1, + .23561326072010832539, -.24472785656448415351, 0.0, .24472785656448415351, + -.23561326072010832539, .24861810572835756217e-1, + .49981120317798783134e-1, .79691635865674781228e-1, + -.95725617891693941833e-1, -.57957553356854386344e-1, + .21164072460540271452, -.27529837844505833514, .21164072460540271452, + -.57957553356854386344e-1, -.95725617891693941833e-1, + .79691635865674781228e-1, + -.10894869830716590913, .20131094491947531782, -.15407672674888869038, + .83385723639789791384e-1, 0.0, -.83385723639789791384e-1, + .15407672674888869038, -.20131094491947531782, .10894869830716590913, + .54581057089643838221e-1, -.10916211417928767644, .10916211417928767644, + -.10916211417928767644, .10916211417928767644, -.10916211417928767644, + .10916211417928767644, -.10916211417928767644, .54581057089643838221e-1 +}; + +static const double V3inv[17 * 17] = { + .27729677693590098996e-2, .26423663180333065153e-1, + .53374068493933898312e-1, .77007854739523195947e-1, + .98257061072911596869e-1, .11538049741786835604, .12832134344120884559, + .13612785914022865001, .13888293186236181317, .13612785914022865001, + .12832134344120884559, .11538049741786835604, .98257061072911596869e-1, + .77007854739523195947e-1, .53374068493933898312e-1, + .26423663180333065153e-1, .27729677693590098996e-2, + -.48029210642807413690e-2, -.44887724635478800254e-1, + -.85409520147301089416e-1, -.11090267822061423050, -.12033983162705862441, + -.11102786862182788886, -.85054870109799336515e-1, + -.45998467987742225160e-1, 0.0, .45998467987742225160e-1, + .85054870109799336515e-1, .11102786862182788886, .12033983162705862441, + .11090267822061423050, .85409520147301089416e-1, .44887724635478800254e-1, + .48029210642807413690e-2, .62758546879582030087e-2, + .55561297093529155869e-1, + .93281491021051539742e-1, .92320151237493695139e-1, + .55077987469605684531e-1, + -.96998141716497488255e-2, -.80285961895427405567e-1, + -.13496839655913850224, + -.15512521776684524331, -.13496839655913850224, -.80285961895427405567e-1, + -.96998141716497488255e-2, .55077987469605684531e-1, + .92320151237493695139e-1, .93281491021051539742e-1, + .55561297093529155869e-1, .62758546879582030087e-2, + -.74850969394858555939e-2, -.61751608943839234096e-1, + -.82974150437304275958e-1, -.38437763431942633378e-1, + .45745502025779701366e-1, .12369235652734542162, .14720439712852868239, + .98768034347019704401e-1, 0.0, + -.98768034347019704401e-1, -.14720439712852868239, -.12369235652734542162, + -.45745502025779701366e-1, .38437763431942633378e-1, + .82974150437304275958e-1, .61751608943839234096e-1, + .74850969394858555939e-2, .86710099994384056338e-2, + .64006230103659573344e-1, .58517426396091675690e-1, + -.29743410528985802680e-1, + -.11934127779157114754, -.12686773515361299409, -.30729137153877447035e-1, + .97307836256600731568e-1, .15635811574451401023, .97307836256600731568e-1, + -.30729137153877447035e-1, -.12686773515361299409, -.11934127779157114754, + -.29743410528985802680e-1, .58517426396091675690e-1, + .64006230103659573344e-1, .86710099994384056338e-2, + -.97486395666294840165e-2, -.62995604908060224672e-1, + -.24373234450275529219e-1, .87760984413626872730e-1, + .12205204576993351394, + .16216004196864002088e-1, -.12422320942156845775, -.13682714580929614678, + 0.0, .13682714580929614678, .12422320942156845775, + -.16216004196864002088e-1, -.12205204576993351394, + -.87760984413626872730e-1, .24373234450275529219e-1, + .62995604908060224672e-1, .97486395666294840165e-2, + .10956271233750488468e-1, .58613204255294358939e-1, + -.13306063940736618859e-1, -.11606666444978454399, + -.52059598001115805639e-1, .10868540217796151849, .12594452879014618005, + -.44678658254872910434e-1, -.15617684362128533405, + -.44678658254872910434e-1, .12594452879014618005, .10868540217796151849, + -.52059598001115805639e-1, -.11606666444978454399, + -.13306063940736618859e-1, .58613204255294358939e-1, + .10956271233750488468e-1, -.12098893000863087230e-1, + -.51626244709126208453e-1, .48919433304746979330e-1, + .10467644465949427090, + -.48729879523084673782e-1, -.13668732103524749234, .28190838706814496438e-1, + .15434223333238741600, 0.0, -.15434223333238741600, + -.28190838706814496438e-1, .13668732103524749234, + .48729879523084673782e-1, -.10467644465949427090, + -.48919433304746979330e-1, .51626244709126208453e-1, + .12098893000863087230e-1, .13542668300437944822e-1, + .41712033418258689308e-1, + -.76190463272803434388e-1, -.58303943170068132010e-1, .12158068748245606853, + .42121099930651007882e-1, -.14684425840766337756, + -.16108203535058647043e-1, .15698075850757976092, + -.16108203535058647043e-1, -.14684425840766337756, + .42121099930651007882e-1, .12158068748245606853, + -.58303943170068132010e-1, -.76190463272803434388e-1, + .41712033418258689308e-1, .13542668300437944822e-1, + -.14939634995117694417e-1, -.30047246373341564039e-1, + .91624635082546425678e-1, -.79133374319110026377e-2, + -.12292558212072233355, .90013382617762643524e-1, + .84013717196539593395e-1, -.14813033309980695856, 0.0, + .14813033309980695856, -.84013717196539593395e-1, + -.90013382617762643524e-1, + .12292558212072233355, .79133374319110026377e-2, -.91624635082546425678e-1, + .30047246373341564039e-1, .14939634995117694417e-1, + .16986031342807474208e-1, + .15760203882617033601e-1, -.91494054040950941996e-1, + .70082459207876130806e-1, + .53390713710144539104e-1, -.14340746778352039430, .84048122493418898508e-1, + .72456667788091316868e-1, -.15564535320096811360, + .72456667788091316868e-1, .84048122493418898508e-1, + -.14340746778352039430, .53390713710144539104e-1, + .70082459207876130806e-1, -.91494054040950941996e-1, + .15760203882617033601e-1, + .16986031342807474208e-1, -.18994065631858742028e-1, + -.82901821370405592927e-3, .77239669773015192888e-1, + -.10850735431039424680, .47524484622086496464e-1, + .69148184871588737021e-1, -.14829314646228194928, .11992057742398672066, + 0.0, -.11992057742398672066, .14829314646228194928, + -.69148184871588737021e-1, -.47524484622086496464e-1, + .10850735431039424680, -.77239669773015192888e-1, + .82901821370405592927e-3, .18994065631858742028e-1, + .22761703826371535132e-1, -.17728848711449643358e-1, + -.47496371572480503788e-1, .10659958402328690063, -.11696013966166296514, + .63073750910894244526e-1, .32928881123602721303e-1, + -.12280950532497593683, .15926189077282729505, -.12280950532497593683, + .32928881123602721303e-1, .63073750910894244526e-1, + -.11696013966166296514, .10659958402328690063, -.47496371572480503788e-1, + -.17728848711449643358e-1, .22761703826371535132e-1, + -.26493215276042203434e-1, .35579780856128386192e-1, + .10447309718398935122e-1, -.68616154085314996709e-1, + .11775363082763954214, -.13918901977011837274, .12312819418827395690, + -.72053565748259077905e-1, 0.0, .72053565748259077905e-1, + -.12312819418827395690, .13918901977011837274, -.11775363082763954214, + .68616154085314996709e-1, -.10447309718398935122e-1, + -.35579780856128386192e-1, + .26493215276042203434e-1, .40742523354399706918e-1, + -.73124912999529117195e-1, .49317266444153837821e-1, + -.13686605413876015320e-1, -.28342624942191100464e-1, + .70371855298258216249e-1, -.10600251632853603875, .12981016288391131812, + -.13817029659318161476, .12981016288391131812, -.10600251632853603875, + .70371855298258216249e-1, -.28342624942191100464e-1, + -.13686605413876015320e-1, + .49317266444153837821e-1, -.73124912999529117195e-1, + .40742523354399706918e-1, -.54944368958699908688e-1, + .10777725663147408190, -.10152395581538265428, .91369146312596428468e-1, + -.77703071757424700773e-1, .61050911730999815031e-1, + -.42052599404498348871e-1, .21438229266251454773e-1, 0.0, + -.21438229266251454773e-1, .42052599404498348871e-1, + -.61050911730999815031e-1, .77703071757424700773e-1, + -.91369146312596428468e-1, + .10152395581538265428, -.10777725663147408190, .54944368958699908688e-1, + .27485608464748840573e-1, -.54971216929497681146e-1, + .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .54971216929497681146e-1, + -.54971216929497681146e-1, .27485608464748840573e-1 +}; + +static const double V4inv[33 * 33] = { + .69120897476690862600e-3, .66419939766331555194e-2, + .13600665164323186111e-1, .20122785860913684493e-1, + .26583214101668429944e-1, .32712713318999268739e-1, + .38576221976287138036e-1, .44033030938268925133e-1, + .49092709529622799673e-1, .53657949874312515646e-1, + .57724533144734311859e-1, .61219564530655179096e-1, + .64138907503837875026e-1, .66427905189318792009e-1, + .68088956652280022887e-1, .69083051391555695878e-1, + .69422738116739271449e-1, .69083051391555695878e-1, + .68088956652280022887e-1, .66427905189318792009e-1, + .64138907503837875026e-1, .61219564530655179096e-1, + .57724533144734311859e-1, .53657949874312515646e-1, + .49092709529622799673e-1, .44033030938268925133e-1, + .38576221976287138036e-1, .32712713318999268739e-1, + .26583214101668429944e-1, .20122785860913684493e-1, + .13600665164323186111e-1, .66419939766331555194e-2, + .69120897476690862600e-3, -.11972090629438798134e-2, + -.11448874821643225573e-1, -.23104401104002905904e-1, + -.33352899418646530133e-1, -.42538626424075425908e-1, + -.49969730733911825941e-1, -.55555454015360728353e-1, + -.58955533624852604918e-1, -.60126044219122513907e-1, + -.58959430451175833624e-1, -.55546925396227130606e-1, + -.49984739749347973762e-1, -.42513009141170294365e-1, + -.33399140950669746346e-1, -.23007690803851790829e-1, + -.11728275717520066169e-1, 0.0, .11728275717520066169e-1, + .23007690803851790829e-1, .33399140950669746346e-1, + .42513009141170294365e-1, .49984739749347973762e-1, + .55546925396227130606e-1, .58959430451175833624e-1, + .60126044219122513907e-1, .58955533624852604918e-1, + .55555454015360728353e-1, .49969730733911825941e-1, + .42538626424075425908e-1, .33352899418646530133e-1, + .23104401104002905904e-1, .11448874821643225573e-1, + .11972090629438798134e-2, .15501585012936019146e-2, + .14628781502199620482e-1, .28684915921474815271e-1, + .39299396074628048026e-1, .46393418975496284204e-1, + .48756902531094699526e-1, .46331333488337494692e-1, + .39012645376980228775e-1, .27452795421085791153e-1, + .12430953621169863781e-1, -.47682978056024928800e-2, + -.22825828045428973853e-1, + -.40195512090720278312e-1, -.55503004262826221955e-1, + -.67424537752827046308e-1, -.75020199300113606452e-1, + -.77607844312483656131e-1, -.75020199300113606452e-1, + -.67424537752827046308e-1, -.55503004262826221955e-1, + -.40195512090720278312e-1, -.22825828045428973853e-1, + -.47682978056024928800e-2, .12430953621169863781e-1, + .27452795421085791153e-1, .39012645376980228775e-1, + .46331333488337494692e-1, .48756902531094699526e-1, + .46393418975496284204e-1, .39299396074628048026e-1, + .28684915921474815271e-1, .14628781502199620482e-1, + .15501585012936019146e-2, -.18377757558949194214e-2, + -.17050470050949761565e-1, -.31952119564923250836e-1, + -.40197423449026348155e-1, + -.41205649520281371624e-1, -.33909965817492272248e-1, + -.19393664422115332144e-1, .56661049630886784692e-3, + .22948272173686561721e-1, .44489719570904738207e-1, + .61790363672287920596e-1, .72121014727028013894e-1, + .73627151185287858579e-1, .65784665375961398923e-1, + .49369676372333667559e-1, .26444326317059715065e-1, 0.0, + -.26444326317059715065e-1, -.49369676372333667559e-1, + -.65784665375961398923e-1, -.73627151185287858579e-1, + -.72121014727028013894e-1, -.61790363672287920596e-1, + -.44489719570904738207e-1, -.22948272173686561721e-1, + -.56661049630886784692e-3, .19393664422115332144e-1, + .33909965817492272248e-1, .41205649520281371624e-1, + .40197423449026348155e-1, .31952119564923250836e-1, + .17050470050949761565e-1, .18377757558949194214e-2, + .20942714740729767769e-2, .18935902405146518232e-1, + .33335840852491735126e-1, .36770680999102286065e-1, + .28873194534132768509e-1, .10267303017729535513e-1, + -.14607738306201572890e-1, -.40139568545572305818e-1, + -.59808326733858291561e-1, -.68528358823372627506e-1, + -.63306535387619244879e-1, -.44508601817574921056e-1, + -.15449116105605395357e-1, .17941083795006546367e-1, + .48747356011657242123e-1, .70329553984201665523e-1, + .78106117292526169663e-1, .70329553984201665523e-1, + .48747356011657242123e-1, .17941083795006546367e-1, + -.15449116105605395357e-1, -.44508601817574921056e-1, + -.63306535387619244879e-1, -.68528358823372627506e-1, + -.59808326733858291561e-1, + -.40139568545572305818e-1, -.14607738306201572890e-1, + .10267303017729535513e-1, .28873194534132768509e-1, + .36770680999102286065e-1, .33335840852491735126e-1, + .18935902405146518232e-1, .20942714740729767769e-2, + -.23245285491878278419e-2, -.20401404737639389919e-1, + -.33019548231022514097e-1, -.29709828426463720091e-1, + -.11760070922697422156e-1, .15987584743850393793e-1, + .43619012891472813485e-1, .61177322409671487721e-1, + .61144030218486655594e-1, + .41895377620089086167e-1, .80232011820644308033e-2, + -.30574701186675900915e-1, + -.62072243008844865848e-1, -.76336186183574765586e-1, + -.68435466095345537115e-1, -.40237669208466966207e-1, 0.0, + .40237669208466966207e-1, .68435466095345537115e-1, + .76336186183574765586e-1, .62072243008844865848e-1, + .30574701186675900915e-1, -.80232011820644308033e-2, + -.41895377620089086167e-1, -.61144030218486655594e-1, + -.61177322409671487721e-1, -.43619012891472813485e-1, + -.15987584743850393793e-1, .11760070922697422156e-1, + .29709828426463720091e-1, .33019548231022514097e-1, + .20401404737639389919e-1, .23245285491878278419e-2, + .25451717261579269307e-2, .21480418595666878775e-1, + .31177212469293007998e-1, .19816333607013379373e-1, + -.72439496274458793681e-2, -.38404203906598342397e-1, + -.57633632255322221046e-1, -.54070547403585392952e-1, + -.26249823354368866005e-1, .15643058212336881516e-1, + .54539832735118677194e-1, .73283028002473989724e-1, + .62835303524135936213e-1, .26175977027801048141e-1, + -.22193636309998606610e-1, -.62597049956093311234e-1, + -.78206986173170212505e-1, -.62597049956093311234e-1, + -.22193636309998606610e-1, .26175977027801048141e-1, + .62835303524135936213e-1, + .73283028002473989724e-1, .54539832735118677194e-1, + .15643058212336881516e-1, + -.26249823354368866005e-1, -.54070547403585392952e-1, + -.57633632255322221046e-1, -.38404203906598342397e-1, + -.72439496274458793681e-2, .19816333607013379373e-1, + .31177212469293007998e-1, .21480418595666878775e-1, + .25451717261579269307e-2, -.27506573922483820005e-2, + -.22224442095099251870e-1, -.27949927254215773020e-1, + -.80918481053370034987e-2, .25121859354449306916e-1, + .51563535009373061074e-1, .51936965107145960512e-1, + .22146626648171527753e-1, + -.24172689882103382748e-1, -.61731229104853568296e-1, + -.68477262429344201201e-1, -.38311232728303704742e-1, + .14160578713659552679e-1, .61248813427564184033e-1, + .77136328841293031805e-1, .52514801765183697988e-1, 0.0, + -.52514801765183697988e-1, -.77136328841293031805e-1, + -.61248813427564184033e-1, -.14160578713659552679e-1, + .38311232728303704742e-1, + .68477262429344201201e-1, .61731229104853568296e-1, + .24172689882103382748e-1, + -.22146626648171527753e-1, -.51936965107145960512e-1, + -.51563535009373061074e-1, -.25121859354449306916e-1, + .80918481053370034987e-2, .27949927254215773020e-1, + .22224442095099251870e-1, .27506573922483820005e-2, + .29562461131654311467e-2, .22630271480554450613e-1, + .23547399831373800971e-1, -.43964593440902476642e-2, + -.39055315767504970597e-1, -.52369643937940066804e-1, + -.28506131614971613422e-1, .19906048093338832322e-1, + .60408880866392420279e-1, .62493397473656883090e-1, + .21391278377641297859e-1, -.37302864786623254746e-1, + -.73665127933539496872e-1, -.61706142476854010202e-1, + -.78065168882546327888e-2, .52335307373945544428e-1, + .78278746279419264777e-1, .52335307373945544428e-1, + -.78065168882546327888e-2, -.61706142476854010202e-1, + -.73665127933539496872e-1, -.37302864786623254746e-1, + .21391278377641297859e-1, .62493397473656883090e-1, + .60408880866392420279e-1, .19906048093338832322e-1, + -.28506131614971613422e-1, -.52369643937940066804e-1, + -.39055315767504970597e-1, -.43964593440902476642e-2, + .23547399831373800971e-1, .22630271480554450613e-1, + .29562461131654311467e-2, -.31515718415504761303e-2, + -.22739451096655080673e-1, -.18157123602272119779e-1, + .16496480897167303621e-1, .46921166788569301124e-1, + .40644395739978416354e-1, -.46275803430732216900e-2, + -.52883375891308909486e-1, -.61116483226324111734e-1, + -.17411698764545629853e-1, .44773430013166822765e-1, + .73441577962383869198e-1, .42127368371995472815e-1, + -.25504645957196772465e-1, -.74126818045972742488e-1, + -.62780077864719287317e-1, 0.0, .62780077864719287317e-1, + .74126818045972742488e-1, .25504645957196772465e-1, + -.42127368371995472815e-1, -.73441577962383869198e-1, + -.44773430013166822765e-1, .17411698764545629853e-1, + .61116483226324111734e-1, .52883375891308909486e-1, + .46275803430732216900e-2, -.40644395739978416354e-1, + -.46921166788569301124e-1, -.16496480897167303621e-1, + .18157123602272119779e-1, .22739451096655080673e-1, + .31515718415504761303e-2, .33536559294882188208e-2, + .22535348942792006185e-1, + .12048629300953560767e-1, -.27166076791299493403e-1, + -.47492745604230978367e-1, -.19246623430993153174e-1, + .36231297307556299322e-1, .61713617181636122004e-1, + .25928029734266134490e-1, -.40478700752883602818e-1, + -.71053889866326412049e-1, -.31870824482961751482e-1, + .41515251100219081281e-1, .76481960760098381651e-1, + .36726509155999912440e-1, -.40090067032627055969e-1, + -.78270742903374539397e-1, -.40090067032627055969e-1, + .36726509155999912440e-1, .76481960760098381651e-1, + .41515251100219081281e-1, -.31870824482961751482e-1, + -.71053889866326412049e-1, -.40478700752883602818e-1, + .25928029734266134490e-1, .61713617181636122004e-1, + .36231297307556299322e-1, -.19246623430993153174e-1, + -.47492745604230978367e-1, -.27166076791299493403e-1, + .12048629300953560767e-1, .22535348942792006185e-1, + .33536559294882188208e-2, + -.35481220456925318865e-2, -.22062913693073191150e-1, + -.54487362861834144999e-2, .35438821865804087489e-1, + .40733077820527411302e-1, -.67403098138950720914e-2, + -.55559584405239171054e-1, -.42417050790865158745e-1, + .24499901971884704925e-1, .68721232891705409302e-1, + .34086082787461126592e-1, -.43441000373118474002e-1, + -.73878085292669148950e-1, -.18846995664706657127e-1, + .59827776178286834498e-1, .70644634584085901794e-1, 0.0, + -.70644634584085901794e-1, -.59827776178286834498e-1, + .18846995664706657127e-1, .73878085292669148950e-1, + .43441000373118474002e-1, -.34086082787461126592e-1, + -.68721232891705409302e-1, -.24499901971884704925e-1, + .42417050790865158745e-1, .55559584405239171054e-1, + .67403098138950720914e-2, -.40733077820527411302e-1, + -.35438821865804087489e-1, .54487362861834144999e-2, + .22062913693073191150e-1, .35481220456925318865e-2, + .37554176816665075631e-2, .21297045781589919482e-1, + -.13327293083183431816e-2, + -.40635299172764596484e-1, -.27659860508374175359e-1, + .31089232744083445986e-1, .56113781541334176109e-1, + .37577840643257763400e-2, -.60511227350664590865e-1, + -.46670556446129053853e-1, .33263195878575888247e-1, + .72757324720645228775e-1, .15011712351692283635e-1, + -.65601212994924119078e-1, -.60016855838843789772e-1, + .26220858553188665966e-1, .78322776605833552980e-1, + .26220858553188665966e-1, -.60016855838843789772e-1, + -.65601212994924119078e-1, + .15011712351692283635e-1, .72757324720645228775e-1, + .33263195878575888247e-1, + -.46670556446129053853e-1, -.60511227350664590865e-1, + .37577840643257763400e-2, .56113781541334176109e-1, + .31089232744083445986e-1, -.27659860508374175359e-1, + -.40635299172764596484e-1, -.13327293083183431816e-2, + .21297045781589919482e-1, .37554176816665075631e-2, + -.39566995305720591229e-2, -.20291873414438919995e-1, + .80617453830770930551e-2, .42270189157016547906e-1, + .10332624526759093004e-1, -.48054759547616142024e-1, + -.37678032941171643972e-1, + .36617192625732482394e-1, .61009425973424865714e-1, + -.95589113168026591466e-2, + -.71023202645076922361e-1, -.25097788086808784456e-1, + .62406621963267050244e-1, .56907293171100693511e-1, + -.36435383083882206257e-1, -.75790105119208756348e-1, 0.0, + .75790105119208756348e-1, .36435383083882206257e-1, + -.56907293171100693511e-1, -.62406621963267050244e-1, + .25097788086808784456e-1, .71023202645076922361e-1, + .95589113168026591466e-2, + -.61009425973424865714e-1, -.36617192625732482394e-1, + .37678032941171643972e-1, .48054759547616142024e-1, + -.10332624526759093004e-1, -.42270189157016547906e-1, + -.80617453830770930551e-2, .20291873414438919995e-1, + .39566995305720591229e-2, .41776092289182138591e-2, + .19013221163904414395e-1, -.14420609729849899876e-1, + -.40259160586844441220e-1, .86327811113710831649e-2, + .53564430703021034399e-1, .65469185402150431933e-2, + -.60383116311280629856e-1, + -.25657793784058876939e-1, .58745680576829226900e-1, + .45649937869034420296e-1, + -.49167932056844167772e-1, -.62696614328552187977e-1, + .32540234556426699997e-1, .74280410383464269758e-1, + -.11425672633410999870e-1, -.78280649404686404903e-1, + -.11425672633410999870e-1, .74280410383464269758e-1, + .32540234556426699997e-1, -.62696614328552187977e-1, + -.49167932056844167772e-1, .45649937869034420296e-1, + .58745680576829226900e-1, -.25657793784058876939e-1, + -.60383116311280629856e-1, .65469185402150431933e-2, + .53564430703021034399e-1, + .86327811113710831649e-2, -.40259160586844441220e-1, + -.14420609729849899876e-1, .19013221163904414395e-1, + .41776092289182138591e-2, -.43935502082478059199e-2, + -.17528761237509401631e-1, .20208915249153872535e-1, + .34734743119040669109e-1, -.26275910172353637955e-1, + -.46368003346018878786e-1, + .26800056330709381025e-1, .56681476464606609921e-1, + -.24749011438127255898e-1, + -.64934612189056658992e-1, .20333742247679279535e-1, + .71429299070059318651e-1, + -.14452513210428671266e-1, -.75793341281736586582e-1, + .74717094137184935270e-2, .78034921554757317374e-1, 0.0, + -.78034921554757317374e-1, -.74717094137184935270e-2, + .75793341281736586582e-1, .14452513210428671266e-1, + -.71429299070059318651e-1, -.20333742247679279535e-1, + .64934612189056658992e-1, .24749011438127255898e-1, + -.56681476464606609921e-1, + -.26800056330709381025e-1, .46368003346018878786e-1, + .26275910172353637955e-1, + -.34734743119040669109e-1, -.20208915249153872535e-1, + .17528761237509401631e-1, .43935502082478059199e-2, + .46379089482818671473e-2, .15791188144791287229e-1, + -.25134290048737455284e-1, -.26249795071946841205e-1, + .39960457575789924651e-1, .28111892450146525404e-1, + -.51026476400767918226e-1, + -.27266747278681831364e-1, .60708796647861610865e-1, + .23532306960642115854e-1, + -.68169639871532441111e-1, -.18204924701958312032e-1, + .73822890510656128485e-1, .11373392486424717019e-1, + -.77133324017644609416e-1, -.39295877480342619961e-2, + .78351902829418987960e-1, -.39295877480342619961e-2, + -.77133324017644609416e-1, .11373392486424717019e-1, + .73822890510656128485e-1, -.18204924701958312032e-1, + -.68169639871532441111e-1, .23532306960642115854e-1, + .60708796647861610865e-1, -.27266747278681831364e-1, + -.51026476400767918226e-1, .28111892450146525404e-1, + .39960457575789924651e-1, -.26249795071946841205e-1, + -.25134290048737455284e-1, .15791188144791287229e-1, + .46379089482818671473e-2, -.48780095920069827068e-2, + -.13886961667516983541e-1, .29071311049368895844e-1, + .15480559452075811600e-1, -.47527977686242313065e-1, + -.31929089844361042178e-2, .58015667638415922967e-1, + -.14547915466597622925e-1, -.61067668299848923244e-1, + .35093678009090186851e-1, .55378399159800654657e-1, + -.54277226474891610385e-1, -.42023830782434076509e-1, + .69197384645944912066e-1, .22610783557709586445e-1, + -.77269275900637030185e-1, 0.0, .77269275900637030185e-1, + -.22610783557709586445e-1, + -.69197384645944912066e-1, .42023830782434076509e-1, + .54277226474891610385e-1, + -.55378399159800654657e-1, -.35093678009090186851e-1, + .61067668299848923244e-1, .14547915466597622925e-1, + -.58015667638415922967e-1, .31929089844361042178e-2, + .47527977686242313065e-1, -.15480559452075811600e-1, + -.29071311049368895844e-1, .13886961667516983541e-1, + .48780095920069827068e-2, .51591759101720291381e-2, + .11747497650231330965e-1, -.31777863364694653331e-1, + -.34555825499804605557e-2, .47914131921157015198e-1, + -.22573685920142225247e-1, -.45320344390022666738e-1, + .49660630547172186418e-1, .25707858143963615736e-1, + -.68132707341917233933e-1, .67534860185243140399e-2, + .69268150370037450063e-1, -.41585011920451477177e-1, + -.51622397460510041271e-1, .68408139576363036148e-1, + .18981259024768933323e-1, -.78265472429342305554e-1, + .18981259024768933323e-1, .68408139576363036148e-1, + -.51622397460510041271e-1, + -.41585011920451477177e-1, .69268150370037450063e-1, + .67534860185243140399e-2, + -.68132707341917233933e-1, .25707858143963615736e-1, + .49660630547172186418e-1, + -.45320344390022666738e-1, -.22573685920142225247e-1, + .47914131921157015198e-1, -.34555825499804605557e-2, + -.31777863364694653331e-1, .11747497650231330965e-1, + .51591759101720291381e-2, -.54365757412741340377e-2, + -.94862516619529080191e-2, .33240472093448190877e-1, + -.88698898099681552229e-2, + -.40973252097216337576e-1, .42995673349795657065e-1, + .17320914507876958783e-1, + -.62201292691914856803e-1, .24726274174637346693e-1, + .51320859246515407288e-1, + -.62882063373810501763e-1, -.11003569131725622672e-1, + .73842261324108943465e-1, -.39240120294802923208e-1, + -.49293966443941122807e-1, .73552644778818223475e-1, 0.0, + -.73552644778818223475e-1, .49293966443941122807e-1, + .39240120294802923208e-1, -.73842261324108943465e-1, + .11003569131725622672e-1, .62882063373810501763e-1, + -.51320859246515407288e-1, + -.24726274174637346693e-1, .62201292691914856803e-1, + -.17320914507876958783e-1, -.42995673349795657065e-1, + .40973252097216337576e-1, .88698898099681552229e-2, + -.33240472093448190877e-1, .94862516619529080191e-2, + .54365757412741340377e-2, .57750194549356126240e-2, + .69981166020044116791e-2, -.33274982140403110792e-1, + .20297071020698356116e-1, .27898517839646066582e-1, + -.53368678853282030262e-1, .16656482990394548343e-1, + .46342901447260614255e-1, + -.60536796508149003365e-1, .29109107483842596340e-2, + .63224486124385124504e-1, + -.59028872851312033411e-1, -.14783105962696191734e-1, + .74269399241069253865e-1, -.49053677339382384625e-1, + -.33525466624811186739e-1, .78397349622515386647e-1, + -.33525466624811186739e-1, -.49053677339382384625e-1, + .74269399241069253865e-1, -.14783105962696191734e-1, + -.59028872851312033411e-1, + .63224486124385124504e-1, .29109107483842596340e-2, + -.60536796508149003365e-1, + .46342901447260614255e-1, .16656482990394548343e-1, + -.53368678853282030262e-1, + .27898517839646066582e-1, .20297071020698356116e-1, + -.33274982140403110792e-1, + .69981166020044116791e-2, .57750194549356126240e-2, + -.61100308370519200637e-2, -.44383614355738148616e-2, + .32011283412619094811e-1, -.29965011866372897633e-1, + -.10560682331349193348e-1, .51110336443392506342e-1, + -.45012284729681775492e-1, -.94236825555873320102e-2, + .60860695783141264746e-1, + -.55014628647083368926e-1, -.73474782382499482121e-2, + .66640148475243034781e-1, -.62533116045749887988e-1, + -.38650525912400102585e-2, .68429769005837003777e-1, + -.66984505412544901945e-1, 0.0, .66984505412544901945e-1, + -.68429769005837003777e-1, .38650525912400102585e-2, + .62533116045749887988e-1, -.66640148475243034781e-1, + .73474782382499482121e-2, + .55014628647083368926e-1, -.60860695783141264746e-1, + .94236825555873320102e-2, + .45012284729681775492e-1, -.51110336443392506342e-1, + .10560682331349193348e-1, + .29965011866372897633e-1, -.32011283412619094811e-1, + .44383614355738148616e-2, + .61100308370519200637e-2, .65409373892036191538e-2, + .16350101107071157065e-2, -.29301957285983144319e-1, + .36838667173388832579e-1, -.81922703976491586393e-2, + -.36955670021050133434e-1, .58374851095540469865e-1, + -.31977016246946181856e-1, -.25311073698658094646e-1, + .66674413950106952577e-1, + -.54865713324521039571e-1, -.39797027891537985440e-2, + .62830285264808449064e-1, -.72226313251296100676e-1, + .22560232697133353980e-1, .46455784709904033738e-1, + -.78200930751070349956e-1, .46455784709904033738e-1, + .22560232697133353980e-1, -.72226313251296100676e-1, + .62830285264808449064e-1, -.39797027891537985440e-2, + -.54865713324521039571e-1, .66674413950106952577e-1, + -.25311073698658094646e-1, -.31977016246946181856e-1, + .58374851095540469865e-1, -.36955670021050133434e-1, + -.81922703976491586393e-2, .36838667173388832579e-1, + -.29301957285983144319e-1, .16350101107071157065e-2, + .65409373892036191538e-2, -.69686180931868703196e-2, + .11849538727632789870e-2, .25452286414610537766e-1, + -.40522480651713943230e-1, .25694679053362813183e-1, + .14057118113748390637e-1, -.52037614725803488893e-1, + .58849342223684035589e-1, + -.25075229077361409271e-1, -.29559771094034181083e-1, + .68296746944165720199e-1, -.62890462146423984955e-1, + .14457636466274596445e-1, .45787612031322361496e-1, + -.77231759014655809742e-1, .57881203613910543657e-1, 0.0, + -.57881203613910543657e-1, .77231759014655809742e-1, + -.45787612031322361496e-1, -.14457636466274596445e-1, + .62890462146423984955e-1, + -.68296746944165720199e-1, .29559771094034181083e-1, + .25075229077361409271e-1, + -.58849342223684035589e-1, .52037614725803488893e-1, + -.14057118113748390637e-1, -.25694679053362813183e-1, + .40522480651713943230e-1, -.25452286414610537766e-1, + -.11849538727632789870e-2, .69686180931868703196e-2, + .75611653617520254845e-2, -.43290610418608409141e-2, + -.20277062025115566914e-1, + .40362947027704828926e-1, -.38938808024132120254e-1, + .11831186195916702262e-1, + .28476667401744525357e-1, -.59320969056617684621e-1, + .61101629747436200186e-1, + -.29514834848355389223e-1, -.20668001885001084821e-1, + .62923592802445122793e-1, -.73558456263588833115e-1, + .45314556330160999776e-1, .79031645918426015574e-2, + -.58136953576334689357e-1, .78538474524006405758e-1, + -.58136953576334689357e-1, .79031645918426015574e-2, + .45314556330160999776e-1, -.73558456263588833115e-1, + .62923592802445122793e-1, -.20668001885001084821e-1, + -.29514834848355389223e-1, .61101629747436200186e-1, + -.59320969056617684621e-1, .28476667401744525357e-1, + .11831186195916702262e-1, -.38938808024132120254e-1, + .40362947027704828926e-1, -.20277062025115566914e-1, + -.43290610418608409141e-2, .75611653617520254845e-2, + -.81505692478987769484e-2, .74297333588288568430e-2, + .14314212513540223314e-1, -.36711242251332751607e-1, + .46240027755503814626e-1, -.34921532671769023773e-1, + .46930051972353714773e-2, + .32842770336385381562e-1, -.61317813706529588466e-1, + .67000809902468893103e-1, + -.45337449655535622885e-1, .35794459576271920867e-2, + .41830061526027213385e-1, + -.72091371931944711708e-1, .74150028530317793195e-1, + -.46487632538609942002e-1, 0.0, .46487632538609942002e-1, + -.74150028530317793195e-1, .72091371931944711708e-1, + -.41830061526027213385e-1, -.35794459576271920867e-2, + .45337449655535622885e-1, -.67000809902468893103e-1, + .61317813706529588466e-1, -.32842770336385381562e-1, + -.46930051972353714773e-2, .34921532671769023773e-1, + -.46240027755503814626e-1, .36711242251332751607e-1, + -.14314212513540223314e-1, -.74297333588288568430e-2, + .81505692478987769484e-2, .90693182942442189743e-2, + -.11121000903959576737e-1, -.71308296141317458546e-2, + .29219439765986671645e-1, -.45820286629778129593e-1, + .49088381175879124421e-1, -.35614888785023038938e-1, + .78906970900092777895e-2, + .26262843038404929480e-1, -.56143674270125757857e-1, + .71700220472378350694e-1, + -.66963544500697307945e-1, .42215091779892228883e-1, + -.41338867413966866997e-2, -.36164891772995367321e-1, + .66584367783847858225e-1, -.77874712365070098328e-1, + .66584367783847858225e-1, -.36164891772995367321e-1, + -.41338867413966866997e-2, .42215091779892228883e-1, + -.66963544500697307945e-1, + .71700220472378350694e-1, -.56143674270125757857e-1, + .26262843038404929480e-1, + .78906970900092777895e-2, -.35614888785023038938e-1, + .49088381175879124421e-1, + -.45820286629778129593e-1, .29219439765986671645e-1, + -.71308296141317458546e-2, -.11121000903959576737e-1, + .90693182942442189743e-2, -.99848472706332791043e-2, + .14701271465939718856e-1, -.32917820356048383366e-3, + -.19201195309873585230e-1, .38409681836626963278e-1, + -.51647324405878909521e-1, .54522171113149311354e-1, + -.45040302741689006270e-1, .24183738595685990149e-1, + .42204134165479735097e-2, -.34317295181348742251e-1, + .59542472465494579941e-1, -.74135115907618101263e-1, + .74491937840566532596e-1, -.60042604725161994304e-1, + .33437677409000083169e-1, 0.0, + -.33437677409000083169e-1, .60042604725161994304e-1, + -.74491937840566532596e-1, .74135115907618101263e-1, + -.59542472465494579941e-1, .34317295181348742251e-1, + -.42204134165479735097e-2, -.24183738595685990149e-1, + .45040302741689006270e-1, -.54522171113149311354e-1, + .51647324405878909521e-1, -.38409681836626963278e-1, + .19201195309873585230e-1, .32917820356048383366e-3, + -.14701271465939718856e-1, .99848472706332791043e-2, + .11775579274769383373e-1, -.19892153937316935880e-1, + .95335114477449041055e-2, .57661528440359081617e-2, + -.23382690532380910781e-1, .40237257037170725321e-1, + -.53280289903551636474e-1, .59974361806023689068e-1, + -.58701684061992853224e-1, .49033407111597129616e-1, + -.31818835267847249219e-1, .90800541261162098886e-2, + .16272906819312603838e-1, -.40863896581186229487e-1, + .61346046297517367703e-1, + -.74896047554167268919e-1, .79632642148310325817e-1, + -.74896047554167268919e-1, .61346046297517367703e-1, + -.40863896581186229487e-1, .16272906819312603838e-1, + .90800541261162098886e-2, -.31818835267847249219e-1, + .49033407111597129616e-1, -.58701684061992853224e-1, + .59974361806023689068e-1, -.53280289903551636474e-1, + .40237257037170725321e-1, -.23382690532380910781e-1, + .57661528440359081617e-2, .95335114477449041055e-2, + -.19892153937316935880e-1, + .11775579274769383373e-1, -.13562702617218467450e-1, + .24885419969649845849e-1, -.18368693901908875583e-1, + .81673147806084084638e-2, .47890591326129587131e-2, + -.19313752945227974024e-1, .34065953398362954708e-1, + -.47667045133463415672e-1, .58820377816690514309e-1, + -.66424139824618415970e-1, + .69667606260856092515e-1, -.68102459384364543253e-1, + .61683024923302547971e-1, + -.50771943476441639136e-1, .36110771847327189215e-1, + -.18758028464284563358e-1, 0.0, .18758028464284563358e-1, + -.36110771847327189215e-1, .50771943476441639136e-1, + -.61683024923302547971e-1, .68102459384364543253e-1, + -.69667606260856092515e-1, .66424139824618415970e-1, + -.58820377816690514309e-1, .47667045133463415672e-1, + -.34065953398362954708e-1, .19313752945227974024e-1, + -.47890591326129587131e-2, -.81673147806084084638e-2, + .18368693901908875583e-1, -.24885419969649845849e-1, + .13562702617218467450e-1, .20576545037980523979e-1, + -.40093155172981004337e-1, .36954083167944054826e-1, + -.31856506837591907746e-1, .24996323181546255126e-1, + -.16637165210473614136e-1, .71002706773325085237e-2, + .32478629093205201133e-2, + -.14009562579050569518e-1, .24771262248780618922e-1, + -.35119395835433647559e-1, .44656290368574753171e-1, + -.53015448339647394161e-1, .59875631995693046782e-1, + -.64973208326045193862e-1, .68112280331082143373e-1, + -.69172215234062186994e-1, .68112280331082143373e-1, + -.64973208326045193862e-1, .59875631995693046782e-1, + -.53015448339647394161e-1, .44656290368574753171e-1, + -.35119395835433647559e-1, .24771262248780618922e-1, + -.14009562579050569518e-1, .32478629093205201133e-2, + .71002706773325085237e-2, -.16637165210473614136e-1, + .24996323181546255126e-1, -.31856506837591907746e-1, + .36954083167944054826e-1, -.40093155172981004337e-1, + .20576545037980523979e-1, -.27584914609096156163e-1, + .54904171411058497973e-1, -.54109756419563083153e-1, + .52794234894345577483e-1, -.50970276026831042415e-1, + .48655445537990983379e-1, + -.45872036510847994332e-1, .42646854695899611372e-1, + -.39010960357087507670e-1, .34999369144476467749e-1, + -.30650714874402762189e-1, .26006877464703437057e-1, + -.21112579608213651273e-1, .16014956068786763273e-1, + -.10763099747751940252e-1, .54075888924374485533e-2, 0.0, + -.54075888924374485533e-2, .10763099747751940252e-1, + -.16014956068786763273e-1, + .21112579608213651273e-1, -.26006877464703437057e-1, + .30650714874402762189e-1, + -.34999369144476467749e-1, .39010960357087507670e-1, + -.42646854695899611372e-1, .45872036510847994332e-1, + -.48655445537990983379e-1, .50970276026831042415e-1, + -.52794234894345577483e-1, .54109756419563083153e-1, + -.54904171411058497973e-1, .27584914609096156163e-1, + .13794141262469565740e-1, -.27588282524939131481e-1, + .27588282524939131481e-1, -.27588282524939131481e-1, + .27588282524939131481e-1, -.27588282524939131481e-1, + .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .27588282524939131481e-1, + -.27588282524939131481e-1, .13794141262469565740e-1 +}; + +static const double Tleft[33 * 33] = { + 1., -.86602540378443864678, 0., .33071891388307382381, 0., + -.20728904939721249057, 0., .15128841196122722208, 0., + -.11918864298744029244, 0., .98352013661686631224e-1, 0., + -.83727065404940845733e-1, 0., .72893399403505841203e-1, 0., + -.64544632643375022436e-1, 0., .57913170372415565639e-1, 0., + -.52518242575729562263e-1, 0., .48043311993977520457e-1, 0., + -.44271433659733990243e-1, 0., .41048928022856771981e-1, 0., + -.38263878662008271459e-1, 0., .35832844026365304501e-1, 0., 0., + .50000000000000000000, -.96824583655185422130, .57282196186948000082, + .21650635094610966169, -.35903516540862679125, -.97578093724974971969e-1, + .26203921611325660506, .55792409597991015609e-1, -.20644078533943456204, + -.36172381205961199479e-1, .17035068468874958194, + .25371838001497225980e-1, -.14501953125000000000, + -.18786835250972344757e-1, .12625507130328301066, + .14473795929590520582e-1, -.11179458309419422675, + -.11494434254897626155e-1, .10030855351241635862, + .93498556820544479096e-2, -.90964264465390582629e-1, + -.77546391824364392762e-2, .83213457337452292745e-1, + .65358085945588638605e-2, -.76680372422574234569e-1, + -.55835321940047427169e-2, .71098828931825789428e-1, + .48253327982967591019e-2, -.66274981937248958553e-1, + -.42118078245337801387e-2, .62064306433355646267e-1, + .37083386598903548973e-2, 0., 0., .25000000000000000000, + -.73950997288745200531, .83852549156242113615, -.23175620272173946716, + -.37791833195149451496, .25710129174850522325, .21608307321780204633, + -.22844049245646009157, -.14009503000335388415, .19897685605518413847, + .98264706042471226893e-1, -.17445445004279014046, + -.72761100054958328401e-1, .15463589893742108388, + .56056770591708784481e-1, -.13855313872640495158, + -.44517752443294564781e-1, .12534277657695128850, + .36211835346039665762e-1, -.11434398255136139683, + -.30033588409423828125e-1, .10506705408753910481, + .25313077840725783008e-1, -.97149327637744872155e-1, + -.21624927200393328444e-1, .90319582367202122625e-1, + .18688433567711780666e-1, -.84372291635345108584e-1, + -.16312261561845420752e-1, .79149526894804751586e-1, + .14362333871852474757e-1, 0., 0., 0., .12500000000000000000, + -.49607837082461073572, .82265291131801144317, -.59621200088559103072, + -.80054302859059362371e-1, .42612156697795759420, + -.90098145270865592887e-1, -.29769623255090078484, .13630307904779758221, + .21638835185708931831, -.14600247270306082052, -.16348801804014290453, + .14340708728599057249, .12755243353979286190, -.13661523715071346961, + -.10215585947881057394, .12864248070157166547, .83592528025348693602e-1, + -.12066728689302565222, -.69633728678718053052e-1, .11314245177331919532, + .58882939251410088028e-1, -.10621835858758221487, + -.50432266865187597572e-1, .99916834723527771581e-1, + .43672094283057258509e-1, -.94206380251950852413e-1, + -.38181356812697746418e-1, .89035739656537771225e-1, + .33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1, + -.31093357409581873586, .67604086414949799246, -.75644205980613611039, + .28990586430124175741, .30648508196770360914, -.35801372616842500052, + -.91326869828709014708e-1, .31127929687500000000, + -.90915752838698393094e-2, -.25637381283965534330, + .57601077850322797594e-1, .21019685709225757945, + -.81244992138514014256e-1, -.17375078516720988858, + .92289437277967051125e-1, .14527351914265391374, + -.96675340792832019889e-1, -.12289485697108543415, + .97448175340011084006e-1, .10511755943298339844, + -.96242247086378239657e-1, -.90822942272780513537e-1, + .93966350452322132384e-1, .79189411876493712558e-1, + -.91139307067989309325e-1, -.69613039934383197265e-1, + .88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0., + .31250000000000000000e-1, -.18684782411095934408, .50176689760410660236, + -.74784031498626095398, .56472001151566251186, .14842464993721351203e-1, + -.41162920273003120936, .20243071230196532282, .23772054897172750436, + -.24963810923972235950, -.12116179938394678936, .24330535483519110663, + .47903849781124471359e-1, -.22133299683101224293, + -.20542915138527200983e-2, .19653465717678146728, + -.26818172626509178444e-1, -.17319122357631210944, + .45065391411065545445e-1, .15253391395444065941, + -.56543897711725408302e-1, -.13469154928743585367, + .63632471400208840155e-1, .11941684923913523817, + -.67828850207933293098e-1, -.10636309084510652670, + .70095786922999181504e-1, .95187373095150709082e-1, 0., 0., 0., 0., 0., + 0., .15625000000000000000e-1, -.10909562534194485289, + .34842348626527747318, -.64461114561628111443, .69382480527334683659, + -.29551102358528827763, -.25527584713978439819, .38878771718544715394, + -.82956185835347407489e-2, -.31183177761966943912, .12831420840372374767, + .22067618205599434368, -.17569196937129496961, -.14598057000132284135, + .18864406621763419484, .89921002550386645767e-1, -.18571835020187122114, + -.48967672227195481777e-1, .17584685670380332798, + .19267984545067426324e-1, -.16335437520503462738, + .22598055455032407594e-2, .15032800884170631129, + -.17883358353754640871e-1, -.13774837869432209951, + .29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0., + 0., .78125000000000000000e-2, -.62377810244809812496e-1, + .23080781467370883845, -.50841310636012325368, .69834547012574056043, + -.52572723156526459672, .11464215704954976471e-1, .38698869011491210342, + -.26125646622255207507, -.16951698812361607510, .29773875898928782269, + .20130501202570367491e-1, -.26332493149159310198, + .67734613690401207009e-1, .21207315477103762715, -.11541543390889415193, + -.16249634759782417533, .13885887405041735068, .11996491328010275427, + -.14810432001630926895, -.85177658352556243411e-1, .14918860659904380587, + .57317789510444151564e-1, -.14569827645586660151, + -.35213090145965327390e-1, .13975998126844578198, 0., 0., 0., 0., 0., 0., + 0., 0., .39062500000000000000e-2, -.35101954600803571207e-1, + .14761284084133737720, -.37655033076080192966, .62410290231517322776, + -.64335622317683389875, .28188168266139524244, .22488495672137010675, + -.39393811089283576186, .75184777995770096714e-1, .28472023119398293003, + -.20410910833705899572, -.15590046962908511750, .23814567544617953125, + .54442805556829031204e-1, -.22855930338589720954, + .16303223615756629897e-1, .20172722433875559213, + -.62723406421217419404e-1, -.17012230831020922010, + .91754642766136561612e-1, .13927644821381121197, -.10886600968068418181, + -.11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0., + 0., 0., .19531250000000000000e-2, -.19506820659607596598e-1, + .91865676095362231937e-1, -.26604607809696493849, .51425874205091288223, + -.66047561132505329292, .48660109511591303851, -.17575661168678285615e-1, + -.36594333408055703366, .29088854695378694533, .11318677346656537927, + -.31110645235730182168, .60733219161008787341e-1, .24333848233620420826, + -.15254312332655419708, -.15995968483455388613, .19010344455215289289, + .86040636766440260000e-1, -.19652589954665259945, + -.27633388517205837713e-1, .18660848552712880387, + -.15942583868416775867e-1, -.16902042462382064786, + .47278526495327740646e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .97656250000000000000e-3, -.10731084460857378207e-1, + .55939644713816406331e-1, -.18118487371914493668, .39914857299829864263, + -.60812322949933902435, .60011887183061967583, -.26002695805835928795, + -.20883922404786010096, .38988130966114638081, -.11797833550782589082, + -.25231824756239520077, .24817859972953934712, .90516417677868996417e-1, + -.26079073291293066798, .30259468817169480161e-1, .22178195264114178432, + -.10569877864302048175, -.16679648389266977455, .14637718550245050850, + .11219272032739559870, -.16359363640525750353, -.64358194509092101393e-1, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3, + -.58542865274813470967e-2, .33461741635290096452e-1, + -.11979993155896201271, .29580223766987206958, -.51874761979436016742, + .62861483498014306968, -.44868895761051453296, .12567502628371529386e-1, + .35040366183235474275, -.30466868455569500886, -.70903913601490112666e-1, + .30822791893032512740, -.11969443264190207736, -.20764760317621313946, + .20629838355452128532, .95269702915334718507e-1, -.22432624768705133300, + -.33103381593477797101e-2, .20570036048155716333, + -.62208282720094518964e-1, -.17095309330441436348, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., .24414062500000000000e-3, + -.31714797501871532475e-2, .19721062526127334100e-1, + -.77311181185536498246e-1, .21124871792841566575, -.41777980401893650886, + .59401977834943551650, -.56132417807488349048, .23433675061367565951, + .20222775295220942126, -.38280372496506190127, .14443804214023095767, + .22268950939178466797, -.27211314150777981984, -.34184876506180717313e-1, + .26006498895669734842, -.97650425186005090107e-1, -.19024527660129101293, + .16789164198044635671, .10875811641651905252, -.19276785058805921298, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3, + -.17078941137247586143e-2, .11477733754843910060e-1, + -.48887017020924625462e-1, .14634927241421789683, -.32156282683019547854, + .52165811920227223937, -.60001958466396926460, .41208501541480733755, + -.11366945503190350975e-2, -.33968093962672089159, .30955190935923386766, + .40657421856578262210e-1, -.29873400409871531764, .16094481791768257440, + .16876122436206497694, -.23650217045022161255, -.33070260090574765012e-1, + .22985258456375907796, -.68645651043827097771e-1, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4, + -.91501857608428649078e-3, .66085179496951987952e-2, + -.30383171695850355404e-1, .98840838845366876117e-1, + -.23855447246420318989, .43322017468145613917, -.58049033744876107191, + .52533893203742699346, -.20681056202371946180, -.20180000924562504384, + .37503922291962681797, -.15988102869837429062, -.19823558102762374094, + .28393023878803799622, -.11188133439357510403e-1, -.24730368377168229255, + .14731529061377942839, .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .30517578125000000000e-4, + -.48804277318479845551e-3, .37696080990601968396e-2, + -.18603912108994738255e-1, .65325006755649582964e-1, + -.17162960707938819795, .34411527956476971322, -.52289350347082497959, + .57319653625674910592, -.37662253421045430413, -.14099055105384663902e-1, + .33265570610216904208, -.30921265572647566661, -.19911390594166455281e-1, + .28738590811031797718, -.18912130469738472647, -.13235936203215819193, + .25076406142356675279, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., .15258789062500000000e-4, -.25928719280954633249e-3, + .21327398937568540428e-2, -.11244626133630732010e-1, + .42375605740664331966e-1, -.12031130345907846211, .26352562258934426830, + -.44590628258512682078, .56682835613700749379, -.49116715128261660395, + .17845943097110339078, .20541650677432497477, -.36739803642257458221, + .16776034069210108273, .17920950989905112908, -.28867732805385066532, + .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., .76293945312500000000e-5, -.13727610943181290891e-3, + .11979683091449349286e-2, -.67195313034570709806e-2, + .27044920779931968175e-1, -.82472196498517457862e-1, + .19570475044896150093, -.36391620788543817693, .52241392782736588032, + -.54727504974907879912, .34211551468813581183, .31580472732719957762e-1, + -.32830006549176759667, .30563797665254420769, .64905014620683140120e-2, + -.27642986248995073032, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., .38146972656250000000e-5, + -.72454147007837596854e-4, .66859847582761390285e-3, + -.39751311980366118437e-2, .17015198650201528366e-1, + -.55443621868993855715e-1, .14157060481641692131, -.28641242619559616836, + .45610665490966615415, -.55262786406029265394, .45818352706035500108, + -.14984403004611673047, -.21163807462970713245, .36007252928843413718, + -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5, + -.38135049864067468562e-4, .37101393638555730015e-3, + -.23305339886279723213e-2, .10569913448297127219e-1, + -.36640175162216897547e-1, .10010476414320235508, -.21860074212675559892, + .38124757096345313719, -.52020999209879669177, .52172632730659212045, + -.30841620620308814614, -.50322546186721500184e-1, .32577618885114899053, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .95367431640625000000e-6, -.20021483206955925244e-4, + .20481807322420625431e-3, -.13553476938058909882e-2, + .64919676350791905019e-2, -.23848725425069251903e-1, + .69384632678886421292e-1, -.16249711393618776934, .30736618106830314788, + -.46399909601971539157, .53765031034002467225, -.42598991476520183929, + .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6, + -.10487707828484902486e-4, .11254146162337528943e-3, + -.78248929534271987118e-3, .39468337145306794566e-2, + -.15313546659475671763e-1, .47249070825218564146e-1, + -.11804374107101480543, .24031796927792491122, -.39629215049166341285, + .51629108968402548545, -.49622372075429782915, 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .23841857910156250000e-6, -.54823314130625337326e-5, + .61575377321535518154e-4, -.44877834366497538134e-3, + .23774612048621955857e-2, -.97136347645161687796e-2, + .31671599547606636717e-1, -.84028665767000747480e-1, + .18298487576742964949, -.32647878537696945218, .46970971486488895077, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., .11920928955078125000e-6, -.28604020001177375838e-5, + .33559227978295551013e-4, -.25583821662860610560e-3, + .14201552747787302339e-2, -.60938046986874414969e-2, + .20930869247951926793e-1, -.58745021125678072911e-1, + .13613725780285953720, -.26083988356030237586, 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .59604644775390625000e-7, -.14898180663526043291e-5, + .18224991282807693921e-4, -.14504433444608833821e-3, + .84184722720281809548e-3, -.37846965430000478789e-2, + .13656355548211376864e-1, -.40409541997718853934e-1, + .99226988101858325902e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .29802322387695312500e-7, -.77471708843445529468e-6, + .98649879372606876995e-5, -.81814934772838523887e-4, + .49554483992403011328e-3, -.23290922072351413938e-2, + .88068134250844034186e-2, -.27393666952485719070e-1, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., .14901161193847656250e-7, -.40226235946098233685e-6, + .53236418690561306700e-5, -.45933829691164002269e-4, + .28982005232838857913e-3, -.14212974043211018374e-2, + .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .74505805969238281250e-8, -.20858299254133430408e-6, + .28648457300134381744e-5, -.25677535898258910850e-4, + .16849420429491355445e-3, -.86062824010315834002e-3, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., .37252902984619140625e-8, -.10801736017613096861e-6, + .15376606719887104015e-5, -.14296523739727437959e-4, + .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .18626451492309570312e-8, -.55871592916438890146e-7, + .82331193828137454068e-6, -.79302250528382787666e-5, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9, + -.28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9, + -.14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .23283064365386962891e-9 +}; + +static const double Tright[33 * 33] = { + 1., .86602540378443864678, 0., -.33071891388307382381, 0., + .20728904939721249057, 0., -.15128841196122722208, 0., + .11918864298744029244, 0., -.98352013661686631224e-1, 0., + .83727065404940845733e-1, 0., -.72893399403505841203e-1, 0., + .64544632643375022436e-1, 0., -.57913170372415565639e-1, 0., + .52518242575729562263e-1, 0., -.48043311993977520457e-1, 0., + .44271433659733990243e-1, 0., -.41048928022856771981e-1, 0., + .38263878662008271459e-1, 0., -.35832844026365304501e-1, 0., 0., + .50000000000000000000, .96824583655185422130, .57282196186948000082, + -.21650635094610966169, -.35903516540862679125, .97578093724974971969e-1, + .26203921611325660506, -.55792409597991015609e-1, -.20644078533943456204, + .36172381205961199479e-1, .17035068468874958194, + -.25371838001497225980e-1, -.14501953125000000000, + .18786835250972344757e-1, .12625507130328301066, + -.14473795929590520582e-1, -.11179458309419422675, + .11494434254897626155e-1, .10030855351241635862, + -.93498556820544479096e-2, -.90964264465390582629e-1, + .77546391824364392762e-2, .83213457337452292745e-1, + -.65358085945588638605e-2, -.76680372422574234569e-1, + .55835321940047427169e-2, .71098828931825789428e-1, + -.48253327982967591019e-2, -.66274981937248958553e-1, + .42118078245337801387e-2, .62064306433355646267e-1, + -.37083386598903548973e-2, 0., 0., .25000000000000000000, + .73950997288745200531, .83852549156242113615, .23175620272173946716, + -.37791833195149451496, -.25710129174850522325, .21608307321780204633, + .22844049245646009157, -.14009503000335388415, -.19897685605518413847, + .98264706042471226893e-1, .17445445004279014046, + -.72761100054958328401e-1, -.15463589893742108388, + .56056770591708784481e-1, .13855313872640495158, + -.44517752443294564781e-1, -.12534277657695128850, + .36211835346039665762e-1, .11434398255136139683, + -.30033588409423828125e-1, -.10506705408753910481, + .25313077840725783008e-1, .97149327637744872155e-1, + -.21624927200393328444e-1, -.90319582367202122625e-1, + .18688433567711780666e-1, .84372291635345108584e-1, + -.16312261561845420752e-1, -.79149526894804751586e-1, + .14362333871852474757e-1, 0., 0., 0., .12500000000000000000, + .49607837082461073572, .82265291131801144317, .59621200088559103072, + -.80054302859059362371e-1, -.42612156697795759420, + -.90098145270865592887e-1, .29769623255090078484, .13630307904779758221, + -.21638835185708931831, -.14600247270306082052, .16348801804014290453, + .14340708728599057249, -.12755243353979286190, -.13661523715071346961, + .10215585947881057394, .12864248070157166547, -.83592528025348693602e-1, + -.12066728689302565222, .69633728678718053052e-1, .11314245177331919532, + -.58882939251410088028e-1, -.10621835858758221487, + .50432266865187597572e-1, .99916834723527771581e-1, + -.43672094283057258509e-1, -.94206380251950852413e-1, + .38181356812697746418e-1, .89035739656537771225e-1, + -.33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1, + .31093357409581873586, .67604086414949799246, .75644205980613611039, + .28990586430124175741, -.30648508196770360914, -.35801372616842500052, + .91326869828709014708e-1, .31127929687500000000, .90915752838698393094e-2, + -.25637381283965534330, -.57601077850322797594e-1, .21019685709225757945, + .81244992138514014256e-1, -.17375078516720988858, + -.92289437277967051125e-1, .14527351914265391374, + .96675340792832019889e-1, -.12289485697108543415, + -.97448175340011084006e-1, .10511755943298339844, + .96242247086378239657e-1, -.90822942272780513537e-1, + -.93966350452322132384e-1, .79189411876493712558e-1, + .91139307067989309325e-1, -.69613039934383197265e-1, + -.88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0., + .31250000000000000000e-1, .18684782411095934408, .50176689760410660236, + .74784031498626095398, .56472001151566251186, -.14842464993721351203e-1, + -.41162920273003120936, -.20243071230196532282, .23772054897172750436, + .24963810923972235950, -.12116179938394678936, -.24330535483519110663, + .47903849781124471359e-1, .22133299683101224293, + -.20542915138527200983e-2, -.19653465717678146728, + -.26818172626509178444e-1, .17319122357631210944, + .45065391411065545445e-1, -.15253391395444065941, + -.56543897711725408302e-1, .13469154928743585367, + .63632471400208840155e-1, -.11941684923913523817, + -.67828850207933293098e-1, .10636309084510652670, + .70095786922999181504e-1, -.95187373095150709082e-1, 0., 0., 0., 0., 0., + 0., .15625000000000000000e-1, .10909562534194485289, + .34842348626527747318, .64461114561628111443, .69382480527334683659, + .29551102358528827763, -.25527584713978439819, -.38878771718544715394, + -.82956185835347407489e-2, .31183177761966943912, .12831420840372374767, + -.22067618205599434368, -.17569196937129496961, .14598057000132284135, + .18864406621763419484, -.89921002550386645767e-1, -.18571835020187122114, + .48967672227195481777e-1, .17584685670380332798, + -.19267984545067426324e-1, -.16335437520503462738, + -.22598055455032407594e-2, .15032800884170631129, + .17883358353754640871e-1, -.13774837869432209951, + -.29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0., + 0., .78125000000000000000e-2, .62377810244809812496e-1, + .23080781467370883845, .50841310636012325368, .69834547012574056043, + .52572723156526459672, .11464215704954976471e-1, -.38698869011491210342, + -.26125646622255207507, .16951698812361607510, .29773875898928782269, + -.20130501202570367491e-1, -.26332493149159310198, + -.67734613690401207009e-1, .21207315477103762715, .11541543390889415193, + -.16249634759782417533, -.13885887405041735068, .11996491328010275427, + .14810432001630926895, -.85177658352556243411e-1, -.14918860659904380587, + .57317789510444151564e-1, .14569827645586660151, + -.35213090145965327390e-1, -.13975998126844578198, 0., 0., 0., 0., 0., 0., + 0., 0., .39062500000000000000e-2, .35101954600803571207e-1, + .14761284084133737720, .37655033076080192966, .62410290231517322776, + .64335622317683389875, .28188168266139524244, -.22488495672137010675, + -.39393811089283576186, -.75184777995770096714e-1, .28472023119398293003, + .20410910833705899572, -.15590046962908511750, -.23814567544617953125, + .54442805556829031204e-1, .22855930338589720954, .16303223615756629897e-1, + -.20172722433875559213, -.62723406421217419404e-1, .17012230831020922010, + .91754642766136561612e-1, -.13927644821381121197, -.10886600968068418181, + .11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0., + 0., 0., .19531250000000000000e-2, .19506820659607596598e-1, + .91865676095362231937e-1, .26604607809696493849, .51425874205091288223, + .66047561132505329292, .48660109511591303851, .17575661168678285615e-1, + -.36594333408055703366, -.29088854695378694533, .11318677346656537927, + .31110645235730182168, .60733219161008787341e-1, -.24333848233620420826, + -.15254312332655419708, .15995968483455388613, .19010344455215289289, + -.86040636766440260000e-1, -.19652589954665259945, + .27633388517205837713e-1, .18660848552712880387, .15942583868416775867e-1, + -.16902042462382064786, -.47278526495327740646e-1, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., .97656250000000000000e-3, .10731084460857378207e-1, + .55939644713816406331e-1, .18118487371914493668, .39914857299829864263, + .60812322949933902435, .60011887183061967583, .26002695805835928795, + -.20883922404786010096, -.38988130966114638081, -.11797833550782589082, + .25231824756239520077, .24817859972953934712, -.90516417677868996417e-1, + -.26079073291293066798, -.30259468817169480161e-1, .22178195264114178432, + .10569877864302048175, -.16679648389266977455, -.14637718550245050850, + .11219272032739559870, .16359363640525750353, -.64358194509092101393e-1, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3, + .58542865274813470967e-2, .33461741635290096452e-1, .11979993155896201271, + .29580223766987206958, .51874761979436016742, .62861483498014306968, + .44868895761051453296, .12567502628371529386e-1, -.35040366183235474275, + -.30466868455569500886, .70903913601490112666e-1, .30822791893032512740, + .11969443264190207736, -.20764760317621313946, -.20629838355452128532, + .95269702915334718507e-1, .22432624768705133300, + -.33103381593477797101e-2, -.20570036048155716333, + -.62208282720094518964e-1, .17095309330441436348, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., .24414062500000000000e-3, + .31714797501871532475e-2, .19721062526127334100e-1, + .77311181185536498246e-1, .21124871792841566575, .41777980401893650886, + .59401977834943551650, .56132417807488349048, .23433675061367565951, + -.20222775295220942126, -.38280372496506190127, -.14443804214023095767, + .22268950939178466797, .27211314150777981984, -.34184876506180717313e-1, + -.26006498895669734842, -.97650425186005090107e-1, .19024527660129101293, + .16789164198044635671, -.10875811641651905252, -.19276785058805921298, 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3, + .17078941137247586143e-2, .11477733754843910060e-1, + .48887017020924625462e-1, .14634927241421789683, .32156282683019547854, + .52165811920227223937, .60001958466396926460, .41208501541480733755, + .11366945503190350975e-2, -.33968093962672089159, -.30955190935923386766, + .40657421856578262210e-1, .29873400409871531764, .16094481791768257440, + -.16876122436206497694, -.23650217045022161255, .33070260090574765012e-1, + .22985258456375907796, .68645651043827097771e-1, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4, + .91501857608428649078e-3, .66085179496951987952e-2, + .30383171695850355404e-1, .98840838845366876117e-1, .23855447246420318989, + .43322017468145613917, .58049033744876107191, .52533893203742699346, + .20681056202371946180, -.20180000924562504384, -.37503922291962681797, + -.15988102869837429062, .19823558102762374094, .28393023878803799622, + .11188133439357510403e-1, -.24730368377168229255, -.14731529061377942839, + .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .30517578125000000000e-4, .48804277318479845551e-3, + .37696080990601968396e-2, .18603912108994738255e-1, + .65325006755649582964e-1, .17162960707938819795, .34411527956476971322, + .52289350347082497959, .57319653625674910592, .37662253421045430413, + -.14099055105384663902e-1, -.33265570610216904208, -.30921265572647566661, + .19911390594166455281e-1, .28738590811031797718, .18912130469738472647, + -.13235936203215819193, -.25076406142356675279, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .15258789062500000000e-4, + .25928719280954633249e-3, .21327398937568540428e-2, + .11244626133630732010e-1, .42375605740664331966e-1, .12031130345907846211, + .26352562258934426830, .44590628258512682078, .56682835613700749379, + .49116715128261660395, .17845943097110339078, -.20541650677432497477, + -.36739803642257458221, -.16776034069210108273, .17920950989905112908, + .28867732805385066532, .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .76293945312500000000e-5, + .13727610943181290891e-3, .11979683091449349286e-2, + .67195313034570709806e-2, .27044920779931968175e-1, + .82472196498517457862e-1, .19570475044896150093, .36391620788543817693, + .52241392782736588032, .54727504974907879912, .34211551468813581183, + -.31580472732719957762e-1, -.32830006549176759667, -.30563797665254420769, + .64905014620683140120e-2, .27642986248995073032, 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .38146972656250000000e-5, + .72454147007837596854e-4, .66859847582761390285e-3, + .39751311980366118437e-2, .17015198650201528366e-1, + .55443621868993855715e-1, .14157060481641692131, .28641242619559616836, + .45610665490966615415, .55262786406029265394, .45818352706035500108, + .14984403004611673047, -.21163807462970713245, -.36007252928843413718, + -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5, + .38135049864067468562e-4, .37101393638555730015e-3, + .23305339886279723213e-2, .10569913448297127219e-1, + .36640175162216897547e-1, .10010476414320235508, .21860074212675559892, + .38124757096345313719, .52020999209879669177, .52172632730659212045, + .30841620620308814614, -.50322546186721500184e-1, -.32577618885114899053, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .95367431640625000000e-6, .20021483206955925244e-4, + .20481807322420625431e-3, .13553476938058909882e-2, + .64919676350791905019e-2, .23848725425069251903e-1, + .69384632678886421292e-1, .16249711393618776934, .30736618106830314788, + .46399909601971539157, .53765031034002467225, .42598991476520183929, + .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6, + .10487707828484902486e-4, .11254146162337528943e-3, + .78248929534271987118e-3, .39468337145306794566e-2, + .15313546659475671763e-1, .47249070825218564146e-1, .11804374107101480543, + .24031796927792491122, .39629215049166341285, .51629108968402548545, + .49622372075429782915, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., .23841857910156250000e-6, + .54823314130625337326e-5, .61575377321535518154e-4, + .44877834366497538134e-3, .23774612048621955857e-2, + .97136347645161687796e-2, .31671599547606636717e-1, + .84028665767000747480e-1, .18298487576742964949, .32647878537696945218, + .46970971486488895077, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .11920928955078125000e-6, + .28604020001177375838e-5, .33559227978295551013e-4, + .25583821662860610560e-3, .14201552747787302339e-2, + .60938046986874414969e-2, .20930869247951926793e-1, + .58745021125678072911e-1, .13613725780285953720, .26083988356030237586, + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., .59604644775390625000e-7, + .14898180663526043291e-5, .18224991282807693921e-4, + .14504433444608833821e-3, .84184722720281809548e-3, + .37846965430000478789e-2, .13656355548211376864e-1, + .40409541997718853934e-1, .99226988101858325902e-1, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .29802322387695312500e-7, .77471708843445529468e-6, + .98649879372606876995e-5, .81814934772838523887e-4, + .49554483992403011328e-3, .23290922072351413938e-2, + .88068134250844034186e-2, .27393666952485719070e-1, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., .14901161193847656250e-7, .40226235946098233685e-6, + .53236418690561306700e-5, .45933829691164002269e-4, + .28982005232838857913e-3, .14212974043211018374e-2, + .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .74505805969238281250e-8, .20858299254133430408e-6, + .28648457300134381744e-5, .25677535898258910850e-4, + .16849420429491355445e-3, .86062824010315834002e-3, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., .37252902984619140625e-8, .10801736017613096861e-6, + .15376606719887104015e-5, .14296523739727437959e-4, + .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + .18626451492309570312e-8, .55871592916438890146e-7, + .82331193828137454068e-6, .79302250528382787666e-5, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9, + .28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9, + .14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., .23283064365386962891e-9 +}; + +/* Allocates a workspace for the given maximum number of intervals. + Note that if the workspace gets filled, the intervals with the + lowest error estimates are dropped. The maximum number of + intervals is therefore not the maximum number of intervals + that will be computed, but merely the size of the buffer. + */ + +/* Compute the product of the fx with one of the inverse + Vandermonde-like matrices. */ + +void +Vinvfx (const double *fx, double *c, const int d) +{ + + int i, j; + + switch (d) + { + case 0: + for (i = 0; i <= 4; i++) + { + c[i] = 0.0; + for (j = 0; j <= 4; j++) + c[i] += V1inv[i * 5 + j] * fx[j * 8]; + } + break; + case 1: + for (i = 0; i <= 8; i++) + { + c[i] = 0.0; + for (j = 0; j <= 8; j++) + c[i] += V2inv[i * 9 + j] * fx[j * 4]; + } + break; + case 2: + for (i = 0; i <= 16; i++) + { + c[i] = 0.0; + for (j = 0; j <= 16; j++) + c[i] += V3inv[i * 17 + j] * fx[j * 2]; + } + break; + case 3: + for (i = 0; i <= 32; i++) + { + c[i] = 0.0; + for (j = 0; j <= 32; j++) + c[i] += V4inv[i * 33 + j] * fx[j]; + } + break; + } + +} + + +/* Downdate the interpolation given by the n coefficients c + by removing the nodes with indices in nans. */ + +void +downdate (double *c, int n, int d, int *nans, int nnans) +{ + + static const int bidx[4] = { 0, 6, 16, 34 }; + double b_new[34], alpha; + int i, j; + + for (i = 0; i <= n + 1; i++) + b_new[i] = bee[bidx[d] + i]; + for (i = 0; i < nnans; i++) + { + b_new[n + 1] = b_new[n + 1] / Lalpha[n]; + b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1]; + for (j = n - 1; j > 0; j--) + b_new[j] = + (b_new[j] + xi[nans[i]] * b_new[j + 1] - + Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1]; + for (j = 0; j <= n; j++) + b_new[j] = b_new[j + 1]; + alpha = c[n] / b_new[n]; + for (j = 0; j < n; j++) + c[j] -= alpha * b_new[j]; + c[n] = 0; + n--; + } + +} + + + +/* The actual integration routine. + */ + +DEFUN_DLD (quadcc, args, nargout, +"-*- texinfo -*-\n\ +@deftypefn {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})\n\ +@deftypefnx {Function File} {[@var{int}, @var{err}, @var{nr_points}] =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})\n\ +Numerically evaluates an integral using the doubly-adaptive\n\ +quadrature described by P. Gonnet in @cite{Increasing the\n\ +Reliability of Adaptive Quadrature Using Explicit Interpolants},\n\ +ACM Transactions on Mathematical Software, in Press, 2010.\n\ +The algorithm uses Clenshaw-Curtis quadrature rules of increasing\n\ +degree in each interval and bisects the interval if either the\n\ +function does not appear to be smooth or a rule of maximum\n\ +degree has been reached. The error estimate is computed from the\n\ +L2-norm of the difference between two successive interpolations\n\ +of the integrand over the nodes of the respective quadrature rules.\n\ +\n\ +For example,\n\ +\n\ +@example\n\ + int = quadcc (f, a, b, 1.0e-6);\n\ +@end example\n\ +\n\ +@noindent\n\ +computes the integral of a function @var{f} in the interval\n\ +[@var{a}, @var{b}] to the relative precision of six\n\ +decimal digits.\n\ +The integrand @var{f} should accept a vector argument and return a vector\n\ +result containing the integrand evaluated at each element of the\n\ +argument, for example:\n\ +\n\ +@example\n\ + f = @@(x) x .* sin (1 ./ x) .* sqrt (abs (1 - x));\n\ +@end example\n\ +\n\ +If the integrand has known singularities or discontinuities\n\ +in any of its derivatives inside the interval,\n\ +as does the above example at x=1, these can be specified in\n\ +the additional argument @var{sing} as follows\n\ +\n\ +@example\n\ + int = quadcc (f, a, b, 1.0e-6, [ 1 ]);\n\ +@end example\n\ +\n\ +The two additional output variables @var{err} and @var{nr_points}\n\ +return an estimate of the absolute integration error and\n\ +the number of points at which the integrand was evaluated\n\ +respectively.\n\ +If the adaptive integration did not converge, the value of\n\ +@var{err} will be larger than the requested tolerance. It is\n\ +therefore recommended to verify this value for difficult\n\ +integrands.\n\ +\n\ +If either @var{a} or @var{b} are @code{+/-Inf}, @code{quadcc}\n\ +integrates @var{f} by substituting the variable of integration\n\ +with @code{x=tan(pi/2*u)}.\n\ +\n\ +@code{quadcc} is capable of dealing with non-numeric\n\ +values of the integrand such as @code{NaN}, @code{Inf}\n\ +or @code{-Inf}, as in the above example at x=0.\n\ +If the integral diverges and @code{quadcc} detects this, \n\ +a warning is issued and @code{Inf} or @code{-Inf} is returned.\n\ +\n\ +Note that @code{quadcc} is a general purpose quadrature algorithm\n\ +and as such may be less efficient for smooth or otherwise\n\ +well-behaved integrand than other methods such as\n\ +@code{quadgk} or @code{trapz}.\n\ +\n\ +@seealso{quad,quadv,quadl,quadgk,trapz,dblquad,triplequad}\n\ +@end deftypefn") +{ + + /* Some constants that we will need. */ + static const int n[4] = { 4, 8, 16, 32 }; + static const int skip[4] = { 8, 4, 2, 1 }; + static const int idx[4] = { 0, 5, 14, 31 }; + static const double w = M_SQRT2 / 2; + static const int ndiv_max = 20; + + /* The interval heap. */ + cquad_ival ivals[cquad_heapsize]; + int heap[cquad_heapsize]; + + /* Arguments left and right */ + int nargin = args.length (); + octave_function *fcn; + double a, b, tol, iivals[cquad_heapsize], *sing; + + /* Variables needed for transforming the integrand. */ + int wrap = 0; + double xw; + + /* Stuff we will need to call the integrand. */ + octave_value_list fargs, retval; + + /* Actual variables (as opposed to constants above). */ + double m, h, ml, hl, mr, hr, temp; + double igral, err, igral_final, err_final, err_excess; + int nivals, neval = 0; + int i, j, d, split, t; + int nnans, nans[33]; + cquad_ival *iv, *ivl, *ivr; + double nc, ncdiff; + + + /* Parse the input arguments. */ + if (nargin < 1) + { + error + ("quadcc: first argument (integrand) of type function handle required."); + return octave_value_list (); + } + else + { + if (args (0).is_function_handle () || args (0).is_inline_function ()) + fcn = args (0).function_value (); + else { + error + ("quadcc: first argument (integrand) must be a function handle or an inline function."); + return octave_value_list(); + } + } + + if (nargin < 2 || !args (1).is_real_scalar ()) + { + error + ("quadcc: second argument (left interval edge) must be a single real scalar."); + return octave_value_list (); + } + else + a = args (1).double_value (); + + if (nargin < 3 || !args (2).is_real_scalar ()) + { + error + ("quadcc: third argument (right interval edge) must be a single real scalar."); + return octave_value_list (); + } + else + b = args (2).double_value (); + + if (nargin < 4) + tol = 1.0e-6; + else if (!args (3).is_real_scalar ()) + { + error + ("quadcc: fourth argument (tolerance) must be a single real scalar."); + return octave_value_list (); + } + else + tol = args (3).double_value (); + + if (nargin < 5) + { + nivals = 1; + iivals[0] = a; + iivals[1] = b; + } + else if (!(args (4).is_real_scalar () || args (4).is_real_matrix ())) + { + error + ("quadcc: fifth argument (singularities) must be a vector of real values."); + return octave_value_list (); + } + else + { + nivals = 1 + args (4).length (); + if ( nivals > cquad_heapsize ) { + error("quadcc: maximum number of singular points is limited to %i.",cquad_heapsize-1); + return octave_value_list(); + } + sing = args (4).array_value ().fortran_vec (); + iivals[0] = a; + for (i = 0; i < nivals - 2; i++) + iivals[i + 1] = sing[i]; + iivals[nivals] = b; + } + + /* If a or b are +/-Inf, transform the integral. */ + if (xisinf (a) || xisinf (b)) + { + wrap = 1; + for (i = 0; i <= nivals; i++) + if (xisinf (iivals[i])) + iivals[i] = copysign (1.0, iivals[i]); + else + iivals[i] = 2.0 * atan (iivals[i]) / M_PI; + } + + + /* Initialize the heaps. */ + for (i = 0; i < cquad_heapsize; i++) + heap[i] = i; + + + /* Create the first interval(s). */ + igral = 0.0; + err = 0.0; + for (j = 0; j < nivals; j++) + { + + /* Initialize the interval. */ + iv = &(ivals[heap[j]]); + m = (iivals[j] + iivals[j + 1]) / 2; + h = (iivals[j + 1] - iivals[j]) / 2; + nnans = 0; + ColumnVector ex (33); + if (wrap) + { + for (i = 0; i <= n[3]; i++) + ex (i) = tan (M_PI / 2 * (m + xi[i] * h)); + } + else + { + for (i = 0; i <= n[3]; i++) + ex (i) = m + xi[i] * h; + } + fargs (0) = ex; + retval = feval (fcn, fargs, 1); + if (retval.length () != 1 || !retval (0).is_real_matrix ()) + { + error + ("quadcc: integrand must return a single, real-valued vector."); + return octave_value_list (); + } + Matrix effex = retval (0).matrix_value (); + if (effex.length () != ex.length ()) + { + error + ("quadcc: integrand must return a single, real-valued vector of the same size as the input"); + return octave_value_list (); + } + for (i = 0; i <= n[3]; i++) + { + iv->fx[i] = effex (i); + if (wrap) + { + xw = ex(i); + iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2; + } + neval++; + if (!xfinite (iv->fx[i])) + { + nans[nnans++] = i; + iv->fx[i] = 0.0; + } + } + Vinvfx (iv->fx, &(iv->c[idx[3]]), 3); + Vinvfx (iv->fx, &(iv->c[idx[2]]), 2); + Vinvfx (iv->fx, &(iv->c[0]), 0); + for (i = 0; i < nnans; i++) + iv->fx[i] = octave_NaN; + iv->a = iivals[j]; + iv->b = iivals[j + 1]; + iv->depth = 3; + iv->rdepth = 1; + iv->ndiv = 0; + iv->igral = 2 * h * iv->c[idx[3]] * w; + nc = 0.0; + for (i = n[2] + 1; i <= n[3]; i++) + { + temp = iv->c[idx[3] + i]; + nc += temp * temp; + } + ncdiff = nc; + for (i = 0; i <= n[2]; i++) + { + temp = iv->c[idx[2] + i] - iv->c[idx[3] + i]; + ncdiff += temp * temp; + nc += iv->c[idx[3] + i] * iv->c[idx[3] + i]; + } + ncdiff = sqrt (ncdiff); + nc = sqrt (nc); + iv->err = ncdiff * 2 * h; + if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc) + iv->err = 2 * h * nc; + + /* Tabulate this interval's data. */ + igral += iv->igral; + err += iv->err; + + /* Sift it up the heap. */ + i = j; + while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err) + { + temp = heap[i]; + heap[i] = heap[i / 2]; + heap[i / 2] = temp; + i /= 2; + } + + } + + + /* Initialize some global values. */ + igral_final = 0.0; + err_final = 0.0; + err_excess = 0.0; + + + /* Main loop. */ + while (nivals > 0 && err > 0.0 && err > fabs (igral) * tol + && !(err_final > fabs (igral) * tol + && err - err_final < fabs (igral) * tol)) + { + + /* Allow the user to interrupt. */ + OCTAVE_QUIT; + + /* Put our finger on the interval with the largest error. */ + iv = &(ivals[heap[0]]); + m = (iv->a + iv->b) / 2; + h = (iv->b - iv->a) / 2; + +/* printf + ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", + heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth); +*/ + /* Should we try to increase the degree? */ + if (iv->depth < 3) + { + + /* Keep tabs on some variables. */ + d = ++iv->depth; + + /* Get the new (missing) function values */ + { + ColumnVector ex (n[d] / 2); + if (wrap) + { + for (i = 0; i < n[d] / 2; i++) + ex (i) = + tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h)); + } + else + { + for (i = 0; i < n[d] / 2; i++) + ex (i) = m + xi[(2 * i + 1) * skip[d]] * h; + } + fargs (0) = ex; + retval = feval (fcn, fargs, 1); + if (retval.length () != 1 || !retval (0).is_real_matrix ()) + { + error + ("quadcc: integrand must return a single, real-valued vector."); + return octave_value_list (); + } + Matrix effex = retval (0).matrix_value (); + if (effex.length () != ex.length ()) + { + error + ("quadcc: integrand must return a single, real-valued vector of the same size as the input."); + return octave_value_list (); + } + neval += effex.length (); + for (i = 0; i < n[d] / 2; i++) + { + j = (2 * i + 1) * skip[d]; + iv->fx[j] = effex (i); + if (wrap) + { + xw = ex(i); + iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2; + } + } + } + nnans = 0; + for (i = 0; i <= 32; i += skip[d]) + { + if (!xfinite (iv->fx[i])) + { + nans[nnans++] = i; + iv->fx[i] = 0.0; + } + } + + /* Compute the new coefficients. */ + Vinvfx (iv->fx, &(iv->c[idx[d]]), d); + /* Downdate any NaNs. */ + if (nnans > 0) + { + downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans); + for (i = 0; i < nnans; i++) + iv->fx[i] = octave_NaN; + } + + /* Compute the error estimate. */ + nc = 0.0; + for (i = n[d - 1] + 1; i <= n[d]; i++) + { + temp = iv->c[idx[d] + i]; + nc += temp * temp; + } + ncdiff = nc; + for (i = 0; i <= n[d - 1]; i++) + { + temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i]; + ncdiff += temp * temp; + nc += iv->c[idx[d] + i] * iv->c[idx[d] + i]; + } + ncdiff = sqrt (ncdiff); + nc = sqrt (nc); + iv->err = ncdiff * 2 * h; + /* Compute the local integral. */ + iv->igral = 2 * h * w * iv->c[idx[d]]; + /* Split the interval prematurely? */ + split = (nc > 0 && ncdiff / nc > 0.1); + } + + /* Maximum degree reached, just split. */ + else + { + split = 1; + } + + + /* Should we drop this interval? */ + if ((m + h * xi[0]) >= (m + h * xi[1]) + || (m + h * xi[31]) >= (m + h * xi[32]) + || iv->err < fabs (iv->igral) * DBL_EPSILON * 10) + { + +/* printf + ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", + heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, + iv->depth); +*/ + /* Keep this interval's contribution */ + err_final += iv->err; + igral_final += iv->igral; + /* Swap with the last element on the heap */ + t = heap[nivals - 1]; + heap[nivals - 1] = heap[0]; + heap[0] = t; + nivals--; + /* Fix up the heap */ + i = 0; + while (2 * i + 1 < nivals) + { + + /* Get the kids */ + j = 2 * i + 1; + /* If the j+1st entry exists and is larger than the jth, + use it instead. */ + if (j + 1 < nivals + && ivals[heap[j + 1]].err >= ivals[heap[j]].err) + j++; + /* Do we need to move the ith entry up? */ + if (ivals[heap[j]].err <= ivals[heap[i]].err) + break; + else + { + t = heap[j]; + heap[j] = heap[i]; + heap[i] = t; + i = j; + } + } + + } + + /* Do we need to split this interval? */ + else if (split) + { + + /* Some values we will need often... */ + d = iv->depth; + /* Generate the interval on the left */ + ivl = &(ivals[heap[nivals++]]); + ivl->a = iv->a; + ivl->b = m; + ml = (ivl->a + ivl->b) / 2; + hl = h / 2; + ivl->depth = 0; + ivl->rdepth = iv->rdepth + 1; + ivl->fx[0] = iv->fx[0]; + ivl->fx[32] = iv->fx[16]; + { + ColumnVector ex (n[0] - 1); + if (wrap) + { + for (i = 0; i < n[0] - 1; i++) + ex (i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl)); + } + else + { + for (i = 0; i < n[0] - 1; i++) + ex (i) = ml + xi[(i + 1) * skip[0]] * hl; + } + fargs (0) = ex; + retval = feval (fcn, fargs, 1); + if (retval.length () != 1 || !retval (0).is_real_matrix ()) + { + error + ("quadcc: integrand must return a single, real-valued vector."); + return octave_value_list (); + } + Matrix effex = retval (0).matrix_value (); + if (effex.length () != ex.length ()) + { + error + ("quadcc: integrand must return a single, real-valued vector of the same size as the input."); + return octave_value_list (); + } + neval += effex.length (); + for (i = 0; i < n[0] - 1; i++) + { + j = (i + 1) * skip[0]; + ivl->fx[j] = effex (i); + if (wrap) + { + xw = ex(i); + ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2; + } + } + } + nnans = 0; + for (i = 0; i <= 32; i += skip[0]) + { + if (!xfinite (ivl->fx[i])) + { + nans[nnans++] = i; + ivl->fx[i] = 0.0; + } + } + Vinvfx (ivl->fx, ivl->c, 0); + if (nnans > 0) + { + downdate (ivl->c, n[0], 0, nans, nnans); + for (i = 0; i < nnans; i++) + ivl->fx[i] = octave_NaN; + } + for (i = 0; i <= n[d]; i++) + { + ivl->c[idx[d] + i] = 0.0; + for (j = i; j <= n[d]; j++) + ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j]; + } + ncdiff = 0.0; + for (i = 0; i <= n[0]; i++) + { + temp = ivl->c[i] - ivl->c[idx[d] + i]; + ncdiff += temp * temp; + } + for (i = n[0] + 1; i <= n[d]; i++) + { + temp = ivl->c[idx[d] + i]; + ncdiff += temp * temp; + } + ncdiff = sqrt (ncdiff); + ivl->err = ncdiff * h; + /* Check for divergence. */ + ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0 + && ivl->c[0] / iv->c[0] > 2); + if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth) + { + igral = copysign (octave_Inf, igral); + warning ("quadcc: divergent integral detected."); + break; + } + + /* Compute the local integral. */ + ivl->igral = h * w * ivl->c[0]; + + + /* Generate the interval on the right */ + ivr = &(ivals[heap[nivals++]]); + ivr->a = m; + ivr->b = iv->b; + mr = (ivr->a + ivr->b) / 2; + hr = h / 2; + ivr->depth = 0; + ivr->rdepth = iv->rdepth + 1; + ivr->fx[0] = iv->fx[16]; + ivr->fx[32] = iv->fx[32]; + { + ColumnVector ex (n[0] - 1); + if (wrap) + { + for (i = 0; i < n[0] - 1; i++) + ex (i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr)); + } + else + { + for (i = 0; i < n[0] - 1; i++) + ex (i) = mr + xi[(i + 1) * skip[0]] * hr; + } + fargs (0) = ex; + retval = feval (fcn, fargs, 1); + if (retval.length () != 1 || !retval (0).is_real_matrix ()) + { + error + ("quadcc: integrand must return a single, real-valued vector."); + return octave_value_list (); + } + Matrix effex = retval (0).matrix_value (); + if (effex.length () != ex.length ()) + { + error + ("quadcc: integrand must return a single, real-valued vector of the same size as the input."); + return octave_value_list (); + } + neval += effex.length (); + for (i = 0; i < n[0] - 1; i++) + { + j = (i + 1) * skip[0]; + ivr->fx[j] = effex (i); + if (wrap) + { + xw = ex(i); + ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2; + } + } + } + nnans = 0; + for (i = 0; i <= 32; i += skip[0]) + { + if (!xfinite (ivr->fx[i])) + { + nans[nnans++] = i; + ivr->fx[i] = 0.0; + } + } + Vinvfx (ivr->fx, ivr->c, 0); + if (nnans > 0) + { + downdate (ivr->c, n[0], 0, nans, nnans); + for (i = 0; i < nnans; i++) + ivr->fx[i] = octave_NaN; + } + for (i = 0; i <= n[d]; i++) + { + ivr->c[idx[d] + i] = 0.0; + for (j = i; j <= n[d]; j++) + ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j]; + } + ncdiff = 0.0; + for (i = 0; i <= n[0]; i++) + { + temp = ivr->c[i] - ivr->c[idx[d] + i]; + ncdiff += temp * temp; + } + for (i = n[0] + 1; i <= n[d]; i++) + { + temp = ivr->c[idx[d] + i]; + ncdiff += temp * temp; + } + ncdiff = sqrt (ncdiff); + ivr->err = ncdiff * h; + /* Check for divergence. */ + ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0 + && ivr->c[0] / iv->c[0] > 2); + if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth) + { + igral = copysign (octave_Inf, igral); + warning ("quadcc: divergent integral detected."); + break; + } + + /* Compute the local integral. */ + ivr->igral = h * w * ivr->c[0]; + + + /* Fix-up the heap: we now have one interval on top + that we don't need any more and two new, unsorted + ones at the bottom. */ + /* Flip the last interval to the top of the heap and + sift down. */ + t = heap[nivals - 1]; + heap[nivals - 1] = heap[0]; + heap[0] = t; + nivals--; + /* Sift this interval back down the heap. */ + i = 0; + while (2 * i + 1 < nivals - 1) + { + j = 2 * i + 1; + if (j + 1 < nivals - 1 + && ivals[heap[j + 1]].err >= ivals[heap[j]].err) + j++; + if (ivals[heap[j]].err <= ivals[heap[i]].err) + break; + else + { + t = heap[j]; + heap[j] = heap[i]; + heap[i] = t; + i = j; + } + } + + /* Now grab the last interval and sift it up the heap. */ + i = nivals - 1; + while (i > 0) + { + j = (i - 1) / 2; + if (ivals[heap[j]].err < ivals[heap[i]].err) + { + t = heap[j]; + heap[j] = heap[i]; + heap[i] = t; + i = j; + } + else + break; + } + + + } + + /* Otherwise, just fix-up the heap. */ + else + { + i = 0; + while (2 * i + 1 < nivals) + { + j = 2 * i + 1; + if (j + 1 < nivals + && ivals[heap[j + 1]].err >= ivals[heap[j]].err) + j++; + if (ivals[heap[j]].err <= ivals[heap[i]].err) + break; + else + { + t = heap[j]; + heap[j] = heap[i]; + heap[i] = t; + i = j; + } + } + + } + + /* If the heap is about to overflow, remove the last two + intervals. */ + while (nivals > cquad_heapsize - 2) + { + iv = &(ivals[heap[nivals - 1]]); +/* printf + ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n", + heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, + iv->depth); +*/ + err_final += iv->err; + igral_final += iv->igral; + nivals--; + } + + /* Collect the value of the integral and error. */ + igral = igral_final; + err = err_final; + for (i = 0; i < nivals; i++) + { + igral += ivals[heap[i]].igral; + err += ivals[heap[i]].err; + } + + } + + /* Dump the contents of the heap. */ +/* for (i = 0; i < nivals; i++) + { + iv = &(ivals[heap[i]]); + printf + ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n", + i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth, + iv->rdepth, iv->ndiv); + } +*/ + /* Clean up and present the results. */ + retval (0) = igral; + if (nargout > 1) + retval (1) = err; + if (nargout > 2) + retval (2) = neval; + /* All is well that ends well. */ + return retval; +}