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