Bug Summary

File:libinterp/corefcn/quadcc.cc
Location:line 1469, column 28
Description:The left operand of '+' is a garbage value

Annotated Source Code

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