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