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