GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
quadcc.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2010-2025 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 "interpreter.h"
41#include "ovl.h"
42#include "utils.h"
43#include "variables.h"
44
45// Extended debugging.
46#define DEBUG_QUADCC 0
47
49
50// Define the minimum size of the interval heap.
51static const int MIN_CQUAD_HEAPSIZE = 200;
52
53// Data of a single interval.
54struct cquad_ival
55{
56 double a, b;
57 double c[64];
58 double fx[33];
59 double igral, err;
60 int depth, rdepth, ndiv;
61};
62
63// Define relative tolerance used when deciding to drop an interval.
64static constexpr 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, nargout,
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. If only a single output is requested then
1541a warning will be emitted when the requested tolerance is not met. If the
1542second output @var{err} is requested then no warning is issued and it is the
1543responsibility of the programmer to inspect and determine whether the results
1544are satisfactory.
1545
1546@code{quadcc} is capable of dealing with non-numeric values of the integrand
1547such as @code{NaN} or @code{Inf}. If the integral diverges, and @code{quadcc}
1548detects this, then a warning is issued and @code{Inf} or @code{-Inf} is
1549returned.
1550
1551Note: @code{quadcc} is a general purpose quadrature algorithm and, as such,
1552may be less efficient for a smooth or otherwise well-behaved integrand than
1553other methods such as @code{quadgk}.
1554
1555The algorithm uses @nospell{Clenshaw-Curtis} quadrature rules of increasing
1556degree in each interval and bisects the interval if either the function does
1557not appear to be smooth or a rule of maximum degree has been reached. The
1558error estimate is computed from the L2-norm of the difference between two
1559successive interpolations of the integrand over the nodes of the respective
1560quadrature rules.
1561
1562Reference: @nospell{P. Gonnet}, @cite{Increasing the Reliability of Adaptive
1563Quadrature Using Explicit Interpolants}, @nospell{ACM} Transactions on
1564Mathematical Software, Vol.@: 37, Issue 3, Article No.@: 3, 2010.
1565@seealso{quad, quadv, quadl, quadgk, trapz, dblquad, triplequad}
1566@end deftypefn */)
1567{
1568 // Some constants that we will need.
1569 static const int n[4] = { 4, 8, 16, 32 };
1570 static const int skip[4] = { 8, 4, 2, 1 };
1571 static const int idx[4] = { 0, 5, 14, 31 };
1572 static const double w = M_SQRT2 / 2;
1573 static const int ndiv_max = 20;
1574
1575 // Arguments left and right.
1576 int nargin = args.length ();
1577 octave_value fcn;
1578 double a, b, abstol, reltol, *sing;
1579 bool issingle;
1580
1581 // Variables needed for transforming the integrand.
1582 bool wrap = false;
1583 double xw;
1584
1585 // Stuff we will need to call the integrand.
1586 octave_value_list fargs, fvals;
1587
1588 // Actual variables (as opposed to constants above).
1589 double m, h, ml, hl, mr, hr, temp;
1590 double igral, err, igral_final, err_final;
1591 int nivals;
1592 int neval = 0;
1593 int i, j, d, split;
1594 int nnans, nans[33];
1595 cquad_ival *iv, *ivl, *ivr;
1596 double nc, ncdiff;
1597
1598 // Parse the input arguments.
1599 if (nargin < 3)
1600 print_usage ();
1601
1602 fcn = get_function_handle (interp, args(0), "x");
1603
1604 if (! args(1).is_real_scalar ())
1605 error ("quadcc: lower limit of integration (A) must be a real scalar");
1606 a = args(1).double_value ();
1607 issingle = args(1).is_single_type ();
1608
1609 if (! args(2).is_real_scalar ())
1610 error ("quadcc: upper limit of integration (B) must be a real scalar");
1611 b = args(2).double_value ();
1612 issingle = (issingle || args(2).is_single_type ());
1613
1614 if (nargin < 4 || args(3).isempty ())
1615 {
1616 if (issingle)
1617 {
1618 abstol = 1.0e-5;
1619 reltol = 1.0e-4;
1620 }
1621 else
1622 {
1623 abstol = 1.0e-10;
1624 reltol = 1.0e-6;
1625 }
1626 }
1627 else if (! args(3).isnumeric () || args(3).numel () > 2)
1628 error ("quadcc: TOL must be a 1- or 2-element vector [AbsTol, RelTol]");
1629 else
1630 {
1631 NDArray tol = args(3).array_value ();
1632
1633 abstol = tol(0);
1634 if (abstol < 0)
1635 error ("quadcc: absolute tolerance must be >=0");
1636
1637 if (tol.numel () == 1)
1638 {
1639 if (issingle)
1640 reltol = 1.0e-4;
1641 else
1642 reltol = 1.0e-6;
1643 }
1644 else
1645 {
1646 reltol = tol(1);
1647 if (reltol < 0)
1648 error ("quadcc: relative tolerance must be >=0");
1649 }
1650 }
1651
1652 if (nargin < 5)
1653 nivals = 1;
1654 else if (! (args(4).is_real_scalar () || args(4).is_real_matrix ()))
1655 error ("quadcc: list of singularities (SING) must be a vector of real values");
1656 else
1657 nivals = 1 + args(4).numel ();
1658
1659 int cquad_heapsize = (nivals >= MIN_CQUAD_HEAPSIZE ? nivals + 1
1660 : MIN_CQUAD_HEAPSIZE);
1661 // The interval heap.
1662 OCTAVE_LOCAL_BUFFER (cquad_ival, ivals, cquad_heapsize);
1663 OCTAVE_LOCAL_BUFFER (double, iivals, cquad_heapsize);
1664 OCTAVE_LOCAL_BUFFER (int, heap, cquad_heapsize);
1665
1666 if (nivals == 1)
1667 {
1668 iivals[0] = a;
1669 iivals[1] = b;
1670 }
1671 else
1672 {
1673 // Intervals around singularities.
1674 NDArray arr = args(4).array_value ();
1675 sing = arr.rwdata ();
1676 iivals[0] = a;
1677 std::copy_n (sing, nivals-1, iivals+1);
1678 iivals[nivals] = b;
1679 }
1680
1681 // If a or b are +/-Inf, transform the integral.
1682 if (math::isinf (a) || math::isinf (b))
1683 {
1684 wrap = true;
1685 for (i = 0; i < nivals + 1; i++)
1686 if (math::isinf (iivals[i]))
1687 iivals[i] = std::copysign (1.0, iivals[i]);
1688 else
1689 iivals[i] = 2.0 * atan (iivals[i]) / M_PI;
1690 }
1691
1692 // Initialize the heaps.
1693 for (i = 0; i < cquad_heapsize; i++)
1694 heap[i] = i;
1695
1696 // Create the first interval(s).
1697 igral = 0.0;
1698 err = 0.0;
1699 for (j = 0; j < nivals; j++)
1700 {
1701 // Initialize the interval.
1702 iv = &(ivals[heap[j]]);
1703 m = (iivals[j] + iivals[j + 1]) / 2;
1704 h = (iivals[j + 1] - iivals[j]) / 2;
1705 nnans = 0;
1706 ColumnVector ex (33);
1707 if (wrap)
1708 {
1709 for (i = 0; i <= n[3]; i++)
1710 ex(i) = tan (M_PI/2 * (m + xi[i]*h));
1711 }
1712 else
1713 {
1714 for (i = 0; i <= n[3]; i++)
1715 ex(i) = m + xi[i]*h;
1716 }
1717 fargs(0) = ex;
1718 fvals = interp.feval (fcn, fargs, 1);
1719 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1720 error ("quadcc: integrand F must return a single, real-valued vector");
1721
1722 Matrix effex = fvals(0).matrix_value ();
1723 if (effex.numel () != ex.numel ())
1724 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1725
1726 for (i = 0; i <= n[3]; i++)
1727 {
1728 iv->fx[i] = effex(i);
1729 if (wrap)
1730 {
1731 xw = ex(i);
1732 iv->fx[i] *= (1.0 + xw*xw) * M_PI/2;
1733 }
1734 neval++;
1735 if (! math::isfinite (iv->fx[i]))
1736 {
1737 nans[nnans++] = i;
1738 iv->fx[i] = 0.0;
1739 }
1740 }
1741 Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
1742 Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
1743 Vinvfx (iv->fx, &(iv->c[0]), 0);
1744 for (i = 0; i < nnans; i++)
1745 iv->fx[nans[i]] = numeric_limits<double>::NaN ();
1746 iv->a = iivals[j];
1747 iv->b = iivals[j + 1];
1748 iv->depth = 3;
1749 iv->rdepth = 1;
1750 iv->ndiv = 0;
1751 iv->igral = 2 * h * iv->c[idx[3]] * w;
1752 nc = 0.0;
1753 for (i = n[2] + 1; i <= n[3]; i++)
1754 {
1755 temp = iv->c[idx[3] + i];
1756 nc += temp * temp;
1757 }
1758 ncdiff = nc;
1759 for (i = 0; i <= n[2]; i++)
1760 {
1761 temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
1762 ncdiff += temp * temp;
1763 temp = iv->c[idx[3] + i];
1764 nc += temp * temp;
1765 }
1766 ncdiff = sqrt (ncdiff);
1767 nc = sqrt (nc);
1768 iv->err = ncdiff * 2 * h;
1769 if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
1770 iv->err = 2 * h * nc;
1771
1772 // Tabulate this interval's data.
1773 igral += iv->igral;
1774 err += iv->err;
1775
1776 // Sift it up the heap.
1777 i = j;
1778 while (i > 0 && ivals[heap[i / 2]].err < ivals[heap[i]].err)
1779 {
1780 std::swap (heap[i], heap[i / 2]);
1781 i /= 2;
1782 }
1783 }
1784
1785 // Initialize some global values.
1786 igral_final = 0.0;
1787 err_final = 0.0;
1788
1789 // Main loop.
1790 while (nivals > 0
1791 && err > std::max (abstol, fabs (igral) * reltol)
1792 && ! (err_final > std::max (abstol, fabs (igral) * reltol)
1793 && err - err_final < std::max (abstol, fabs (igral) * reltol)))
1794 {
1795 // Allow the user to interrupt.
1796 octave_quit ();
1797
1798 // Put our finger on the interval with the largest error.
1799 iv = &(ivals[heap[0]]);
1800 m = (iv->a + iv->b) / 2;
1801 h = (iv->b - iv->a) / 2;
1802
1803#if (DEBUG_QUADCC)
1804 printf ("quadcc: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1805 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1806#endif
1807
1808 // Should we try to increase the degree?
1809 if (iv->depth < 3)
1810 {
1811 // Keep tabs on some variables.
1812 d = ++iv->depth;
1813
1814 // Get the new (missing) function values.
1815 {
1816 ColumnVector ex (n[d] / 2);
1817 if (wrap)
1818 {
1819 for (i = 0; i < n[d] / 2; i++)
1820 ex(i) = tan (M_PI/2 * (m + xi[(2*i + 1) * skip[d]] * h));
1821 }
1822 else
1823 {
1824 for (i = 0; i < n[d] / 2; i++)
1825 ex(i) = m + xi[(2*i + 1) * skip[d]] * h;
1826 }
1827 fargs(0) = ex;
1828 fvals = interp.feval (fcn, fargs, 1);
1829 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1830 error ("quadcc: integrand F must return a single, real-valued vector");
1831
1832 Matrix effex = fvals(0).matrix_value ();
1833 if (effex.numel () != ex.numel ())
1834 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1835
1836 neval += effex.numel ();
1837 for (i = 0; i < n[d] / 2; i++)
1838 {
1839 j = (2*i + 1) * skip[d];
1840 iv->fx[j] = effex(i);
1841 if (wrap)
1842 {
1843 xw = ex(i);
1844 iv->fx[j] *= (1.0 + xw*xw) * M_PI/2;
1845 }
1846 }
1847 }
1848 nnans = 0;
1849 for (i = 0; i <= 32; i += skip[d])
1850 {
1851 if (! math::isfinite (iv->fx[i]))
1852 {
1853 nans[nnans++] = i;
1854 iv->fx[i] = 0.0;
1855 }
1856 }
1857
1858 // Compute the new coefficients.
1859 Vinvfx (iv->fx, &(iv->c[idx[d]]), d);
1860 // Downdate any NaNs.
1861 if (nnans > 0)
1862 {
1863 downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
1864 for (i = 0; i < nnans; i++)
1865 iv->fx[nans[i]] = numeric_limits<double>::NaN ();
1866 }
1867
1868 // Compute the error estimate.
1869 nc = 0.0;
1870 for (i = n[d - 1] + 1; i <= n[d]; i++)
1871 {
1872 temp = iv->c[idx[d] + i];
1873 nc += temp * temp;
1874 }
1875 ncdiff = nc;
1876 for (i = 0; i <= n[d - 1]; i++)
1877 {
1878 temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
1879 ncdiff += temp * temp;
1880 temp = iv->c[idx[d] + i];
1881 nc += temp * temp;
1882 }
1883 ncdiff = sqrt (ncdiff);
1884 nc = sqrt (nc);
1885 iv->err = ncdiff * 2 * h;
1886 // Compute the local integral.
1887 iv->igral = 2 * h * w * iv->c[idx[d]];
1888 // Split the interval prematurely?
1889 split = (nc > 0 && ncdiff / nc > 0.1);
1890 }
1891 else
1892 {
1893 // Maximum degree reached, just split.
1894 split = 1;
1895 }
1896
1897 // Should we drop this interval?
1898 if ((m + h*xi[0]) >= (m + h*xi[1])
1899 || (m + h*xi[31]) >= (m + h*xi[32])
1900 || iv->err < fabs (iv->igral) * DROP_RELTOL)
1901 {
1902#if (DEBUG_QUADCC)
1903 printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
1904 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
1905#endif
1906
1907 // Keep this interval's contribution.
1908 err_final += iv->err;
1909 igral_final += iv->igral;
1910 // Swap with the last element on the heap.
1911 std::swap (heap[nivals - 1], heap[0]);
1912 nivals--;
1913 // Fix up the heap.
1914 i = 0;
1915 while (2*i + 1 < nivals)
1916 {
1917 // Get the kids.
1918 j = 2*i + 1;
1919 // If the j+1st entry exists and is larger than the jth,
1920 // use it instead.
1921 if (j + 1 < nivals
1922 && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
1923 j++;
1924 // Do we need to move the ith entry up?
1925 if (ivals[heap[j]].err <= ivals[heap[i]].err)
1926 break;
1927 else
1928 {
1929 std::swap (heap[j], heap[i]);
1930 i = j;
1931 }
1932 }
1933 }
1934 else if (split)
1935 {
1936 // Some values we will need often...
1937 d = iv->depth;
1938 // Generate the interval on the left.
1939 ivl = &(ivals[heap[nivals++]]);
1940 ivl->a = iv->a;
1941 ivl->b = m;
1942 ml = (ivl->a + ivl->b) / 2;
1943 hl = h / 2;
1944 ivl->depth = 0;
1945 ivl->rdepth = iv->rdepth + 1;
1946 ivl->fx[0] = iv->fx[0];
1947 ivl->fx[32] = iv->fx[16];
1948 {
1949 ColumnVector ex (n[0] - 1);
1950 if (wrap)
1951 {
1952 for (i = 0; i < n[0] - 1; i++)
1953 ex(i) = tan (M_PI/2 * (ml + xi[(i + 1) * skip[0]] * hl));
1954 }
1955 else
1956 {
1957 for (i = 0; i < n[0] - 1; i++)
1958 ex(i) = ml + xi[(i + 1) * skip[0]] * hl;
1959 }
1960 fargs(0) = ex;
1961 fvals = interp.feval (fcn, fargs, 1);
1962 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
1963 error ("quadcc: integrand F must return a single, real-valued vector");
1964
1965 Matrix effex = fvals(0).matrix_value ();
1966 if (effex.numel () != ex.numel ())
1967 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
1968
1969 neval += effex.numel ();
1970 for (i = 0; i < n[0] - 1; i++)
1971 {
1972 j = (i + 1) * skip[0];
1973 ivl->fx[j] = effex(i);
1974 if (wrap)
1975 {
1976 xw = ex(i);
1977 ivl->fx[j] *= (1.0 + xw*xw) * M_PI/2;
1978 }
1979 }
1980 }
1981 nnans = 0;
1982 for (i = 0; i <= 32; i += skip[0])
1983 {
1984 if (! math::isfinite (ivl->fx[i]))
1985 {
1986 nans[nnans++] = i;
1987 ivl->fx[i] = 0.0;
1988 }
1989 }
1990 Vinvfx (ivl->fx, ivl->c, 0);
1991 if (nnans > 0)
1992 {
1993 downdate (ivl->c, n[0], 0, nans, nnans);
1994 for (i = 0; i < nnans; i++)
1995 ivl->fx[nans[i]] = numeric_limits<double>::NaN ();
1996 }
1997 for (i = 0; i <= n[d]; i++)
1998 {
1999 ivl->c[idx[d] + i] = 0.0;
2000 for (j = i; j <= n[d]; j++)
2001 ivl->c[idx[d] + i] += Tleft[i*33 + j] * iv->c[idx[d] + j];
2002 }
2003 ncdiff = 0.0;
2004 for (i = 0; i <= n[0]; i++)
2005 {
2006 temp = ivl->c[i] - ivl->c[idx[d] + i];
2007 ncdiff += temp * temp;
2008 }
2009 for (i = n[0] + 1; i <= n[d]; i++)
2010 {
2011 temp = ivl->c[idx[d] + i];
2012 ncdiff += temp * temp;
2013 }
2014 ncdiff = sqrt (ncdiff);
2015 ivl->err = ncdiff * h;
2016 // Check for divergence.
2017 ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2018 && ivl->c[0] / iv->c[0] > 2);
2019 if (ivl->ndiv > ndiv_max && 2*ivl->ndiv > ivl->rdepth)
2020 {
2021 igral = std::copysign (numeric_limits<double>::Inf (), igral);
2022 warning ("quadcc: divergent integral detected");
2023 break;
2024 }
2025
2026 // Compute the local integral.
2027 ivl->igral = h * w * ivl->c[0];
2028
2029 // Generate the interval on the right.
2030 ivr = &(ivals[heap[nivals++]]);
2031 ivr->a = m;
2032 ivr->b = iv->b;
2033 mr = (ivr->a + ivr->b) / 2;
2034 hr = h / 2;
2035 ivr->depth = 0;
2036 ivr->rdepth = iv->rdepth + 1;
2037 ivr->fx[0] = iv->fx[16];
2038 ivr->fx[32] = iv->fx[32];
2039 {
2040 ColumnVector ex (n[0] - 1);
2041 if (wrap)
2042 {
2043 for (i = 0; i < n[0] - 1; i++)
2044 ex(i) = tan (M_PI/2 * (mr + xi[(i + 1) * skip[0]] * hr));
2045 }
2046 else
2047 {
2048 for (i = 0; i < n[0] - 1; i++)
2049 ex(i) = mr + xi[(i + 1) * skip[0]] * hr;
2050 }
2051 fargs(0) = ex;
2052 fvals = interp.feval (fcn, fargs, 1);
2053 if (fvals.length () != 1 || ! fvals(0).is_real_matrix ())
2054 error ("quadcc: integrand F must return a single, real-valued vector");
2055
2056 Matrix effex = fvals(0).matrix_value ();
2057 if (effex.numel () != ex.numel ())
2058 error ("quadcc: integrand F must return a single, real-valued vector of the same size as the input");
2059
2060 neval += effex.numel ();
2061 for (i = 0; i < n[0] - 1; i++)
2062 {
2063 j = (i + 1) * skip[0];
2064 ivr->fx[j] = effex(i);
2065 if (wrap)
2066 {
2067 xw = ex(i);
2068 ivr->fx[j] *= (1.0 + xw*xw) * M_PI/2;
2069 }
2070 }
2071 }
2072 nnans = 0;
2073 for (i = 0; i <= 32; i += skip[0])
2074 {
2075 if (! math::isfinite (ivr->fx[i]))
2076 {
2077 nans[nnans++] = i;
2078 ivr->fx[i] = 0.0;
2079 }
2080 }
2081 Vinvfx (ivr->fx, ivr->c, 0);
2082 if (nnans > 0)
2083 {
2084 downdate (ivr->c, n[0], 0, nans, nnans);
2085 for (i = 0; i < nnans; i++)
2086 ivr->fx[nans[i]] = numeric_limits<double>::NaN ();
2087 }
2088 for (i = 0; i <= n[d]; i++)
2089 {
2090 ivr->c[idx[d] + i] = 0.0;
2091 for (j = i; j <= n[d]; j++)
2092 ivr->c[idx[d] + i] += Tright[i*33 + j] * iv->c[idx[d] + j];
2093 }
2094 ncdiff = 0.0;
2095 for (i = 0; i <= n[0]; i++)
2096 {
2097 temp = ivr->c[i] - ivr->c[idx[d] + i];
2098 ncdiff += temp * temp;
2099 }
2100 for (i = n[0] + 1; i <= n[d]; i++)
2101 {
2102 temp = ivr->c[idx[d] + i];
2103 ncdiff += temp * temp;
2104 }
2105 ncdiff = sqrt (ncdiff);
2106 ivr->err = ncdiff * h;
2107 // Check for divergence.
2108 ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
2109 && ivr->c[0] / iv->c[0] > 2);
2110 if (ivr->ndiv > ndiv_max && 2*ivr->ndiv > ivr->rdepth)
2111 {
2112 igral = std::copysign (numeric_limits<double>::Inf (), igral);
2113 warning ("quadcc: divergent integral detected");
2114 break;
2115 }
2116
2117 // Compute the local integral.
2118 ivr->igral = h * w * ivr->c[0];
2119
2120 // Fix-up the heap: we now have one interval on top that we
2121 // don't need any more and two new, unsorted ones at the bottom.
2122
2123 // Flip the last interval to the top of the heap and sift down.
2124 std::swap (heap[nivals - 1], heap[0]);
2125 nivals--;
2126 // Sift this interval back down the heap.
2127 i = 0;
2128 while (2*i + 1 < nivals - 1)
2129 {
2130 j = 2*i + 1;
2131 if (j + 1 < nivals - 1
2132 && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2133 j++;
2134 if (ivals[heap[j]].err <= ivals[heap[i]].err)
2135 break;
2136 else
2137 {
2138 std::swap (heap[j], heap[i]);
2139 i = j;
2140 }
2141 }
2142
2143 // Now grab the last interval and sift it up the heap.
2144 i = nivals - 1;
2145 while (i > 0)
2146 {
2147 j = (i - 1) / 2;
2148 if (ivals[heap[j]].err < ivals[heap[i]].err)
2149 {
2150 std::swap (heap[j], heap[i]);
2151 i = j;
2152 }
2153 else
2154 break;
2155 }
2156 }
2157 else
2158 {
2159 // Otherwise, just fix-up the heap.
2160 i = 0;
2161 while (2*i + 1 < nivals)
2162 {
2163 j = 2*i + 1;
2164 if (j + 1 < nivals
2165 && ivals[heap[j + 1]].err >= ivals[heap[j]].err)
2166 j++;
2167 if (ivals[heap[j]].err <= ivals[heap[i]].err)
2168 break;
2169 else
2170 {
2171 std::swap (heap[j], heap[i]);
2172 i = j;
2173 }
2174 }
2175 }
2176
2177 // If the heap is about to overflow, remove the last two intervals.
2178 while (nivals > cquad_heapsize - 2)
2179 {
2180 iv = &(ivals[heap[nivals - 1]]);
2181#if (DEBUG_QUADCC)
2182 printf ("quadcc: dropping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
2183 heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
2184#endif
2185 err_final += iv->err;
2186 igral_final += iv->igral;
2187 nivals--;
2188 }
2189
2190 // Collect the value of the integral and error.
2191 igral = igral_final;
2192 err = err_final;
2193 for (i = 0; i < nivals; i++)
2194 {
2195 igral += ivals[heap[i]].igral;
2196 err += ivals[heap[i]].err;
2197 }
2198 }
2199
2200#if (DEBUG_QUADCC)
2201 // Dump the contents of the heap.
2202 for (i = 0; i < nivals; i++)
2203 {
2204 iv = &(ivals[heap[i]]);
2205 printf ("quadcc: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i, ndiv=%i\n",
2206 i, heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
2207 iv->rdepth, iv->ndiv);
2208 }
2209#endif
2210
2211 if (nargout < 2 && err > std::max (abstol, reltol * std::abs (igral)))
2212 warning ("quadcc: Error tolerance not met. Estimated error: %g\n", err);
2213
2214 if (issingle)
2215 return ovl (static_cast<float> (igral), err, neval);
2216 else
2217 return ovl (igral, err, neval);
2218}
2219
2220/*
2221%!assert (quadcc (@sin, -pi, pi), 0, 1e-10)
2222%!test
2223%! warning ("off", "Octave:legacy-function", "local");
2224%! assert (quadcc (inline ("sin"), -pi, pi), 0, 1e-10);
2225%!assert (quadcc ("sin", -pi, pi), 0, 1e-10)
2226
2227%!assert (quadcc (@sin, -pi, 0), -2, 1e-10)
2228%!assert (quadcc (@sin, 0, pi), 2, 1e-10)
2229%!assert (quadcc (@(x) 1./sqrt (x), 0, 1), 2, -1e-6)
2230%!assert (quadcc (@(x) 1./(sqrt (x).*(x+1)), 0, Inf), pi, -1e-6)
2231%!assert (quadcc (@(x) 1./(sqrt (x).*(x+1)), 0, Inf, [0, 2e-7]), pi, -2e-7)
2232
2233%!assert (quadcc (@(x) exp (-x .^ 2), -Inf, Inf), sqrt (pi), 1e-10)
2234%!assert (quadcc (@(x) exp (-x .^ 2), -Inf, 0), sqrt (pi)/2, 1e-10)
2235
2236## Test function with NaNs in interval
2237%!function y = __nansin (x)
2238%! nan_locs = [-3*pi/4, -pi/4, 0, pi/3, pi/2, pi];
2239%! y = sin (x);
2240%! idx = min (abs (bsxfun (@minus, x(:), nan_locs)), [], 2);
2241%! y(idx < 1e-10) = NaN;
2242%!endfunction
2243
2244%!test
2245%! [q, err, npoints] = quadcc ("__nansin", -pi, pi, [0, 1e-6]);
2246%! assert (q, 0, -1e-6);
2247%! assert (err, 0, 15*eps);
2248
2249%!test
2250%! assert (class (quadcc (@sin, 0, 1)), "double");
2251%! assert (class (quadcc (@sin, single (0), 1)), "single");
2252%! assert (class (quadcc (@sin, 0, single (1))), "single");
2253%! assert (class (quadcc (@sin, single (0), single (1))), "single");
2254
2255%!test <*62412>
2256%! f = @(t) -1 ./ t.^1.1;
2257%! fail ("quadcc (f, 1, Inf)", "warning", "Error tolerance not met");
2258%! [q, err] = quadcc (f, 1, Inf);
2259%! assert (err > 1e-5);
2260
2261## Test input validation
2262%!error quadcc ()
2263%!error quadcc (@sin)
2264%!error quadcc (@sin, 0)
2265%!error <lower limit .* must be a .* scalar> (quadcc (@sin, ones (2), pi))
2266%!error <lower limit .* must be a real scalar> (quadcc (@sin, -i, pi))
2267%!error <upper limit .* must be a .* scalar> (quadcc (@sin, 0, ones (2)))
2268%!error <upper limit .* must be a real scalar> (quadcc (@sin, 0, i))
2269%!error <TOL must be a 1- or 2-element vector> (quadcc (@sin, 0, pi, {1}))
2270%!error <TOL must be a 1- or 2-element vector> (quadcc (@sin, 0, pi, [1,2,3]))
2271%!error <absolute tolerance must be .=0> (quadcc (@sin, 0, pi, -1))
2272%!error <relative tolerance must be .=0> (quadcc (@sin, 0, pi, [1, -1]))
2273%!error <SING.* must be .* real values> (quadcc (@sin, 0, pi, 1e-6, [ i ]))
2274*/
2275
2276OCTAVE_END_NAMESPACE(octave)
T * rwdata()
Size of the specified dimension.
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
octave_idx_type length() const
Definition ovl.h:111
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
void print_usage()
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:1078
void error(const char *fmt,...)
Definition error.cc:1003
octave_value get_function_handle(interpreter &interp, const octave_value &arg, const std::string &parameter_name)
Complex atan(const Complex &x)
Definition lo-mappers.h:71
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition ovl.h:217