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