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