GNU Octave  4.0.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-2015 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 "oct-locbuf.h"
29 
30 #include "defun.h"
31 #include "error.h"
32 #include "oct-obj.h"
33 #include "parse.h"
34 #include "utils.h"
35 #include "variables.h"
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 doubly-adaptive @nospell{Clenshaw-Curtis} quadrature.\n\
1496 \n\
1497 @var{f} is a function handle, inline function, or string containing the name\n\
1498 of the function to evaluate. The function @var{f} must be vectorized and\n\
1499 must return a vector of output values if given a vector of input values.\n\
1500 For example,\n\
1501 \n\
1502 @example\n\
1503 f = @@(x) x .* sin (1./x) .* sqrt (abs (1 - x));\n\
1504 @end example\n\
1505 \n\
1506 @noindent\n\
1507 which uses the element-by-element ``dot'' form for all operators.\n\
1508 \n\
1509 @var{a} and @var{b} are the lower and upper limits of integration. Either\n\
1510 or both limits may be infinite. @code{quadcc} handles an inifinite limit\n\
1511 by substituting the variable of integration with @code{x = tan (pi/2*u)}.\n\
1512 \n\
1513 The optional argument @var{tol} defines the relative tolerance used to stop\n\
1514 the integration procedure. The default value is @math{1e^{-6}}.\n\
1515 \n\
1516 The optional argument @var{sing} contains a list of points where the\n\
1517 integrand has known singularities, or discontinuities\n\
1518 in any of its derivatives, inside the integration interval.\n\
1519 For the example above, which has a discontinuity at x=1, the call to\n\
1520 @code{quadcc} would be as follows\n\
1521 \n\
1522 @example\n\
1523 int = quadcc (f, a, b, 1.0e-6, [ 1 ]);\n\
1524 @end example\n\
1525 \n\
1526 The result of the integration is returned in @var{q}.\n\
1527 \n\
1528 @var{err} is an estimate of the absolute integration error.\n\
1529 \n\
1530 @var{nr_points} is the number of points at which the integrand was evaluated.\n\
1531 \n\
1532 If the adaptive integration did not converge, the value of @var{err} will be\n\
1533 larger than the requested tolerance. Therefore, it is recommended to verify\n\
1534 this value for difficult integrands.\n\
1535 \n\
1536 @code{quadcc} is capable of dealing with non-numeric values of the integrand\n\
1537 such as @code{NaN} or @code{Inf}. If the integral diverges, and\n\
1538 @code{quadcc} detects this, then a warning is issued and @code{Inf} or\n\
1539 @code{-Inf} is returned.\n\
1540 \n\
1541 Note: @code{quadcc} is a general purpose quadrature algorithm and, as such,\n\
1542 may be less efficient for a smooth or otherwise well-behaved integrand than\n\
1543 other methods such as @code{quadgk}.\n\
1544 \n\
1545 The algorithm uses @nospell{Clenshaw-Curtis} quadrature rules of increasing\n\
1546 degree in each interval and bisects the interval if either the function does\n\
1547 not appear to be smooth or a rule of maximum degree has been reached. The\n\
1548 error estimate is computed from the L2-norm of the difference between two\n\
1549 successive interpolations of the integrand over the nodes of the respective\n\
1550 quadrature rules.\n\
1551 \n\
1552 Reference: @nospell{P. Gonnet}, @cite{Increasing the Reliability of Adaptive\n\
1553 Quadrature Using Explicit Interpolants}, ACM Transactions on\n\
1554 Mathematical Software, Vol. 37, Issue 3, Article No. 3, 2010.\n\
1555 @seealso{quad, quadv, quadl, quadgk, trapz, dblquad, triplequad}\n\
1556 @end deftypefn")
1557 {
1558  octave_value_list retval;
1559 
1560  /* Some constants that we will need. */
1561  static const int n[4] = { 4, 8, 16, 32 };
1562  static const int skip[4] = { 8, 4, 2, 1 };
1563  static const int idx[4] = { 0, 5, 14, 31 };
1564  static const double w = M_SQRT2 / 2;
1565  static const int ndiv_max = 20;
1566 
1567  /* Arguments left and right */
1568  int nargin = args.length ();
1569  octave_function *fcn;
1570  double a, b, tol, *sing;
1571 
1572  /* Variables needed for transforming the integrand. */
1573  bool wrap = false;
1574  double xw;
1575 
1576  /* Stuff we will need to call the integrand. */
1577  octave_value_list fargs, fvals;
1578 
1579  /* Actual variables (as opposed to constants above). */
1580  double m, h, ml, hl, mr, hr, temp;
1581  double igral, err, igral_final, err_final;
1582  int nivals, neval = 0;
1583  int i, j, d, split, t;
1584  int nnans, nans[33];
1585  cquad_ival *iv, *ivl, *ivr;
1586  double nc, ncdiff;
1587 
1588 
1589  /* Parse the input arguments. */
1590  if (nargin < 3)
1591  {
1592  print_usage ();
1593  return retval;
1594  }
1595 
1596  if (args(0).is_function_handle () || args(0).is_inline_function ())
1597  fcn = args(0).function_value ();
1598  else
1599  {
1600  std::string fcn_name = unique_symbol_name ("__quadcc_fcn__");
1601  std::string fname = "function y = ";
1602  fname.append (fcn_name);
1603  fname.append ("(x) y = ");
1604  fcn = extract_function (args(0), "quadcc", fcn_name, fname,
1605  "; endfunction");
1606  }
1607 
1608  if (! args(1).is_real_scalar ())
1609  {
1610  error ("quadcc: lower limit of integration (A) must be a single real scalar");
1611  return retval;
1612  }
1613  else
1614  a = args(1).double_value ();
1615 
1616  if (! args(2).is_real_scalar ())
1617  {
1618  error ("quadcc: upper limit of integration (B) must be a single real scalar");
1619  return retval;
1620  }
1621  else
1622  b = args(2).double_value ();
1623 
1624  if (nargin < 4 || args(3).is_empty ())
1625  tol = 1.0e-6;
1626  else if (! args(3).is_real_scalar () || args(3).double_value () <= 0)
1627  {
1628  error ("quadcc: tolerance (TOL) must be a single real scalar > 0");
1629  return retval;
1630  }
1631  else
1632  tol = args(3).double_value ();
1633 
1634  if (nargin < 5)
1635  {
1636  nivals = 1;
1637  }
1638  else if (!(args(4).is_real_scalar () || args(4).is_real_matrix ()))
1639  {
1640  error ("quadcc: list of singularities (SING) must be a vector of real values");
1641  return retval;
1642  }
1643  else
1644  {
1645  nivals = 1 + args(4).numel ();
1646  }
1647 
1648  int cquad_heapsize = (nivals >= min_cquad_heapsize ? nivals + 1
1649  : min_cquad_heapsize);
1650  /* The interval heap. */
1651  OCTAVE_LOCAL_BUFFER (cquad_ival, ivals, cquad_heapsize);
1652  OCTAVE_LOCAL_BUFFER (double, iivals, cquad_heapsize);
1653  OCTAVE_LOCAL_BUFFER (int, heap, cquad_heapsize);
1654 
1655  if (nivals == 1)
1656  {
1657  iivals[0] = a;
1658  iivals[1] = b;
1659  }
1660  else
1661  {
1662  // Intervals around singularities
1663  sing = args(4).array_value ().fortran_vec ();
1664  iivals[0] = a;
1665  for (i = 0; i < nivals - 1; i++)
1666  iivals[i + 1] = sing[i];
1667  iivals[nivals] = b;
1668  }
1669 
1670  /* If a or b are +/-Inf, transform the integral. */
1671  if (xisinf (a) || xisinf (b))
1672  {
1673  wrap = true;
1674  for (i = 0; i < nivals + 1; i++)
1675  if (xisinf (iivals[i]))
1676  iivals[i] = gnulib::copysign (1.0, iivals[i]);
1677  else
1678  iivals[i] = 2.0 * atan (iivals[i]) / M_PI;
1679  }
1680 
1681 
1682  /* Initialize the heaps. */
1683  for (i = 0; i < cquad_heapsize; i++)
1684  heap[i] = i;
1685 
1686  /* Create the first interval(s). */
1687  igral = 0.0;
1688  err = 0.0;
1689  for (j = 0; j < nivals; j++)
1690  {
1691 
1692  /* Initialize the interval. */
1693  iv = &(ivals[heap[j]]);
1694  m = (iivals[j] + iivals[j + 1]) / 2;
1695  h = (iivals[j + 1] - iivals[j]) / 2;
1696  nnans = 0;
1697  ColumnVector ex (33);
1698  if (wrap)
1699  {
1700  for (i = 0; i <= n[3]; i++)
1701  ex(i) = tan (M_PI / 2 * (m + xi[i] * h));
1702  }
1703  else
1704  {
1705  for (i = 0; i <= n[3]; i++)
1706  ex(i) = m + xi[i] * h;
1707  }
1708  fargs(0) = ex;
1709  fvals = feval (fcn, fargs, 1);
1710  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1711  {
1712  error ("quadcc: integrand F must return a single, real-valued vector");
1713  return retval;
1714  }
1715  Matrix effex = fvals(0).matrix_value ();
1716  if (effex.length () != ex.length ())
1717  {
1718  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1719  return retval;
1720  }
1721  for (i = 0; i <= n[3]; i++)
1722  {
1723  iv->fx[i] = effex(i);
1724  if (wrap)
1725  {
1726  xw = ex(i);
1727  iv->fx[i] *= (1.0 + xw * xw) * M_PI / 2;
1728  }
1729  neval++;
1730  if (! xfinite (iv->fx[i]))
1731  {
1732  nans[nnans++] = i;
1733  iv->fx[i] = 0.0;
1734  }
1735  }
1736  Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
1737  Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
1738  Vinvfx (iv->fx, &(iv->c[0]), 0);
1739  for (i = 0; i < nnans; i++)
1740  iv->fx[nans[i]] = octave_NaN;
1741  iv->a = iivals[j];
1742  iv->b = iivals[j + 1];
1743  iv->depth = 3;
1744  iv->rdepth = 1;
1745  iv->ndiv = 0;
1746  iv->igral = 2 * h * iv->c[idx[3]] * w;
1747  nc = 0.0;
1748  for (i = n[2] + 1; i <= n[3]; i++)
1749  {
1750  temp = iv->c[idx[3] + i];
1751  nc += temp * temp;
1752  }
1753  ncdiff = nc;
1754  for (i = 0; i <= n[2]; i++)
1755  {
1756  temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
1757  ncdiff += temp * temp;
1758  nc += iv->c[idx[3] + i] * iv->c[idx[3] + i];
1759  }
1760  ncdiff = sqrt (ncdiff);
1761  nc = sqrt (nc);
1762  iv->err = ncdiff * 2 * h;
1763  if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
1764  iv->err = 2 * h * nc;
1765 
1766  /* Tabulate this interval's data. */
1767  igral += iv->igral;
1768  err += iv->err;
1769 
1770  /* Sift it up the heap. */
1771  i = j;
1772  while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err)
1773  {
1774  temp = heap[i];
1775  heap[i] = heap[i / 2];
1776  heap[i / 2] = temp;
1777  i /= 2;
1778  }
1779 
1780  }
1781 
1782 
1783  /* Initialize some global values. */
1784  igral_final = 0.0;
1785  err_final = 0.0;
1786 
1787 
1788  /* Main loop. */
1789  while (nivals > 0 && err > 0.0 && err > fabs (igral) * tol
1790  && !(err_final > fabs (igral) * tol
1791  && err - err_final < fabs (igral) * tol))
1792  {
1793 
1794  /* Allow the user to interrupt. */
1795  OCTAVE_QUIT;
1796 
1797  /* Put our finger on the interval with the largest error. */
1798  iv = &(ivals[heap[0]]);
1799  m = (iv->a + iv->b) / 2;
1800  h = (iv->b - iv->a) / 2;
1801 
1802 #if (DEBUG_QUADCC)
1803  printf ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1804  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1805 #endif
1806 
1807  /* Should we try to increase the degree? */
1808  if (iv->depth < 3)
1809  {
1810 
1811  /* Keep tabs on some variables. */
1812  d = ++iv->depth;
1813 
1814  /* Get the new (missing) function values */
1815  {
1816  ColumnVector ex (n[d] / 2);
1817  if (wrap)
1818  {
1819  for (i = 0; i < n[d] / 2; i++)
1820  ex(i) = tan (M_PI / 2 * (m + xi[(2 * i + 1) * skip[d]] * h));
1821  }
1822  else
1823  {
1824  for (i = 0; i < n[d] / 2; i++)
1825  ex(i) = m + xi[(2 * i + 1) * skip[d]] * h;
1826  }
1827  fargs(0) = ex;
1828  fvals = feval (fcn, fargs, 1);
1829  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1830  {
1831  error ("quadcc: integrand F must return a single, real-valued vector");
1832  return retval;
1833  }
1834  Matrix effex = fvals(0).matrix_value ();
1835  if (effex.length () != ex.length ())
1836  {
1837  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1838  return retval;
1839  }
1840  neval += effex.length ();
1841  for (i = 0; i < n[d] / 2; i++)
1842  {
1843  j = (2 * i + 1) * skip[d];
1844  iv->fx[j] = effex(i);
1845  if (wrap)
1846  {
1847  xw = ex(i);
1848  iv->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
1849  }
1850  }
1851  }
1852  nnans = 0;
1853  for (i = 0; i <= 32; i += skip[d])
1854  {
1855  if (! xfinite (iv->fx[i]))
1856  {
1857  nans[nnans++] = i;
1858  iv->fx[i] = 0.0;
1859  }
1860  }
1861 
1862  /* Compute the new coefficients. */
1863  Vinvfx (iv->fx, &(iv->c[idx[d]]), d);
1864  /* Downdate any NaNs. */
1865  if (nnans > 0)
1866  {
1867  downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
1868  for (i = 0; i < nnans; i++)
1869  iv->fx[nans[i]] = octave_NaN;
1870  }
1871 
1872  /* Compute the error estimate. */
1873  nc = 0.0;
1874  for (i = n[d - 1] + 1; i <= n[d]; i++)
1875  {
1876  temp = iv->c[idx[d] + i];
1877  nc += temp * temp;
1878  }
1879  ncdiff = nc;
1880  for (i = 0; i <= n[d - 1]; i++)
1881  {
1882  temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
1883  ncdiff += temp * temp;
1884  nc += iv->c[idx[d] + i] * iv->c[idx[d] + i];
1885  }
1886  ncdiff = sqrt (ncdiff);
1887  nc = sqrt (nc);
1888  iv->err = ncdiff * 2 * h;
1889  /* Compute the local integral. */
1890  iv->igral = 2 * h * w * iv->c[idx[d]];
1891  /* Split the interval prematurely? */
1892  split = (nc > 0 && ncdiff / nc > 0.1);
1893  }
1894 
1895  /* Maximum degree reached, just split. */
1896  else
1897  {
1898  split = 1;
1899  }
1900 
1901 
1902  /* Should we drop this interval? */
1903  if ((m + h * xi[0]) >= (m + h * xi[1])
1904  || (m + h * xi[31]) >= (m + h * xi[32])
1905  || iv->err < fabs (iv->igral)
1906  * std::numeric_limits<double>::epsilon () * 10)
1907  {
1908 
1909 #if (DEBUG_QUADCC)
1910  printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1911  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1912 #endif
1913 
1914  /* Keep this interval's contribution */
1915  err_final += iv->err;
1916  igral_final += iv->igral;
1917  /* Swap with the last element on the heap */
1918  t = heap[nivals - 1];
1919  heap[nivals - 1] = heap[0];
1920  heap[0] = t;
1921  nivals--;
1922  /* Fix up the heap */
1923  i = 0;
1924  while (2 * i + 1 < nivals)
1925  {
1926  /* Get the kids */
1927  j = 2 * i + 1;
1928  /* If the j+1st entry exists and is larger than the jth,
1929  use it instead. */
1930  if (j + 1 < nivals
1931  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
1932  j++;
1933  /* Do we need to move the ith entry up? */
1934  if (ivals[heap[j]].err <= ivals[heap[i]].err)
1935  break;
1936  else
1937  {
1938  t = heap[j];
1939  heap[j] = heap[i];
1940  heap[i] = t;
1941  i = j;
1942  }
1943  }
1944 
1945  }
1946 
1947  /* Do we need to split this interval? */
1948  else if (split)
1949  {
1950 
1951  /* Some values we will need often... */
1952  d = iv->depth;
1953  /* Generate the interval on the left */
1954  ivl = &(ivals[heap[nivals++]]);
1955  ivl->a = iv->a;
1956  ivl->b = m;
1957  ml = (ivl->a + ivl->b) / 2;
1958  hl = h / 2;
1959  ivl->depth = 0;
1960  ivl->rdepth = iv->rdepth + 1;
1961  ivl->fx[0] = iv->fx[0];
1962  ivl->fx[32] = iv->fx[16];
1963  {
1964  ColumnVector ex (n[0] - 1);
1965  if (wrap)
1966  {
1967  for (i = 0; i < n[0] - 1; i++)
1968  ex(i) = tan (M_PI / 2 * (ml + xi[(i + 1) * skip[0]] * hl));
1969  }
1970  else
1971  {
1972  for (i = 0; i < n[0] - 1; i++)
1973  ex(i) = ml + xi[(i + 1) * skip[0]] * hl;
1974  }
1975  fargs(0) = ex;
1976  fvals = feval (fcn, fargs, 1);
1977  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1978  {
1979  error ("quadcc: integrand F must return a single, real-valued vector");
1980  return retval;
1981  }
1982  Matrix effex = fvals(0).matrix_value ();
1983  if (effex.length () != ex.length ())
1984  {
1985  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1986  return retval;
1987  }
1988  neval += effex.length ();
1989  for (i = 0; i < n[0] - 1; i++)
1990  {
1991  j = (i + 1) * skip[0];
1992  ivl->fx[j] = effex(i);
1993  if (wrap)
1994  {
1995  xw = ex(i);
1996  ivl->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
1997  }
1998  }
1999  }
2000  nnans = 0;
2001  for (i = 0; i <= 32; i += skip[0])
2002  {
2003  if (! xfinite (ivl->fx[i]))
2004  {
2005  nans[nnans++] = i;
2006  ivl->fx[i] = 0.0;
2007  }
2008  }
2009  Vinvfx (ivl->fx, ivl->c, 0);
2010  if (nnans > 0)
2011  {
2012  downdate (ivl->c, n[0], 0, nans, nnans);
2013  for (i = 0; i < nnans; i++)
2014  ivl->fx[nans[i]] = octave_NaN;
2015  }
2016  for (i = 0; i <= n[d]; i++)
2017  {
2018  ivl->c[idx[d] + i] = 0.0;
2019  for (j = i; j <= n[d]; j++)
2020  ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j];
2021  }
2022  ncdiff = 0.0;
2023  for (i = 0; i <= n[0]; i++)
2024  {
2025  temp = ivl->c[i] - ivl->c[idx[d] + i];
2026  ncdiff += temp * temp;
2027  }
2028  for (i = n[0] + 1; i <= n[d]; i++)
2029  {
2030  temp = ivl->c[idx[d] + i];
2031  ncdiff += temp * temp;
2032  }
2033  ncdiff = sqrt (ncdiff);
2034  ivl->err = ncdiff * h;
2035  /* Check for divergence. */
2036  ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2037  && ivl->c[0] / iv->c[0] > 2);
2038  if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth)
2039  {
2040  igral = gnulib::copysign (octave_Inf, igral);
2041  warning ("quadcc: divergent integral detected");
2042  break;
2043  }
2044 
2045  /* Compute the local integral. */
2046  ivl->igral = h * w * ivl->c[0];
2047 
2048 
2049  /* Generate the interval on the right */
2050  ivr = &(ivals[heap[nivals++]]);
2051  ivr->a = m;
2052  ivr->b = iv->b;
2053  mr = (ivr->a + ivr->b) / 2;
2054  hr = h / 2;
2055  ivr->depth = 0;
2056  ivr->rdepth = iv->rdepth + 1;
2057  ivr->fx[0] = iv->fx[16];
2058  ivr->fx[32] = iv->fx[32];
2059  {
2060  ColumnVector ex (n[0] - 1);
2061  if (wrap)
2062  {
2063  for (i = 0; i < n[0] - 1; i++)
2064  ex(i) = tan (M_PI / 2 * (mr + xi[(i + 1) * skip[0]] * hr));
2065  }
2066  else
2067  {
2068  for (i = 0; i < n[0] - 1; i++)
2069  ex(i) = mr + xi[(i + 1) * skip[0]] * hr;
2070  }
2071  fargs(0) = ex;
2072  fvals = feval (fcn, fargs, 1);
2073  if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
2074  {
2075  error ("quadcc: integrand F must return a single, real-valued vector");
2076  return retval;
2077  }
2078  Matrix effex = fvals(0).matrix_value ();
2079  if (effex.length () != ex.length ())
2080  {
2081  error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
2082  return retval;
2083  }
2084  neval += effex.length ();
2085  for (i = 0; i < n[0] - 1; i++)
2086  {
2087  j = (i + 1) * skip[0];
2088  ivr->fx[j] = effex(i);
2089  if (wrap)
2090  {
2091  xw = ex(i);
2092  ivr->fx[j] *= (1.0 + xw * xw) * M_PI / 2;
2093  }
2094  }
2095  }
2096  nnans = 0;
2097  for (i = 0; i <= 32; i += skip[0])
2098  {
2099  if (! xfinite (ivr->fx[i]))
2100  {
2101  nans[nnans++] = i;
2102  ivr->fx[i] = 0.0;
2103  }
2104  }
2105  Vinvfx (ivr->fx, ivr->c, 0);
2106  if (nnans > 0)
2107  {
2108  downdate (ivr->c, n[0], 0, nans, nnans);
2109  for (i = 0; i < nnans; i++)
2110  ivr->fx[nans[i]] = octave_NaN;
2111  }
2112  for (i = 0; i <= n[d]; i++)
2113  {
2114  ivr->c[idx[d] + i] = 0.0;
2115  for (j = i; j <= n[d]; j++)
2116  ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j];
2117  }
2118  ncdiff = 0.0;
2119  for (i = 0; i <= n[0]; i++)
2120  {
2121  temp = ivr->c[i] - ivr->c[idx[d] + i];
2122  ncdiff += temp * temp;
2123  }
2124  for (i = n[0] + 1; i <= n[d]; i++)
2125  {
2126  temp = ivr->c[idx[d] + i];
2127  ncdiff += temp * temp;
2128  }
2129  ncdiff = sqrt (ncdiff);
2130  ivr->err = ncdiff * h;
2131  /* Check for divergence. */
2132  ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2133  && ivr->c[0] / iv->c[0] > 2);
2134  if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth)
2135  {
2136  igral = gnulib::copysign (octave_Inf, igral);
2137  warning ("quadcc: divergent integral detected");
2138  break;
2139  }
2140 
2141  /* Compute the local integral. */
2142  ivr->igral = h * w * ivr->c[0];
2143 
2144 
2145  /* Fix-up the heap: we now have one interval on top
2146  that we don't need any more and two new, unsorted
2147  ones at the bottom. */
2148  /* Flip the last interval to the top of the heap and
2149  sift down. */
2150  t = heap[nivals - 1];
2151  heap[nivals - 1] = heap[0];
2152  heap[0] = t;
2153  nivals--;
2154  /* Sift this interval back down the heap. */
2155  i = 0;
2156  while (2 * i + 1 < nivals - 1)
2157  {
2158  j = 2 * i + 1;
2159  if (j + 1 < nivals - 1
2160  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2161  j++;
2162  if (ivals[heap[j]].err <= ivals[heap[i]].err)
2163  break;
2164  else
2165  {
2166  t = heap[j];
2167  heap[j] = heap[i];
2168  heap[i] = t;
2169  i = j;
2170  }
2171  }
2172 
2173  /* Now grab the last interval and sift it up the heap. */
2174  i = nivals - 1;
2175  while (i > 0)
2176  {
2177  j = (i - 1) / 2;
2178  if (ivals[heap[j]].err < ivals[heap[i]].err)
2179  {
2180  t = heap[j];
2181  heap[j] = heap[i];
2182  heap[i] = t;
2183  i = j;
2184  }
2185  else
2186  break;
2187  }
2188 
2189 
2190  }
2191 
2192  /* Otherwise, just fix-up the heap. */
2193  else
2194  {
2195  i = 0;
2196  while (2 * i + 1 < nivals)
2197  {
2198  j = 2 * i + 1;
2199  if (j + 1 < nivals
2200  && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2201  j++;
2202  if (ivals[heap[j]].err <= ivals[heap[i]].err)
2203  break;
2204  else
2205  {
2206  t = heap[j];
2207  heap[j] = heap[i];
2208  heap[i] = t;
2209  i = j;
2210  }
2211  }
2212 
2213  }
2214 
2215  /* If the heap is about to overflow, remove the last two intervals. */
2216  while (nivals > cquad_heapsize - 2)
2217  {
2218  iv = &(ivals[heap[nivals - 1]]);
2219 #if (DEBUG_QUADCC)
2220  printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
2221  heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
2222 #endif
2223  err_final += iv->err;
2224  igral_final += iv->igral;
2225  nivals--;
2226  }
2227 
2228  /* Collect the value of the integral and error. */
2229  igral = igral_final;
2230  err = err_final;
2231  for (i = 0; i < nivals; i++)
2232  {
2233  igral += ivals[heap[i]].igral;
2234  err += ivals[heap[i]].err;
2235  }
2236 
2237  }
2238 
2239  /* Dump the contents of the heap. */
2240 #if (DEBUG_QUADCC)
2241  for (i = 0; i < nivals; i++)
2242  {
2243  iv = &(ivals[heap[i]]);
2244  printf ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n",
2245  i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
2246  iv->rdepth, iv->ndiv);
2247  }
2248 #endif
2249 
2250  /* Clean up and present the results. */
2251  if (nargout > 2)
2252  retval(2) = neval;
2253  if (nargout > 1)
2254  retval(1) = err;
2255  retval(0) = igral;
2256  /* All is well that ends well. */
2257  return retval;
2258 }
2259 
2260 
2261 /*
2262 %!assert (quadcc (@sin, -pi, pi), 0, 1e-6)
2263 %!assert (quadcc (inline ("sin"),- pi, pi), 0, 1e-6)
2264 %!assert (quadcc ("sin", -pi, pi), 0, 1e-6)
2265 
2266 %!assert (quadcc (@sin, -pi, 0), -2, 1e-6)
2267 %!assert (quadcc (@sin, 0, pi), 2, 1e-6)
2268 %!assert (quadcc (@(x) 1./sqrt (x), 0, 1), 2, 1e-6)
2269 %!assert (quadcc (@(x) 1./(sqrt (x).*(x+1)), 0, Inf), pi, 1e-6)
2270 
2271 %!assert (quadcc (@(x) exp (-x .^ 2), -Inf, Inf), sqrt (pi), 1e-6)
2272 %!assert (quadcc (@(x) exp (-x .^ 2), -Inf, 0), sqrt (pi)/2, 1e-6)
2273 
2274 ## Test function with NaNs in interval
2275 %!function y = __nansin (x)
2276 %! nan_locs = [-3*pi/4, -pi/4, 0, pi/3, pi/2, pi];
2277 %! y = sin (x);
2278 %! idx = min (abs (bsxfun (@minus, x(:), nan_locs)), [], 2);
2279 %! y(idx < 1e-10) = NaN;
2280 %!endfunction
2281 
2282 %!test
2283 %! [q, err, npoints] = quadcc ("__nansin", -pi, pi);
2284 %! assert (q, 0, 1e-6);
2285 %! assert (err, 0, 15*eps);
2286 
2287 %% Test input validation
2288 %!error (quadcc ())
2289 %!error (quadcc (@sin))
2290 %!error (quadcc (@sin, 0))
2291 %!error (quadcc (@sin, ones (2), pi))
2292 %!error (quadcc (@sin, -i, pi))
2293 %!error (quadcc (@sin, 0, ones (2)))
2294 %!error (quadcc (@sin, 0, i))
2295 %!error (quadcc (@sin, 0, pi, 0))
2296 %!error (quadcc (@sin, 0, pi, 1e-6, [ i ]))
2297 */
static const double Lgamma[33]
Definition: quadcc.cc:115
double err
Definition: quadcc.cc:50
static const double Tright[33 *33]
Definition: quadcc.cc:1148
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
double b
Definition: quadcc.cc:47
double fx[33]
Definition: quadcc.cc:49
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
octave_value_list feval(const std::string &name, const octave_value_list &args, int nargout)
Definition: oct-parse.cc:8625
bool xisinf(double x)
Definition: lo-mappers.cc:160
#define min_cquad_heapsize
Definition: quadcc.cc:41
int ndiv
Definition: quadcc.cc:51
F77_RET_T const double const double double * d
static const double V2inv[9 *9]
Definition: quadcc.cc:143
int depth
Definition: quadcc.cc:51
std::string unique_symbol_name(const std::string &basename)
Definition: variables.cc:495
static const double bee[68]
Definition: quadcc.cc:73
static const double Tleft[33 *33]
Definition: quadcc.cc:893
double a
Definition: quadcc.cc:47
#define octave_Inf
Definition: lo-ieee.h:31
static const double V3inv[17 *17]
Definition: quadcc.cc:178
int rdepth
Definition: quadcc.cc:51
std::complex< double > w(std::complex< double > z, double relerr=0)
double c[64]
Definition: quadcc.cc:48
Definition: dMatrix.h:35
static const double V1inv[5 *5]
Definition: quadcc.cc:130
void downdate(double *c, int n, int d, int *nans, int nnans)
Definition: quadcc.cc:1457
static const double V4inv[33 *33]
Definition: quadcc.cc:316
void warning(const char *fmt,...)
Definition: error.cc:681
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
bool xfinite(double x)
Definition: lo-mappers.cc:152
#define octave_NaN
Definition: lo-ieee.h:37
double igral
Definition: quadcc.cc:50
static const double Lalpha[33]
Definition: quadcc.cc:100
octave_function * extract_function(const octave_value &arg, const std::string &warn_for, const std::string &fname, const std::string &header, const std::string &trailer)
Definition: variables.cc:140
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
virtual octave_function * function_value(bool silent=false)
Definition: ov-base.cc:992
#define OCTAVE_QUIT
Definition: quit.h:130
Complex atan(const Complex &x)
Definition: lo-mappers.cc:231
static const double xi[33]
Definition: quadcc.cc:56
void Vinvfx(const double *fx, double *c, const int d)
Definition: quadcc.cc:1409