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