GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
quadcc.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2010-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include <cmath>
31 
32 #include <algorithm>
33 
34 #include "lo-ieee.h"
35 #include "oct-locbuf.h"
36 
37 #include "defun.h"
38 #include "error.h"
39 #include "interpreter-private.h"
40 #include "ovl.h"
41 #include "parse.h"
42 #include "utils.h"
43 #include "variables.h"
44 
45 // Extended debugging.
46 #define DEBUG_QUADCC 0
47 
48 // Define the minimum size of the interval heap.
49 static const int MIN_CQUAD_HEAPSIZE = 200;
50 
51 // Data of a single interval.
52 typedef struct
53 {
54  double a, b;
55  double c[64];
56  double fx[33];
57  double igral, err;
58  int depth, rdepth, ndiv;
59 } cquad_ival;
60 
61 // Define relative tolerance used when deciding to drop an interval.
62 static const double DROP_RELTOL = std::numeric_limits<double>::epsilon () * 10;
63 
64 // Some constants and matrices that we'll need.
65 
66 static const double xi[33] =
67 {
68  -1., -0.99518472667219688624, -0.98078528040323044912,
69  -0.95694033573220886493, -0.92387953251128675612,
70  -0.88192126434835502970, -0.83146961230254523708,
71  -0.77301045336273696082, -0.70710678118654752440,
72  -0.63439328416364549822, -0.55557023301960222475,
73  -0.47139673682599764857, -0.38268343236508977173,
74  -0.29028467725446236764, -0.19509032201612826785,
75  -0.098017140329560601995, 0., 0.098017140329560601995,
76  0.19509032201612826785, 0.29028467725446236764, 0.38268343236508977173,
77  0.47139673682599764857, 0.55557023301960222475, 0.63439328416364549822,
78  0.70710678118654752440, 0.77301045336273696082, 0.83146961230254523708,
79  0.88192126434835502970, 0.92387953251128675612, 0.95694033573220886493,
80  0.98078528040323044912, 0.99518472667219688624, 1.
81 };
82 
83 static const double bee[68] =
84 {
85  0.00000000000000e+00, 2.28868854108532e-01, 0.00000000000000e+00,
86  -8.15740215243451e-01, 0.00000000000000e+00, 5.31212715259731e-01,
87  0.00000000000000e+00, 1.38538036812454e-02, 0.00000000000000e+00,
88  3.74405228908818e-02, 0.00000000000000e+00, 2.12224115039342e-01,
89  0.00000000000000e+00, -8.16362644507898e-01, 0.00000000000000e+00,
90  5.35648426691481e-01, 0.00000000000000e+00, 1.52417902753662e-03,
91  0.00000000000000e+00, 2.63058840550873e-03, 0.00000000000000e+00,
92  4.15292106318904e-03, 0.00000000000000e+00, 6.97106011119775e-03,
93  0.00000000000000e+00, 1.35535708431058e-02, 0.00000000000000e+00,
94  3.52132898424856e-02, 0.00000000000000e+00, 2.06946714741884e-01,
95  0.00000000000000e+00, -8.15674251283876e-01, 0.00000000000000e+00,
96  5.38841175520580e-01, 0.00000000000000e+00, 1.84909689577590e-04,
97  0.00000000000000e+00, 2.90936325007499e-04, 0.00000000000000e+00,
98  3.84877750950089e-04, 0.00000000000000e+00, 4.86436656735046e-04,
99  0.00000000000000e+00, 6.08688640346879e-04, 0.00000000000000e+00,
100  7.66732830740331e-04, 0.00000000000000e+00, 9.82753336104205e-04,
101  0.00000000000000e+00, 1.29359957505615e-03, 0.00000000000000e+00,
102  1.76616363801885e-03, 0.00000000000000e+00, 2.53323433039089e-03,
103  0.00000000000000e+00, 3.88872172121956e-03, 0.00000000000000e+00,
104  6.58635106468291e-03, 0.00000000000000e+00, 1.30326736343254e-02,
105  0.00000000000000e+00, 3.44353850696714e-02, 0.00000000000000e+00,
106  2.05025409531915e-01, 0.00000000000000e+00, -8.14985893995401e-01,
107  0.00000000000000e+00, 5.40679930965238e-01
108 };
109 
110 static const double Lalpha[33] =
111 {
112  5.77350269189626e-01, 5.16397779494322e-01, 5.07092552837110e-01,
113  5.03952630678970e-01, 5.02518907629606e-01, 5.01745206004255e-01,
114  5.01280411827603e-01, 5.00979432868120e-01, 5.00773395667191e-01,
115  5.00626174321759e-01, 5.00517330712619e-01, 5.00434593736979e-01,
116  5.00370233297676e-01, 5.00319182924304e-01, 5.00278009473803e-01,
117  5.00244319584578e-01, 5.00216403386025e-01, 5.00193012939056e-01,
118  5.00173220168024e-01, 5.00156323280355e-01, 5.00141783641018e-01,
119  5.00129182278347e-01, 5.00118189340972e-01, 5.00108542278496e-01,
120  5.00100030010004e-01, 5.00092481273333e-01, 5.00085755939229e-01,
121  5.00079738458365e-01, 5.00074332862969e-01, 5.00069458915387e-01,
122  5.00065049112355e-01, 5.00061046334395e-01, 5.00057401986298e-01
123 };
124 
125 static const double Lgamma[33] =
126 {
127  0.0, 0.0, 5.16397779494322e-01, 5.07092552837110e-01, 5.03952630678970e-01,
128  5.02518907629606e-01, 5.01745206004255e-01, 5.01280411827603e-01,
129  5.00979432868120e-01, 5.00773395667191e-01, 5.00626174321759e-01,
130  5.00517330712619e-01, 5.00434593736979e-01, 5.00370233297676e-01,
131  5.00319182924304e-01, 5.00278009473803e-01, 5.00244319584578e-01,
132  5.00216403386025e-01, 5.00193012939056e-01, 5.00173220168024e-01,
133  5.00156323280355e-01, 5.00141783641018e-01, 5.00129182278347e-01,
134  5.00118189340972e-01, 5.00108542278496e-01, 5.00100030010003e-01,
135  5.00092481273333e-01, 5.00085755939229e-01, 5.00079738458365e-01,
136  5.00074332862969e-01, 5.00069458915387e-01, 5.00065049112355e-01,
137  5.00061046334395e-01
138 };
139 
140 static const double V1inv[5 * 5] =
141 {
142  .47140452079103168293e-1, .37712361663282534635, .56568542494923801952,
143  .37712361663282534635, .47140452079103168293e-1,
144  -.81649658092772603273e-1, -.46188021535170061160, 0,
145  .46188021535170061160, .81649658092772603273e-1, .15058465048420853962,
146  .12046772038736683169, -.54210474174315074262, .12046772038736683169,
147  .15058465048420853962, -.21380899352993950775, .30237157840738178177, -0.,
148  -.30237157840738178177, .21380899352993950775, .10774960475223581324,
149  -.21549920950447162648, .21549920950447162648, -.21549920950447162648,
150  .10774960475223581324
151 };
152 
153 static const double V2inv[9 * 9] =
154 {
155  .11223917161691230546e-1, .10339219839658349826, .19754094204576565761,
156  .25577315077753587922, .27835314560994251755, .25577315077753587922,
157  .19754094204576565761, .10339219839658349826, .11223917161691230546e-1,
158  -.19440394783993476970e-1, -.16544884625069155470, -.24193725566041460608,
159  -.16953338808305493604, 0.0, .16953338808305493604, .24193725566041460608,
160  .16544884625069155470, .19440394783993476970e-1, .26466393115406349388e-1,
161  .17766815796285469394, .11316664642449611462, -.16306601003711325980,
162  -.30847037493128779631, -.16306601003711325980, .11316664642449611462,
163  .17766815796285469394, .26466393115406349388e-1,
164  -.32395302049990834508e-1, -.15521142532414866547,
165  .88573492664788602740e-1, .29570405784974857322, 0.0,
166  -.29570405784974857322, -.88573492664788602740e-1, .15521142532414866547,
167  .32395302049990834508e-1, .41442155673936851246e-1,
168  .98186757907405608245e-1, -.23056908429499411784,
169  -.68047008326360625520e-1, .31797435808002456774,
170  -.68047008326360625520e-1, -.23056908429499411784,
171  .98186757907405608245e-1, .41442155673936851246e-1,
172  -.49981120317798783134e-1, -.24861810572835756217e-1,
173  .23561326072010832539, -.24472785656448415351, 0.0, .24472785656448415351,
174  -.23561326072010832539, .24861810572835756217e-1,
175  .49981120317798783134e-1, .79691635865674781228e-1,
176  -.95725617891693941833e-1, -.57957553356854386344e-1,
177  .21164072460540271452, -.27529837844505833514, .21164072460540271452,
178  -.57957553356854386344e-1, -.95725617891693941833e-1,
179  .79691635865674781228e-1,
180  -.10894869830716590913, .20131094491947531782, -.15407672674888869038,
181  .83385723639789791384e-1, 0.0, -.83385723639789791384e-1,
182  .15407672674888869038, -.20131094491947531782, .10894869830716590913,
183  .54581057089643838221e-1, -.10916211417928767644, .10916211417928767644,
184  -.10916211417928767644, .10916211417928767644, -.10916211417928767644,
185  .10916211417928767644, -.10916211417928767644, .54581057089643838221e-1
186 };
187 
188 static const double V3inv[17 * 17] =
189 {
190  .27729677693590098996e-2, .26423663180333065153e-1,
191  .53374068493933898312e-1, .77007854739523195947e-1,
192  .98257061072911596869e-1, .11538049741786835604, .12832134344120884559,
193  .13612785914022865001, .13888293186236181317, .13612785914022865001,
194  .12832134344120884559, .11538049741786835604, .98257061072911596869e-1,
195  .77007854739523195947e-1, .53374068493933898312e-1,
196  .26423663180333065153e-1, .27729677693590098996e-2,
197  -.48029210642807413690e-2, -.44887724635478800254e-1,
198  -.85409520147301089416e-1, -.11090267822061423050, -.12033983162705862441,
199  -.11102786862182788886, -.85054870109799336515e-1,
200  -.45998467987742225160e-1, 0.0, .45998467987742225160e-1,
201  .85054870109799336515e-1, .11102786862182788886, .12033983162705862441,
202  .11090267822061423050, .85409520147301089416e-1, .44887724635478800254e-1,
203  .48029210642807413690e-2, .62758546879582030087e-2,
204  .55561297093529155869e-1,
205  .93281491021051539742e-1, .92320151237493695139e-1,
206  .55077987469605684531e-1,
207  -.96998141716497488255e-2, -.80285961895427405567e-1,
208  -.13496839655913850224,
209  -.15512521776684524331, -.13496839655913850224, -.80285961895427405567e-1,
210  -.96998141716497488255e-2, .55077987469605684531e-1,
211  .92320151237493695139e-1, .93281491021051539742e-1,
212  .55561297093529155869e-1, .62758546879582030087e-2,
213  -.74850969394858555939e-2, -.61751608943839234096e-1,
214  -.82974150437304275958e-1, -.38437763431942633378e-1,
215  .45745502025779701366e-1, .12369235652734542162, .14720439712852868239,
216  .98768034347019704401e-1, 0.0,
217  -.98768034347019704401e-1, -.14720439712852868239, -.12369235652734542162,
218  -.45745502025779701366e-1, .38437763431942633378e-1,
219  .82974150437304275958e-1, .61751608943839234096e-1,
220  .74850969394858555939e-2, .86710099994384056338e-2,
221  .64006230103659573344e-1, .58517426396091675690e-1,
222  -.29743410528985802680e-1,
223  -.11934127779157114754, -.12686773515361299409, -.30729137153877447035e-1,
224  .97307836256600731568e-1, .15635811574451401023, .97307836256600731568e-1,
225  -.30729137153877447035e-1, -.12686773515361299409, -.11934127779157114754,
226  -.29743410528985802680e-1, .58517426396091675690e-1,
227  .64006230103659573344e-1, .86710099994384056338e-2,
228  -.97486395666294840165e-2, -.62995604908060224672e-1,
229  -.24373234450275529219e-1, .87760984413626872730e-1,
230  .12205204576993351394,
231  .16216004196864002088e-1, -.12422320942156845775, -.13682714580929614678,
232  0.0, .13682714580929614678, .12422320942156845775,
233  -.16216004196864002088e-1, -.12205204576993351394,
234  -.87760984413626872730e-1, .24373234450275529219e-1,
235  .62995604908060224672e-1, .97486395666294840165e-2,
236  .10956271233750488468e-1, .58613204255294358939e-1,
237  -.13306063940736618859e-1, -.11606666444978454399,
238  -.52059598001115805639e-1, .10868540217796151849, .12594452879014618005,
239  -.44678658254872910434e-1, -.15617684362128533405,
240  -.44678658254872910434e-1, .12594452879014618005, .10868540217796151849,
241  -.52059598001115805639e-1, -.11606666444978454399,
242  -.13306063940736618859e-1, .58613204255294358939e-1,
243  .10956271233750488468e-1, -.12098893000863087230e-1,
244  -.51626244709126208453e-1, .48919433304746979330e-1,
245  .10467644465949427090,
246  -.48729879523084673782e-1, -.13668732103524749234, .28190838706814496438e-1,
247  .15434223333238741600, 0.0, -.15434223333238741600,
248  -.28190838706814496438e-1, .13668732103524749234,
249  .48729879523084673782e-1, -.10467644465949427090,
250  -.48919433304746979330e-1, .51626244709126208453e-1,
251  .12098893000863087230e-1, .13542668300437944822e-1,
252  .41712033418258689308e-1,
253  -.76190463272803434388e-1, -.58303943170068132010e-1, .12158068748245606853,
254  .42121099930651007882e-1, -.14684425840766337756,
255  -.16108203535058647043e-1, .15698075850757976092,
256  -.16108203535058647043e-1, -.14684425840766337756,
257  .42121099930651007882e-1, .12158068748245606853,
258  -.58303943170068132010e-1, -.76190463272803434388e-1,
259  .41712033418258689308e-1, .13542668300437944822e-1,
260  -.14939634995117694417e-1, -.30047246373341564039e-1,
261  .91624635082546425678e-1, -.79133374319110026377e-2,
262  -.12292558212072233355, .90013382617762643524e-1,
263  .84013717196539593395e-1, -.14813033309980695856, 0.0,
264  .14813033309980695856, -.84013717196539593395e-1,
265  -.90013382617762643524e-1,
266  .12292558212072233355, .79133374319110026377e-2, -.91624635082546425678e-1,
267  .30047246373341564039e-1, .14939634995117694417e-1,
268  .16986031342807474208e-1,
269  .15760203882617033601e-1, -.91494054040950941996e-1,
270  .70082459207876130806e-1,
271  .53390713710144539104e-1, -.14340746778352039430, .84048122493418898508e-1,
272  .72456667788091316868e-1, -.15564535320096811360,
273  .72456667788091316868e-1, .84048122493418898508e-1,
274  -.14340746778352039430, .53390713710144539104e-1,
275  .70082459207876130806e-1, -.91494054040950941996e-1,
276  .15760203882617033601e-1,
277  .16986031342807474208e-1, -.18994065631858742028e-1,
278  -.82901821370405592927e-3, .77239669773015192888e-1,
279  -.10850735431039424680, .47524484622086496464e-1,
280  .69148184871588737021e-1, -.14829314646228194928, .11992057742398672066,
281  0.0, -.11992057742398672066, .14829314646228194928,
282  -.69148184871588737021e-1, -.47524484622086496464e-1,
283  .10850735431039424680, -.77239669773015192888e-1,
284  .82901821370405592927e-3, .18994065631858742028e-1,
285  .22761703826371535132e-1, -.17728848711449643358e-1,
286  -.47496371572480503788e-1, .10659958402328690063, -.11696013966166296514,
287  .63073750910894244526e-1, .32928881123602721303e-1,
288  -.12280950532497593683, .15926189077282729505, -.12280950532497593683,
289  .32928881123602721303e-1, .63073750910894244526e-1,
290  -.11696013966166296514, .10659958402328690063, -.47496371572480503788e-1,
291  -.17728848711449643358e-1, .22761703826371535132e-1,
292  -.26493215276042203434e-1, .35579780856128386192e-1,
293  .10447309718398935122e-1, -.68616154085314996709e-1,
294  .11775363082763954214, -.13918901977011837274, .12312819418827395690,
295  -.72053565748259077905e-1, 0.0, .72053565748259077905e-1,
296  -.12312819418827395690, .13918901977011837274, -.11775363082763954214,
297  .68616154085314996709e-1, -.10447309718398935122e-1,
298  -.35579780856128386192e-1,
299  .26493215276042203434e-1, .40742523354399706918e-1,
300  -.73124912999529117195e-1, .49317266444153837821e-1,
301  -.13686605413876015320e-1, -.28342624942191100464e-1,
302  .70371855298258216249e-1, -.10600251632853603875, .12981016288391131812,
303  -.13817029659318161476, .12981016288391131812, -.10600251632853603875,
304  .70371855298258216249e-1, -.28342624942191100464e-1,
305  -.13686605413876015320e-1,
306  .49317266444153837821e-1, -.73124912999529117195e-1,
307  .40742523354399706918e-1, -.54944368958699908688e-1,
308  .10777725663147408190, -.10152395581538265428, .91369146312596428468e-1,
309  -.77703071757424700773e-1, .61050911730999815031e-1,
310  -.42052599404498348871e-1, .21438229266251454773e-1, 0.0,
311  -.21438229266251454773e-1, .42052599404498348871e-1,
312  -.61050911730999815031e-1, .77703071757424700773e-1,
313  -.91369146312596428468e-1,
314  .10152395581538265428, -.10777725663147408190, .54944368958699908688e-1,
315  .27485608464748840573e-1, -.54971216929497681146e-1,
316  .54971216929497681146e-1,
317  -.54971216929497681146e-1, .54971216929497681146e-1,
318  -.54971216929497681146e-1, .54971216929497681146e-1,
319  -.54971216929497681146e-1, .54971216929497681146e-1,
320  -.54971216929497681146e-1, .54971216929497681146e-1,
321  -.54971216929497681146e-1, .54971216929497681146e-1,
322  -.54971216929497681146e-1, .54971216929497681146e-1,
323  -.54971216929497681146e-1, .27485608464748840573e-1
324 };
325 
326 static const double V4inv[33 * 33] =
327 {
328  .69120897476690862600e-3, .66419939766331555194e-2,
329  .13600665164323186111e-1, .20122785860913684493e-1,
330  .26583214101668429944e-1, .32712713318999268739e-1,
331  .38576221976287138036e-1, .44033030938268925133e-1,
332  .49092709529622799673e-1, .53657949874312515646e-1,
333  .57724533144734311859e-1, .61219564530655179096e-1,
334  .64138907503837875026e-1, .66427905189318792009e-1,
335  .68088956652280022887e-1, .69083051391555695878e-1,
336  .69422738116739271449e-1, .69083051391555695878e-1,
337  .68088956652280022887e-1, .66427905189318792009e-1,
338  .64138907503837875026e-1, .61219564530655179096e-1,
339  .57724533144734311859e-1, .53657949874312515646e-1,
340  .49092709529622799673e-1, .44033030938268925133e-1,
341  .38576221976287138036e-1, .32712713318999268739e-1,
342  .26583214101668429944e-1, .20122785860913684493e-1,
343  .13600665164323186111e-1, .66419939766331555194e-2,
344  .69120897476690862600e-3, -.11972090629438798134e-2,
345  -.11448874821643225573e-1, -.23104401104002905904e-1,
346  -.33352899418646530133e-1, -.42538626424075425908e-1,
347  -.49969730733911825941e-1, -.55555454015360728353e-1,
348  -.58955533624852604918e-1, -.60126044219122513907e-1,
349  -.58959430451175833624e-1, -.55546925396227130606e-1,
350  -.49984739749347973762e-1, -.42513009141170294365e-1,
351  -.33399140950669746346e-1, -.23007690803851790829e-1,
352  -.11728275717520066169e-1, 0.0, .11728275717520066169e-1,
353  .23007690803851790829e-1, .33399140950669746346e-1,
354  .42513009141170294365e-1, .49984739749347973762e-1,
355  .55546925396227130606e-1, .58959430451175833624e-1,
356  .60126044219122513907e-1, .58955533624852604918e-1,
357  .55555454015360728353e-1, .49969730733911825941e-1,
358  .42538626424075425908e-1, .33352899418646530133e-1,
359  .23104401104002905904e-1, .11448874821643225573e-1,
360  .11972090629438798134e-2, .15501585012936019146e-2,
361  .14628781502199620482e-1, .28684915921474815271e-1,
362  .39299396074628048026e-1, .46393418975496284204e-1,
363  .48756902531094699526e-1, .46331333488337494692e-1,
364  .39012645376980228775e-1, .27452795421085791153e-1,
365  .12430953621169863781e-1, -.47682978056024928800e-2,
366  -.22825828045428973853e-1,
367  -.40195512090720278312e-1, -.55503004262826221955e-1,
368  -.67424537752827046308e-1, -.75020199300113606452e-1,
369  -.77607844312483656131e-1, -.75020199300113606452e-1,
370  -.67424537752827046308e-1, -.55503004262826221955e-1,
371  -.40195512090720278312e-1, -.22825828045428973853e-1,
372  -.47682978056024928800e-2, .12430953621169863781e-1,
373  .27452795421085791153e-1, .39012645376980228775e-1,
374  .46331333488337494692e-1, .48756902531094699526e-1,
375  .46393418975496284204e-1, .39299396074628048026e-1,
376  .28684915921474815271e-1, .14628781502199620482e-1,
377  .15501585012936019146e-2, -.18377757558949194214e-2,
378  -.17050470050949761565e-1, -.31952119564923250836e-1,
379  -.40197423449026348155e-1,
380  -.41205649520281371624e-1, -.33909965817492272248e-1,
381  -.19393664422115332144e-1, .56661049630886784692e-3,
382  .22948272173686561721e-1, .44489719570904738207e-1,
383  .61790363672287920596e-1, .72121014727028013894e-1,
384  .73627151185287858579e-1, .65784665375961398923e-1,
385  .49369676372333667559e-1, .26444326317059715065e-1, 0.0,
386  -.26444326317059715065e-1, -.49369676372333667559e-1,
387  -.65784665375961398923e-1, -.73627151185287858579e-1,
388  -.72121014727028013894e-1, -.61790363672287920596e-1,
389  -.44489719570904738207e-1, -.22948272173686561721e-1,
390  -.56661049630886784692e-3, .19393664422115332144e-1,
391  .33909965817492272248e-1, .41205649520281371624e-1,
392  .40197423449026348155e-1, .31952119564923250836e-1,
393  .17050470050949761565e-1, .18377757558949194214e-2,
394  .20942714740729767769e-2, .18935902405146518232e-1,
395  .33335840852491735126e-1, .36770680999102286065e-1,
396  .28873194534132768509e-1, .10267303017729535513e-1,
397  -.14607738306201572890e-1, -.40139568545572305818e-1,
398  -.59808326733858291561e-1, -.68528358823372627506e-1,
399  -.63306535387619244879e-1, -.44508601817574921056e-1,
400  -.15449116105605395357e-1, .17941083795006546367e-1,
401  .48747356011657242123e-1, .70329553984201665523e-1,
402  .78106117292526169663e-1, .70329553984201665523e-1,
403  .48747356011657242123e-1, .17941083795006546367e-1,
404  -.15449116105605395357e-1, -.44508601817574921056e-1,
405  -.63306535387619244879e-1, -.68528358823372627506e-1,
406  -.59808326733858291561e-1,
407  -.40139568545572305818e-1, -.14607738306201572890e-1,
408  .10267303017729535513e-1, .28873194534132768509e-1,
409  .36770680999102286065e-1, .33335840852491735126e-1,
410  .18935902405146518232e-1, .20942714740729767769e-2,
411  -.23245285491878278419e-2, -.20401404737639389919e-1,
412  -.33019548231022514097e-1, -.29709828426463720091e-1,
413  -.11760070922697422156e-1, .15987584743850393793e-1,
414  .43619012891472813485e-1, .61177322409671487721e-1,
415  .61144030218486655594e-1,
416  .41895377620089086167e-1, .80232011820644308033e-2,
417  -.30574701186675900915e-1,
418  -.62072243008844865848e-1, -.76336186183574765586e-1,
419  -.68435466095345537115e-1, -.40237669208466966207e-1, 0.0,
420  .40237669208466966207e-1, .68435466095345537115e-1,
421  .76336186183574765586e-1, .62072243008844865848e-1,
422  .30574701186675900915e-1, -.80232011820644308033e-2,
423  -.41895377620089086167e-1, -.61144030218486655594e-1,
424  -.61177322409671487721e-1, -.43619012891472813485e-1,
425  -.15987584743850393793e-1, .11760070922697422156e-1,
426  .29709828426463720091e-1, .33019548231022514097e-1,
427  .20401404737639389919e-1, .23245285491878278419e-2,
428  .25451717261579269307e-2, .21480418595666878775e-1,
429  .31177212469293007998e-1, .19816333607013379373e-1,
430  -.72439496274458793681e-2, -.38404203906598342397e-1,
431  -.57633632255322221046e-1, -.54070547403585392952e-1,
432  -.26249823354368866005e-1, .15643058212336881516e-1,
433  .54539832735118677194e-1, .73283028002473989724e-1,
434  .62835303524135936213e-1, .26175977027801048141e-1,
435  -.22193636309998606610e-1, -.62597049956093311234e-1,
436  -.78206986173170212505e-1, -.62597049956093311234e-1,
437  -.22193636309998606610e-1, .26175977027801048141e-1,
438  .62835303524135936213e-1,
439  .73283028002473989724e-1, .54539832735118677194e-1,
440  .15643058212336881516e-1,
441  -.26249823354368866005e-1, -.54070547403585392952e-1,
442  -.57633632255322221046e-1, -.38404203906598342397e-1,
443  -.72439496274458793681e-2, .19816333607013379373e-1,
444  .31177212469293007998e-1, .21480418595666878775e-1,
445  .25451717261579269307e-2, -.27506573922483820005e-2,
446  -.22224442095099251870e-1, -.27949927254215773020e-1,
447  -.80918481053370034987e-2, .25121859354449306916e-1,
448  .51563535009373061074e-1, .51936965107145960512e-1,
449  .22146626648171527753e-1,
450  -.24172689882103382748e-1, -.61731229104853568296e-1,
451  -.68477262429344201201e-1, -.38311232728303704742e-1,
452  .14160578713659552679e-1, .61248813427564184033e-1,
453  .77136328841293031805e-1, .52514801765183697988e-1, 0.0,
454  -.52514801765183697988e-1, -.77136328841293031805e-1,
455  -.61248813427564184033e-1, -.14160578713659552679e-1,
456  .38311232728303704742e-1,
457  .68477262429344201201e-1, .61731229104853568296e-1,
458  .24172689882103382748e-1,
459  -.22146626648171527753e-1, -.51936965107145960512e-1,
460  -.51563535009373061074e-1, -.25121859354449306916e-1,
461  .80918481053370034987e-2, .27949927254215773020e-1,
462  .22224442095099251870e-1, .27506573922483820005e-2,
463  .29562461131654311467e-2, .22630271480554450613e-1,
464  .23547399831373800971e-1, -.43964593440902476642e-2,
465  -.39055315767504970597e-1, -.52369643937940066804e-1,
466  -.28506131614971613422e-1, .19906048093338832322e-1,
467  .60408880866392420279e-1, .62493397473656883090e-1,
468  .21391278377641297859e-1, -.37302864786623254746e-1,
469  -.73665127933539496872e-1, -.61706142476854010202e-1,
470  -.78065168882546327888e-2, .52335307373945544428e-1,
471  .78278746279419264777e-1, .52335307373945544428e-1,
472  -.78065168882546327888e-2, -.61706142476854010202e-1,
473  -.73665127933539496872e-1, -.37302864786623254746e-1,
474  .21391278377641297859e-1, .62493397473656883090e-1,
475  .60408880866392420279e-1, .19906048093338832322e-1,
476  -.28506131614971613422e-1, -.52369643937940066804e-1,
477  -.39055315767504970597e-1, -.43964593440902476642e-2,
478  .23547399831373800971e-1, .22630271480554450613e-1,
479  .29562461131654311467e-2, -.31515718415504761303e-2,
480  -.22739451096655080673e-1, -.18157123602272119779e-1,
481  .16496480897167303621e-1, .46921166788569301124e-1,
482  .40644395739978416354e-1, -.46275803430732216900e-2,
483  -.52883375891308909486e-1, -.61116483226324111734e-1,
484  -.17411698764545629853e-1, .44773430013166822765e-1,
485  .73441577962383869198e-1, .42127368371995472815e-1,
486  -.25504645957196772465e-1, -.74126818045972742488e-1,
487  -.62780077864719287317e-1, 0.0, .62780077864719287317e-1,
488  .74126818045972742488e-1, .25504645957196772465e-1,
489  -.42127368371995472815e-1, -.73441577962383869198e-1,
490  -.44773430013166822765e-1, .17411698764545629853e-1,
491  .61116483226324111734e-1, .52883375891308909486e-1,
492  .46275803430732216900e-2, -.40644395739978416354e-1,
493  -.46921166788569301124e-1, -.16496480897167303621e-1,
494  .18157123602272119779e-1, .22739451096655080673e-1,
495  .31515718415504761303e-2, .33536559294882188208e-2,
496  .22535348942792006185e-1,
497  .12048629300953560767e-1, -.27166076791299493403e-1,
498  -.47492745604230978367e-1, -.19246623430993153174e-1,
499  .36231297307556299322e-1, .61713617181636122004e-1,
500  .25928029734266134490e-1, -.40478700752883602818e-1,
501  -.71053889866326412049e-1, -.31870824482961751482e-1,
502  .41515251100219081281e-1, .76481960760098381651e-1,
503  .36726509155999912440e-1, -.40090067032627055969e-1,
504  -.78270742903374539397e-1, -.40090067032627055969e-1,
505  .36726509155999912440e-1, .76481960760098381651e-1,
506  .41515251100219081281e-1, -.31870824482961751482e-1,
507  -.71053889866326412049e-1, -.40478700752883602818e-1,
508  .25928029734266134490e-1, .61713617181636122004e-1,
509  .36231297307556299322e-1, -.19246623430993153174e-1,
510  -.47492745604230978367e-1, -.27166076791299493403e-1,
511  .12048629300953560767e-1, .22535348942792006185e-1,
512  .33536559294882188208e-2,
513  -.35481220456925318865e-2, -.22062913693073191150e-1,
514  -.54487362861834144999e-2, .35438821865804087489e-1,
515  .40733077820527411302e-1, -.67403098138950720914e-2,
516  -.55559584405239171054e-1, -.42417050790865158745e-1,
517  .24499901971884704925e-1, .68721232891705409302e-1,
518  .34086082787461126592e-1, -.43441000373118474002e-1,
519  -.73878085292669148950e-1, -.18846995664706657127e-1,
520  .59827776178286834498e-1, .70644634584085901794e-1, 0.0,
521  -.70644634584085901794e-1, -.59827776178286834498e-1,
522  .18846995664706657127e-1, .73878085292669148950e-1,
523  .43441000373118474002e-1, -.34086082787461126592e-1,
524  -.68721232891705409302e-1, -.24499901971884704925e-1,
525  .42417050790865158745e-1, .55559584405239171054e-1,
526  .67403098138950720914e-2, -.40733077820527411302e-1,
527  -.35438821865804087489e-1, .54487362861834144999e-2,
528  .22062913693073191150e-1, .35481220456925318865e-2,
529  .37554176816665075631e-2, .21297045781589919482e-1,
530  -.13327293083183431816e-2,
531  -.40635299172764596484e-1, -.27659860508374175359e-1,
532  .31089232744083445986e-1, .56113781541334176109e-1,
533  .37577840643257763400e-2, -.60511227350664590865e-1,
534  -.46670556446129053853e-1, .33263195878575888247e-1,
535  .72757324720645228775e-1, .15011712351692283635e-1,
536  -.65601212994924119078e-1, -.60016855838843789772e-1,
537  .26220858553188665966e-1, .78322776605833552980e-1,
538  .26220858553188665966e-1, -.60016855838843789772e-1,
539  -.65601212994924119078e-1,
540  .15011712351692283635e-1, .72757324720645228775e-1,
541  .33263195878575888247e-1,
542  -.46670556446129053853e-1, -.60511227350664590865e-1,
543  .37577840643257763400e-2, .56113781541334176109e-1,
544  .31089232744083445986e-1, -.27659860508374175359e-1,
545  -.40635299172764596484e-1, -.13327293083183431816e-2,
546  .21297045781589919482e-1, .37554176816665075631e-2,
547  -.39566995305720591229e-2, -.20291873414438919995e-1,
548  .80617453830770930551e-2, .42270189157016547906e-1,
549  .10332624526759093004e-1, -.48054759547616142024e-1,
550  -.37678032941171643972e-1,
551  .36617192625732482394e-1, .61009425973424865714e-1,
552  -.95589113168026591466e-2,
553  -.71023202645076922361e-1, -.25097788086808784456e-1,
554  .62406621963267050244e-1, .56907293171100693511e-1,
555  -.36435383083882206257e-1, -.75790105119208756348e-1, 0.0,
556  .75790105119208756348e-1, .36435383083882206257e-1,
557  -.56907293171100693511e-1, -.62406621963267050244e-1,
558  .25097788086808784456e-1, .71023202645076922361e-1,
559  .95589113168026591466e-2,
560  -.61009425973424865714e-1, -.36617192625732482394e-1,
561  .37678032941171643972e-1, .48054759547616142024e-1,
562  -.10332624526759093004e-1, -.42270189157016547906e-1,
563  -.80617453830770930551e-2, .20291873414438919995e-1,
564  .39566995305720591229e-2, .41776092289182138591e-2,
565  .19013221163904414395e-1, -.14420609729849899876e-1,
566  -.40259160586844441220e-1, .86327811113710831649e-2,
567  .53564430703021034399e-1, .65469185402150431933e-2,
568  -.60383116311280629856e-1,
569  -.25657793784058876939e-1, .58745680576829226900e-1,
570  .45649937869034420296e-1,
571  -.49167932056844167772e-1, -.62696614328552187977e-1,
572  .32540234556426699997e-1, .74280410383464269758e-1,
573  -.11425672633410999870e-1, -.78280649404686404903e-1,
574  -.11425672633410999870e-1, .74280410383464269758e-1,
575  .32540234556426699997e-1, -.62696614328552187977e-1,
576  -.49167932056844167772e-1, .45649937869034420296e-1,
577  .58745680576829226900e-1, -.25657793784058876939e-1,
578  -.60383116311280629856e-1, .65469185402150431933e-2,
579  .53564430703021034399e-1,
580  .86327811113710831649e-2, -.40259160586844441220e-1,
581  -.14420609729849899876e-1, .19013221163904414395e-1,
582  .41776092289182138591e-2, -.43935502082478059199e-2,
583  -.17528761237509401631e-1, .20208915249153872535e-1,
584  .34734743119040669109e-1, -.26275910172353637955e-1,
585  -.46368003346018878786e-1,
586  .26800056330709381025e-1, .56681476464606609921e-1,
587  -.24749011438127255898e-1,
588  -.64934612189056658992e-1, .20333742247679279535e-1,
589  .71429299070059318651e-1,
590  -.14452513210428671266e-1, -.75793341281736586582e-1,
591  .74717094137184935270e-2, .78034921554757317374e-1, 0.0,
592  -.78034921554757317374e-1, -.74717094137184935270e-2,
593  .75793341281736586582e-1, .14452513210428671266e-1,
594  -.71429299070059318651e-1, -.20333742247679279535e-1,
595  .64934612189056658992e-1, .24749011438127255898e-1,
596  -.56681476464606609921e-1,
597  -.26800056330709381025e-1, .46368003346018878786e-1,
598  .26275910172353637955e-1,
599  -.34734743119040669109e-1, -.20208915249153872535e-1,
600  .17528761237509401631e-1, .43935502082478059199e-2,
601  .46379089482818671473e-2, .15791188144791287229e-1,
602  -.25134290048737455284e-1, -.26249795071946841205e-1,
603  .39960457575789924651e-1, .28111892450146525404e-1,
604  -.51026476400767918226e-1,
605  -.27266747278681831364e-1, .60708796647861610865e-1,
606  .23532306960642115854e-1,
607  -.68169639871532441111e-1, -.18204924701958312032e-1,
608  .73822890510656128485e-1, .11373392486424717019e-1,
609  -.77133324017644609416e-1, -.39295877480342619961e-2,
610  .78351902829418987960e-1, -.39295877480342619961e-2,
611  -.77133324017644609416e-1, .11373392486424717019e-1,
612  .73822890510656128485e-1, -.18204924701958312032e-1,
613  -.68169639871532441111e-1, .23532306960642115854e-1,
614  .60708796647861610865e-1, -.27266747278681831364e-1,
615  -.51026476400767918226e-1, .28111892450146525404e-1,
616  .39960457575789924651e-1, -.26249795071946841205e-1,
617  -.25134290048737455284e-1, .15791188144791287229e-1,
618  .46379089482818671473e-2, -.48780095920069827068e-2,
619  -.13886961667516983541e-1, .29071311049368895844e-1,
620  .15480559452075811600e-1, -.47527977686242313065e-1,
621  -.31929089844361042178e-2, .58015667638415922967e-1,
622  -.14547915466597622925e-1, -.61067668299848923244e-1,
623  .35093678009090186851e-1, .55378399159800654657e-1,
624  -.54277226474891610385e-1, -.42023830782434076509e-1,
625  .69197384645944912066e-1, .22610783557709586445e-1,
626  -.77269275900637030185e-1, 0.0, .77269275900637030185e-1,
627  -.22610783557709586445e-1,
628  -.69197384645944912066e-1, .42023830782434076509e-1,
629  .54277226474891610385e-1,
630  -.55378399159800654657e-1, -.35093678009090186851e-1,
631  .61067668299848923244e-1, .14547915466597622925e-1,
632  -.58015667638415922967e-1, .31929089844361042178e-2,
633  .47527977686242313065e-1, -.15480559452075811600e-1,
634  -.29071311049368895844e-1, .13886961667516983541e-1,
635  .48780095920069827068e-2, .51591759101720291381e-2,
636  .11747497650231330965e-1, -.31777863364694653331e-1,
637  -.34555825499804605557e-2, .47914131921157015198e-1,
638  -.22573685920142225247e-1, -.45320344390022666738e-1,
639  .49660630547172186418e-1, .25707858143963615736e-1,
640  -.68132707341917233933e-1, .67534860185243140399e-2,
641  .69268150370037450063e-1, -.41585011920451477177e-1,
642  -.51622397460510041271e-1, .68408139576363036148e-1,
643  .18981259024768933323e-1, -.78265472429342305554e-1,
644  .18981259024768933323e-1, .68408139576363036148e-1,
645  -.51622397460510041271e-1,
646  -.41585011920451477177e-1, .69268150370037450063e-1,
647  .67534860185243140399e-2,
648  -.68132707341917233933e-1, .25707858143963615736e-1,
649  .49660630547172186418e-1,
650  -.45320344390022666738e-1, -.22573685920142225247e-1,
651  .47914131921157015198e-1, -.34555825499804605557e-2,
652  -.31777863364694653331e-1, .11747497650231330965e-1,
653  .51591759101720291381e-2, -.54365757412741340377e-2,
654  -.94862516619529080191e-2, .33240472093448190877e-1,
655  -.88698898099681552229e-2,
656  -.40973252097216337576e-1, .42995673349795657065e-1,
657  .17320914507876958783e-1,
658  -.62201292691914856803e-1, .24726274174637346693e-1,
659  .51320859246515407288e-1,
660  -.62882063373810501763e-1, -.11003569131725622672e-1,
661  .73842261324108943465e-1, -.39240120294802923208e-1,
662  -.49293966443941122807e-1, .73552644778818223475e-1, 0.0,
663  -.73552644778818223475e-1, .49293966443941122807e-1,
664  .39240120294802923208e-1, -.73842261324108943465e-1,
665  .11003569131725622672e-1, .62882063373810501763e-1,
666  -.51320859246515407288e-1,
667  -.24726274174637346693e-1, .62201292691914856803e-1,
668  -.17320914507876958783e-1, -.42995673349795657065e-1,
669  .40973252097216337576e-1, .88698898099681552229e-2,
670  -.33240472093448190877e-1, .94862516619529080191e-2,
671  .54365757412741340377e-2, .57750194549356126240e-2,
672  .69981166020044116791e-2, -.33274982140403110792e-1,
673  .20297071020698356116e-1, .27898517839646066582e-1,
674  -.53368678853282030262e-1, .16656482990394548343e-1,
675  .46342901447260614255e-1,
676  -.60536796508149003365e-1, .29109107483842596340e-2,
677  .63224486124385124504e-1,
678  -.59028872851312033411e-1, -.14783105962696191734e-1,
679  .74269399241069253865e-1, -.49053677339382384625e-1,
680  -.33525466624811186739e-1, .78397349622515386647e-1,
681  -.33525466624811186739e-1, -.49053677339382384625e-1,
682  .74269399241069253865e-1, -.14783105962696191734e-1,
683  -.59028872851312033411e-1,
684  .63224486124385124504e-1, .29109107483842596340e-2,
685  -.60536796508149003365e-1,
686  .46342901447260614255e-1, .16656482990394548343e-1,
687  -.53368678853282030262e-1,
688  .27898517839646066582e-1, .20297071020698356116e-1,
689  -.33274982140403110792e-1,
690  .69981166020044116791e-2, .57750194549356126240e-2,
691  -.61100308370519200637e-2, -.44383614355738148616e-2,
692  .32011283412619094811e-1, -.29965011866372897633e-1,
693  -.10560682331349193348e-1, .51110336443392506342e-1,
694  -.45012284729681775492e-1, -.94236825555873320102e-2,
695  .60860695783141264746e-1,
696  -.55014628647083368926e-1, -.73474782382499482121e-2,
697  .66640148475243034781e-1, -.62533116045749887988e-1,
698  -.38650525912400102585e-2, .68429769005837003777e-1,
699  -.66984505412544901945e-1, 0.0, .66984505412544901945e-1,
700  -.68429769005837003777e-1, .38650525912400102585e-2,
701  .62533116045749887988e-1, -.66640148475243034781e-1,
702  .73474782382499482121e-2,
703  .55014628647083368926e-1, -.60860695783141264746e-1,
704  .94236825555873320102e-2,
705  .45012284729681775492e-1, -.51110336443392506342e-1,
706  .10560682331349193348e-1,
707  .29965011866372897633e-1, -.32011283412619094811e-1,
708  .44383614355738148616e-2,
709  .61100308370519200637e-2, .65409373892036191538e-2,
710  .16350101107071157065e-2, -.29301957285983144319e-1,
711  .36838667173388832579e-1, -.81922703976491586393e-2,
712  -.36955670021050133434e-1, .58374851095540469865e-1,
713  -.31977016246946181856e-1, -.25311073698658094646e-1,
714  .66674413950106952577e-1,
715  -.54865713324521039571e-1, -.39797027891537985440e-2,
716  .62830285264808449064e-1, -.72226313251296100676e-1,
717  .22560232697133353980e-1, .46455784709904033738e-1,
718  -.78200930751070349956e-1, .46455784709904033738e-1,
719  .22560232697133353980e-1, -.72226313251296100676e-1,
720  .62830285264808449064e-1, -.39797027891537985440e-2,
721  -.54865713324521039571e-1, .66674413950106952577e-1,
722  -.25311073698658094646e-1, -.31977016246946181856e-1,
723  .58374851095540469865e-1, -.36955670021050133434e-1,
724  -.81922703976491586393e-2, .36838667173388832579e-1,
725  -.29301957285983144319e-1, .16350101107071157065e-2,
726  .65409373892036191538e-2, -.69686180931868703196e-2,
727  .11849538727632789870e-2, .25452286414610537766e-1,
728  -.40522480651713943230e-1, .25694679053362813183e-1,
729  .14057118113748390637e-1, -.52037614725803488893e-1,
730  .58849342223684035589e-1,
731  -.25075229077361409271e-1, -.29559771094034181083e-1,
732  .68296746944165720199e-1, -.62890462146423984955e-1,
733  .14457636466274596445e-1, .45787612031322361496e-1,
734  -.77231759014655809742e-1, .57881203613910543657e-1, 0.0,
735  -.57881203613910543657e-1, .77231759014655809742e-1,
736  -.45787612031322361496e-1, -.14457636466274596445e-1,
737  .62890462146423984955e-1,
738  -.68296746944165720199e-1, .29559771094034181083e-1,
739  .25075229077361409271e-1,
740  -.58849342223684035589e-1, .52037614725803488893e-1,
741  -.14057118113748390637e-1, -.25694679053362813183e-1,
742  .40522480651713943230e-1, -.25452286414610537766e-1,
743  -.11849538727632789870e-2, .69686180931868703196e-2,
744  .75611653617520254845e-2, -.43290610418608409141e-2,
745  -.20277062025115566914e-1,
746  .40362947027704828926e-1, -.38938808024132120254e-1,
747  .11831186195916702262e-1,
748  .28476667401744525357e-1, -.59320969056617684621e-1,
749  .61101629747436200186e-1,
750  -.29514834848355389223e-1, -.20668001885001084821e-1,
751  .62923592802445122793e-1, -.73558456263588833115e-1,
752  .45314556330160999776e-1, .79031645918426015574e-2,
753  -.58136953576334689357e-1, .78538474524006405758e-1,
754  -.58136953576334689357e-1, .79031645918426015574e-2,
755  .45314556330160999776e-1, -.73558456263588833115e-1,
756  .62923592802445122793e-1, -.20668001885001084821e-1,
757  -.29514834848355389223e-1, .61101629747436200186e-1,
758  -.59320969056617684621e-1, .28476667401744525357e-1,
759  .11831186195916702262e-1, -.38938808024132120254e-1,
760  .40362947027704828926e-1, -.20277062025115566914e-1,
761  -.43290610418608409141e-2, .75611653617520254845e-2,
762  -.81505692478987769484e-2, .74297333588288568430e-2,
763  .14314212513540223314e-1, -.36711242251332751607e-1,
764  .46240027755503814626e-1, -.34921532671769023773e-1,
765  .46930051972353714773e-2,
766  .32842770336385381562e-1, -.61317813706529588466e-1,
767  .67000809902468893103e-1,
768  -.45337449655535622885e-1, .35794459576271920867e-2,
769  .41830061526027213385e-1,
770  -.72091371931944711708e-1, .74150028530317793195e-1,
771  -.46487632538609942002e-1, 0.0, .46487632538609942002e-1,
772  -.74150028530317793195e-1, .72091371931944711708e-1,
773  -.41830061526027213385e-1, -.35794459576271920867e-2,
774  .45337449655535622885e-1, -.67000809902468893103e-1,
775  .61317813706529588466e-1, -.32842770336385381562e-1,
776  -.46930051972353714773e-2, .34921532671769023773e-1,
777  -.46240027755503814626e-1, .36711242251332751607e-1,
778  -.14314212513540223314e-1, -.74297333588288568430e-2,
779  .81505692478987769484e-2, .90693182942442189743e-2,
780  -.11121000903959576737e-1, -.71308296141317458546e-2,
781  .29219439765986671645e-1, -.45820286629778129593e-1,
782  .49088381175879124421e-1, -.35614888785023038938e-1,
783  .78906970900092777895e-2,
784  .26262843038404929480e-1, -.56143674270125757857e-1,
785  .71700220472378350694e-1,
786  -.66963544500697307945e-1, .42215091779892228883e-1,
787  -.41338867413966866997e-2, -.36164891772995367321e-1,
788  .66584367783847858225e-1, -.77874712365070098328e-1,
789  .66584367783847858225e-1, -.36164891772995367321e-1,
790  -.41338867413966866997e-2, .42215091779892228883e-1,
791  -.66963544500697307945e-1,
792  .71700220472378350694e-1, -.56143674270125757857e-1,
793  .26262843038404929480e-1,
794  .78906970900092777895e-2, -.35614888785023038938e-1,
795  .49088381175879124421e-1,
796  -.45820286629778129593e-1, .29219439765986671645e-1,
797  -.71308296141317458546e-2, -.11121000903959576737e-1,
798  .90693182942442189743e-2, -.99848472706332791043e-2,
799  .14701271465939718856e-1, -.32917820356048383366e-3,
800  -.19201195309873585230e-1, .38409681836626963278e-1,
801  -.51647324405878909521e-1, .54522171113149311354e-1,
802  -.45040302741689006270e-1, .24183738595685990149e-1,
803  .42204134165479735097e-2, -.34317295181348742251e-1,
804  .59542472465494579941e-1, -.74135115907618101263e-1,
805  .74491937840566532596e-1, -.60042604725161994304e-1,
806  .33437677409000083169e-1, 0.0,
807  -.33437677409000083169e-1, .60042604725161994304e-1,
808  -.74491937840566532596e-1, .74135115907618101263e-1,
809  -.59542472465494579941e-1, .34317295181348742251e-1,
810  -.42204134165479735097e-2, -.24183738595685990149e-1,
811  .45040302741689006270e-1, -.54522171113149311354e-1,
812  .51647324405878909521e-1, -.38409681836626963278e-1,
813  .19201195309873585230e-1, .32917820356048383366e-3,
814  -.14701271465939718856e-1, .99848472706332791043e-2,
815  .11775579274769383373e-1, -.19892153937316935880e-1,
816  .95335114477449041055e-2, .57661528440359081617e-2,
817  -.23382690532380910781e-1, .40237257037170725321e-1,
818  -.53280289903551636474e-1, .59974361806023689068e-1,
819  -.58701684061992853224e-1, .49033407111597129616e-1,
820  -.31818835267847249219e-1, .90800541261162098886e-2,
821  .16272906819312603838e-1, -.40863896581186229487e-1,
822  .61346046297517367703e-1,
823  -.74896047554167268919e-1, .79632642148310325817e-1,
824  -.74896047554167268919e-1, .61346046297517367703e-1,
825  -.40863896581186229487e-1, .16272906819312603838e-1,
826  .90800541261162098886e-2, -.31818835267847249219e-1,
827  .49033407111597129616e-1, -.58701684061992853224e-1,
828  .59974361806023689068e-1, -.53280289903551636474e-1,
829  .40237257037170725321e-1, -.23382690532380910781e-1,
830  .57661528440359081617e-2, .95335114477449041055e-2,
831  -.19892153937316935880e-1,
832  .11775579274769383373e-1, -.13562702617218467450e-1,
833  .24885419969649845849e-1, -.18368693901908875583e-1,
834  .81673147806084084638e-2, .47890591326129587131e-2,
835  -.19313752945227974024e-1, .34065953398362954708e-1,
836  -.47667045133463415672e-1, .58820377816690514309e-1,
837  -.66424139824618415970e-1,
838  .69667606260856092515e-1, -.68102459384364543253e-1,
839  .61683024923302547971e-1,
840  -.50771943476441639136e-1, .36110771847327189215e-1,
841  -.18758028464284563358e-1, 0.0, .18758028464284563358e-1,
842  -.36110771847327189215e-1, .50771943476441639136e-1,
843  -.61683024923302547971e-1, .68102459384364543253e-1,
844  -.69667606260856092515e-1, .66424139824618415970e-1,
845  -.58820377816690514309e-1, .47667045133463415672e-1,
846  -.34065953398362954708e-1, .19313752945227974024e-1,
847  -.47890591326129587131e-2, -.81673147806084084638e-2,
848  .18368693901908875583e-1, -.24885419969649845849e-1,
849  .13562702617218467450e-1, .20576545037980523979e-1,
850  -.40093155172981004337e-1, .36954083167944054826e-1,
851  -.31856506837591907746e-1, .24996323181546255126e-1,
852  -.16637165210473614136e-1, .71002706773325085237e-2,
853  .32478629093205201133e-2,
854  -.14009562579050569518e-1, .24771262248780618922e-1,
855  -.35119395835433647559e-1, .44656290368574753171e-1,
856  -.53015448339647394161e-1, .59875631995693046782e-1,
857  -.64973208326045193862e-1, .68112280331082143373e-1,
858  -.69172215234062186994e-1, .68112280331082143373e-1,
859  -.64973208326045193862e-1, .59875631995693046782e-1,
860  -.53015448339647394161e-1, .44656290368574753171e-1,
861  -.35119395835433647559e-1, .24771262248780618922e-1,
862  -.14009562579050569518e-1, .32478629093205201133e-2,
863  .71002706773325085237e-2, -.16637165210473614136e-1,
864  .24996323181546255126e-1, -.31856506837591907746e-1,
865  .36954083167944054826e-1, -.40093155172981004337e-1,
866  .20576545037980523979e-1, -.27584914609096156163e-1,
867  .54904171411058497973e-1, -.54109756419563083153e-1,
868  .52794234894345577483e-1, -.50970276026831042415e-1,
869  .48655445537990983379e-1,
870  -.45872036510847994332e-1, .42646854695899611372e-1,
871  -.39010960357087507670e-1, .34999369144476467749e-1,
872  -.30650714874402762189e-1, .26006877464703437057e-1,
873  -.21112579608213651273e-1, .16014956068786763273e-1,
874  -.10763099747751940252e-1, .54075888924374485533e-2, 0.0,
875  -.54075888924374485533e-2, .10763099747751940252e-1,
876  -.16014956068786763273e-1,
877  .21112579608213651273e-1, -.26006877464703437057e-1,
878  .30650714874402762189e-1,
879  -.34999369144476467749e-1, .39010960357087507670e-1,
880  -.42646854695899611372e-1, .45872036510847994332e-1,
881  -.48655445537990983379e-1, .50970276026831042415e-1,
882  -.52794234894345577483e-1, .54109756419563083153e-1,
883  -.54904171411058497973e-1, .27584914609096156163e-1,
884  .13794141262469565740e-1, -.27588282524939131481e-1,
885  .27588282524939131481e-1, -.27588282524939131481e-1,
886  .27588282524939131481e-1, -.27588282524939131481e-1,
887  .27588282524939131481e-1,
888  -.27588282524939131481e-1, .27588282524939131481e-1,
889  -.27588282524939131481e-1, .27588282524939131481e-1,
890  -.27588282524939131481e-1, .27588282524939131481e-1,
891  -.27588282524939131481e-1, .27588282524939131481e-1,
892  -.27588282524939131481e-1, .27588282524939131481e-1,
893  -.27588282524939131481e-1, .27588282524939131481e-1,
894  -.27588282524939131481e-1, .27588282524939131481e-1,
895  -.27588282524939131481e-1, .27588282524939131481e-1,
896  -.27588282524939131481e-1, .27588282524939131481e-1,
897  -.27588282524939131481e-1, .27588282524939131481e-1,
898  -.27588282524939131481e-1, .27588282524939131481e-1,
899  -.27588282524939131481e-1, .27588282524939131481e-1,
900  -.27588282524939131481e-1, .13794141262469565740e-1
901 };
902 
903 static const double Tleft[33 * 33] =
904 {
905  1., -.86602540378443864678, 0., .33071891388307382381, 0.,
906  -.20728904939721249057, 0., .15128841196122722208, 0.,
907  -.11918864298744029244, 0., .98352013661686631224e-1, 0.,
908  -.83727065404940845733e-1, 0., .72893399403505841203e-1, 0.,
909  -.64544632643375022436e-1, 0., .57913170372415565639e-1, 0.,
910  -.52518242575729562263e-1, 0., .48043311993977520457e-1, 0.,
911  -.44271433659733990243e-1, 0., .41048928022856771981e-1, 0.,
912  -.38263878662008271459e-1, 0., .35832844026365304501e-1, 0., 0.,
913  .50000000000000000000, -.96824583655185422130, .57282196186948000082,
914  .21650635094610966169, -.35903516540862679125, -.97578093724974971969e-1,
915  .26203921611325660506, .55792409597991015609e-1, -.20644078533943456204,
916  -.36172381205961199479e-1, .17035068468874958194,
917  .25371838001497225980e-1, -.14501953125000000000,
918  -.18786835250972344757e-1, .12625507130328301066,
919  .14473795929590520582e-1, -.11179458309419422675,
920  -.11494434254897626155e-1, .10030855351241635862,
921  .93498556820544479096e-2, -.90964264465390582629e-1,
922  -.77546391824364392762e-2, .83213457337452292745e-1,
923  .65358085945588638605e-2, -.76680372422574234569e-1,
924  -.55835321940047427169e-2, .71098828931825789428e-1,
925  .48253327982967591019e-2, -.66274981937248958553e-1,
926  -.42118078245337801387e-2, .62064306433355646267e-1,
927  .37083386598903548973e-2, 0., 0., .25000000000000000000,
928  -.73950997288745200531, .83852549156242113615, -.23175620272173946716,
929  -.37791833195149451496, .25710129174850522325, .21608307321780204633,
930  -.22844049245646009157, -.14009503000335388415, .19897685605518413847,
931  .98264706042471226893e-1, -.17445445004279014046,
932  -.72761100054958328401e-1, .15463589893742108388,
933  .56056770591708784481e-1, -.13855313872640495158,
934  -.44517752443294564781e-1, .12534277657695128850,
935  .36211835346039665762e-1, -.11434398255136139683,
936  -.30033588409423828125e-1, .10506705408753910481,
937  .25313077840725783008e-1, -.97149327637744872155e-1,
938  -.21624927200393328444e-1, .90319582367202122625e-1,
939  .18688433567711780666e-1, -.84372291635345108584e-1,
940  -.16312261561845420752e-1, .79149526894804751586e-1,
941  .14362333871852474757e-1, 0., 0., 0., .12500000000000000000,
942  -.49607837082461073572, .82265291131801144317, -.59621200088559103072,
943  -.80054302859059362371e-1, .42612156697795759420,
944  -.90098145270865592887e-1, -.29769623255090078484, .13630307904779758221,
945  .21638835185708931831, -.14600247270306082052, -.16348801804014290453,
946  .14340708728599057249, .12755243353979286190, -.13661523715071346961,
947  -.10215585947881057394, .12864248070157166547, .83592528025348693602e-1,
948  -.12066728689302565222, -.69633728678718053052e-1, .11314245177331919532,
949  .58882939251410088028e-1, -.10621835858758221487,
950  -.50432266865187597572e-1, .99916834723527771581e-1,
951  .43672094283057258509e-1, -.94206380251950852413e-1,
952  -.38181356812697746418e-1, .89035739656537771225e-1,
953  .33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1,
954  -.31093357409581873586, .67604086414949799246, -.75644205980613611039,
955  .28990586430124175741, .30648508196770360914, -.35801372616842500052,
956  -.91326869828709014708e-1, .31127929687500000000,
957  -.90915752838698393094e-2, -.25637381283965534330,
958  .57601077850322797594e-1, .21019685709225757945,
959  -.81244992138514014256e-1, -.17375078516720988858,
960  .92289437277967051125e-1, .14527351914265391374,
961  -.96675340792832019889e-1, -.12289485697108543415,
962  .97448175340011084006e-1, .10511755943298339844,
963  -.96242247086378239657e-1, -.90822942272780513537e-1,
964  .93966350452322132384e-1, .79189411876493712558e-1,
965  -.91139307067989309325e-1, -.69613039934383197265e-1,
966  .88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0.,
967  .31250000000000000000e-1, -.18684782411095934408, .50176689760410660236,
968  -.74784031498626095398, .56472001151566251186, .14842464993721351203e-1,
969  -.41162920273003120936, .20243071230196532282, .23772054897172750436,
970  -.24963810923972235950, -.12116179938394678936, .24330535483519110663,
971  .47903849781124471359e-1, -.22133299683101224293,
972  -.20542915138527200983e-2, .19653465717678146728,
973  -.26818172626509178444e-1, -.17319122357631210944,
974  .45065391411065545445e-1, .15253391395444065941,
975  -.56543897711725408302e-1, -.13469154928743585367,
976  .63632471400208840155e-1, .11941684923913523817,
977  -.67828850207933293098e-1, -.10636309084510652670,
978  .70095786922999181504e-1, .95187373095150709082e-1, 0., 0., 0., 0., 0.,
979  0., .15625000000000000000e-1, -.10909562534194485289,
980  .34842348626527747318, -.64461114561628111443, .69382480527334683659,
981  -.29551102358528827763, -.25527584713978439819, .38878771718544715394,
982  -.82956185835347407489e-2, -.31183177761966943912, .12831420840372374767,
983  .22067618205599434368, -.17569196937129496961, -.14598057000132284135,
984  .18864406621763419484, .89921002550386645767e-1, -.18571835020187122114,
985  -.48967672227195481777e-1, .17584685670380332798,
986  .19267984545067426324e-1, -.16335437520503462738,
987  .22598055455032407594e-2, .15032800884170631129,
988  -.17883358353754640871e-1, -.13774837869432209951,
989  .29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0.,
990  0., .78125000000000000000e-2, -.62377810244809812496e-1,
991  .23080781467370883845, -.50841310636012325368, .69834547012574056043,
992  -.52572723156526459672, .11464215704954976471e-1, .38698869011491210342,
993  -.26125646622255207507, -.16951698812361607510, .29773875898928782269,
994  .20130501202570367491e-1, -.26332493149159310198,
995  .67734613690401207009e-1, .21207315477103762715, -.11541543390889415193,
996  -.16249634759782417533, .13885887405041735068, .11996491328010275427,
997  -.14810432001630926895, -.85177658352556243411e-1, .14918860659904380587,
998  .57317789510444151564e-1, -.14569827645586660151,
999  -.35213090145965327390e-1, .13975998126844578198, 0., 0., 0., 0., 0., 0.,
1000  0., 0., .39062500000000000000e-2, -.35101954600803571207e-1,
1001  .14761284084133737720, -.37655033076080192966, .62410290231517322776,
1002  -.64335622317683389875, .28188168266139524244, .22488495672137010675,
1003  -.39393811089283576186, .75184777995770096714e-1, .28472023119398293003,
1004  -.20410910833705899572, -.15590046962908511750, .23814567544617953125,
1005  .54442805556829031204e-1, -.22855930338589720954,
1006  .16303223615756629897e-1, .20172722433875559213,
1007  -.62723406421217419404e-1, -.17012230831020922010,
1008  .91754642766136561612e-1, .13927644821381121197, -.10886600968068418181,
1009  -.11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0.,
1010  0., 0., .19531250000000000000e-2, -.19506820659607596598e-1,
1011  .91865676095362231937e-1, -.26604607809696493849, .51425874205091288223,
1012  -.66047561132505329292, .48660109511591303851, -.17575661168678285615e-1,
1013  -.36594333408055703366, .29088854695378694533, .11318677346656537927,
1014  -.31110645235730182168, .60733219161008787341e-1, .24333848233620420826,
1015  -.15254312332655419708, -.15995968483455388613, .19010344455215289289,
1016  .86040636766440260000e-1, -.19652589954665259945,
1017  -.27633388517205837713e-1, .18660848552712880387,
1018  -.15942583868416775867e-1, -.16902042462382064786,
1019  .47278526495327740646e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1020  .97656250000000000000e-3, -.10731084460857378207e-1,
1021  .55939644713816406331e-1, -.18118487371914493668, .39914857299829864263,
1022  -.60812322949933902435, .60011887183061967583, -.26002695805835928795,
1023  -.20883922404786010096, .38988130966114638081, -.11797833550782589082,
1024  -.25231824756239520077, .24817859972953934712, .90516417677868996417e-1,
1025  -.26079073291293066798, .30259468817169480161e-1, .22178195264114178432,
1026  -.10569877864302048175, -.16679648389266977455, .14637718550245050850,
1027  .11219272032739559870, -.16359363640525750353, -.64358194509092101393e-1,
1028  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3,
1029  -.58542865274813470967e-2, .33461741635290096452e-1,
1030  -.11979993155896201271, .29580223766987206958, -.51874761979436016742,
1031  .62861483498014306968, -.44868895761051453296, .12567502628371529386e-1,
1032  .35040366183235474275, -.30466868455569500886, -.70903913601490112666e-1,
1033  .30822791893032512740, -.11969443264190207736, -.20764760317621313946,
1034  .20629838355452128532, .95269702915334718507e-1, -.22432624768705133300,
1035  -.33103381593477797101e-2, .20570036048155716333,
1036  -.62208282720094518964e-1, -.17095309330441436348, 0., 0., 0., 0., 0., 0.,
1037  0., 0., 0., 0., 0., 0., .24414062500000000000e-3,
1038  -.31714797501871532475e-2, .19721062526127334100e-1,
1039  -.77311181185536498246e-1, .21124871792841566575, -.41777980401893650886,
1040  .59401977834943551650, -.56132417807488349048, .23433675061367565951,
1041  .20222775295220942126, -.38280372496506190127, .14443804214023095767,
1042  .22268950939178466797, -.27211314150777981984, -.34184876506180717313e-1,
1043  .26006498895669734842, -.97650425186005090107e-1, -.19024527660129101293,
1044  .16789164198044635671, .10875811641651905252, -.19276785058805921298, 0.,
1045  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3,
1046  -.17078941137247586143e-2, .11477733754843910060e-1,
1047  -.48887017020924625462e-1, .14634927241421789683, -.32156282683019547854,
1048  .52165811920227223937, -.60001958466396926460, .41208501541480733755,
1049  -.11366945503190350975e-2, -.33968093962672089159, .30955190935923386766,
1050  .40657421856578262210e-1, -.29873400409871531764, .16094481791768257440,
1051  .16876122436206497694, -.23650217045022161255, -.33070260090574765012e-1,
1052  .22985258456375907796, -.68645651043827097771e-1, 0., 0., 0., 0., 0., 0.,
1053  0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4,
1054  -.91501857608428649078e-3, .66085179496951987952e-2,
1055  -.30383171695850355404e-1, .98840838845366876117e-1,
1056  -.23855447246420318989, .43322017468145613917, -.58049033744876107191,
1057  .52533893203742699346, -.20681056202371946180, -.20180000924562504384,
1058  .37503922291962681797, -.15988102869837429062, -.19823558102762374094,
1059  .28393023878803799622, -.11188133439357510403e-1, -.24730368377168229255,
1060  .14731529061377942839, .14878558042884266021, 0., 0., 0., 0., 0., 0., 0.,
1061  0., 0., 0., 0., 0., 0., 0., 0., .30517578125000000000e-4,
1062  -.48804277318479845551e-3, .37696080990601968396e-2,
1063  -.18603912108994738255e-1, .65325006755649582964e-1,
1064  -.17162960707938819795, .34411527956476971322, -.52289350347082497959,
1065  .57319653625674910592, -.37662253421045430413, -.14099055105384663902e-1,
1066  .33265570610216904208, -.30921265572647566661, -.19911390594166455281e-1,
1067  .28738590811031797718, -.18912130469738472647, -.13235936203215819193,
1068  .25076406142356675279, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1069  0., 0., 0., .15258789062500000000e-4, -.25928719280954633249e-3,
1070  .21327398937568540428e-2, -.11244626133630732010e-1,
1071  .42375605740664331966e-1, -.12031130345907846211, .26352562258934426830,
1072  -.44590628258512682078, .56682835613700749379, -.49116715128261660395,
1073  .17845943097110339078, .20541650677432497477, -.36739803642257458221,
1074  .16776034069210108273, .17920950989905112908, -.28867732805385066532,
1075  .46473465543376206337e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1076  0., 0., 0., 0., 0., .76293945312500000000e-5, -.13727610943181290891e-3,
1077  .11979683091449349286e-2, -.67195313034570709806e-2,
1078  .27044920779931968175e-1, -.82472196498517457862e-1,
1079  .19570475044896150093, -.36391620788543817693, .52241392782736588032,
1080  -.54727504974907879912, .34211551468813581183, .31580472732719957762e-1,
1081  -.32830006549176759667, .30563797665254420769, .64905014620683140120e-2,
1082  -.27642986248995073032, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1083  0., 0., 0., 0., 0., 0., .38146972656250000000e-5,
1084  -.72454147007837596854e-4, .66859847582761390285e-3,
1085  -.39751311980366118437e-2, .17015198650201528366e-1,
1086  -.55443621868993855715e-1, .14157060481641692131, -.28641242619559616836,
1087  .45610665490966615415, -.55262786406029265394, .45818352706035500108,
1088  -.14984403004611673047, -.21163807462970713245, .36007252928843413718,
1089  -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1090  0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5,
1091  -.38135049864067468562e-4, .37101393638555730015e-3,
1092  -.23305339886279723213e-2, .10569913448297127219e-1,
1093  -.36640175162216897547e-1, .10010476414320235508, -.21860074212675559892,
1094  .38124757096345313719, -.52020999209879669177, .52172632730659212045,
1095  -.30841620620308814614, -.50322546186721500184e-1, .32577618885114899053,
1096  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1097  0., 0., .95367431640625000000e-6, -.20021483206955925244e-4,
1098  .20481807322420625431e-3, -.13553476938058909882e-2,
1099  .64919676350791905019e-2, -.23848725425069251903e-1,
1100  .69384632678886421292e-1, -.16249711393618776934, .30736618106830314788,
1101  -.46399909601971539157, .53765031034002467225, -.42598991476520183929,
1102  .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1103  0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6,
1104  -.10487707828484902486e-4, .11254146162337528943e-3,
1105  -.78248929534271987118e-3, .39468337145306794566e-2,
1106  -.15313546659475671763e-1, .47249070825218564146e-1,
1107  -.11804374107101480543, .24031796927792491122, -.39629215049166341285,
1108  .51629108968402548545, -.49622372075429782915, 0., 0., 0., 0., 0., 0., 0.,
1109  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1110  .23841857910156250000e-6, -.54823314130625337326e-5,
1111  .61575377321535518154e-4, -.44877834366497538134e-3,
1112  .23774612048621955857e-2, -.97136347645161687796e-2,
1113  .31671599547606636717e-1, -.84028665767000747480e-1,
1114  .18298487576742964949, -.32647878537696945218, .46970971486488895077, 0.,
1115  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1116  0., 0., 0., 0., .11920928955078125000e-6, -.28604020001177375838e-5,
1117  .33559227978295551013e-4, -.25583821662860610560e-3,
1118  .14201552747787302339e-2, -.60938046986874414969e-2,
1119  .20930869247951926793e-1, -.58745021125678072911e-1,
1120  .13613725780285953720, -.26083988356030237586, 0., 0., 0., 0., 0., 0., 0.,
1121  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1122  .59604644775390625000e-7, -.14898180663526043291e-5,
1123  .18224991282807693921e-4, -.14504433444608833821e-3,
1124  .84184722720281809548e-3, -.37846965430000478789e-2,
1125  .13656355548211376864e-1, -.40409541997718853934e-1,
1126  .99226988101858325902e-1, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1127  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1128  .29802322387695312500e-7, -.77471708843445529468e-6,
1129  .98649879372606876995e-5, -.81814934772838523887e-4,
1130  .49554483992403011328e-3, -.23290922072351413938e-2,
1131  .88068134250844034186e-2, -.27393666952485719070e-1, 0., 0., 0., 0., 0.,
1132  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1133  0., 0., 0., .14901161193847656250e-7, -.40226235946098233685e-6,
1134  .53236418690561306700e-5, -.45933829691164002269e-4,
1135  .28982005232838857913e-3, -.14212974043211018374e-2,
1136  .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1137  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1138  .74505805969238281250e-8, -.20858299254133430408e-6,
1139  .28648457300134381744e-5, -.25677535898258910850e-4,
1140  .16849420429491355445e-3, -.86062824010315834002e-3, 0., 0., 0., 0., 0.,
1141  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1142  0., 0., 0., 0., 0., .37252902984619140625e-8, -.10801736017613096861e-6,
1143  .15376606719887104015e-5, -.14296523739727437959e-4,
1144  .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1145  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1146  .18626451492309570312e-8, -.55871592916438890146e-7,
1147  .82331193828137454068e-6, -.79302250528382787666e-5, 0., 0., 0., 0., 0.,
1148  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1149  0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9,
1150  -.28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0.,
1151  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1152  0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9,
1153  -.14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1154  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1155  0., 0., .23283064365386962891e-9
1156 };
1157 
1158 static const double Tright[33 * 33] =
1159 {
1160  1., .86602540378443864678, 0., -.33071891388307382381, 0.,
1161  .20728904939721249057, 0., -.15128841196122722208, 0.,
1162  .11918864298744029244, 0., -.98352013661686631224e-1, 0.,
1163  .83727065404940845733e-1, 0., -.72893399403505841203e-1, 0.,
1164  .64544632643375022436e-1, 0., -.57913170372415565639e-1, 0.,
1165  .52518242575729562263e-1, 0., -.48043311993977520457e-1, 0.,
1166  .44271433659733990243e-1, 0., -.41048928022856771981e-1, 0.,
1167  .38263878662008271459e-1, 0., -.35832844026365304501e-1, 0., 0.,
1168  .50000000000000000000, .96824583655185422130, .57282196186948000082,
1169  -.21650635094610966169, -.35903516540862679125, .97578093724974971969e-1,
1170  .26203921611325660506, -.55792409597991015609e-1, -.20644078533943456204,
1171  .36172381205961199479e-1, .17035068468874958194,
1172  -.25371838001497225980e-1, -.14501953125000000000,
1173  .18786835250972344757e-1, .12625507130328301066,
1174  -.14473795929590520582e-1, -.11179458309419422675,
1175  .11494434254897626155e-1, .10030855351241635862,
1176  -.93498556820544479096e-2, -.90964264465390582629e-1,
1177  .77546391824364392762e-2, .83213457337452292745e-1,
1178  -.65358085945588638605e-2, -.76680372422574234569e-1,
1179  .55835321940047427169e-2, .71098828931825789428e-1,
1180  -.48253327982967591019e-2, -.66274981937248958553e-1,
1181  .42118078245337801387e-2, .62064306433355646267e-1,
1182  -.37083386598903548973e-2, 0., 0., .25000000000000000000,
1183  .73950997288745200531, .83852549156242113615, .23175620272173946716,
1184  -.37791833195149451496, -.25710129174850522325, .21608307321780204633,
1185  .22844049245646009157, -.14009503000335388415, -.19897685605518413847,
1186  .98264706042471226893e-1, .17445445004279014046,
1187  -.72761100054958328401e-1, -.15463589893742108388,
1188  .56056770591708784481e-1, .13855313872640495158,
1189  -.44517752443294564781e-1, -.12534277657695128850,
1190  .36211835346039665762e-1, .11434398255136139683,
1191  -.30033588409423828125e-1, -.10506705408753910481,
1192  .25313077840725783008e-1, .97149327637744872155e-1,
1193  -.21624927200393328444e-1, -.90319582367202122625e-1,
1194  .18688433567711780666e-1, .84372291635345108584e-1,
1195  -.16312261561845420752e-1, -.79149526894804751586e-1,
1196  .14362333871852474757e-1, 0., 0., 0., .12500000000000000000,
1197  .49607837082461073572, .82265291131801144317, .59621200088559103072,
1198  -.80054302859059362371e-1, -.42612156697795759420,
1199  -.90098145270865592887e-1, .29769623255090078484, .13630307904779758221,
1200  -.21638835185708931831, -.14600247270306082052, .16348801804014290453,
1201  .14340708728599057249, -.12755243353979286190, -.13661523715071346961,
1202  .10215585947881057394, .12864248070157166547, -.83592528025348693602e-1,
1203  -.12066728689302565222, .69633728678718053052e-1, .11314245177331919532,
1204  -.58882939251410088028e-1, -.10621835858758221487,
1205  .50432266865187597572e-1, .99916834723527771581e-1,
1206  -.43672094283057258509e-1, -.94206380251950852413e-1,
1207  .38181356812697746418e-1, .89035739656537771225e-1,
1208  -.33661934598216332678e-1, 0., 0., 0., 0., .62500000000000000000e-1,
1209  .31093357409581873586, .67604086414949799246, .75644205980613611039,
1210  .28990586430124175741, -.30648508196770360914, -.35801372616842500052,
1211  .91326869828709014708e-1, .31127929687500000000, .90915752838698393094e-2,
1212  -.25637381283965534330, -.57601077850322797594e-1, .21019685709225757945,
1213  .81244992138514014256e-1, -.17375078516720988858,
1214  -.92289437277967051125e-1, .14527351914265391374,
1215  .96675340792832019889e-1, -.12289485697108543415,
1216  -.97448175340011084006e-1, .10511755943298339844,
1217  .96242247086378239657e-1, -.90822942272780513537e-1,
1218  -.93966350452322132384e-1, .79189411876493712558e-1,
1219  .91139307067989309325e-1, -.69613039934383197265e-1,
1220  -.88062491671135767870e-1, .61646331729340817494e-1, 0., 0., 0., 0., 0.,
1221  .31250000000000000000e-1, .18684782411095934408, .50176689760410660236,
1222  .74784031498626095398, .56472001151566251186, -.14842464993721351203e-1,
1223  -.41162920273003120936, -.20243071230196532282, .23772054897172750436,
1224  .24963810923972235950, -.12116179938394678936, -.24330535483519110663,
1225  .47903849781124471359e-1, .22133299683101224293,
1226  -.20542915138527200983e-2, -.19653465717678146728,
1227  -.26818172626509178444e-1, .17319122357631210944,
1228  .45065391411065545445e-1, -.15253391395444065941,
1229  -.56543897711725408302e-1, .13469154928743585367,
1230  .63632471400208840155e-1, -.11941684923913523817,
1231  -.67828850207933293098e-1, .10636309084510652670,
1232  .70095786922999181504e-1, -.95187373095150709082e-1, 0., 0., 0., 0., 0.,
1233  0., .15625000000000000000e-1, .10909562534194485289,
1234  .34842348626527747318, .64461114561628111443, .69382480527334683659,
1235  .29551102358528827763, -.25527584713978439819, -.38878771718544715394,
1236  -.82956185835347407489e-2, .31183177761966943912, .12831420840372374767,
1237  -.22067618205599434368, -.17569196937129496961, .14598057000132284135,
1238  .18864406621763419484, -.89921002550386645767e-1, -.18571835020187122114,
1239  .48967672227195481777e-1, .17584685670380332798,
1240  -.19267984545067426324e-1, -.16335437520503462738,
1241  -.22598055455032407594e-2, .15032800884170631129,
1242  .17883358353754640871e-1, -.13774837869432209951,
1243  -.29227555960587143675e-1, .12604194747513151053, 0., 0., 0., 0., 0., 0.,
1244  0., .78125000000000000000e-2, .62377810244809812496e-1,
1245  .23080781467370883845, .50841310636012325368, .69834547012574056043,
1246  .52572723156526459672, .11464215704954976471e-1, -.38698869011491210342,
1247  -.26125646622255207507, .16951698812361607510, .29773875898928782269,
1248  -.20130501202570367491e-1, -.26332493149159310198,
1249  -.67734613690401207009e-1, .21207315477103762715, .11541543390889415193,
1250  -.16249634759782417533, -.13885887405041735068, .11996491328010275427,
1251  .14810432001630926895, -.85177658352556243411e-1, -.14918860659904380587,
1252  .57317789510444151564e-1, .14569827645586660151,
1253  -.35213090145965327390e-1, -.13975998126844578198, 0., 0., 0., 0., 0., 0.,
1254  0., 0., .39062500000000000000e-2, .35101954600803571207e-1,
1255  .14761284084133737720, .37655033076080192966, .62410290231517322776,
1256  .64335622317683389875, .28188168266139524244, -.22488495672137010675,
1257  -.39393811089283576186, -.75184777995770096714e-1, .28472023119398293003,
1258  .20410910833705899572, -.15590046962908511750, -.23814567544617953125,
1259  .54442805556829031204e-1, .22855930338589720954, .16303223615756629897e-1,
1260  -.20172722433875559213, -.62723406421217419404e-1, .17012230831020922010,
1261  .91754642766136561612e-1, -.13927644821381121197, -.10886600968068418181,
1262  .11139075654373395292, .11797455976331702879, 0., 0., 0., 0., 0., 0., 0.,
1263  0., 0., .19531250000000000000e-2, .19506820659607596598e-1,
1264  .91865676095362231937e-1, .26604607809696493849, .51425874205091288223,
1265  .66047561132505329292, .48660109511591303851, .17575661168678285615e-1,
1266  -.36594333408055703366, -.29088854695378694533, .11318677346656537927,
1267  .31110645235730182168, .60733219161008787341e-1, -.24333848233620420826,
1268  -.15254312332655419708, .15995968483455388613, .19010344455215289289,
1269  -.86040636766440260000e-1, -.19652589954665259945,
1270  .27633388517205837713e-1, .18660848552712880387, .15942583868416775867e-1,
1271  -.16902042462382064786, -.47278526495327740646e-1, 0., 0., 0., 0., 0., 0.,
1272  0., 0., 0., 0., .97656250000000000000e-3, .10731084460857378207e-1,
1273  .55939644713816406331e-1, .18118487371914493668, .39914857299829864263,
1274  .60812322949933902435, .60011887183061967583, .26002695805835928795,
1275  -.20883922404786010096, -.38988130966114638081, -.11797833550782589082,
1276  .25231824756239520077, .24817859972953934712, -.90516417677868996417e-1,
1277  -.26079073291293066798, -.30259468817169480161e-1, .22178195264114178432,
1278  .10569877864302048175, -.16679648389266977455, -.14637718550245050850,
1279  .11219272032739559870, .16359363640525750353, -.64358194509092101393e-1,
1280  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .48828125000000000000e-3,
1281  .58542865274813470967e-2, .33461741635290096452e-1, .11979993155896201271,
1282  .29580223766987206958, .51874761979436016742, .62861483498014306968,
1283  .44868895761051453296, .12567502628371529386e-1, -.35040366183235474275,
1284  -.30466868455569500886, .70903913601490112666e-1, .30822791893032512740,
1285  .11969443264190207736, -.20764760317621313946, -.20629838355452128532,
1286  .95269702915334718507e-1, .22432624768705133300,
1287  -.33103381593477797101e-2, -.20570036048155716333,
1288  -.62208282720094518964e-1, .17095309330441436348, 0., 0., 0., 0., 0., 0.,
1289  0., 0., 0., 0., 0., 0., .24414062500000000000e-3,
1290  .31714797501871532475e-2, .19721062526127334100e-1,
1291  .77311181185536498246e-1, .21124871792841566575, .41777980401893650886,
1292  .59401977834943551650, .56132417807488349048, .23433675061367565951,
1293  -.20222775295220942126, -.38280372496506190127, -.14443804214023095767,
1294  .22268950939178466797, .27211314150777981984, -.34184876506180717313e-1,
1295  -.26006498895669734842, -.97650425186005090107e-1, .19024527660129101293,
1296  .16789164198044635671, -.10875811641651905252, -.19276785058805921298, 0.,
1297  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .12207031250000000000e-3,
1298  .17078941137247586143e-2, .11477733754843910060e-1,
1299  .48887017020924625462e-1, .14634927241421789683, .32156282683019547854,
1300  .52165811920227223937, .60001958466396926460, .41208501541480733755,
1301  .11366945503190350975e-2, -.33968093962672089159, -.30955190935923386766,
1302  .40657421856578262210e-1, .29873400409871531764, .16094481791768257440,
1303  -.16876122436206497694, -.23650217045022161255, .33070260090574765012e-1,
1304  .22985258456375907796, .68645651043827097771e-1, 0., 0., 0., 0., 0., 0.,
1305  0., 0., 0., 0., 0., 0., 0., 0., .61035156250000000000e-4,
1306  .91501857608428649078e-3, .66085179496951987952e-2,
1307  .30383171695850355404e-1, .98840838845366876117e-1, .23855447246420318989,
1308  .43322017468145613917, .58049033744876107191, .52533893203742699346,
1309  .20681056202371946180, -.20180000924562504384, -.37503922291962681797,
1310  -.15988102869837429062, .19823558102762374094, .28393023878803799622,
1311  .11188133439357510403e-1, -.24730368377168229255, -.14731529061377942839,
1312  .14878558042884266021, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1313  0., 0., .30517578125000000000e-4, .48804277318479845551e-3,
1314  .37696080990601968396e-2, .18603912108994738255e-1,
1315  .65325006755649582964e-1, .17162960707938819795, .34411527956476971322,
1316  .52289350347082497959, .57319653625674910592, .37662253421045430413,
1317  -.14099055105384663902e-1, -.33265570610216904208, -.30921265572647566661,
1318  .19911390594166455281e-1, .28738590811031797718, .18912130469738472647,
1319  -.13235936203215819193, -.25076406142356675279, 0., 0., 0., 0., 0., 0.,
1320  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .15258789062500000000e-4,
1321  .25928719280954633249e-3, .21327398937568540428e-2,
1322  .11244626133630732010e-1, .42375605740664331966e-1, .12031130345907846211,
1323  .26352562258934426830, .44590628258512682078, .56682835613700749379,
1324  .49116715128261660395, .17845943097110339078, -.20541650677432497477,
1325  -.36739803642257458221, -.16776034069210108273, .17920950989905112908,
1326  .28867732805385066532, .46473465543376206337e-1, 0., 0., 0., 0., 0., 0.,
1327  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .76293945312500000000e-5,
1328  .13727610943181290891e-3, .11979683091449349286e-2,
1329  .67195313034570709806e-2, .27044920779931968175e-1,
1330  .82472196498517457862e-1, .19570475044896150093, .36391620788543817693,
1331  .52241392782736588032, .54727504974907879912, .34211551468813581183,
1332  -.31580472732719957762e-1, -.32830006549176759667, -.30563797665254420769,
1333  .64905014620683140120e-2, .27642986248995073032, 0., 0., 0., 0., 0., 0.,
1334  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .38146972656250000000e-5,
1335  .72454147007837596854e-4, .66859847582761390285e-3,
1336  .39751311980366118437e-2, .17015198650201528366e-1,
1337  .55443621868993855715e-1, .14157060481641692131, .28641242619559616836,
1338  .45610665490966615415, .55262786406029265394, .45818352706035500108,
1339  .14984403004611673047, -.21163807462970713245, -.36007252928843413718,
1340  -.17030961385712954159, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1341  0., 0., 0., 0., 0., 0., 0., .19073486328125000000e-5,
1342  .38135049864067468562e-4, .37101393638555730015e-3,
1343  .23305339886279723213e-2, .10569913448297127219e-1,
1344  .36640175162216897547e-1, .10010476414320235508, .21860074212675559892,
1345  .38124757096345313719, .52020999209879669177, .52172632730659212045,
1346  .30841620620308814614, -.50322546186721500184e-1, -.32577618885114899053,
1347  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1348  0., 0., .95367431640625000000e-6, .20021483206955925244e-4,
1349  .20481807322420625431e-3, .13553476938058909882e-2,
1350  .64919676350791905019e-2, .23848725425069251903e-1,
1351  .69384632678886421292e-1, .16249711393618776934, .30736618106830314788,
1352  .46399909601971539157, .53765031034002467225, .42598991476520183929,
1353  .12130445348350215652, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1354  0., 0., 0., 0., 0., 0., 0., 0., .47683715820312500000e-6,
1355  .10487707828484902486e-4, .11254146162337528943e-3,
1356  .78248929534271987118e-3, .39468337145306794566e-2,
1357  .15313546659475671763e-1, .47249070825218564146e-1, .11804374107101480543,
1358  .24031796927792491122, .39629215049166341285, .51629108968402548545,
1359  .49622372075429782915, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1360  0., 0., 0., 0., 0., 0., 0., 0., 0., .23841857910156250000e-6,
1361  .54823314130625337326e-5, .61575377321535518154e-4,
1362  .44877834366497538134e-3, .23774612048621955857e-2,
1363  .97136347645161687796e-2, .31671599547606636717e-1,
1364  .84028665767000747480e-1, .18298487576742964949, .32647878537696945218,
1365  .46970971486488895077, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1366  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., .11920928955078125000e-6,
1367  .28604020001177375838e-5, .33559227978295551013e-4,
1368  .25583821662860610560e-3, .14201552747787302339e-2,
1369  .60938046986874414969e-2, .20930869247951926793e-1,
1370  .58745021125678072911e-1, .13613725780285953720, .26083988356030237586,
1371  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1372  0., 0., 0., 0., 0., 0., .59604644775390625000e-7,
1373  .14898180663526043291e-5, .18224991282807693921e-4,
1374  .14504433444608833821e-3, .84184722720281809548e-3,
1375  .37846965430000478789e-2, .13656355548211376864e-1,
1376  .40409541997718853934e-1, .99226988101858325902e-1, 0., 0., 0., 0., 0.,
1377  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1378  0., 0., .29802322387695312500e-7, .77471708843445529468e-6,
1379  .98649879372606876995e-5, .81814934772838523887e-4,
1380  .49554483992403011328e-3, .23290922072351413938e-2,
1381  .88068134250844034186e-2, .27393666952485719070e-1, 0., 0., 0., 0., 0.,
1382  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1383  0., 0., 0., .14901161193847656250e-7, .40226235946098233685e-6,
1384  .53236418690561306700e-5, .45933829691164002269e-4,
1385  .28982005232838857913e-3, .14212974043211018374e-2,
1386  .56192363087488842264e-2, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1387  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1388  .74505805969238281250e-8, .20858299254133430408e-6,
1389  .28648457300134381744e-5, .25677535898258910850e-4,
1390  .16849420429491355445e-3, .86062824010315834002e-3, 0., 0., 0., 0., 0.,
1391  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1392  0., 0., 0., 0., 0., .37252902984619140625e-8, .10801736017613096861e-6,
1393  .15376606719887104015e-5, .14296523739727437959e-4,
1394  .97419023656050887203e-4, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1395  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1396  .18626451492309570312e-8, .55871592916438890146e-7,
1397  .82331193828137454068e-6, .79302250528382787666e-5, 0., 0., 0., 0., 0.,
1398  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1399  0., 0., 0., 0., 0., 0., 0., .93132257461547851562e-9,
1400  .28867244235852488244e-7, .43982811713864556957e-6, 0., 0., 0., 0., 0.,
1401  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1402  0., 0., 0., 0., 0., 0., 0., 0., .46566128730773925781e-9,
1403  .14899342093408253335e-7, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1404  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
1405  0., 0., .23283064365386962891e-9
1406 };
1407 
1408 // Allocates a workspace for the given maximum number of intervals.
1409 // Note that if the workspace gets filled, the intervals with the lowest
1410 // error estimates are dropped. The maximum number of intervals is
1411 // therefore not the maximum number of intervals that will be computed,
1412 // but merely the size of the buffer.
1413 
1414 // Compute the product of fx with one of the inverse Vandermonde-like matrices.
1415 
1416 static void
1417 Vinvfx (const double *fx, double *c, const int d)
1418 {
1419  int i, j;
1420 
1421  switch (d)
1422  {
1423  case 0:
1424  for (i = 0; i <= 4; i++)
1425  {
1426  c[i] = 0.0;
1427  for (j = 0; j <= 4; j++)
1428  c[i] += V1inv[i*5 + j] * fx[j * 8];
1429  }
1430  break;
1431  case 1:
1432  for (i = 0; i <= 8; i++)
1433  {
1434  c[i] = 0.0;
1435  for (j = 0; j <= 8; j++)
1436  c[i] += V2inv[i*9 + j] * fx[j * 4];
1437  }
1438  break;
1439  case 2:
1440  for (i = 0; i <= 16; i++)
1441  {
1442  c[i] = 0.0;
1443  for (j = 0; j <= 16; j++)
1444  c[i] += V3inv[i*17 + j] * fx[j * 2];
1445  }
1446  break;
1447  case 3:
1448  for (i = 0; i <= 32; i++)
1449  {
1450  c[i] = 0.0;
1451  for (j = 0; j <= 32; j++)
1452  c[i] += V4inv[i*33 + j] * fx[j];
1453  }
1454  break;
1455  }
1456 }
1457 
1458 // Downdate the interpolation given by the N coefficients C by removing
1459 // the nodes with indices in NaNs.
1460 
1461 static void
1462 downdate (double *c, int n, int d, int *nans, int nnans)
1463 {
1464  static const int bidx[4] = { 0, 6, 16, 34 };
1465  double b_new[34], alpha;
1466  int i, j;
1467 
1468  for (i = 0; i <= n + 1; i++)
1469  b_new[i] = bee[bidx[d] + i];
1470  for (i = 0; i < nnans; i++)
1471  {
1472  b_new[n + 1] = b_new[n + 1] / Lalpha[n];
1473  b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1];
1474  for (j = n - 1; j > 0; j--)
1475  b_new[j] = (b_new[j] + xi[nans[i]] * b_new[j + 1]
1476  - Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1];
1477  for (j = 0; j <= n; j++)
1478  b_new[j] = b_new[j + 1];
1479  alpha = c[n] / b_new[n];
1480  for (j = 0; j < n; j++)
1481  c[j] -= alpha * b_new[j];
1482  c[n] = 0;
1483  n--;
1484  }
1485 }
1486 
1487 // The actual integration routine.
1488 
1489 DEFMETHOD (quadcc, interp, args, ,
1490  doc: /* -*- texinfo -*-
1491 @deftypefn {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b})
1492 @deftypefnx {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol})
1493 @deftypefnx {} {@var{q} =} quadcc (@var{f}, @var{a}, @var{b}, @var{tol}, @var{sing})
1494 @deftypefnx {} {[@var{q}, @var{err}, @var{nr_points}] =} quadcc (@dots{})
1495 Numerically evaluate the integral of @var{f} from @var{a} to @var{b} using
1496 doubly-adaptive @nospell{Clenshaw-Curtis} quadrature.
1497 
1498 @var{f} is a function handle, inline function, or string containing the name
1499 of the function to evaluate. The function @var{f} must be vectorized and
1500 must return a vector of output values if given a vector of input values.
1501 For example,
1502 
1503 @example
1504 f = @@(x) x .* sin (1./x) .* sqrt (abs (1 - x));
1505 @end example
1506 
1507 @noindent
1508 which uses the element-by-element ``dot'' form for all operators.
1509 
1510 @var{a} and @var{b} are the lower and upper limits of integration. Either or
1511 both limits may be infinite. @code{quadcc} handles an infinite limit by
1512 substituting the variable of integration with @code{x = tan (pi/2*u)}.
1513 
1514 The optional argument @var{tol} is a 1- or 2-element vector that specifies the
1515 desired accuracy of the result. The first element of the vector is the desired
1516 absolute tolerance, and the second element is the desired relative tolerance.
1517 To choose a relative test only, set the absolute tolerance to zero. To choose
1518 an absolute test only, set the relative tolerance to zero. The default
1519 absolute tolerance is 1e-10 (1e-5 for single), and the default relative
1520 tolerance is 1e-6 (1e-4 for single).
1521 
1522 The optional argument @var{sing} contains a list of points where the integrand
1523 has known singularities, or discontinuities in any of its derivatives, inside
1524 the integration interval. For the example above, which has a discontinuity at
1525 x=1, the call to @code{quadcc} would be as follows
1526 
1527 @example
1528 int = quadcc (f, a, b, [], [ 1 ]);
1529 @end example
1530 
1531 The result of the integration is returned in @var{q}.
1532 
1533 @var{err} is an estimate of the absolute integration error.
1534 
1535 @var{nr_points} is the number of points at which the integrand was evaluated.
1536 
1537 If the adaptive integration did not converge, the value of @var{err} will be
1538 larger than the requested tolerance. Therefore, it is recommended to verify
1539 this value for difficult integrands.
1540 
1541 @code{quadcc} is capable of dealing with non-numeric values of the integrand
1542 such as @code{NaN} or @code{Inf}. If the integral diverges, and @code{quadcc}
1543 detects this, then a warning is issued and @code{Inf} or @code{-Inf} is
1544 returned.
1545 
1546 Note: @code{quadcc} is a general purpose quadrature algorithm and, as such,
1547 may be less efficient for a smooth or otherwise well-behaved integrand than
1548 other methods such as @code{quadgk}.
1549 
1550 The algorithm uses @nospell{Clenshaw-Curtis} quadrature rules of increasing
1551 degree in each interval and bisects the interval if either the function does
1552 not appear to be smooth or a rule of maximum degree has been reached. The
1553 error estimate is computed from the L2-norm of the difference between two
1554 successive interpolations of the integrand over the nodes of the respective
1555 quadrature rules.
1556 
1557 Reference: @nospell{P. Gonnet}, @cite{Increasing the Reliability of Adaptive
1558 Quadrature Using Explicit Interpolants}, @nospell{ACM} Transactions on
1559 Mathematical Software, Vol.@: 37, Issue 3, Article No.@: 3, 2010.
1560 @seealso{quad, quadv, quadl, quadgk, trapz, dblquad, triplequad}
1561 @end deftypefn */)
1562 {
1563  // Some constants that we will need.
1564  static const int n[4] = { 4, 8, 16, 32 };
1565  static const int skip[4] = { 8, 4, 2, 1 };
1566  static const int idx[4] = { 0, 5, 14, 31 };
1567  static const double w = M_SQRT2 / 2;
1568  static const int ndiv_max = 20;
1569 
1570  // Arguments left and right.
1571  int nargin = args.length ();
1572  octave_value fcn;
1573  double a, b, abstol, reltol, *sing;
1574  bool issingle;
1575 
1576  // Variables needed for transforming the integrand.
1577  bool wrap = false;
1578  double xw;
1579 
1580  // Stuff we will need to call the integrand.
1581  octave_value_list fargs, fvals;
1582 
1583  // Actual variables (as opposed to constants above).
1584  double m, h, ml, hl, mr, hr, temp;
1585  double igral, err, igral_final, err_final;
1586  int nivals;
1587  int neval = 0;
1588  int i, j, d, split;
1589  int nnans, nans[33];
1590  cquad_ival *iv, *ivl, *ivr;
1591  double nc, ncdiff;
1592 
1593  // Parse the input arguments.
1594  if (nargin < 3)
1595  print_usage ();
1596 
1597  fcn = octave::get_function_handle (interp, args(0), "x");
1598 
1599  if (! args(1).is_real_scalar ())
1600  error ("quadcc: lower limit of integration (A) must be a real scalar");
1601  a = args(1).double_value ();
1602  issingle = args(1).is_single_type ();
1603 
1604  if (! args(2).is_real_scalar ())
1605  error ("quadcc: upper limit of integration (B) must be a real scalar");
1606  b = args(2).double_value ();
1607  issingle = (issingle || args(2).is_single_type ());
1608 
1609  if (nargin < 4 || args(3).isempty ())
1610  {
1611  if (issingle)
1612  {
1613  abstol = 1.0e-5;
1614  reltol = 1.0e-4;
1615  }
1616  else
1617  {
1618  abstol = 1.0e-10;
1619  reltol = 1.0e-6;
1620  }
1621  }
1622  else if (! args(3).isnumeric () || args(3).numel () > 2)
1623  error ("quadcc: TOL must be a 1- or 2-element vector [AbsTol, RelTol]");
1624  else
1625  {
1626  NDArray tol = args(3).array_value ();
1627 
1628  abstol = tol(0);
1629  if (abstol < 0)
1630  error ("quadcc: absolute tolerance must be >=0");
1631 
1632  if (tol.numel () == 1)
1633  {
1634  if (issingle)
1635  reltol = 1.0e-4;
1636  else
1637  reltol = 1.0e-6;
1638  }
1639  else
1640  {
1641  reltol = tol(1);
1642  if (reltol < 0)
1643  error ("quadcc: relative tolerance must be >=0");
1644  }
1645  }
1646 
1647  if (nargin < 5)
1648  nivals = 1;
1649  else if (! (args(4).is_real_scalar () || args(4).is_real_matrix ()))
1650  error ("quadcc: list of singularities (SING) must be a vector of real values");
1651  else
1652  nivals = 1 + args(4).numel ();
1653 
1654  int cquad_heapsize = (nivals >= MIN_CQUAD_HEAPSIZE ? nivals + 1
1655  : MIN_CQUAD_HEAPSIZE);
1656  // The interval heap.
1657  OCTAVE_LOCAL_BUFFER (cquad_ival, ivals, cquad_heapsize);
1658  OCTAVE_LOCAL_BUFFER (double, iivals, cquad_heapsize);
1659  OCTAVE_LOCAL_BUFFER (int, heap, cquad_heapsize);
1660 
1661  if (nivals == 1)
1662  {
1663  iivals[0] = a;
1664  iivals[1] = b;
1665  }
1666  else
1667  {
1668  // Intervals around singularities.
1669  sing = args(4).array_value ().fortran_vec ();
1670  iivals[0] = a;
1671  std::copy_n (sing, nivals-1, iivals+1);
1672  iivals[nivals] = b;
1673  }
1674 
1675  // If a or b are +/-Inf, transform the integral.
1677  {
1678  wrap = true;
1679  for (i = 0; i < nivals + 1; i++)
1680  if (octave::math::isinf (iivals[i]))
1681  iivals[i] = std::copysign (1.0, iivals[i]);
1682  else
1683  iivals[i] = 2.0 * atan (iivals[i]) / M_PI;
1684  }
1685 
1686  // Initialize the heaps.
1687  for (i = 0; i < cquad_heapsize; i++)
1688  heap[i] = i;
1689 
1690  // Create the first interval(s).
1691  igral = 0.0;
1692  err = 0.0;
1693  for (j = 0; j < nivals; j++)
1694  {
1695  // Initialize the interval.
1696  iv = &(ivals[heap[j]]);
1697  m = (iivals[j] + iivals[j + 1]) / 2;
1698  h = (iivals[j + 1] - iivals[j]) / 2;
1699  nnans = 0;
1700  ColumnVector ex (33);
1701  if (wrap)
1702  {
1703  for (i = 0; i <= n[3]; i++)
1704  ex(i) = tan (M_PI/2 * (m + xi[i]*h));
1705  }
1706  else
1707  {
1708  for (i = 0; i <= n[3]; i++)
1709  ex(i) = m + xi[i]*h;
1710  }
1711  fargs(0) = ex;
1712  fvals = octave::feval (fcn, fargs, 1);
1713  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1714  error ("quadcc: integrand F must return a single, real-valued vector");
1715 
1716  Matrix effex = fvals(0).matrix_value ();
1717  if (effex.numel () != ex.numel ())
1718  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1719 
1720  for (i = 0; i <= n[3]; i++)
1721  {
1722  iv->fx[i] = effex(i);
1723  if (wrap)
1724  {
1725  xw = ex(i);
1726  iv->fx[i] *= (1.0 + xw*xw) * M_PI/2;
1727  }
1728  neval++;
1729  if (! octave::math::isfinite (iv->fx[i]))
1730  {
1731  nans[nnans++] = i;
1732  iv->fx[i] = 0.0;
1733  }
1734  }
1735  Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
1736  Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
1737  Vinvfx (iv->fx, &(iv->c[0]), 0);
1738  for (i = 0; i < nnans; i++)
1739  iv->fx[nans[i]] = octave::numeric_limits<double>::NaN ();
1740  iv->a = iivals[j];
1741  iv->b = iivals[j + 1];
1742  iv->depth = 3;
1743  iv->rdepth = 1;
1744  iv->ndiv = 0;
1745  iv->igral = 2 * h * iv->c[idx[3]] * w;
1746  nc = 0.0;
1747  for (i = n[2] + 1; i <= n[3]; i++)
1748  {
1749  temp = iv->c[idx[3] + i];
1750  nc += temp * temp;
1751  }
1752  ncdiff = nc;
1753  for (i = 0; i <= n[2]; i++)
1754  {
1755  temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
1756  ncdiff += temp * temp;
1757  temp = iv->c[idx[3] + i];
1758  nc += temp * temp;
1759  }
1760  ncdiff = sqrt (ncdiff);
1761  nc = sqrt (nc);
1762  iv->err = ncdiff * 2 * h;
1763  if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
1764  iv->err = 2 * h * nc;
1765 
1766  // Tabulate this interval's data.
1767  igral += iv->igral;
1768  err += iv->err;
1769 
1770  // Sift it up the heap.
1771  i = j;
1772  while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err)
1773  {
1774  std::swap (heap[i], heap[i / 2]);
1775  i /= 2;
1776  }
1777  }
1778 
1779  // Initialize some global values.
1780  igral_final = 0.0;
1781  err_final = 0.0;
1782 
1783  // Main loop.
1784  while (nivals > 0
1785  && err > std::max (abstol, fabs (igral) * reltol)
1786  && ! (err_final > std::max (abstol, fabs (igral) * reltol)
1787  && err - err_final < std::max (abstol, fabs (igral) * reltol)))
1788  {
1789  // Allow the user to interrupt.
1790  octave_quit ();
1791 
1792  // Put our finger on the interval with the largest error.
1793  iv = &(ivals[heap[0]]);
1794  m = (iv->a + iv->b) / 2;
1795  h = (iv->b - iv->a) / 2;
1796 
1797 #if (DEBUG_QUADCC)
1798  printf ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1799  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1800 #endif
1801 
1802  // Should we try to increase the degree?
1803  if (iv->depth < 3)
1804  {
1805  // Keep tabs on some variables.
1806  d = ++iv->depth;
1807 
1808  // Get the new (missing) function values.
1809  {
1810  ColumnVector ex (n[d] / 2);
1811  if (wrap)
1812  {
1813  for (i = 0; i < n[d] / 2; i++)
1814  ex(i) = tan (M_PI/2 * (m + xi[(2*i + 1) * skip[d]] * h));
1815  }
1816  else
1817  {
1818  for (i = 0; i < n[d] / 2; i++)
1819  ex(i) = m + xi[(2*i + 1) * skip[d]] * h;
1820  }
1821  fargs(0) = ex;
1822  fvals = octave::feval (fcn, fargs, 1);
1823  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1824  error ("quadcc: integrand F must return a single, real-valued vector");
1825 
1826  Matrix effex = fvals(0).matrix_value ();
1827  if (effex.numel () != ex.numel ())
1828  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1829 
1830  neval += effex.numel ();
1831  for (i = 0; i < n[d] / 2; i++)
1832  {
1833  j = (2*i + 1) * skip[d];
1834  iv->fx[j] = effex(i);
1835  if (wrap)
1836  {
1837  xw = ex(i);
1838  iv->fx[j] *= (1.0 + xw*xw) * M_PI/2;
1839  }
1840  }
1841  }
1842  nnans = 0;
1843  for (i = 0; i <= 32; i += skip[d])
1844  {
1845  if (! octave::math::isfinite (iv->fx[i]))
1846  {
1847  nans[nnans++] = i;
1848  iv->fx[i] = 0.0;
1849  }
1850  }
1851 
1852  // Compute the new coefficients.
1853  Vinvfx (iv->fx, &(iv->c[idx[d]]), d);
1854  // Downdate any NaNs.
1855  if (nnans > 0)
1856  {
1857  downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
1858  for (i = 0; i < nnans; i++)
1859  iv->fx[nans[i]] = octave::numeric_limits<double>::NaN ();
1860  }
1861 
1862  // Compute the error estimate.
1863  nc = 0.0;
1864  for (i = n[d - 1] + 1; i <= n[d]; i++)
1865  {
1866  temp = iv->c[idx[d] + i];
1867  nc += temp * temp;
1868  }
1869  ncdiff = nc;
1870  for (i = 0; i <= n[d - 1]; i++)
1871  {
1872  temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
1873  ncdiff += temp * temp;
1874  temp = iv->c[idx[d] + i];
1875  nc += temp * temp;
1876  }
1877  ncdiff = sqrt (ncdiff);
1878  nc = sqrt (nc);
1879  iv->err = ncdiff * 2 * h;
1880  // Compute the local integral.
1881  iv->igral = 2 * h * w * iv->c[idx[d]];
1882  // Split the interval prematurely?
1883  split = (nc > 0 && ncdiff / nc > 0.1);
1884  }
1885  else
1886  {
1887  // Maximum degree reached, just split.
1888  split = 1;
1889  }
1890 
1891  // Should we drop this interval?
1892  if ((m + h*xi[0]) >= (m + h*xi[1])
1893  || (m + h*xi[31]) >= (m + h*xi[32])
1894  || iv->err < fabs (iv->igral) * DROP_RELTOL)
1895  {
1896 #if (DEBUG_QUADCC)
1897  printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1898  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1899 #endif
1900 
1901  // Keep this interval's contribution.
1902  err_final += iv->err;
1903  igral_final += iv->igral;
1904  // Swap with the last element on the heap.
1905  std::swap (heap[nivals - 1], heap[0]);
1906  nivals--;
1907  // Fix up the heap.
1908  i = 0;
1909  while (2*i + 1 < nivals)
1910  {
1911  // Get the kids.
1912  j = 2*i + 1;
1913  // If the j+1st entry exists and is larger than the jth,
1914  // use it instead.
1915  if (j + 1 < nivals
1916  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
1917  j++;
1918  // Do we need to move the ith entry up?
1919  if (ivals[heap[j]].err <= ivals[heap[i]].err)
1920  break;
1921  else
1922  {
1923  std::swap (heap[j], heap[i]);
1924  i = j;
1925  }
1926  }
1927  }
1928  else if (split)
1929  {
1930  // Some values we will need often...
1931  d = iv->depth;
1932  // Generate the interval on the left.
1933  ivl = &(ivals[heap[nivals++]]);
1934  ivl->a = iv->a;
1935  ivl->b = m;
1936  ml = (ivl->a + ivl->b) / 2;
1937  hl = h / 2;
1938  ivl->depth = 0;
1939  ivl->rdepth = iv->rdepth + 1;
1940  ivl->fx[0] = iv->fx[0];
1941  ivl->fx[32] = iv->fx[16];
1942  {
1943  ColumnVector ex (n[0] - 1);
1944  if (wrap)
1945  {
1946  for (i = 0; i < n[0] - 1; i++)
1947  ex(i) = tan (M_PI/2 * (ml + xi[(i + 1) * skip[0]] * hl));
1948  }
1949  else
1950  {
1951  for (i = 0; i < n[0] - 1; i++)
1952  ex(i) = ml + xi[(i + 1) * skip[0]] * hl;
1953  }
1954  fargs(0) = ex;
1955  fvals = octave::feval (fcn, fargs, 1);
1956  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1957  error ("quadcc: integrand F must return a single, real-valued vector");
1958 
1959  Matrix effex = fvals(0).matrix_value ();
1960  if (effex.numel () != ex.numel ())
1961  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1962 
1963  neval += effex.numel ();
1964  for (i = 0; i < n[0] - 1; i++)
1965  {
1966  j = (i + 1) * skip[0];
1967  ivl->fx[j] = effex(i);
1968  if (wrap)
1969  {
1970  xw = ex(i);
1971  ivl->fx[j] *= (1.0 + xw*xw) * M_PI/2;
1972  }
1973  }
1974  }
1975  nnans = 0;
1976  for (i = 0; i <= 32; i += skip[0])
1977  {
1978  if (! octave::math::isfinite (ivl->fx[i]))
1979  {
1980  nans[nnans++] = i;
1981  ivl->fx[i] = 0.0;
1982  }
1983  }
1984  Vinvfx (ivl->fx, ivl->c, 0);
1985  if (nnans > 0)
1986  {
1987  downdate (ivl->c, n[0], 0, nans, nnans);
1988  for (i = 0; i < nnans; i++)
1989  ivl->fx[nans[i]] = octave::numeric_limits<double>::NaN ();
1990  }
1991  for (i = 0; i <= n[d]; i++)
1992  {
1993  ivl->c[idx[d] + i] = 0.0;
1994  for (j = i; j <= n[d]; j++)
1995  ivl->c[idx[d] + i] += Tleft[i*33 + j] * iv->c[idx[d] + j];
1996  }
1997  ncdiff = 0.0;
1998  for (i = 0; i <= n[0]; i++)
1999  {
2000  temp = ivl->c[i] - ivl->c[idx[d] + i];
2001  ncdiff += temp * temp;
2002  }
2003  for (i = n[0] + 1; i <= n[d]; i++)
2004  {
2005  temp = ivl->c[idx[d] + i];
2006  ncdiff += temp * temp;
2007  }
2008  ncdiff = sqrt (ncdiff);
2009  ivl->err = ncdiff * h;
2010  // Check for divergence.
2011  ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2012  && ivl->c[0] / iv->c[0] > 2);
2013  if (ivl->ndiv > ndiv_max && 2*ivl->ndiv > ivl->rdepth)
2014  {
2016  warning ("quadcc: divergent integral detected");
2017  break;
2018  }
2019 
2020  // Compute the local integral.
2021  ivl->igral = h * w * ivl->c[0];
2022 
2023  // Generate the interval on the right.
2024  ivr = &(ivals[heap[nivals++]]);
2025  ivr->a = m;
2026  ivr->b = iv->b;
2027  mr = (ivr->a + ivr->b) / 2;
2028  hr = h / 2;
2029  ivr->depth = 0;
2030  ivr->rdepth = iv->rdepth + 1;
2031  ivr->fx[0] = iv->fx[16];
2032  ivr->fx[32] = iv->fx[32];
2033  {
2034  ColumnVector ex (n[0] - 1);
2035  if (wrap)
2036  {
2037  for (i = 0; i < n[0] - 1; i++)
2038  ex(i) = tan (M_PI/2 * (mr + xi[(i + 1) * skip[0]] * hr));
2039  }
2040  else
2041  {
2042  for (i = 0; i < n[0] - 1; i++)
2043  ex(i) = mr + xi[(i + 1) * skip[0]] * hr;
2044  }
2045  fargs(0) = ex;
2046  fvals = octave::feval (fcn, fargs, 1);
2047  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
2048  error ("quadcc: integrand F must return a single, real-valued vector");
2049 
2050  Matrix effex = fvals(0).matrix_value ();
2051  if (effex.numel () != ex.numel ())
2052  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
2053 
2054  neval += effex.numel ();
2055  for (i = 0; i < n[0] - 1; i++)
2056  {
2057  j = (i + 1) * skip[0];
2058  ivr->fx[j] = effex(i);
2059  if (wrap)
2060  {
2061  xw = ex(i);
2062  ivr->fx[j] *= (1.0 + xw*xw) * M_PI/2;
2063  }
2064  }
2065  }
2066  nnans = 0;
2067  for (i = 0; i <= 32; i += skip[0])
2068  {
2069  if (! octave::math::isfinite (ivr->fx[i]))
2070  {
2071  nans[nnans++] = i;
2072  ivr->fx[i] = 0.0;
2073  }
2074  }
2075  Vinvfx (ivr->fx, ivr->c, 0);
2076  if (nnans > 0)
2077  {
2078  downdate (ivr->c, n[0], 0, nans, nnans);
2079  for (i = 0; i < nnans; i++)
2080  ivr->fx[nans[i]] = octave::numeric_limits<double>::NaN ();
2081  }
2082  for (i = 0; i <= n[d]; i++)
2083  {
2084  ivr->c[idx[d] + i] = 0.0;
2085  for (j = i; j <= n[d]; j++)
2086  ivr->c[idx[d] + i] += Tright[i*33 + j] * iv->c[idx[d] + j];
2087  }
2088  ncdiff = 0.0;
2089  for (i = 0; i <= n[0]; i++)
2090  {
2091  temp = ivr->c[i] - ivr->c[idx[d] + i];
2092  ncdiff += temp * temp;
2093  }
2094  for (i = n[0] + 1; i <= n[d]; i++)
2095  {
2096  temp = ivr->c[idx[d] + i];
2097  ncdiff += temp * temp;
2098  }
2099  ncdiff = sqrt (ncdiff);
2100  ivr->err = ncdiff * h;
2101  // Check for divergence.
2102  ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2103  && ivr->c[0] / iv->c[0] > 2);
2104  if (ivr->ndiv > ndiv_max && 2*ivr->ndiv > ivr->rdepth)
2105  {
2107  warning ("quadcc: divergent integral detected");
2108  break;
2109  }
2110 
2111  // Compute the local integral.
2112  ivr->igral = h * w * ivr->c[0];
2113 
2114  // Fix-up the heap: we now have one interval on top that we
2115  // don't need any more and two new, unsorted ones at the bottom.
2116 
2117  // Flip the last interval to the top of the heap and sift down.
2118  std::swap (heap[nivals - 1], heap[0]);
2119  nivals--;
2120  // Sift this interval back down the heap.
2121  i = 0;
2122  while (2*i + 1 < nivals - 1)
2123  {
2124  j = 2*i + 1;
2125  if (j + 1 < nivals - 1
2126  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2127  j++;
2128  if (ivals[heap[j]].err <= ivals[heap[i]].err)
2129  break;
2130  else
2131  {
2132  std::swap (heap[j], heap[i]);
2133  i = j;
2134  }
2135  }
2136 
2137  // Now grab the last interval and sift it up the heap.
2138  i = nivals - 1;
2139  while (i > 0)
2140  {
2141  j = (i - 1) / 2;
2142  if (ivals[heap[j]].err < ivals[heap[i]].err)
2143  {
2144  std::swap (heap[j], heap[i]);
2145  i = j;
2146  }
2147  else
2148  break;
2149  }
2150  }
2151  else
2152  {
2153  // Otherwise, just fix-up the heap.
2154  i = 0;
2155  while (2*i + 1 < nivals)
2156  {
2157  j = 2*i + 1;
2158  if (j + 1 < nivals
2159  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2160  j++;
2161  if (ivals[heap[j]].err <= ivals[heap[i]].err)
2162  break;
2163  else
2164  {
2165  std::swap (heap[j], heap[i]);
2166  i = j;
2167  }
2168  }
2169  }
2170 
2171  // If the heap is about to overflow, remove the last two intervals.
2172  while (nivals > cquad_heapsize - 2)
2173  {
2174  iv = &(ivals[heap[nivals - 1]]);
2175 #if (DEBUG_QUADCC)
2176  printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
2177  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
2178 #endif
2179  err_final += iv->err;
2180  igral_final += iv->igral;
2181  nivals--;
2182  }
2183 
2184  // Collect the value of the integral and error.
2185  igral = igral_final;
2186  err = err_final;
2187  for (i = 0; i < nivals; i++)
2188  {
2189  igral += ivals[heap[i]].igral;
2190  err += ivals[heap[i]].err;
2191  }
2192  }
2193 
2194 #if (DEBUG_QUADCC)
2195  // Dump the contents of the heap.
2196  for (i = 0; i < nivals; i++)
2197  {
2198  iv = &(ivals[heap[i]]);
2199  printf ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n",
2200  i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
2201  iv->rdepth, iv->ndiv);
2202  }
2203 #endif
2204 
2205  if (issingle)
2206  return ovl (static_cast<float> (igral), err, neval);
2207  else
2208  return ovl (igral, err, neval);
2209 }
2210 
2211 /*
2212 %!assert (quadcc (@sin, -pi, pi), 0, 1e-10)
2213 %!assert (quadcc (inline ("sin"), -pi, pi), 0, 1e-10)
2214 %!assert (quadcc ("sin", -pi, pi), 0, 1e-10)
2215 
2216 %!assert (quadcc (@sin, -pi, 0), -2, 1e-10)
2217 %!assert (quadcc (@sin, 0, pi), 2, 1e-10)
2218 %!assert (quadcc (@(x) 1./sqrt (x), 0, 1), 2, -1e-6)
2219 %!assert (quadcc (@(x) 1./(sqrt (x).*(x+1)), 0, Inf), pi, -1e-6)
2220 %!assert (quadcc (@(x) 1./(sqrt (x).*(x+1)), 0, Inf, [0, 1e-8]), pi, -1e-8)
2221 
2222 %!assert (quadcc (@(x) exp (-x .^ 2), -Inf, Inf), sqrt (pi), 1e-10)
2223 %!assert (quadcc (@(x) exp (-x .^ 2), -Inf, 0), sqrt (pi)/2, 1e-10)
2224 
2225 ## Test function with NaNs in interval
2226 %!function y = __nansin (x)
2227 %! nan_locs = [-3*pi/4, -pi/4, 0, pi/3, pi/2, pi];
2228 %! y = sin (x);
2229 %! idx = min (abs (bsxfun (@minus, x(:), nan_locs)), [], 2);
2230 %! y(idx < 1e-10) = NaN;
2231 %!endfunction
2232 
2233 %!test
2234 %! [q, err, npoints] = quadcc ("__nansin", -pi, pi, [0, 1e-6]);
2235 %! assert (q, 0, -1e-6);
2236 %! assert (err, 0, 15*eps);
2237 
2238 %!test
2239 %! assert (class (quadcc (@sin, 0, 1)), "double");
2240 %! assert (class (quadcc (@sin, single (0), 1)), "single");
2241 %! assert (class (quadcc (@sin, 0, single (1))), "single");
2242 %! assert (class (quadcc (@sin, single (0), single (1))), "single");
2243 
2244 ## Test input validation
2245 %!error (quadcc ())
2246 %!error (quadcc (@sin))
2247 %!error (quadcc (@sin, 0))
2248 %!error <lower limit .* must be a .* scalar> (quadcc (@sin, ones (2), pi))
2249 %!error <lower limit .* must be a real scalar> (quadcc (@sin, -i, pi))
2250 %!error <upper limit .* must be a .* scalar> (quadcc (@sin, 0, ones (2)))
2251 %!error <upper limit .* must be a real scalar> (quadcc (@sin, 0, i))
2252 %!error <TOL must be a 1- or 2-element vector> (quadcc (@sin, 0, pi, {1}))
2253 %!error <TOL must be a 1- or 2-element vector> (quadcc (@sin, 0, pi, [1,2,3]))
2254 %!error <absolute tolerance must be .=0> (quadcc (@sin, 0, pi, -1))
2255 %!error <relative tolerance must be .=0> (quadcc (@sin, 0, pi, [1, -1]))
2256 %!error <SING.* must be .* real values> (quadcc (@sin, 0, pi, 1e-6, [ i ]))
2257 */
#define Inf
Definition: Faddeeva.cc:247
#define NaN
Definition: Faddeeva.cc:248
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:230
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
Definition: dMatrix.h:42
octave_idx_type length(void) const
Definition: ovl.h:113
OCTINTERP_API void print_usage(void)
Definition: defun.cc:53
#define DEFMETHOD(name, interp_name, args_name, nargout_name, doc)
Macro to define a builtin method.
Definition: defun.h:138
void warning(const char *fmt,...)
Definition: error.cc:1050
void error(const char *fmt,...)
Definition: error.cc:968
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
T octave_idx_type m
Definition: mx-inlines.cc:773
octave_idx_type n
Definition: mx-inlines.cc:753
std::complex< double > w(std::complex< double > z, double relerr=0)
Complex atan(const Complex &x)
Definition: lo-mappers.h:71
bool isfinite(double x)
Definition: lo-mappers.h:192
double copysign(double x, double y)
Definition: lo-mappers.h:51
bool isinf(double x)
Definition: lo-mappers.h:203
octave_value get_function_handle(interpreter &interp, const octave_value &arg, const std::string &parameter_name)
octave_value_list feval(const char *name, const octave_value_list &args, int nargout)
Evaluate an Octave function (built-in or interpreted) and return the list of result values.
Definition: oct-parse.cc:9580
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition: ovl.h:211
static const double V3inv[17 *17]
Definition: quadcc.cc:188
static const double V4inv[33 *33]
Definition: quadcc.cc:326
static const double V2inv[9 *9]
Definition: quadcc.cc:153
static const double DROP_RELTOL
Definition: quadcc.cc:62
static const double Lalpha[33]
Definition: quadcc.cc:110
static const double Tleft[33 *33]
Definition: quadcc.cc:903
static const double xi[33]
Definition: quadcc.cc:66
static const double bee[68]
Definition: quadcc.cc:83
static const double V1inv[5 *5]
Definition: quadcc.cc:140
static void downdate(double *c, int n, int d, int *nans, int nnans)
Definition: quadcc.cc:1462
static const int MIN_CQUAD_HEAPSIZE
Definition: quadcc.cc:49
static const double Lgamma[33]
Definition: quadcc.cc:125
static void Vinvfx(const double *fx, double *c, const int d)
Definition: quadcc.cc:1417
static const double Tright[33 *33]
Definition: quadcc.cc:1158
int ndiv
Definition: quadcc.cc:58
double fx[33]
Definition: quadcc.cc:56
double err
Definition: quadcc.cc:57
int depth
Definition: quadcc.cc:58
double b
Definition: quadcc.cc:54
int rdepth
Definition: quadcc.cc:58
double c[64]
Definition: quadcc.cc:55
double igral
Definition: quadcc.cc:57
double a
Definition: quadcc.cc:54