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