Line data Source code
1 : /*************************************************************************
2 : ALGLIB 3.17.0 (source code generated 2020-12-27)
3 : Copyright (c) Sergey Bochkanov (ALGLIB project).
4 :
5 : >>> SOURCE LICENSE >>>
6 : This program is free software; you can redistribute it and/or modify
7 : it under the terms of the GNU General Public License as published by
8 : the Free Software Foundation (www.fsf.org); either version 2 of the
9 : License, or (at your option) any later version.
10 :
11 : This program is distributed in the hope that it will be useful,
12 : but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : GNU General Public License for more details.
15 :
16 : A copy of the GNU General Public License is available at
17 : http://www.fsf.org/licensing/licenses
18 : >>> END OF LICENSE >>>
19 : *************************************************************************/
20 : #ifdef _MSC_VER
21 : #define _CRT_SECURE_NO_WARNINGS
22 : #endif
23 : #include "stdafx.h"
24 : #include "integration.h"
25 :
26 : // disable some irrelevant warnings
27 : #if (AE_COMPILER==AE_MSVC) && !defined(AE_ALL_WARNINGS)
28 : #pragma warning(disable:4100)
29 : #pragma warning(disable:4127)
30 : #pragma warning(disable:4611)
31 : #pragma warning(disable:4702)
32 : #pragma warning(disable:4996)
33 : #endif
34 :
35 : /////////////////////////////////////////////////////////////////////////
36 : //
37 : // THIS SECTION CONTAINS IMPLEMENTATION OF C++ INTERFACE
38 : //
39 : /////////////////////////////////////////////////////////////////////////
40 : namespace alglib
41 : {
42 :
43 : #if defined(AE_COMPILE_GQ) || !defined(AE_PARTIAL_BUILD)
44 :
45 : #endif
46 :
47 : #if defined(AE_COMPILE_GKQ) || !defined(AE_PARTIAL_BUILD)
48 :
49 : #endif
50 :
51 : #if defined(AE_COMPILE_AUTOGK) || !defined(AE_PARTIAL_BUILD)
52 :
53 : #endif
54 :
55 : #if defined(AE_COMPILE_GQ) || !defined(AE_PARTIAL_BUILD)
56 : /*************************************************************************
57 : Computation of nodes and weights for a Gauss quadrature formula
58 :
59 : The algorithm generates the N-point Gauss quadrature formula with weight
60 : function given by coefficients alpha and beta of a recurrence relation
61 : which generates a system of orthogonal polynomials:
62 :
63 : P-1(x) = 0
64 : P0(x) = 1
65 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
66 :
67 : and zeroth moment Mu0
68 :
69 : Mu0 = integral(W(x)dx,a,b)
70 :
71 : INPUT PARAMETERS:
72 : Alpha - array[0..N-1], alpha coefficients
73 : Beta - array[0..N-1], beta coefficients
74 : Zero-indexed element is not used and may be arbitrary.
75 : Beta[I]>0.
76 : Mu0 - zeroth moment of the weight function.
77 : N - number of nodes of the quadrature formula, N>=1
78 :
79 : OUTPUT PARAMETERS:
80 : Info - error code:
81 : * -3 internal eigenproblem solver hasn't converged
82 : * -2 Beta[i]<=0
83 : * -1 incorrect N was passed
84 : * 1 OK
85 : X - array[0..N-1] - array of quadrature nodes,
86 : in ascending order.
87 : W - array[0..N-1] - array of quadrature weights.
88 :
89 : -- ALGLIB --
90 : Copyright 2005-2009 by Bochkanov Sergey
91 : *************************************************************************/
92 0 : void gqgeneraterec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
93 : {
94 : jmp_buf _break_jump;
95 : alglib_impl::ae_state _alglib_env_state;
96 0 : alglib_impl::ae_state_init(&_alglib_env_state);
97 0 : if( setjmp(_break_jump) )
98 : {
99 : #if !defined(AE_NO_EXCEPTIONS)
100 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
101 : #else
102 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
103 : return;
104 : #endif
105 : }
106 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
107 0 : if( _xparams.flags!=0x0 )
108 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
109 0 : alglib_impl::gqgeneraterec(const_cast<alglib_impl::ae_vector*>(alpha.c_ptr()), const_cast<alglib_impl::ae_vector*>(beta.c_ptr()), mu0, n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
110 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
111 0 : return;
112 : }
113 :
114 : /*************************************************************************
115 : Computation of nodes and weights for a Gauss-Lobatto quadrature formula
116 :
117 : The algorithm generates the N-point Gauss-Lobatto quadrature formula with
118 : weight function given by coefficients alpha and beta of a recurrence which
119 : generates a system of orthogonal polynomials.
120 :
121 : P-1(x) = 0
122 : P0(x) = 1
123 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
124 :
125 : and zeroth moment Mu0
126 :
127 : Mu0 = integral(W(x)dx,a,b)
128 :
129 : INPUT PARAMETERS:
130 : Alpha - array[0..N-2], alpha coefficients
131 : Beta - array[0..N-2], beta coefficients.
132 : Zero-indexed element is not used, may be arbitrary.
133 : Beta[I]>0
134 : Mu0 - zeroth moment of the weighting function.
135 : A - left boundary of the integration interval.
136 : B - right boundary of the integration interval.
137 : N - number of nodes of the quadrature formula, N>=3
138 : (including the left and right boundary nodes).
139 :
140 : OUTPUT PARAMETERS:
141 : Info - error code:
142 : * -3 internal eigenproblem solver hasn't converged
143 : * -2 Beta[i]<=0
144 : * -1 incorrect N was passed
145 : * 1 OK
146 : X - array[0..N-1] - array of quadrature nodes,
147 : in ascending order.
148 : W - array[0..N-1] - array of quadrature weights.
149 :
150 : -- ALGLIB --
151 : Copyright 2005-2009 by Bochkanov Sergey
152 : *************************************************************************/
153 0 : void gqgenerategausslobattorec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const double b, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
154 : {
155 : jmp_buf _break_jump;
156 : alglib_impl::ae_state _alglib_env_state;
157 0 : alglib_impl::ae_state_init(&_alglib_env_state);
158 0 : if( setjmp(_break_jump) )
159 : {
160 : #if !defined(AE_NO_EXCEPTIONS)
161 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
162 : #else
163 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
164 : return;
165 : #endif
166 : }
167 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
168 0 : if( _xparams.flags!=0x0 )
169 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
170 0 : alglib_impl::gqgenerategausslobattorec(const_cast<alglib_impl::ae_vector*>(alpha.c_ptr()), const_cast<alglib_impl::ae_vector*>(beta.c_ptr()), mu0, a, b, n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
171 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
172 0 : return;
173 : }
174 :
175 : /*************************************************************************
176 : Computation of nodes and weights for a Gauss-Radau quadrature formula
177 :
178 : The algorithm generates the N-point Gauss-Radau quadrature formula with
179 : weight function given by the coefficients alpha and beta of a recurrence
180 : which generates a system of orthogonal polynomials.
181 :
182 : P-1(x) = 0
183 : P0(x) = 1
184 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
185 :
186 : and zeroth moment Mu0
187 :
188 : Mu0 = integral(W(x)dx,a,b)
189 :
190 : INPUT PARAMETERS:
191 : Alpha - array[0..N-2], alpha coefficients.
192 : Beta - array[0..N-1], beta coefficients
193 : Zero-indexed element is not used.
194 : Beta[I]>0
195 : Mu0 - zeroth moment of the weighting function.
196 : A - left boundary of the integration interval.
197 : N - number of nodes of the quadrature formula, N>=2
198 : (including the left boundary node).
199 :
200 : OUTPUT PARAMETERS:
201 : Info - error code:
202 : * -3 internal eigenproblem solver hasn't converged
203 : * -2 Beta[i]<=0
204 : * -1 incorrect N was passed
205 : * 1 OK
206 : X - array[0..N-1] - array of quadrature nodes,
207 : in ascending order.
208 : W - array[0..N-1] - array of quadrature weights.
209 :
210 :
211 : -- ALGLIB --
212 : Copyright 2005-2009 by Bochkanov Sergey
213 : *************************************************************************/
214 0 : void gqgenerategaussradaurec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const double a, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
215 : {
216 : jmp_buf _break_jump;
217 : alglib_impl::ae_state _alglib_env_state;
218 0 : alglib_impl::ae_state_init(&_alglib_env_state);
219 0 : if( setjmp(_break_jump) )
220 : {
221 : #if !defined(AE_NO_EXCEPTIONS)
222 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
223 : #else
224 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
225 : return;
226 : #endif
227 : }
228 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
229 0 : if( _xparams.flags!=0x0 )
230 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
231 0 : alglib_impl::gqgenerategaussradaurec(const_cast<alglib_impl::ae_vector*>(alpha.c_ptr()), const_cast<alglib_impl::ae_vector*>(beta.c_ptr()), mu0, a, n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
232 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
233 0 : return;
234 : }
235 :
236 : /*************************************************************************
237 : Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
238 : nodes.
239 :
240 : INPUT PARAMETERS:
241 : N - number of nodes, >=1
242 :
243 : OUTPUT PARAMETERS:
244 : Info - error code:
245 : * -4 an error was detected when calculating
246 : weights/nodes. N is too large to obtain
247 : weights/nodes with high enough accuracy.
248 : Try to use multiple precision version.
249 : * -3 internal eigenproblem solver hasn't converged
250 : * -1 incorrect N was passed
251 : * +1 OK
252 : X - array[0..N-1] - array of quadrature nodes,
253 : in ascending order.
254 : W - array[0..N-1] - array of quadrature weights.
255 :
256 :
257 : -- ALGLIB --
258 : Copyright 12.05.2009 by Bochkanov Sergey
259 : *************************************************************************/
260 0 : void gqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
261 : {
262 : jmp_buf _break_jump;
263 : alglib_impl::ae_state _alglib_env_state;
264 0 : alglib_impl::ae_state_init(&_alglib_env_state);
265 0 : if( setjmp(_break_jump) )
266 : {
267 : #if !defined(AE_NO_EXCEPTIONS)
268 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
269 : #else
270 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
271 : return;
272 : #endif
273 : }
274 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
275 0 : if( _xparams.flags!=0x0 )
276 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
277 0 : alglib_impl::gqgenerategausslegendre(n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
278 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
279 0 : return;
280 : }
281 :
282 : /*************************************************************************
283 : Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] with weight
284 : function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
285 :
286 : INPUT PARAMETERS:
287 : N - number of nodes, >=1
288 : Alpha - power-law coefficient, Alpha>-1
289 : Beta - power-law coefficient, Beta>-1
290 :
291 : OUTPUT PARAMETERS:
292 : Info - error code:
293 : * -4 an error was detected when calculating
294 : weights/nodes. Alpha or Beta are too close
295 : to -1 to obtain weights/nodes with high enough
296 : accuracy, or, may be, N is too large. Try to
297 : use multiple precision version.
298 : * -3 internal eigenproblem solver hasn't converged
299 : * -1 incorrect N/Alpha/Beta was passed
300 : * +1 OK
301 : X - array[0..N-1] - array of quadrature nodes,
302 : in ascending order.
303 : W - array[0..N-1] - array of quadrature weights.
304 :
305 :
306 : -- ALGLIB --
307 : Copyright 12.05.2009 by Bochkanov Sergey
308 : *************************************************************************/
309 0 : void gqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
310 : {
311 : jmp_buf _break_jump;
312 : alglib_impl::ae_state _alglib_env_state;
313 0 : alglib_impl::ae_state_init(&_alglib_env_state);
314 0 : if( setjmp(_break_jump) )
315 : {
316 : #if !defined(AE_NO_EXCEPTIONS)
317 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
318 : #else
319 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
320 : return;
321 : #endif
322 : }
323 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
324 0 : if( _xparams.flags!=0x0 )
325 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
326 0 : alglib_impl::gqgenerategaussjacobi(n, alpha, beta, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
327 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
328 0 : return;
329 : }
330 :
331 : /*************************************************************************
332 : Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) with
333 : weight function W(x)=Power(x,Alpha)*Exp(-x)
334 :
335 : INPUT PARAMETERS:
336 : N - number of nodes, >=1
337 : Alpha - power-law coefficient, Alpha>-1
338 :
339 : OUTPUT PARAMETERS:
340 : Info - error code:
341 : * -4 an error was detected when calculating
342 : weights/nodes. Alpha is too close to -1 to
343 : obtain weights/nodes with high enough accuracy
344 : or, may be, N is too large. Try to use
345 : multiple precision version.
346 : * -3 internal eigenproblem solver hasn't converged
347 : * -1 incorrect N/Alpha was passed
348 : * +1 OK
349 : X - array[0..N-1] - array of quadrature nodes,
350 : in ascending order.
351 : W - array[0..N-1] - array of quadrature weights.
352 :
353 :
354 : -- ALGLIB --
355 : Copyright 12.05.2009 by Bochkanov Sergey
356 : *************************************************************************/
357 0 : void gqgenerategausslaguerre(const ae_int_t n, const double alpha, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
358 : {
359 : jmp_buf _break_jump;
360 : alglib_impl::ae_state _alglib_env_state;
361 0 : alglib_impl::ae_state_init(&_alglib_env_state);
362 0 : if( setjmp(_break_jump) )
363 : {
364 : #if !defined(AE_NO_EXCEPTIONS)
365 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
366 : #else
367 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
368 : return;
369 : #endif
370 : }
371 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
372 0 : if( _xparams.flags!=0x0 )
373 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
374 0 : alglib_impl::gqgenerategausslaguerre(n, alpha, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
375 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
376 0 : return;
377 : }
378 :
379 : /*************************************************************************
380 : Returns nodes/weights for Gauss-Hermite quadrature on (-inf,+inf) with
381 : weight function W(x)=Exp(-x*x)
382 :
383 : INPUT PARAMETERS:
384 : N - number of nodes, >=1
385 :
386 : OUTPUT PARAMETERS:
387 : Info - error code:
388 : * -4 an error was detected when calculating
389 : weights/nodes. May be, N is too large. Try to
390 : use multiple precision version.
391 : * -3 internal eigenproblem solver hasn't converged
392 : * -1 incorrect N/Alpha was passed
393 : * +1 OK
394 : X - array[0..N-1] - array of quadrature nodes,
395 : in ascending order.
396 : W - array[0..N-1] - array of quadrature weights.
397 :
398 :
399 : -- ALGLIB --
400 : Copyright 12.05.2009 by Bochkanov Sergey
401 : *************************************************************************/
402 0 : void gqgenerategausshermite(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &w, const xparams _xparams)
403 : {
404 : jmp_buf _break_jump;
405 : alglib_impl::ae_state _alglib_env_state;
406 0 : alglib_impl::ae_state_init(&_alglib_env_state);
407 0 : if( setjmp(_break_jump) )
408 : {
409 : #if !defined(AE_NO_EXCEPTIONS)
410 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
411 : #else
412 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
413 : return;
414 : #endif
415 : }
416 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
417 0 : if( _xparams.flags!=0x0 )
418 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
419 0 : alglib_impl::gqgenerategausshermite(n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(w.c_ptr()), &_alglib_env_state);
420 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
421 0 : return;
422 : }
423 : #endif
424 :
425 : #if defined(AE_COMPILE_GKQ) || !defined(AE_PARTIAL_BUILD)
426 : /*************************************************************************
427 : Computation of nodes and weights of a Gauss-Kronrod quadrature formula
428 :
429 : The algorithm generates the N-point Gauss-Kronrod quadrature formula with
430 : weight function given by coefficients alpha and beta of a recurrence
431 : relation which generates a system of orthogonal polynomials:
432 :
433 : P-1(x) = 0
434 : P0(x) = 1
435 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
436 :
437 : and zero moment Mu0
438 :
439 : Mu0 = integral(W(x)dx,a,b)
440 :
441 :
442 : INPUT PARAMETERS:
443 : Alpha - alpha coefficients, array[0..floor(3*K/2)].
444 : Beta - beta coefficients, array[0..ceil(3*K/2)].
445 : Beta[0] is not used and may be arbitrary.
446 : Beta[I]>0.
447 : Mu0 - zeroth moment of the weight function.
448 : N - number of nodes of the Gauss-Kronrod quadrature formula,
449 : N >= 3,
450 : N = 2*K+1.
451 :
452 : OUTPUT PARAMETERS:
453 : Info - error code:
454 : * -5 no real and positive Gauss-Kronrod formula can
455 : be created for such a weight function with a
456 : given number of nodes.
457 : * -4 N is too large, task may be ill conditioned -
458 : x[i]=x[i+1] found.
459 : * -3 internal eigenproblem solver hasn't converged
460 : * -2 Beta[i]<=0
461 : * -1 incorrect N was passed
462 : * +1 OK
463 : X - array[0..N-1] - array of quadrature nodes,
464 : in ascending order.
465 : WKronrod - array[0..N-1] - Kronrod weights
466 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
467 : corresponding to extended Kronrod nodes).
468 :
469 : -- ALGLIB --
470 : Copyright 08.05.2009 by Bochkanov Sergey
471 : *************************************************************************/
472 0 : void gkqgeneraterec(const real_1d_array &alpha, const real_1d_array &beta, const double mu0, const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, const xparams _xparams)
473 : {
474 : jmp_buf _break_jump;
475 : alglib_impl::ae_state _alglib_env_state;
476 0 : alglib_impl::ae_state_init(&_alglib_env_state);
477 0 : if( setjmp(_break_jump) )
478 : {
479 : #if !defined(AE_NO_EXCEPTIONS)
480 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
481 : #else
482 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
483 : return;
484 : #endif
485 : }
486 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
487 0 : if( _xparams.flags!=0x0 )
488 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
489 0 : alglib_impl::gkqgeneraterec(const_cast<alglib_impl::ae_vector*>(alpha.c_ptr()), const_cast<alglib_impl::ae_vector*>(beta.c_ptr()), mu0, n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(wkronrod.c_ptr()), const_cast<alglib_impl::ae_vector*>(wgauss.c_ptr()), &_alglib_env_state);
490 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
491 0 : return;
492 : }
493 :
494 : /*************************************************************************
495 : Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre
496 : quadrature with N points.
497 :
498 : GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed table) is
499 : used depending on machine precision and number of nodes.
500 :
501 : INPUT PARAMETERS:
502 : N - number of Kronrod nodes, must be odd number, >=3.
503 :
504 : OUTPUT PARAMETERS:
505 : Info - error code:
506 : * -4 an error was detected when calculating
507 : weights/nodes. N is too large to obtain
508 : weights/nodes with high enough accuracy.
509 : Try to use multiple precision version.
510 : * -3 internal eigenproblem solver hasn't converged
511 : * -1 incorrect N was passed
512 : * +1 OK
513 : X - array[0..N-1] - array of quadrature nodes, ordered in
514 : ascending order.
515 : WKronrod - array[0..N-1] - Kronrod weights
516 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
517 : corresponding to extended Kronrod nodes).
518 :
519 :
520 : -- ALGLIB --
521 : Copyright 12.05.2009 by Bochkanov Sergey
522 : *************************************************************************/
523 0 : void gkqgenerategausslegendre(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, const xparams _xparams)
524 : {
525 : jmp_buf _break_jump;
526 : alglib_impl::ae_state _alglib_env_state;
527 0 : alglib_impl::ae_state_init(&_alglib_env_state);
528 0 : if( setjmp(_break_jump) )
529 : {
530 : #if !defined(AE_NO_EXCEPTIONS)
531 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
532 : #else
533 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
534 : return;
535 : #endif
536 : }
537 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
538 0 : if( _xparams.flags!=0x0 )
539 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
540 0 : alglib_impl::gkqgenerategausslegendre(n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(wkronrod.c_ptr()), const_cast<alglib_impl::ae_vector*>(wgauss.c_ptr()), &_alglib_env_state);
541 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
542 0 : return;
543 : }
544 :
545 : /*************************************************************************
546 : Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi
547 : quadrature on [-1,1] with weight function
548 :
549 : W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
550 :
551 : INPUT PARAMETERS:
552 : N - number of Kronrod nodes, must be odd number, >=3.
553 : Alpha - power-law coefficient, Alpha>-1
554 : Beta - power-law coefficient, Beta>-1
555 :
556 : OUTPUT PARAMETERS:
557 : Info - error code:
558 : * -5 no real and positive Gauss-Kronrod formula can
559 : be created for such a weight function with a
560 : given number of nodes.
561 : * -4 an error was detected when calculating
562 : weights/nodes. Alpha or Beta are too close
563 : to -1 to obtain weights/nodes with high enough
564 : accuracy, or, may be, N is too large. Try to
565 : use multiple precision version.
566 : * -3 internal eigenproblem solver hasn't converged
567 : * -1 incorrect N was passed
568 : * +1 OK
569 : * +2 OK, but quadrature rule have exterior nodes,
570 : x[0]<-1 or x[n-1]>+1
571 : X - array[0..N-1] - array of quadrature nodes, ordered in
572 : ascending order.
573 : WKronrod - array[0..N-1] - Kronrod weights
574 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
575 : corresponding to extended Kronrod nodes).
576 :
577 :
578 : -- ALGLIB --
579 : Copyright 12.05.2009 by Bochkanov Sergey
580 : *************************************************************************/
581 0 : void gkqgenerategaussjacobi(const ae_int_t n, const double alpha, const double beta, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, const xparams _xparams)
582 : {
583 : jmp_buf _break_jump;
584 : alglib_impl::ae_state _alglib_env_state;
585 0 : alglib_impl::ae_state_init(&_alglib_env_state);
586 0 : if( setjmp(_break_jump) )
587 : {
588 : #if !defined(AE_NO_EXCEPTIONS)
589 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
590 : #else
591 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
592 : return;
593 : #endif
594 : }
595 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
596 0 : if( _xparams.flags!=0x0 )
597 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
598 0 : alglib_impl::gkqgenerategaussjacobi(n, alpha, beta, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(wkronrod.c_ptr()), const_cast<alglib_impl::ae_vector*>(wgauss.c_ptr()), &_alglib_env_state);
599 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
600 0 : return;
601 : }
602 :
603 : /*************************************************************************
604 : Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
605 :
606 : Reduction to tridiagonal eigenproblem is used.
607 :
608 : INPUT PARAMETERS:
609 : N - number of Kronrod nodes, must be odd number, >=3.
610 :
611 : OUTPUT PARAMETERS:
612 : Info - error code:
613 : * -4 an error was detected when calculating
614 : weights/nodes. N is too large to obtain
615 : weights/nodes with high enough accuracy.
616 : Try to use multiple precision version.
617 : * -3 internal eigenproblem solver hasn't converged
618 : * -1 incorrect N was passed
619 : * +1 OK
620 : X - array[0..N-1] - array of quadrature nodes, ordered in
621 : ascending order.
622 : WKronrod - array[0..N-1] - Kronrod weights
623 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
624 : corresponding to extended Kronrod nodes).
625 :
626 : -- ALGLIB --
627 : Copyright 12.05.2009 by Bochkanov Sergey
628 : *************************************************************************/
629 0 : void gkqlegendrecalc(const ae_int_t n, ae_int_t &info, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, const xparams _xparams)
630 : {
631 : jmp_buf _break_jump;
632 : alglib_impl::ae_state _alglib_env_state;
633 0 : alglib_impl::ae_state_init(&_alglib_env_state);
634 0 : if( setjmp(_break_jump) )
635 : {
636 : #if !defined(AE_NO_EXCEPTIONS)
637 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
638 : #else
639 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
640 : return;
641 : #endif
642 : }
643 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
644 0 : if( _xparams.flags!=0x0 )
645 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
646 0 : alglib_impl::gkqlegendrecalc(n, &info, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(wkronrod.c_ptr()), const_cast<alglib_impl::ae_vector*>(wgauss.c_ptr()), &_alglib_env_state);
647 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
648 0 : return;
649 : }
650 :
651 : /*************************************************************************
652 : Returns Gauss and Gauss-Kronrod nodes for quadrature with N points using
653 : pre-calculated table. Nodes/weights were computed with accuracy up to
654 : 1.0E-32 (if MPFR version of ALGLIB is used). In standard double precision
655 : accuracy reduces to something about 2.0E-16 (depending on your compiler's
656 : handling of long floating point constants).
657 :
658 : INPUT PARAMETERS:
659 : N - number of Kronrod nodes.
660 : N can be 15, 21, 31, 41, 51, 61.
661 :
662 : OUTPUT PARAMETERS:
663 : X - array[0..N-1] - array of quadrature nodes, ordered in
664 : ascending order.
665 : WKronrod - array[0..N-1] - Kronrod weights
666 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
667 : corresponding to extended Kronrod nodes).
668 :
669 :
670 : -- ALGLIB --
671 : Copyright 12.05.2009 by Bochkanov Sergey
672 : *************************************************************************/
673 0 : void gkqlegendretbl(const ae_int_t n, real_1d_array &x, real_1d_array &wkronrod, real_1d_array &wgauss, double &eps, const xparams _xparams)
674 : {
675 : jmp_buf _break_jump;
676 : alglib_impl::ae_state _alglib_env_state;
677 0 : alglib_impl::ae_state_init(&_alglib_env_state);
678 0 : if( setjmp(_break_jump) )
679 : {
680 : #if !defined(AE_NO_EXCEPTIONS)
681 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
682 : #else
683 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
684 : return;
685 : #endif
686 : }
687 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
688 0 : if( _xparams.flags!=0x0 )
689 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
690 0 : alglib_impl::gkqlegendretbl(n, const_cast<alglib_impl::ae_vector*>(x.c_ptr()), const_cast<alglib_impl::ae_vector*>(wkronrod.c_ptr()), const_cast<alglib_impl::ae_vector*>(wgauss.c_ptr()), &eps, &_alglib_env_state);
691 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
692 0 : return;
693 : }
694 : #endif
695 :
696 : #if defined(AE_COMPILE_AUTOGK) || !defined(AE_PARTIAL_BUILD)
697 : /*************************************************************************
698 : Integration report:
699 : * TerminationType = completetion code:
700 : * -5 non-convergence of Gauss-Kronrod nodes
701 : calculation subroutine.
702 : * -1 incorrect parameters were specified
703 : * 1 OK
704 : * Rep.NFEV countains number of function calculations
705 : * Rep.NIntervals contains number of intervals [a,b]
706 : was partitioned into.
707 : *************************************************************************/
708 0 : _autogkreport_owner::_autogkreport_owner()
709 : {
710 : jmp_buf _break_jump;
711 : alglib_impl::ae_state _state;
712 :
713 0 : alglib_impl::ae_state_init(&_state);
714 0 : if( setjmp(_break_jump) )
715 : {
716 0 : if( p_struct!=NULL )
717 : {
718 0 : alglib_impl::_autogkreport_destroy(p_struct);
719 0 : alglib_impl::ae_free(p_struct);
720 : }
721 0 : p_struct = NULL;
722 : #if !defined(AE_NO_EXCEPTIONS)
723 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
724 : #else
725 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
726 : return;
727 : #endif
728 : }
729 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
730 0 : p_struct = NULL;
731 0 : p_struct = (alglib_impl::autogkreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::autogkreport), &_state);
732 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkreport));
733 0 : alglib_impl::_autogkreport_init(p_struct, &_state, ae_false);
734 0 : ae_state_clear(&_state);
735 0 : }
736 :
737 0 : _autogkreport_owner::_autogkreport_owner(const _autogkreport_owner &rhs)
738 : {
739 : jmp_buf _break_jump;
740 : alglib_impl::ae_state _state;
741 :
742 0 : alglib_impl::ae_state_init(&_state);
743 0 : if( setjmp(_break_jump) )
744 : {
745 0 : if( p_struct!=NULL )
746 : {
747 0 : alglib_impl::_autogkreport_destroy(p_struct);
748 0 : alglib_impl::ae_free(p_struct);
749 : }
750 0 : p_struct = NULL;
751 : #if !defined(AE_NO_EXCEPTIONS)
752 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
753 : #else
754 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
755 : return;
756 : #endif
757 : }
758 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
759 0 : p_struct = NULL;
760 0 : alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: autogkreport copy constructor failure (source is not initialized)", &_state);
761 0 : p_struct = (alglib_impl::autogkreport*)alglib_impl::ae_malloc(sizeof(alglib_impl::autogkreport), &_state);
762 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkreport));
763 0 : alglib_impl::_autogkreport_init_copy(p_struct, const_cast<alglib_impl::autogkreport*>(rhs.p_struct), &_state, ae_false);
764 0 : ae_state_clear(&_state);
765 0 : }
766 :
767 0 : _autogkreport_owner& _autogkreport_owner::operator=(const _autogkreport_owner &rhs)
768 : {
769 0 : if( this==&rhs )
770 0 : return *this;
771 : jmp_buf _break_jump;
772 : alglib_impl::ae_state _state;
773 :
774 0 : alglib_impl::ae_state_init(&_state);
775 0 : if( setjmp(_break_jump) )
776 : {
777 : #if !defined(AE_NO_EXCEPTIONS)
778 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
779 : #else
780 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
781 : return *this;
782 : #endif
783 : }
784 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
785 0 : alglib_impl::ae_assert(p_struct!=NULL, "ALGLIB: autogkreport assignment constructor failure (destination is not initialized)", &_state);
786 0 : alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: autogkreport assignment constructor failure (source is not initialized)", &_state);
787 0 : alglib_impl::_autogkreport_destroy(p_struct);
788 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkreport));
789 0 : alglib_impl::_autogkreport_init_copy(p_struct, const_cast<alglib_impl::autogkreport*>(rhs.p_struct), &_state, ae_false);
790 0 : ae_state_clear(&_state);
791 0 : return *this;
792 : }
793 :
794 0 : _autogkreport_owner::~_autogkreport_owner()
795 : {
796 0 : if( p_struct!=NULL )
797 : {
798 0 : alglib_impl::_autogkreport_destroy(p_struct);
799 0 : ae_free(p_struct);
800 : }
801 0 : }
802 :
803 0 : alglib_impl::autogkreport* _autogkreport_owner::c_ptr()
804 : {
805 0 : return p_struct;
806 : }
807 :
808 0 : alglib_impl::autogkreport* _autogkreport_owner::c_ptr() const
809 : {
810 0 : return const_cast<alglib_impl::autogkreport*>(p_struct);
811 : }
812 0 : autogkreport::autogkreport() : _autogkreport_owner() ,terminationtype(p_struct->terminationtype),nfev(p_struct->nfev),nintervals(p_struct->nintervals)
813 : {
814 0 : }
815 :
816 0 : autogkreport::autogkreport(const autogkreport &rhs):_autogkreport_owner(rhs) ,terminationtype(p_struct->terminationtype),nfev(p_struct->nfev),nintervals(p_struct->nintervals)
817 : {
818 0 : }
819 :
820 0 : autogkreport& autogkreport::operator=(const autogkreport &rhs)
821 : {
822 0 : if( this==&rhs )
823 0 : return *this;
824 0 : _autogkreport_owner::operator=(rhs);
825 0 : return *this;
826 : }
827 :
828 0 : autogkreport::~autogkreport()
829 : {
830 0 : }
831 :
832 :
833 : /*************************************************************************
834 : This structure stores state of the integration algorithm.
835 :
836 : Although this class has public fields, they are not intended for external
837 : use. You should use ALGLIB functions to work with this class:
838 : * autogksmooth()/AutoGKSmoothW()/... to create objects
839 : * autogkintegrate() to begin integration
840 : * autogkresults() to get results
841 : *************************************************************************/
842 0 : _autogkstate_owner::_autogkstate_owner()
843 : {
844 : jmp_buf _break_jump;
845 : alglib_impl::ae_state _state;
846 :
847 0 : alglib_impl::ae_state_init(&_state);
848 0 : if( setjmp(_break_jump) )
849 : {
850 0 : if( p_struct!=NULL )
851 : {
852 0 : alglib_impl::_autogkstate_destroy(p_struct);
853 0 : alglib_impl::ae_free(p_struct);
854 : }
855 0 : p_struct = NULL;
856 : #if !defined(AE_NO_EXCEPTIONS)
857 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
858 : #else
859 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
860 : return;
861 : #endif
862 : }
863 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
864 0 : p_struct = NULL;
865 0 : p_struct = (alglib_impl::autogkstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::autogkstate), &_state);
866 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkstate));
867 0 : alglib_impl::_autogkstate_init(p_struct, &_state, ae_false);
868 0 : ae_state_clear(&_state);
869 0 : }
870 :
871 0 : _autogkstate_owner::_autogkstate_owner(const _autogkstate_owner &rhs)
872 : {
873 : jmp_buf _break_jump;
874 : alglib_impl::ae_state _state;
875 :
876 0 : alglib_impl::ae_state_init(&_state);
877 0 : if( setjmp(_break_jump) )
878 : {
879 0 : if( p_struct!=NULL )
880 : {
881 0 : alglib_impl::_autogkstate_destroy(p_struct);
882 0 : alglib_impl::ae_free(p_struct);
883 : }
884 0 : p_struct = NULL;
885 : #if !defined(AE_NO_EXCEPTIONS)
886 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
887 : #else
888 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
889 : return;
890 : #endif
891 : }
892 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
893 0 : p_struct = NULL;
894 0 : alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: autogkstate copy constructor failure (source is not initialized)", &_state);
895 0 : p_struct = (alglib_impl::autogkstate*)alglib_impl::ae_malloc(sizeof(alglib_impl::autogkstate), &_state);
896 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkstate));
897 0 : alglib_impl::_autogkstate_init_copy(p_struct, const_cast<alglib_impl::autogkstate*>(rhs.p_struct), &_state, ae_false);
898 0 : ae_state_clear(&_state);
899 0 : }
900 :
901 0 : _autogkstate_owner& _autogkstate_owner::operator=(const _autogkstate_owner &rhs)
902 : {
903 0 : if( this==&rhs )
904 0 : return *this;
905 : jmp_buf _break_jump;
906 : alglib_impl::ae_state _state;
907 :
908 0 : alglib_impl::ae_state_init(&_state);
909 0 : if( setjmp(_break_jump) )
910 : {
911 : #if !defined(AE_NO_EXCEPTIONS)
912 0 : _ALGLIB_CPP_EXCEPTION(_state.error_msg);
913 : #else
914 : _ALGLIB_SET_ERROR_FLAG(_state.error_msg);
915 : return *this;
916 : #endif
917 : }
918 0 : alglib_impl::ae_state_set_break_jump(&_state, &_break_jump);
919 0 : alglib_impl::ae_assert(p_struct!=NULL, "ALGLIB: autogkstate assignment constructor failure (destination is not initialized)", &_state);
920 0 : alglib_impl::ae_assert(rhs.p_struct!=NULL, "ALGLIB: autogkstate assignment constructor failure (source is not initialized)", &_state);
921 0 : alglib_impl::_autogkstate_destroy(p_struct);
922 0 : memset(p_struct, 0, sizeof(alglib_impl::autogkstate));
923 0 : alglib_impl::_autogkstate_init_copy(p_struct, const_cast<alglib_impl::autogkstate*>(rhs.p_struct), &_state, ae_false);
924 0 : ae_state_clear(&_state);
925 0 : return *this;
926 : }
927 :
928 0 : _autogkstate_owner::~_autogkstate_owner()
929 : {
930 0 : if( p_struct!=NULL )
931 : {
932 0 : alglib_impl::_autogkstate_destroy(p_struct);
933 0 : ae_free(p_struct);
934 : }
935 0 : }
936 :
937 0 : alglib_impl::autogkstate* _autogkstate_owner::c_ptr()
938 : {
939 0 : return p_struct;
940 : }
941 :
942 0 : alglib_impl::autogkstate* _autogkstate_owner::c_ptr() const
943 : {
944 0 : return const_cast<alglib_impl::autogkstate*>(p_struct);
945 : }
946 0 : autogkstate::autogkstate() : _autogkstate_owner() ,needf(p_struct->needf),x(p_struct->x),xminusa(p_struct->xminusa),bminusx(p_struct->bminusx),f(p_struct->f)
947 : {
948 0 : }
949 :
950 0 : autogkstate::autogkstate(const autogkstate &rhs):_autogkstate_owner(rhs) ,needf(p_struct->needf),x(p_struct->x),xminusa(p_struct->xminusa),bminusx(p_struct->bminusx),f(p_struct->f)
951 : {
952 0 : }
953 :
954 0 : autogkstate& autogkstate::operator=(const autogkstate &rhs)
955 : {
956 0 : if( this==&rhs )
957 0 : return *this;
958 0 : _autogkstate_owner::operator=(rhs);
959 0 : return *this;
960 : }
961 :
962 0 : autogkstate::~autogkstate()
963 : {
964 0 : }
965 :
966 : /*************************************************************************
967 : Integration of a smooth function F(x) on a finite interval [a,b].
968 :
969 : Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
970 : is calculated with accuracy close to the machine precision.
971 :
972 : Algorithm works well only with smooth integrands. It may be used with
973 : continuous non-smooth integrands, but with less performance.
974 :
975 : It should never be used with integrands which have integrable singularities
976 : at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
977 : cases.
978 :
979 : INPUT PARAMETERS:
980 : A, B - interval boundaries (A<B, A=B or A>B)
981 :
982 : OUTPUT PARAMETERS
983 : State - structure which stores algorithm state
984 :
985 : SEE ALSO
986 : AutoGKSmoothW, AutoGKSingular, AutoGKResults.
987 :
988 :
989 : -- ALGLIB --
990 : Copyright 06.05.2009 by Bochkanov Sergey
991 : *************************************************************************/
992 0 : void autogksmooth(const double a, const double b, autogkstate &state, const xparams _xparams)
993 : {
994 : jmp_buf _break_jump;
995 : alglib_impl::ae_state _alglib_env_state;
996 0 : alglib_impl::ae_state_init(&_alglib_env_state);
997 0 : if( setjmp(_break_jump) )
998 : {
999 : #if !defined(AE_NO_EXCEPTIONS)
1000 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1001 : #else
1002 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1003 : return;
1004 : #endif
1005 : }
1006 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1007 0 : if( _xparams.flags!=0x0 )
1008 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1009 0 : alglib_impl::autogksmooth(a, b, const_cast<alglib_impl::autogkstate*>(state.c_ptr()), &_alglib_env_state);
1010 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1011 0 : return;
1012 : }
1013 :
1014 : /*************************************************************************
1015 : Integration of a smooth function F(x) on a finite interval [a,b].
1016 :
1017 : This subroutine is same as AutoGKSmooth(), but it guarantees that interval
1018 : [a,b] is partitioned into subintervals which have width at most XWidth.
1019 :
1020 : Subroutine can be used when integrating nearly-constant function with
1021 : narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
1022 : subroutine can overlook them.
1023 :
1024 : INPUT PARAMETERS:
1025 : A, B - interval boundaries (A<B, A=B or A>B)
1026 :
1027 : OUTPUT PARAMETERS
1028 : State - structure which stores algorithm state
1029 :
1030 : SEE ALSO
1031 : AutoGKSmooth, AutoGKSingular, AutoGKResults.
1032 :
1033 :
1034 : -- ALGLIB --
1035 : Copyright 06.05.2009 by Bochkanov Sergey
1036 : *************************************************************************/
1037 0 : void autogksmoothw(const double a, const double b, const double xwidth, autogkstate &state, const xparams _xparams)
1038 : {
1039 : jmp_buf _break_jump;
1040 : alglib_impl::ae_state _alglib_env_state;
1041 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1042 0 : if( setjmp(_break_jump) )
1043 : {
1044 : #if !defined(AE_NO_EXCEPTIONS)
1045 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1046 : #else
1047 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1048 : return;
1049 : #endif
1050 : }
1051 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1052 0 : if( _xparams.flags!=0x0 )
1053 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1054 0 : alglib_impl::autogksmoothw(a, b, xwidth, const_cast<alglib_impl::autogkstate*>(state.c_ptr()), &_alglib_env_state);
1055 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1056 0 : return;
1057 : }
1058 :
1059 : /*************************************************************************
1060 : Integration on a finite interval [A,B].
1061 : Integrand have integrable singularities at A/B.
1062 :
1063 : F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known
1064 : alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates
1065 : from below can be used (but these estimates should be greater than -1 too).
1066 :
1067 : One of alpha/beta variables (or even both alpha/beta) may be equal to 0,
1068 : which means than function F(x) is non-singular at A/B. Anyway (singular at
1069 : bounds or not), function F(x) is supposed to be continuous on (A,B).
1070 :
1071 : Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
1072 : is calculated with accuracy close to the machine precision.
1073 :
1074 : INPUT PARAMETERS:
1075 : A, B - interval boundaries (A<B, A=B or A>B)
1076 : Alpha - power-law coefficient of the F(x) at A,
1077 : Alpha>-1
1078 : Beta - power-law coefficient of the F(x) at B,
1079 : Beta>-1
1080 :
1081 : OUTPUT PARAMETERS
1082 : State - structure which stores algorithm state
1083 :
1084 : SEE ALSO
1085 : AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
1086 :
1087 :
1088 : -- ALGLIB --
1089 : Copyright 06.05.2009 by Bochkanov Sergey
1090 : *************************************************************************/
1091 0 : void autogksingular(const double a, const double b, const double alpha, const double beta, autogkstate &state, const xparams _xparams)
1092 : {
1093 : jmp_buf _break_jump;
1094 : alglib_impl::ae_state _alglib_env_state;
1095 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1096 0 : if( setjmp(_break_jump) )
1097 : {
1098 : #if !defined(AE_NO_EXCEPTIONS)
1099 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1100 : #else
1101 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1102 : return;
1103 : #endif
1104 : }
1105 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1106 0 : if( _xparams.flags!=0x0 )
1107 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1108 0 : alglib_impl::autogksingular(a, b, alpha, beta, const_cast<alglib_impl::autogkstate*>(state.c_ptr()), &_alglib_env_state);
1109 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1110 0 : return;
1111 : }
1112 :
1113 : /*************************************************************************
1114 : This function provides reverse communication interface
1115 : Reverse communication interface is not documented or recommended to use.
1116 : See below for functions which provide better documented API
1117 : *************************************************************************/
1118 0 : bool autogkiteration(const autogkstate &state, const xparams _xparams)
1119 : {
1120 : jmp_buf _break_jump;
1121 : alglib_impl::ae_state _alglib_env_state;
1122 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1123 0 : if( setjmp(_break_jump) )
1124 : {
1125 : #if !defined(AE_NO_EXCEPTIONS)
1126 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1127 : #else
1128 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1129 : return 0;
1130 : #endif
1131 : }
1132 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1133 0 : if( _xparams.flags!=0x0 )
1134 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1135 0 : ae_bool result = alglib_impl::autogkiteration(const_cast<alglib_impl::autogkstate*>(state.c_ptr()), &_alglib_env_state);
1136 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1137 0 : return *(reinterpret_cast<bool*>(&result));
1138 : }
1139 :
1140 :
1141 0 : void autogkintegrate(autogkstate &state,
1142 : void (*func)(double x, double xminusa, double bminusx, double &y, void *ptr),
1143 : void *ptr, const xparams _xparams){
1144 : jmp_buf _break_jump;
1145 : alglib_impl::ae_state _alglib_env_state;
1146 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1147 0 : if( setjmp(_break_jump) )
1148 : {
1149 : #if !defined(AE_NO_EXCEPTIONS)
1150 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1151 : #else
1152 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1153 : return;
1154 : #endif
1155 : }
1156 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1157 0 : if( _xparams.flags!=0x0 )
1158 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1159 0 : alglib_impl::ae_assert(func!=NULL, "ALGLIB: error in 'autogkintegrate()' (func is NULL)", &_alglib_env_state);
1160 0 : while( alglib_impl::autogkiteration(state.c_ptr(), &_alglib_env_state) )
1161 : {
1162 : _ALGLIB_CALLBACK_EXCEPTION_GUARD_BEGIN
1163 0 : if( state.needf )
1164 : {
1165 0 : func(state.x, state.xminusa, state.bminusx, state.f, ptr);
1166 0 : continue;
1167 : }
1168 0 : goto lbl_no_callback;
1169 0 : _ALGLIB_CALLBACK_EXCEPTION_GUARD_END
1170 0 : lbl_no_callback:
1171 0 : alglib_impl::ae_assert(ae_false, "ALGLIB: unexpected error in 'autogkintegrate()'", &_alglib_env_state);
1172 : }
1173 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1174 0 : }
1175 :
1176 :
1177 :
1178 : /*************************************************************************
1179 : Adaptive integration results
1180 :
1181 : Called after AutoGKIteration returned False.
1182 :
1183 : Input parameters:
1184 : State - algorithm state (used by AutoGKIteration).
1185 :
1186 : Output parameters:
1187 : V - integral(f(x)dx,a,b)
1188 : Rep - optimization report (see AutoGKReport description)
1189 :
1190 : -- ALGLIB --
1191 : Copyright 14.11.2007 by Bochkanov Sergey
1192 : *************************************************************************/
1193 0 : void autogkresults(const autogkstate &state, double &v, autogkreport &rep, const xparams _xparams)
1194 : {
1195 : jmp_buf _break_jump;
1196 : alglib_impl::ae_state _alglib_env_state;
1197 0 : alglib_impl::ae_state_init(&_alglib_env_state);
1198 0 : if( setjmp(_break_jump) )
1199 : {
1200 : #if !defined(AE_NO_EXCEPTIONS)
1201 0 : _ALGLIB_CPP_EXCEPTION(_alglib_env_state.error_msg);
1202 : #else
1203 : _ALGLIB_SET_ERROR_FLAG(_alglib_env_state.error_msg);
1204 : return;
1205 : #endif
1206 : }
1207 0 : ae_state_set_break_jump(&_alglib_env_state, &_break_jump);
1208 0 : if( _xparams.flags!=0x0 )
1209 0 : ae_state_set_flags(&_alglib_env_state, _xparams.flags);
1210 0 : alglib_impl::autogkresults(const_cast<alglib_impl::autogkstate*>(state.c_ptr()), &v, const_cast<alglib_impl::autogkreport*>(rep.c_ptr()), &_alglib_env_state);
1211 0 : alglib_impl::ae_state_clear(&_alglib_env_state);
1212 0 : return;
1213 : }
1214 : #endif
1215 : }
1216 :
1217 : /////////////////////////////////////////////////////////////////////////
1218 : //
1219 : // THIS SECTION CONTAINS IMPLEMENTATION OF COMPUTATIONAL CORE
1220 : //
1221 : /////////////////////////////////////////////////////////////////////////
1222 : namespace alglib_impl
1223 : {
1224 : #if defined(AE_COMPILE_GQ) || !defined(AE_PARTIAL_BUILD)
1225 :
1226 :
1227 : #endif
1228 : #if defined(AE_COMPILE_GKQ) || !defined(AE_PARTIAL_BUILD)
1229 :
1230 :
1231 : #endif
1232 : #if defined(AE_COMPILE_AUTOGK) || !defined(AE_PARTIAL_BUILD)
1233 : static ae_int_t autogk_maxsubintervals = 10000;
1234 : static void autogk_autogkinternalprepare(double a,
1235 : double b,
1236 : double eps,
1237 : double xwidth,
1238 : autogkinternalstate* state,
1239 : ae_state *_state);
1240 : static ae_bool autogk_autogkinternaliteration(autogkinternalstate* state,
1241 : ae_state *_state);
1242 : static void autogk_mheappop(/* Real */ ae_matrix* heap,
1243 : ae_int_t heapsize,
1244 : ae_int_t heapwidth,
1245 : ae_state *_state);
1246 : static void autogk_mheappush(/* Real */ ae_matrix* heap,
1247 : ae_int_t heapsize,
1248 : ae_int_t heapwidth,
1249 : ae_state *_state);
1250 : static void autogk_mheapresize(/* Real */ ae_matrix* heap,
1251 : ae_int_t* heapsize,
1252 : ae_int_t newheapsize,
1253 : ae_int_t heapwidth,
1254 : ae_state *_state);
1255 :
1256 :
1257 : #endif
1258 :
1259 : #if defined(AE_COMPILE_GQ) || !defined(AE_PARTIAL_BUILD)
1260 :
1261 :
1262 : /*************************************************************************
1263 : Computation of nodes and weights for a Gauss quadrature formula
1264 :
1265 : The algorithm generates the N-point Gauss quadrature formula with weight
1266 : function given by coefficients alpha and beta of a recurrence relation
1267 : which generates a system of orthogonal polynomials:
1268 :
1269 : P-1(x) = 0
1270 : P0(x) = 1
1271 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
1272 :
1273 : and zeroth moment Mu0
1274 :
1275 : Mu0 = integral(W(x)dx,a,b)
1276 :
1277 : INPUT PARAMETERS:
1278 : Alpha - array[0..N-1], alpha coefficients
1279 : Beta - array[0..N-1], beta coefficients
1280 : Zero-indexed element is not used and may be arbitrary.
1281 : Beta[I]>0.
1282 : Mu0 - zeroth moment of the weight function.
1283 : N - number of nodes of the quadrature formula, N>=1
1284 :
1285 : OUTPUT PARAMETERS:
1286 : Info - error code:
1287 : * -3 internal eigenproblem solver hasn't converged
1288 : * -2 Beta[i]<=0
1289 : * -1 incorrect N was passed
1290 : * 1 OK
1291 : X - array[0..N-1] - array of quadrature nodes,
1292 : in ascending order.
1293 : W - array[0..N-1] - array of quadrature weights.
1294 :
1295 : -- ALGLIB --
1296 : Copyright 2005-2009 by Bochkanov Sergey
1297 : *************************************************************************/
1298 0 : void gqgeneraterec(/* Real */ ae_vector* alpha,
1299 : /* Real */ ae_vector* beta,
1300 : double mu0,
1301 : ae_int_t n,
1302 : ae_int_t* info,
1303 : /* Real */ ae_vector* x,
1304 : /* Real */ ae_vector* w,
1305 : ae_state *_state)
1306 : {
1307 : ae_frame _frame_block;
1308 : ae_int_t i;
1309 : ae_vector d;
1310 : ae_vector e;
1311 : ae_matrix z;
1312 :
1313 0 : ae_frame_make(_state, &_frame_block);
1314 0 : memset(&d, 0, sizeof(d));
1315 0 : memset(&e, 0, sizeof(e));
1316 0 : memset(&z, 0, sizeof(z));
1317 0 : *info = 0;
1318 0 : ae_vector_clear(x);
1319 0 : ae_vector_clear(w);
1320 0 : ae_vector_init(&d, 0, DT_REAL, _state, ae_true);
1321 0 : ae_vector_init(&e, 0, DT_REAL, _state, ae_true);
1322 0 : ae_matrix_init(&z, 0, 0, DT_REAL, _state, ae_true);
1323 :
1324 0 : if( n<1 )
1325 : {
1326 0 : *info = -1;
1327 0 : ae_frame_leave(_state);
1328 0 : return;
1329 : }
1330 0 : *info = 1;
1331 :
1332 : /*
1333 : * Initialize
1334 : */
1335 0 : ae_vector_set_length(&d, n, _state);
1336 0 : ae_vector_set_length(&e, n, _state);
1337 0 : for(i=1; i<=n-1; i++)
1338 : {
1339 0 : d.ptr.p_double[i-1] = alpha->ptr.p_double[i-1];
1340 0 : if( ae_fp_less_eq(beta->ptr.p_double[i],(double)(0)) )
1341 : {
1342 0 : *info = -2;
1343 0 : ae_frame_leave(_state);
1344 0 : return;
1345 : }
1346 0 : e.ptr.p_double[i-1] = ae_sqrt(beta->ptr.p_double[i], _state);
1347 : }
1348 0 : d.ptr.p_double[n-1] = alpha->ptr.p_double[n-1];
1349 :
1350 : /*
1351 : * EVD
1352 : */
1353 0 : if( !smatrixtdevd(&d, &e, n, 3, &z, _state) )
1354 : {
1355 0 : *info = -3;
1356 0 : ae_frame_leave(_state);
1357 0 : return;
1358 : }
1359 :
1360 : /*
1361 : * Generate
1362 : */
1363 0 : ae_vector_set_length(x, n, _state);
1364 0 : ae_vector_set_length(w, n, _state);
1365 0 : for(i=1; i<=n; i++)
1366 : {
1367 0 : x->ptr.p_double[i-1] = d.ptr.p_double[i-1];
1368 0 : w->ptr.p_double[i-1] = mu0*ae_sqr(z.ptr.pp_double[0][i-1], _state);
1369 : }
1370 0 : ae_frame_leave(_state);
1371 : }
1372 :
1373 :
1374 : /*************************************************************************
1375 : Computation of nodes and weights for a Gauss-Lobatto quadrature formula
1376 :
1377 : The algorithm generates the N-point Gauss-Lobatto quadrature formula with
1378 : weight function given by coefficients alpha and beta of a recurrence which
1379 : generates a system of orthogonal polynomials.
1380 :
1381 : P-1(x) = 0
1382 : P0(x) = 1
1383 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
1384 :
1385 : and zeroth moment Mu0
1386 :
1387 : Mu0 = integral(W(x)dx,a,b)
1388 :
1389 : INPUT PARAMETERS:
1390 : Alpha - array[0..N-2], alpha coefficients
1391 : Beta - array[0..N-2], beta coefficients.
1392 : Zero-indexed element is not used, may be arbitrary.
1393 : Beta[I]>0
1394 : Mu0 - zeroth moment of the weighting function.
1395 : A - left boundary of the integration interval.
1396 : B - right boundary of the integration interval.
1397 : N - number of nodes of the quadrature formula, N>=3
1398 : (including the left and right boundary nodes).
1399 :
1400 : OUTPUT PARAMETERS:
1401 : Info - error code:
1402 : * -3 internal eigenproblem solver hasn't converged
1403 : * -2 Beta[i]<=0
1404 : * -1 incorrect N was passed
1405 : * 1 OK
1406 : X - array[0..N-1] - array of quadrature nodes,
1407 : in ascending order.
1408 : W - array[0..N-1] - array of quadrature weights.
1409 :
1410 : -- ALGLIB --
1411 : Copyright 2005-2009 by Bochkanov Sergey
1412 : *************************************************************************/
1413 0 : void gqgenerategausslobattorec(/* Real */ ae_vector* alpha,
1414 : /* Real */ ae_vector* beta,
1415 : double mu0,
1416 : double a,
1417 : double b,
1418 : ae_int_t n,
1419 : ae_int_t* info,
1420 : /* Real */ ae_vector* x,
1421 : /* Real */ ae_vector* w,
1422 : ae_state *_state)
1423 : {
1424 : ae_frame _frame_block;
1425 : ae_vector _alpha;
1426 : ae_vector _beta;
1427 : ae_int_t i;
1428 : ae_vector d;
1429 : ae_vector e;
1430 : ae_matrix z;
1431 : double pim1a;
1432 : double pia;
1433 : double pim1b;
1434 : double pib;
1435 : double t;
1436 : double a11;
1437 : double a12;
1438 : double a21;
1439 : double a22;
1440 : double b1;
1441 : double b2;
1442 : double alph;
1443 : double bet;
1444 :
1445 0 : ae_frame_make(_state, &_frame_block);
1446 0 : memset(&_alpha, 0, sizeof(_alpha));
1447 0 : memset(&_beta, 0, sizeof(_beta));
1448 0 : memset(&d, 0, sizeof(d));
1449 0 : memset(&e, 0, sizeof(e));
1450 0 : memset(&z, 0, sizeof(z));
1451 0 : ae_vector_init_copy(&_alpha, alpha, _state, ae_true);
1452 0 : alpha = &_alpha;
1453 0 : ae_vector_init_copy(&_beta, beta, _state, ae_true);
1454 0 : beta = &_beta;
1455 0 : *info = 0;
1456 0 : ae_vector_clear(x);
1457 0 : ae_vector_clear(w);
1458 0 : ae_vector_init(&d, 0, DT_REAL, _state, ae_true);
1459 0 : ae_vector_init(&e, 0, DT_REAL, _state, ae_true);
1460 0 : ae_matrix_init(&z, 0, 0, DT_REAL, _state, ae_true);
1461 :
1462 0 : if( n<=2 )
1463 : {
1464 0 : *info = -1;
1465 0 : ae_frame_leave(_state);
1466 0 : return;
1467 : }
1468 0 : *info = 1;
1469 :
1470 : /*
1471 : * Initialize, D[1:N+1], E[1:N]
1472 : */
1473 0 : n = n-2;
1474 0 : ae_vector_set_length(&d, n+2, _state);
1475 0 : ae_vector_set_length(&e, n+1, _state);
1476 0 : for(i=1; i<=n+1; i++)
1477 : {
1478 0 : d.ptr.p_double[i-1] = alpha->ptr.p_double[i-1];
1479 : }
1480 0 : for(i=1; i<=n; i++)
1481 : {
1482 0 : if( ae_fp_less_eq(beta->ptr.p_double[i],(double)(0)) )
1483 : {
1484 0 : *info = -2;
1485 0 : ae_frame_leave(_state);
1486 0 : return;
1487 : }
1488 0 : e.ptr.p_double[i-1] = ae_sqrt(beta->ptr.p_double[i], _state);
1489 : }
1490 :
1491 : /*
1492 : * Caclulate Pn(a), Pn+1(a), Pn(b), Pn+1(b)
1493 : */
1494 0 : beta->ptr.p_double[0] = (double)(0);
1495 0 : pim1a = (double)(0);
1496 0 : pia = (double)(1);
1497 0 : pim1b = (double)(0);
1498 0 : pib = (double)(1);
1499 0 : for(i=1; i<=n+1; i++)
1500 : {
1501 :
1502 : /*
1503 : * Pi(a)
1504 : */
1505 0 : t = (a-alpha->ptr.p_double[i-1])*pia-beta->ptr.p_double[i-1]*pim1a;
1506 0 : pim1a = pia;
1507 0 : pia = t;
1508 :
1509 : /*
1510 : * Pi(b)
1511 : */
1512 0 : t = (b-alpha->ptr.p_double[i-1])*pib-beta->ptr.p_double[i-1]*pim1b;
1513 0 : pim1b = pib;
1514 0 : pib = t;
1515 : }
1516 :
1517 : /*
1518 : * Calculate alpha'(n+1), beta'(n+1)
1519 : */
1520 0 : a11 = pia;
1521 0 : a12 = pim1a;
1522 0 : a21 = pib;
1523 0 : a22 = pim1b;
1524 0 : b1 = a*pia;
1525 0 : b2 = b*pib;
1526 0 : if( ae_fp_greater(ae_fabs(a11, _state),ae_fabs(a21, _state)) )
1527 : {
1528 0 : a22 = a22-a12*a21/a11;
1529 0 : b2 = b2-b1*a21/a11;
1530 0 : bet = b2/a22;
1531 0 : alph = (b1-bet*a12)/a11;
1532 : }
1533 : else
1534 : {
1535 0 : a12 = a12-a22*a11/a21;
1536 0 : b1 = b1-b2*a11/a21;
1537 0 : bet = b1/a12;
1538 0 : alph = (b2-bet*a22)/a21;
1539 : }
1540 0 : if( ae_fp_less(bet,(double)(0)) )
1541 : {
1542 0 : *info = -3;
1543 0 : ae_frame_leave(_state);
1544 0 : return;
1545 : }
1546 0 : d.ptr.p_double[n+1] = alph;
1547 0 : e.ptr.p_double[n] = ae_sqrt(bet, _state);
1548 :
1549 : /*
1550 : * EVD
1551 : */
1552 0 : if( !smatrixtdevd(&d, &e, n+2, 3, &z, _state) )
1553 : {
1554 0 : *info = -3;
1555 0 : ae_frame_leave(_state);
1556 0 : return;
1557 : }
1558 :
1559 : /*
1560 : * Generate
1561 : */
1562 0 : ae_vector_set_length(x, n+2, _state);
1563 0 : ae_vector_set_length(w, n+2, _state);
1564 0 : for(i=1; i<=n+2; i++)
1565 : {
1566 0 : x->ptr.p_double[i-1] = d.ptr.p_double[i-1];
1567 0 : w->ptr.p_double[i-1] = mu0*ae_sqr(z.ptr.pp_double[0][i-1], _state);
1568 : }
1569 0 : ae_frame_leave(_state);
1570 : }
1571 :
1572 :
1573 : /*************************************************************************
1574 : Computation of nodes and weights for a Gauss-Radau quadrature formula
1575 :
1576 : The algorithm generates the N-point Gauss-Radau quadrature formula with
1577 : weight function given by the coefficients alpha and beta of a recurrence
1578 : which generates a system of orthogonal polynomials.
1579 :
1580 : P-1(x) = 0
1581 : P0(x) = 1
1582 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
1583 :
1584 : and zeroth moment Mu0
1585 :
1586 : Mu0 = integral(W(x)dx,a,b)
1587 :
1588 : INPUT PARAMETERS:
1589 : Alpha - array[0..N-2], alpha coefficients.
1590 : Beta - array[0..N-1], beta coefficients
1591 : Zero-indexed element is not used.
1592 : Beta[I]>0
1593 : Mu0 - zeroth moment of the weighting function.
1594 : A - left boundary of the integration interval.
1595 : N - number of nodes of the quadrature formula, N>=2
1596 : (including the left boundary node).
1597 :
1598 : OUTPUT PARAMETERS:
1599 : Info - error code:
1600 : * -3 internal eigenproblem solver hasn't converged
1601 : * -2 Beta[i]<=0
1602 : * -1 incorrect N was passed
1603 : * 1 OK
1604 : X - array[0..N-1] - array of quadrature nodes,
1605 : in ascending order.
1606 : W - array[0..N-1] - array of quadrature weights.
1607 :
1608 :
1609 : -- ALGLIB --
1610 : Copyright 2005-2009 by Bochkanov Sergey
1611 : *************************************************************************/
1612 0 : void gqgenerategaussradaurec(/* Real */ ae_vector* alpha,
1613 : /* Real */ ae_vector* beta,
1614 : double mu0,
1615 : double a,
1616 : ae_int_t n,
1617 : ae_int_t* info,
1618 : /* Real */ ae_vector* x,
1619 : /* Real */ ae_vector* w,
1620 : ae_state *_state)
1621 : {
1622 : ae_frame _frame_block;
1623 : ae_vector _alpha;
1624 : ae_vector _beta;
1625 : ae_int_t i;
1626 : ae_vector d;
1627 : ae_vector e;
1628 : ae_matrix z;
1629 : double polim1;
1630 : double poli;
1631 : double t;
1632 :
1633 0 : ae_frame_make(_state, &_frame_block);
1634 0 : memset(&_alpha, 0, sizeof(_alpha));
1635 0 : memset(&_beta, 0, sizeof(_beta));
1636 0 : memset(&d, 0, sizeof(d));
1637 0 : memset(&e, 0, sizeof(e));
1638 0 : memset(&z, 0, sizeof(z));
1639 0 : ae_vector_init_copy(&_alpha, alpha, _state, ae_true);
1640 0 : alpha = &_alpha;
1641 0 : ae_vector_init_copy(&_beta, beta, _state, ae_true);
1642 0 : beta = &_beta;
1643 0 : *info = 0;
1644 0 : ae_vector_clear(x);
1645 0 : ae_vector_clear(w);
1646 0 : ae_vector_init(&d, 0, DT_REAL, _state, ae_true);
1647 0 : ae_vector_init(&e, 0, DT_REAL, _state, ae_true);
1648 0 : ae_matrix_init(&z, 0, 0, DT_REAL, _state, ae_true);
1649 :
1650 0 : if( n<2 )
1651 : {
1652 0 : *info = -1;
1653 0 : ae_frame_leave(_state);
1654 0 : return;
1655 : }
1656 0 : *info = 1;
1657 :
1658 : /*
1659 : * Initialize, D[1:N], E[1:N]
1660 : */
1661 0 : n = n-1;
1662 0 : ae_vector_set_length(&d, n+1, _state);
1663 0 : ae_vector_set_length(&e, n, _state);
1664 0 : for(i=1; i<=n; i++)
1665 : {
1666 0 : d.ptr.p_double[i-1] = alpha->ptr.p_double[i-1];
1667 0 : if( ae_fp_less_eq(beta->ptr.p_double[i],(double)(0)) )
1668 : {
1669 0 : *info = -2;
1670 0 : ae_frame_leave(_state);
1671 0 : return;
1672 : }
1673 0 : e.ptr.p_double[i-1] = ae_sqrt(beta->ptr.p_double[i], _state);
1674 : }
1675 :
1676 : /*
1677 : * Caclulate Pn(a), Pn-1(a), and D[N+1]
1678 : */
1679 0 : beta->ptr.p_double[0] = (double)(0);
1680 0 : polim1 = (double)(0);
1681 0 : poli = (double)(1);
1682 0 : for(i=1; i<=n; i++)
1683 : {
1684 0 : t = (a-alpha->ptr.p_double[i-1])*poli-beta->ptr.p_double[i-1]*polim1;
1685 0 : polim1 = poli;
1686 0 : poli = t;
1687 : }
1688 0 : d.ptr.p_double[n] = a-beta->ptr.p_double[n]*polim1/poli;
1689 :
1690 : /*
1691 : * EVD
1692 : */
1693 0 : if( !smatrixtdevd(&d, &e, n+1, 3, &z, _state) )
1694 : {
1695 0 : *info = -3;
1696 0 : ae_frame_leave(_state);
1697 0 : return;
1698 : }
1699 :
1700 : /*
1701 : * Generate
1702 : */
1703 0 : ae_vector_set_length(x, n+1, _state);
1704 0 : ae_vector_set_length(w, n+1, _state);
1705 0 : for(i=1; i<=n+1; i++)
1706 : {
1707 0 : x->ptr.p_double[i-1] = d.ptr.p_double[i-1];
1708 0 : w->ptr.p_double[i-1] = mu0*ae_sqr(z.ptr.pp_double[0][i-1], _state);
1709 : }
1710 0 : ae_frame_leave(_state);
1711 : }
1712 :
1713 :
1714 : /*************************************************************************
1715 : Returns nodes/weights for Gauss-Legendre quadrature on [-1,1] with N
1716 : nodes.
1717 :
1718 : INPUT PARAMETERS:
1719 : N - number of nodes, >=1
1720 :
1721 : OUTPUT PARAMETERS:
1722 : Info - error code:
1723 : * -4 an error was detected when calculating
1724 : weights/nodes. N is too large to obtain
1725 : weights/nodes with high enough accuracy.
1726 : Try to use multiple precision version.
1727 : * -3 internal eigenproblem solver hasn't converged
1728 : * -1 incorrect N was passed
1729 : * +1 OK
1730 : X - array[0..N-1] - array of quadrature nodes,
1731 : in ascending order.
1732 : W - array[0..N-1] - array of quadrature weights.
1733 :
1734 :
1735 : -- ALGLIB --
1736 : Copyright 12.05.2009 by Bochkanov Sergey
1737 : *************************************************************************/
1738 0 : void gqgenerategausslegendre(ae_int_t n,
1739 : ae_int_t* info,
1740 : /* Real */ ae_vector* x,
1741 : /* Real */ ae_vector* w,
1742 : ae_state *_state)
1743 : {
1744 : ae_frame _frame_block;
1745 : ae_vector alpha;
1746 : ae_vector beta;
1747 : ae_int_t i;
1748 :
1749 0 : ae_frame_make(_state, &_frame_block);
1750 0 : memset(&alpha, 0, sizeof(alpha));
1751 0 : memset(&beta, 0, sizeof(beta));
1752 0 : *info = 0;
1753 0 : ae_vector_clear(x);
1754 0 : ae_vector_clear(w);
1755 0 : ae_vector_init(&alpha, 0, DT_REAL, _state, ae_true);
1756 0 : ae_vector_init(&beta, 0, DT_REAL, _state, ae_true);
1757 :
1758 0 : if( n<1 )
1759 : {
1760 0 : *info = -1;
1761 0 : ae_frame_leave(_state);
1762 0 : return;
1763 : }
1764 0 : ae_vector_set_length(&alpha, n, _state);
1765 0 : ae_vector_set_length(&beta, n, _state);
1766 0 : for(i=0; i<=n-1; i++)
1767 : {
1768 0 : alpha.ptr.p_double[i] = (double)(0);
1769 : }
1770 0 : beta.ptr.p_double[0] = (double)(2);
1771 0 : for(i=1; i<=n-1; i++)
1772 : {
1773 0 : beta.ptr.p_double[i] = 1/(4-1/ae_sqr((double)(i), _state));
1774 : }
1775 0 : gqgeneraterec(&alpha, &beta, beta.ptr.p_double[0], n, info, x, w, _state);
1776 :
1777 : /*
1778 : * test basic properties to detect errors
1779 : */
1780 0 : if( *info>0 )
1781 : {
1782 0 : if( ae_fp_less(x->ptr.p_double[0],(double)(-1))||ae_fp_greater(x->ptr.p_double[n-1],(double)(1)) )
1783 : {
1784 0 : *info = -4;
1785 : }
1786 0 : for(i=0; i<=n-2; i++)
1787 : {
1788 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
1789 : {
1790 0 : *info = -4;
1791 : }
1792 : }
1793 : }
1794 0 : ae_frame_leave(_state);
1795 : }
1796 :
1797 :
1798 : /*************************************************************************
1799 : Returns nodes/weights for Gauss-Jacobi quadrature on [-1,1] with weight
1800 : function W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
1801 :
1802 : INPUT PARAMETERS:
1803 : N - number of nodes, >=1
1804 : Alpha - power-law coefficient, Alpha>-1
1805 : Beta - power-law coefficient, Beta>-1
1806 :
1807 : OUTPUT PARAMETERS:
1808 : Info - error code:
1809 : * -4 an error was detected when calculating
1810 : weights/nodes. Alpha or Beta are too close
1811 : to -1 to obtain weights/nodes with high enough
1812 : accuracy, or, may be, N is too large. Try to
1813 : use multiple precision version.
1814 : * -3 internal eigenproblem solver hasn't converged
1815 : * -1 incorrect N/Alpha/Beta was passed
1816 : * +1 OK
1817 : X - array[0..N-1] - array of quadrature nodes,
1818 : in ascending order.
1819 : W - array[0..N-1] - array of quadrature weights.
1820 :
1821 :
1822 : -- ALGLIB --
1823 : Copyright 12.05.2009 by Bochkanov Sergey
1824 : *************************************************************************/
1825 0 : void gqgenerategaussjacobi(ae_int_t n,
1826 : double alpha,
1827 : double beta,
1828 : ae_int_t* info,
1829 : /* Real */ ae_vector* x,
1830 : /* Real */ ae_vector* w,
1831 : ae_state *_state)
1832 : {
1833 : ae_frame _frame_block;
1834 : ae_vector a;
1835 : ae_vector b;
1836 : double alpha2;
1837 : double beta2;
1838 : double apb;
1839 : double t;
1840 : ae_int_t i;
1841 : double s;
1842 :
1843 0 : ae_frame_make(_state, &_frame_block);
1844 0 : memset(&a, 0, sizeof(a));
1845 0 : memset(&b, 0, sizeof(b));
1846 0 : *info = 0;
1847 0 : ae_vector_clear(x);
1848 0 : ae_vector_clear(w);
1849 0 : ae_vector_init(&a, 0, DT_REAL, _state, ae_true);
1850 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
1851 :
1852 0 : if( (n<1||ae_fp_less_eq(alpha,(double)(-1)))||ae_fp_less_eq(beta,(double)(-1)) )
1853 : {
1854 0 : *info = -1;
1855 0 : ae_frame_leave(_state);
1856 0 : return;
1857 : }
1858 0 : ae_vector_set_length(&a, n, _state);
1859 0 : ae_vector_set_length(&b, n, _state);
1860 0 : apb = alpha+beta;
1861 0 : a.ptr.p_double[0] = (beta-alpha)/(apb+2);
1862 0 : t = (apb+1)*ae_log((double)(2), _state)+lngamma(alpha+1, &s, _state)+lngamma(beta+1, &s, _state)-lngamma(apb+2, &s, _state);
1863 0 : if( ae_fp_greater(t,ae_log(ae_maxrealnumber, _state)) )
1864 : {
1865 0 : *info = -4;
1866 0 : ae_frame_leave(_state);
1867 0 : return;
1868 : }
1869 0 : b.ptr.p_double[0] = ae_exp(t, _state);
1870 0 : if( n>1 )
1871 : {
1872 0 : alpha2 = ae_sqr(alpha, _state);
1873 0 : beta2 = ae_sqr(beta, _state);
1874 0 : a.ptr.p_double[1] = (beta2-alpha2)/((apb+2)*(apb+4));
1875 0 : b.ptr.p_double[1] = 4*(alpha+1)*(beta+1)/((apb+3)*ae_sqr(apb+2, _state));
1876 0 : for(i=2; i<=n-1; i++)
1877 : {
1878 0 : a.ptr.p_double[i] = 0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
1879 0 : b.ptr.p_double[i] = 0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*ae_sqr(1+0.5*apb/i, _state));
1880 : }
1881 : }
1882 0 : gqgeneraterec(&a, &b, b.ptr.p_double[0], n, info, x, w, _state);
1883 :
1884 : /*
1885 : * test basic properties to detect errors
1886 : */
1887 0 : if( *info>0 )
1888 : {
1889 0 : if( ae_fp_less(x->ptr.p_double[0],(double)(-1))||ae_fp_greater(x->ptr.p_double[n-1],(double)(1)) )
1890 : {
1891 0 : *info = -4;
1892 : }
1893 0 : for(i=0; i<=n-2; i++)
1894 : {
1895 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
1896 : {
1897 0 : *info = -4;
1898 : }
1899 : }
1900 : }
1901 0 : ae_frame_leave(_state);
1902 : }
1903 :
1904 :
1905 : /*************************************************************************
1906 : Returns nodes/weights for Gauss-Laguerre quadrature on [0,+inf) with
1907 : weight function W(x)=Power(x,Alpha)*Exp(-x)
1908 :
1909 : INPUT PARAMETERS:
1910 : N - number of nodes, >=1
1911 : Alpha - power-law coefficient, Alpha>-1
1912 :
1913 : OUTPUT PARAMETERS:
1914 : Info - error code:
1915 : * -4 an error was detected when calculating
1916 : weights/nodes. Alpha is too close to -1 to
1917 : obtain weights/nodes with high enough accuracy
1918 : or, may be, N is too large. Try to use
1919 : multiple precision version.
1920 : * -3 internal eigenproblem solver hasn't converged
1921 : * -1 incorrect N/Alpha was passed
1922 : * +1 OK
1923 : X - array[0..N-1] - array of quadrature nodes,
1924 : in ascending order.
1925 : W - array[0..N-1] - array of quadrature weights.
1926 :
1927 :
1928 : -- ALGLIB --
1929 : Copyright 12.05.2009 by Bochkanov Sergey
1930 : *************************************************************************/
1931 0 : void gqgenerategausslaguerre(ae_int_t n,
1932 : double alpha,
1933 : ae_int_t* info,
1934 : /* Real */ ae_vector* x,
1935 : /* Real */ ae_vector* w,
1936 : ae_state *_state)
1937 : {
1938 : ae_frame _frame_block;
1939 : ae_vector a;
1940 : ae_vector b;
1941 : double t;
1942 : ae_int_t i;
1943 : double s;
1944 :
1945 0 : ae_frame_make(_state, &_frame_block);
1946 0 : memset(&a, 0, sizeof(a));
1947 0 : memset(&b, 0, sizeof(b));
1948 0 : *info = 0;
1949 0 : ae_vector_clear(x);
1950 0 : ae_vector_clear(w);
1951 0 : ae_vector_init(&a, 0, DT_REAL, _state, ae_true);
1952 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
1953 :
1954 0 : if( n<1||ae_fp_less_eq(alpha,(double)(-1)) )
1955 : {
1956 0 : *info = -1;
1957 0 : ae_frame_leave(_state);
1958 0 : return;
1959 : }
1960 0 : ae_vector_set_length(&a, n, _state);
1961 0 : ae_vector_set_length(&b, n, _state);
1962 0 : a.ptr.p_double[0] = alpha+1;
1963 0 : t = lngamma(alpha+1, &s, _state);
1964 0 : if( ae_fp_greater_eq(t,ae_log(ae_maxrealnumber, _state)) )
1965 : {
1966 0 : *info = -4;
1967 0 : ae_frame_leave(_state);
1968 0 : return;
1969 : }
1970 0 : b.ptr.p_double[0] = ae_exp(t, _state);
1971 0 : if( n>1 )
1972 : {
1973 0 : for(i=1; i<=n-1; i++)
1974 : {
1975 0 : a.ptr.p_double[i] = 2*i+alpha+1;
1976 0 : b.ptr.p_double[i] = i*(i+alpha);
1977 : }
1978 : }
1979 0 : gqgeneraterec(&a, &b, b.ptr.p_double[0], n, info, x, w, _state);
1980 :
1981 : /*
1982 : * test basic properties to detect errors
1983 : */
1984 0 : if( *info>0 )
1985 : {
1986 0 : if( ae_fp_less(x->ptr.p_double[0],(double)(0)) )
1987 : {
1988 0 : *info = -4;
1989 : }
1990 0 : for(i=0; i<=n-2; i++)
1991 : {
1992 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
1993 : {
1994 0 : *info = -4;
1995 : }
1996 : }
1997 : }
1998 0 : ae_frame_leave(_state);
1999 : }
2000 :
2001 :
2002 : /*************************************************************************
2003 : Returns nodes/weights for Gauss-Hermite quadrature on (-inf,+inf) with
2004 : weight function W(x)=Exp(-x*x)
2005 :
2006 : INPUT PARAMETERS:
2007 : N - number of nodes, >=1
2008 :
2009 : OUTPUT PARAMETERS:
2010 : Info - error code:
2011 : * -4 an error was detected when calculating
2012 : weights/nodes. May be, N is too large. Try to
2013 : use multiple precision version.
2014 : * -3 internal eigenproblem solver hasn't converged
2015 : * -1 incorrect N/Alpha was passed
2016 : * +1 OK
2017 : X - array[0..N-1] - array of quadrature nodes,
2018 : in ascending order.
2019 : W - array[0..N-1] - array of quadrature weights.
2020 :
2021 :
2022 : -- ALGLIB --
2023 : Copyright 12.05.2009 by Bochkanov Sergey
2024 : *************************************************************************/
2025 0 : void gqgenerategausshermite(ae_int_t n,
2026 : ae_int_t* info,
2027 : /* Real */ ae_vector* x,
2028 : /* Real */ ae_vector* w,
2029 : ae_state *_state)
2030 : {
2031 : ae_frame _frame_block;
2032 : ae_vector a;
2033 : ae_vector b;
2034 : ae_int_t i;
2035 :
2036 0 : ae_frame_make(_state, &_frame_block);
2037 0 : memset(&a, 0, sizeof(a));
2038 0 : memset(&b, 0, sizeof(b));
2039 0 : *info = 0;
2040 0 : ae_vector_clear(x);
2041 0 : ae_vector_clear(w);
2042 0 : ae_vector_init(&a, 0, DT_REAL, _state, ae_true);
2043 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
2044 :
2045 0 : if( n<1 )
2046 : {
2047 0 : *info = -1;
2048 0 : ae_frame_leave(_state);
2049 0 : return;
2050 : }
2051 0 : ae_vector_set_length(&a, n, _state);
2052 0 : ae_vector_set_length(&b, n, _state);
2053 0 : for(i=0; i<=n-1; i++)
2054 : {
2055 0 : a.ptr.p_double[i] = (double)(0);
2056 : }
2057 0 : b.ptr.p_double[0] = ae_sqrt(4*ae_atan((double)(1), _state), _state);
2058 0 : if( n>1 )
2059 : {
2060 0 : for(i=1; i<=n-1; i++)
2061 : {
2062 0 : b.ptr.p_double[i] = 0.5*i;
2063 : }
2064 : }
2065 0 : gqgeneraterec(&a, &b, b.ptr.p_double[0], n, info, x, w, _state);
2066 :
2067 : /*
2068 : * test basic properties to detect errors
2069 : */
2070 0 : if( *info>0 )
2071 : {
2072 0 : for(i=0; i<=n-2; i++)
2073 : {
2074 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
2075 : {
2076 0 : *info = -4;
2077 : }
2078 : }
2079 : }
2080 0 : ae_frame_leave(_state);
2081 : }
2082 :
2083 :
2084 : #endif
2085 : #if defined(AE_COMPILE_GKQ) || !defined(AE_PARTIAL_BUILD)
2086 :
2087 :
2088 : /*************************************************************************
2089 : Computation of nodes and weights of a Gauss-Kronrod quadrature formula
2090 :
2091 : The algorithm generates the N-point Gauss-Kronrod quadrature formula with
2092 : weight function given by coefficients alpha and beta of a recurrence
2093 : relation which generates a system of orthogonal polynomials:
2094 :
2095 : P-1(x) = 0
2096 : P0(x) = 1
2097 : Pn+1(x) = (x-alpha(n))*Pn(x) - beta(n)*Pn-1(x)
2098 :
2099 : and zero moment Mu0
2100 :
2101 : Mu0 = integral(W(x)dx,a,b)
2102 :
2103 :
2104 : INPUT PARAMETERS:
2105 : Alpha - alpha coefficients, array[0..floor(3*K/2)].
2106 : Beta - beta coefficients, array[0..ceil(3*K/2)].
2107 : Beta[0] is not used and may be arbitrary.
2108 : Beta[I]>0.
2109 : Mu0 - zeroth moment of the weight function.
2110 : N - number of nodes of the Gauss-Kronrod quadrature formula,
2111 : N >= 3,
2112 : N = 2*K+1.
2113 :
2114 : OUTPUT PARAMETERS:
2115 : Info - error code:
2116 : * -5 no real and positive Gauss-Kronrod formula can
2117 : be created for such a weight function with a
2118 : given number of nodes.
2119 : * -4 N is too large, task may be ill conditioned -
2120 : x[i]=x[i+1] found.
2121 : * -3 internal eigenproblem solver hasn't converged
2122 : * -2 Beta[i]<=0
2123 : * -1 incorrect N was passed
2124 : * +1 OK
2125 : X - array[0..N-1] - array of quadrature nodes,
2126 : in ascending order.
2127 : WKronrod - array[0..N-1] - Kronrod weights
2128 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
2129 : corresponding to extended Kronrod nodes).
2130 :
2131 : -- ALGLIB --
2132 : Copyright 08.05.2009 by Bochkanov Sergey
2133 : *************************************************************************/
2134 0 : void gkqgeneraterec(/* Real */ ae_vector* alpha,
2135 : /* Real */ ae_vector* beta,
2136 : double mu0,
2137 : ae_int_t n,
2138 : ae_int_t* info,
2139 : /* Real */ ae_vector* x,
2140 : /* Real */ ae_vector* wkronrod,
2141 : /* Real */ ae_vector* wgauss,
2142 : ae_state *_state)
2143 : {
2144 : ae_frame _frame_block;
2145 : ae_vector _alpha;
2146 : ae_vector _beta;
2147 : ae_vector ta;
2148 : ae_int_t i;
2149 : ae_int_t j;
2150 : ae_vector t;
2151 : ae_vector s;
2152 : ae_int_t wlen;
2153 : ae_int_t woffs;
2154 : double u;
2155 : ae_int_t m;
2156 : ae_int_t l;
2157 : ae_int_t k;
2158 : ae_vector xgtmp;
2159 : ae_vector wgtmp;
2160 :
2161 0 : ae_frame_make(_state, &_frame_block);
2162 0 : memset(&_alpha, 0, sizeof(_alpha));
2163 0 : memset(&_beta, 0, sizeof(_beta));
2164 0 : memset(&ta, 0, sizeof(ta));
2165 0 : memset(&t, 0, sizeof(t));
2166 0 : memset(&s, 0, sizeof(s));
2167 0 : memset(&xgtmp, 0, sizeof(xgtmp));
2168 0 : memset(&wgtmp, 0, sizeof(wgtmp));
2169 0 : ae_vector_init_copy(&_alpha, alpha, _state, ae_true);
2170 0 : alpha = &_alpha;
2171 0 : ae_vector_init_copy(&_beta, beta, _state, ae_true);
2172 0 : beta = &_beta;
2173 0 : *info = 0;
2174 0 : ae_vector_clear(x);
2175 0 : ae_vector_clear(wkronrod);
2176 0 : ae_vector_clear(wgauss);
2177 0 : ae_vector_init(&ta, 0, DT_REAL, _state, ae_true);
2178 0 : ae_vector_init(&t, 0, DT_REAL, _state, ae_true);
2179 0 : ae_vector_init(&s, 0, DT_REAL, _state, ae_true);
2180 0 : ae_vector_init(&xgtmp, 0, DT_REAL, _state, ae_true);
2181 0 : ae_vector_init(&wgtmp, 0, DT_REAL, _state, ae_true);
2182 :
2183 0 : if( n%2!=1||n<3 )
2184 : {
2185 0 : *info = -1;
2186 0 : ae_frame_leave(_state);
2187 0 : return;
2188 : }
2189 0 : for(i=0; i<=ae_iceil((double)(3*(n/2))/(double)2, _state); i++)
2190 : {
2191 0 : if( ae_fp_less_eq(beta->ptr.p_double[i],(double)(0)) )
2192 : {
2193 0 : *info = -2;
2194 0 : ae_frame_leave(_state);
2195 0 : return;
2196 : }
2197 : }
2198 0 : *info = 1;
2199 :
2200 : /*
2201 : * from external conventions about N/Beta/Mu0 to internal
2202 : */
2203 0 : n = n/2;
2204 0 : beta->ptr.p_double[0] = mu0;
2205 :
2206 : /*
2207 : * Calculate Gauss nodes/weights, save them for later processing
2208 : */
2209 0 : gqgeneraterec(alpha, beta, mu0, n, info, &xgtmp, &wgtmp, _state);
2210 0 : if( *info<0 )
2211 : {
2212 0 : ae_frame_leave(_state);
2213 0 : return;
2214 : }
2215 :
2216 : /*
2217 : * Resize:
2218 : * * A from 0..floor(3*n/2) to 0..2*n
2219 : * * B from 0..ceil(3*n/2) to 0..2*n
2220 : */
2221 0 : ae_vector_set_length(&ta, ae_ifloor((double)(3*n)/(double)2, _state)+1, _state);
2222 0 : ae_v_move(&ta.ptr.p_double[0], 1, &alpha->ptr.p_double[0], 1, ae_v_len(0,ae_ifloor((double)(3*n)/(double)2, _state)));
2223 0 : ae_vector_set_length(alpha, 2*n+1, _state);
2224 0 : ae_v_move(&alpha->ptr.p_double[0], 1, &ta.ptr.p_double[0], 1, ae_v_len(0,ae_ifloor((double)(3*n)/(double)2, _state)));
2225 0 : for(i=ae_ifloor((double)(3*n)/(double)2, _state)+1; i<=2*n; i++)
2226 : {
2227 0 : alpha->ptr.p_double[i] = (double)(0);
2228 : }
2229 0 : ae_vector_set_length(&ta, ae_iceil((double)(3*n)/(double)2, _state)+1, _state);
2230 0 : ae_v_move(&ta.ptr.p_double[0], 1, &beta->ptr.p_double[0], 1, ae_v_len(0,ae_iceil((double)(3*n)/(double)2, _state)));
2231 0 : ae_vector_set_length(beta, 2*n+1, _state);
2232 0 : ae_v_move(&beta->ptr.p_double[0], 1, &ta.ptr.p_double[0], 1, ae_v_len(0,ae_iceil((double)(3*n)/(double)2, _state)));
2233 0 : for(i=ae_iceil((double)(3*n)/(double)2, _state)+1; i<=2*n; i++)
2234 : {
2235 0 : beta->ptr.p_double[i] = (double)(0);
2236 : }
2237 :
2238 : /*
2239 : * Initialize T, S
2240 : */
2241 0 : wlen = 2+n/2;
2242 0 : ae_vector_set_length(&t, wlen, _state);
2243 0 : ae_vector_set_length(&s, wlen, _state);
2244 0 : ae_vector_set_length(&ta, wlen, _state);
2245 0 : woffs = 1;
2246 0 : for(i=0; i<=wlen-1; i++)
2247 : {
2248 0 : t.ptr.p_double[i] = (double)(0);
2249 0 : s.ptr.p_double[i] = (double)(0);
2250 : }
2251 :
2252 : /*
2253 : * Algorithm from Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules", 1997.
2254 : */
2255 0 : t.ptr.p_double[woffs+0] = beta->ptr.p_double[n+1];
2256 0 : for(m=0; m<=n-2; m++)
2257 : {
2258 0 : u = (double)(0);
2259 0 : for(k=(m+1)/2; k>=0; k--)
2260 : {
2261 0 : l = m-k;
2262 0 : u = u+(alpha->ptr.p_double[k+n+1]-alpha->ptr.p_double[l])*t.ptr.p_double[woffs+k]+beta->ptr.p_double[k+n+1]*s.ptr.p_double[woffs+k-1]-beta->ptr.p_double[l]*s.ptr.p_double[woffs+k];
2263 0 : s.ptr.p_double[woffs+k] = u;
2264 : }
2265 0 : ae_v_move(&ta.ptr.p_double[0], 1, &t.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2266 0 : ae_v_move(&t.ptr.p_double[0], 1, &s.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2267 0 : ae_v_move(&s.ptr.p_double[0], 1, &ta.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2268 : }
2269 0 : for(j=n/2; j>=0; j--)
2270 : {
2271 0 : s.ptr.p_double[woffs+j] = s.ptr.p_double[woffs+j-1];
2272 : }
2273 0 : for(m=n-1; m<=2*n-3; m++)
2274 : {
2275 0 : u = (double)(0);
2276 0 : for(k=m+1-n; k<=(m-1)/2; k++)
2277 : {
2278 0 : l = m-k;
2279 0 : j = n-1-l;
2280 0 : u = u-(alpha->ptr.p_double[k+n+1]-alpha->ptr.p_double[l])*t.ptr.p_double[woffs+j]-beta->ptr.p_double[k+n+1]*s.ptr.p_double[woffs+j]+beta->ptr.p_double[l]*s.ptr.p_double[woffs+j+1];
2281 0 : s.ptr.p_double[woffs+j] = u;
2282 : }
2283 0 : if( m%2==0 )
2284 : {
2285 0 : k = m/2;
2286 0 : alpha->ptr.p_double[k+n+1] = alpha->ptr.p_double[k]+(s.ptr.p_double[woffs+j]-beta->ptr.p_double[k+n+1]*s.ptr.p_double[woffs+j+1])/t.ptr.p_double[woffs+j+1];
2287 : }
2288 : else
2289 : {
2290 0 : k = (m+1)/2;
2291 0 : beta->ptr.p_double[k+n+1] = s.ptr.p_double[woffs+j]/s.ptr.p_double[woffs+j+1];
2292 : }
2293 0 : ae_v_move(&ta.ptr.p_double[0], 1, &t.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2294 0 : ae_v_move(&t.ptr.p_double[0], 1, &s.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2295 0 : ae_v_move(&s.ptr.p_double[0], 1, &ta.ptr.p_double[0], 1, ae_v_len(0,wlen-1));
2296 : }
2297 0 : alpha->ptr.p_double[2*n] = alpha->ptr.p_double[n-1]-beta->ptr.p_double[2*n]*s.ptr.p_double[woffs+0]/t.ptr.p_double[woffs+0];
2298 :
2299 : /*
2300 : * calculation of Kronrod nodes and weights, unpacking of Gauss weights
2301 : */
2302 0 : gqgeneraterec(alpha, beta, mu0, 2*n+1, info, x, wkronrod, _state);
2303 0 : if( *info==-2 )
2304 : {
2305 0 : *info = -5;
2306 : }
2307 0 : if( *info<0 )
2308 : {
2309 0 : ae_frame_leave(_state);
2310 0 : return;
2311 : }
2312 0 : for(i=0; i<=2*n-1; i++)
2313 : {
2314 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
2315 : {
2316 0 : *info = -4;
2317 : }
2318 : }
2319 0 : if( *info<0 )
2320 : {
2321 0 : ae_frame_leave(_state);
2322 0 : return;
2323 : }
2324 0 : ae_vector_set_length(wgauss, 2*n+1, _state);
2325 0 : for(i=0; i<=2*n; i++)
2326 : {
2327 0 : wgauss->ptr.p_double[i] = (double)(0);
2328 : }
2329 0 : for(i=0; i<=n-1; i++)
2330 : {
2331 0 : wgauss->ptr.p_double[2*i+1] = wgtmp.ptr.p_double[i];
2332 : }
2333 0 : ae_frame_leave(_state);
2334 : }
2335 :
2336 :
2337 : /*************************************************************************
2338 : Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Legendre
2339 : quadrature with N points.
2340 :
2341 : GKQLegendreCalc (calculation) or GKQLegendreTbl (precomputed table) is
2342 : used depending on machine precision and number of nodes.
2343 :
2344 : INPUT PARAMETERS:
2345 : N - number of Kronrod nodes, must be odd number, >=3.
2346 :
2347 : OUTPUT PARAMETERS:
2348 : Info - error code:
2349 : * -4 an error was detected when calculating
2350 : weights/nodes. N is too large to obtain
2351 : weights/nodes with high enough accuracy.
2352 : Try to use multiple precision version.
2353 : * -3 internal eigenproblem solver hasn't converged
2354 : * -1 incorrect N was passed
2355 : * +1 OK
2356 : X - array[0..N-1] - array of quadrature nodes, ordered in
2357 : ascending order.
2358 : WKronrod - array[0..N-1] - Kronrod weights
2359 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
2360 : corresponding to extended Kronrod nodes).
2361 :
2362 :
2363 : -- ALGLIB --
2364 : Copyright 12.05.2009 by Bochkanov Sergey
2365 : *************************************************************************/
2366 0 : void gkqgenerategausslegendre(ae_int_t n,
2367 : ae_int_t* info,
2368 : /* Real */ ae_vector* x,
2369 : /* Real */ ae_vector* wkronrod,
2370 : /* Real */ ae_vector* wgauss,
2371 : ae_state *_state)
2372 : {
2373 : double eps;
2374 :
2375 0 : *info = 0;
2376 0 : ae_vector_clear(x);
2377 0 : ae_vector_clear(wkronrod);
2378 0 : ae_vector_clear(wgauss);
2379 :
2380 0 : if( ae_fp_greater(ae_machineepsilon,1.0E-32)&&(((((n==15||n==21)||n==31)||n==41)||n==51)||n==61) )
2381 : {
2382 0 : *info = 1;
2383 0 : gkqlegendretbl(n, x, wkronrod, wgauss, &eps, _state);
2384 : }
2385 : else
2386 : {
2387 0 : gkqlegendrecalc(n, info, x, wkronrod, wgauss, _state);
2388 : }
2389 0 : }
2390 :
2391 :
2392 : /*************************************************************************
2393 : Returns Gauss and Gauss-Kronrod nodes/weights for Gauss-Jacobi
2394 : quadrature on [-1,1] with weight function
2395 :
2396 : W(x)=Power(1-x,Alpha)*Power(1+x,Beta).
2397 :
2398 : INPUT PARAMETERS:
2399 : N - number of Kronrod nodes, must be odd number, >=3.
2400 : Alpha - power-law coefficient, Alpha>-1
2401 : Beta - power-law coefficient, Beta>-1
2402 :
2403 : OUTPUT PARAMETERS:
2404 : Info - error code:
2405 : * -5 no real and positive Gauss-Kronrod formula can
2406 : be created for such a weight function with a
2407 : given number of nodes.
2408 : * -4 an error was detected when calculating
2409 : weights/nodes. Alpha or Beta are too close
2410 : to -1 to obtain weights/nodes with high enough
2411 : accuracy, or, may be, N is too large. Try to
2412 : use multiple precision version.
2413 : * -3 internal eigenproblem solver hasn't converged
2414 : * -1 incorrect N was passed
2415 : * +1 OK
2416 : * +2 OK, but quadrature rule have exterior nodes,
2417 : x[0]<-1 or x[n-1]>+1
2418 : X - array[0..N-1] - array of quadrature nodes, ordered in
2419 : ascending order.
2420 : WKronrod - array[0..N-1] - Kronrod weights
2421 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
2422 : corresponding to extended Kronrod nodes).
2423 :
2424 :
2425 : -- ALGLIB --
2426 : Copyright 12.05.2009 by Bochkanov Sergey
2427 : *************************************************************************/
2428 0 : void gkqgenerategaussjacobi(ae_int_t n,
2429 : double alpha,
2430 : double beta,
2431 : ae_int_t* info,
2432 : /* Real */ ae_vector* x,
2433 : /* Real */ ae_vector* wkronrod,
2434 : /* Real */ ae_vector* wgauss,
2435 : ae_state *_state)
2436 : {
2437 : ae_frame _frame_block;
2438 : ae_int_t clen;
2439 : ae_vector a;
2440 : ae_vector b;
2441 : double alpha2;
2442 : double beta2;
2443 : double apb;
2444 : double t;
2445 : ae_int_t i;
2446 : double s;
2447 :
2448 0 : ae_frame_make(_state, &_frame_block);
2449 0 : memset(&a, 0, sizeof(a));
2450 0 : memset(&b, 0, sizeof(b));
2451 0 : *info = 0;
2452 0 : ae_vector_clear(x);
2453 0 : ae_vector_clear(wkronrod);
2454 0 : ae_vector_clear(wgauss);
2455 0 : ae_vector_init(&a, 0, DT_REAL, _state, ae_true);
2456 0 : ae_vector_init(&b, 0, DT_REAL, _state, ae_true);
2457 :
2458 0 : if( n%2!=1||n<3 )
2459 : {
2460 0 : *info = -1;
2461 0 : ae_frame_leave(_state);
2462 0 : return;
2463 : }
2464 0 : if( ae_fp_less_eq(alpha,(double)(-1))||ae_fp_less_eq(beta,(double)(-1)) )
2465 : {
2466 0 : *info = -1;
2467 0 : ae_frame_leave(_state);
2468 0 : return;
2469 : }
2470 0 : clen = ae_iceil((double)(3*(n/2))/(double)2, _state)+1;
2471 0 : ae_vector_set_length(&a, clen, _state);
2472 0 : ae_vector_set_length(&b, clen, _state);
2473 0 : for(i=0; i<=clen-1; i++)
2474 : {
2475 0 : a.ptr.p_double[i] = (double)(0);
2476 : }
2477 0 : apb = alpha+beta;
2478 0 : a.ptr.p_double[0] = (beta-alpha)/(apb+2);
2479 0 : t = (apb+1)*ae_log((double)(2), _state)+lngamma(alpha+1, &s, _state)+lngamma(beta+1, &s, _state)-lngamma(apb+2, &s, _state);
2480 0 : if( ae_fp_greater(t,ae_log(ae_maxrealnumber, _state)) )
2481 : {
2482 0 : *info = -4;
2483 0 : ae_frame_leave(_state);
2484 0 : return;
2485 : }
2486 0 : b.ptr.p_double[0] = ae_exp(t, _state);
2487 0 : if( clen>1 )
2488 : {
2489 0 : alpha2 = ae_sqr(alpha, _state);
2490 0 : beta2 = ae_sqr(beta, _state);
2491 0 : a.ptr.p_double[1] = (beta2-alpha2)/((apb+2)*(apb+4));
2492 0 : b.ptr.p_double[1] = 4*(alpha+1)*(beta+1)/((apb+3)*ae_sqr(apb+2, _state));
2493 0 : for(i=2; i<=clen-1; i++)
2494 : {
2495 0 : a.ptr.p_double[i] = 0.25*(beta2-alpha2)/(i*i*(1+0.5*apb/i)*(1+0.5*(apb+2)/i));
2496 0 : b.ptr.p_double[i] = 0.25*(1+alpha/i)*(1+beta/i)*(1+apb/i)/((1+0.5*(apb+1)/i)*(1+0.5*(apb-1)/i)*ae_sqr(1+0.5*apb/i, _state));
2497 : }
2498 : }
2499 0 : gkqgeneraterec(&a, &b, b.ptr.p_double[0], n, info, x, wkronrod, wgauss, _state);
2500 :
2501 : /*
2502 : * test basic properties to detect errors
2503 : */
2504 0 : if( *info>0 )
2505 : {
2506 0 : if( ae_fp_less(x->ptr.p_double[0],(double)(-1))||ae_fp_greater(x->ptr.p_double[n-1],(double)(1)) )
2507 : {
2508 0 : *info = 2;
2509 : }
2510 0 : for(i=0; i<=n-2; i++)
2511 : {
2512 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
2513 : {
2514 0 : *info = -4;
2515 : }
2516 : }
2517 : }
2518 0 : ae_frame_leave(_state);
2519 : }
2520 :
2521 :
2522 : /*************************************************************************
2523 : Returns Gauss and Gauss-Kronrod nodes for quadrature with N points.
2524 :
2525 : Reduction to tridiagonal eigenproblem is used.
2526 :
2527 : INPUT PARAMETERS:
2528 : N - number of Kronrod nodes, must be odd number, >=3.
2529 :
2530 : OUTPUT PARAMETERS:
2531 : Info - error code:
2532 : * -4 an error was detected when calculating
2533 : weights/nodes. N is too large to obtain
2534 : weights/nodes with high enough accuracy.
2535 : Try to use multiple precision version.
2536 : * -3 internal eigenproblem solver hasn't converged
2537 : * -1 incorrect N was passed
2538 : * +1 OK
2539 : X - array[0..N-1] - array of quadrature nodes, ordered in
2540 : ascending order.
2541 : WKronrod - array[0..N-1] - Kronrod weights
2542 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
2543 : corresponding to extended Kronrod nodes).
2544 :
2545 : -- ALGLIB --
2546 : Copyright 12.05.2009 by Bochkanov Sergey
2547 : *************************************************************************/
2548 0 : void gkqlegendrecalc(ae_int_t n,
2549 : ae_int_t* info,
2550 : /* Real */ ae_vector* x,
2551 : /* Real */ ae_vector* wkronrod,
2552 : /* Real */ ae_vector* wgauss,
2553 : ae_state *_state)
2554 : {
2555 : ae_frame _frame_block;
2556 : ae_vector alpha;
2557 : ae_vector beta;
2558 : ae_int_t alen;
2559 : ae_int_t blen;
2560 : double mu0;
2561 : ae_int_t k;
2562 : ae_int_t i;
2563 :
2564 0 : ae_frame_make(_state, &_frame_block);
2565 0 : memset(&alpha, 0, sizeof(alpha));
2566 0 : memset(&beta, 0, sizeof(beta));
2567 0 : *info = 0;
2568 0 : ae_vector_clear(x);
2569 0 : ae_vector_clear(wkronrod);
2570 0 : ae_vector_clear(wgauss);
2571 0 : ae_vector_init(&alpha, 0, DT_REAL, _state, ae_true);
2572 0 : ae_vector_init(&beta, 0, DT_REAL, _state, ae_true);
2573 :
2574 0 : if( n%2!=1||n<3 )
2575 : {
2576 0 : *info = -1;
2577 0 : ae_frame_leave(_state);
2578 0 : return;
2579 : }
2580 0 : mu0 = (double)(2);
2581 0 : alen = ae_ifloor((double)(3*(n/2))/(double)2, _state)+1;
2582 0 : blen = ae_iceil((double)(3*(n/2))/(double)2, _state)+1;
2583 0 : ae_vector_set_length(&alpha, alen, _state);
2584 0 : ae_vector_set_length(&beta, blen, _state);
2585 0 : for(k=0; k<=alen-1; k++)
2586 : {
2587 0 : alpha.ptr.p_double[k] = (double)(0);
2588 : }
2589 0 : beta.ptr.p_double[0] = (double)(2);
2590 0 : for(k=1; k<=blen-1; k++)
2591 : {
2592 0 : beta.ptr.p_double[k] = 1/(4-1/ae_sqr((double)(k), _state));
2593 : }
2594 0 : gkqgeneraterec(&alpha, &beta, mu0, n, info, x, wkronrod, wgauss, _state);
2595 :
2596 : /*
2597 : * test basic properties to detect errors
2598 : */
2599 0 : if( *info>0 )
2600 : {
2601 0 : if( ae_fp_less(x->ptr.p_double[0],(double)(-1))||ae_fp_greater(x->ptr.p_double[n-1],(double)(1)) )
2602 : {
2603 0 : *info = -4;
2604 : }
2605 0 : for(i=0; i<=n-2; i++)
2606 : {
2607 0 : if( ae_fp_greater_eq(x->ptr.p_double[i],x->ptr.p_double[i+1]) )
2608 : {
2609 0 : *info = -4;
2610 : }
2611 : }
2612 : }
2613 0 : ae_frame_leave(_state);
2614 : }
2615 :
2616 :
2617 : /*************************************************************************
2618 : Returns Gauss and Gauss-Kronrod nodes for quadrature with N points using
2619 : pre-calculated table. Nodes/weights were computed with accuracy up to
2620 : 1.0E-32 (if MPFR version of ALGLIB is used). In standard double precision
2621 : accuracy reduces to something about 2.0E-16 (depending on your compiler's
2622 : handling of long floating point constants).
2623 :
2624 : INPUT PARAMETERS:
2625 : N - number of Kronrod nodes.
2626 : N can be 15, 21, 31, 41, 51, 61.
2627 :
2628 : OUTPUT PARAMETERS:
2629 : X - array[0..N-1] - array of quadrature nodes, ordered in
2630 : ascending order.
2631 : WKronrod - array[0..N-1] - Kronrod weights
2632 : WGauss - array[0..N-1] - Gauss weights (interleaved with zeros
2633 : corresponding to extended Kronrod nodes).
2634 :
2635 :
2636 : -- ALGLIB --
2637 : Copyright 12.05.2009 by Bochkanov Sergey
2638 : *************************************************************************/
2639 0 : void gkqlegendretbl(ae_int_t n,
2640 : /* Real */ ae_vector* x,
2641 : /* Real */ ae_vector* wkronrod,
2642 : /* Real */ ae_vector* wgauss,
2643 : double* eps,
2644 : ae_state *_state)
2645 : {
2646 : ae_frame _frame_block;
2647 : ae_int_t i;
2648 : ae_int_t ng;
2649 : ae_vector p1;
2650 : ae_vector p2;
2651 : double tmp;
2652 :
2653 0 : ae_frame_make(_state, &_frame_block);
2654 0 : memset(&p1, 0, sizeof(p1));
2655 0 : memset(&p2, 0, sizeof(p2));
2656 0 : ae_vector_clear(x);
2657 0 : ae_vector_clear(wkronrod);
2658 0 : ae_vector_clear(wgauss);
2659 0 : *eps = 0;
2660 0 : ae_vector_init(&p1, 0, DT_INT, _state, ae_true);
2661 0 : ae_vector_init(&p2, 0, DT_INT, _state, ae_true);
2662 :
2663 :
2664 : /*
2665 : * these initializers are not really necessary,
2666 : * but without them compiler complains about uninitialized locals
2667 : */
2668 0 : ng = 0;
2669 :
2670 : /*
2671 : * Process
2672 : */
2673 0 : ae_assert(((((n==15||n==21)||n==31)||n==41)||n==51)||n==61, "GKQNodesTbl: incorrect N!", _state);
2674 0 : ae_vector_set_length(x, n, _state);
2675 0 : ae_vector_set_length(wkronrod, n, _state);
2676 0 : ae_vector_set_length(wgauss, n, _state);
2677 0 : for(i=0; i<=n-1; i++)
2678 : {
2679 0 : x->ptr.p_double[i] = (double)(0);
2680 0 : wkronrod->ptr.p_double[i] = (double)(0);
2681 0 : wgauss->ptr.p_double[i] = (double)(0);
2682 : }
2683 0 : *eps = ae_maxreal(ae_machineepsilon, 1.0E-32, _state);
2684 0 : if( n==15 )
2685 : {
2686 0 : ng = 4;
2687 0 : wgauss->ptr.p_double[0] = 0.129484966168869693270611432679082;
2688 0 : wgauss->ptr.p_double[1] = 0.279705391489276667901467771423780;
2689 0 : wgauss->ptr.p_double[2] = 0.381830050505118944950369775488975;
2690 0 : wgauss->ptr.p_double[3] = 0.417959183673469387755102040816327;
2691 0 : x->ptr.p_double[0] = 0.991455371120812639206854697526329;
2692 0 : x->ptr.p_double[1] = 0.949107912342758524526189684047851;
2693 0 : x->ptr.p_double[2] = 0.864864423359769072789712788640926;
2694 0 : x->ptr.p_double[3] = 0.741531185599394439863864773280788;
2695 0 : x->ptr.p_double[4] = 0.586087235467691130294144838258730;
2696 0 : x->ptr.p_double[5] = 0.405845151377397166906606412076961;
2697 0 : x->ptr.p_double[6] = 0.207784955007898467600689403773245;
2698 0 : x->ptr.p_double[7] = 0.000000000000000000000000000000000;
2699 0 : wkronrod->ptr.p_double[0] = 0.022935322010529224963732008058970;
2700 0 : wkronrod->ptr.p_double[1] = 0.063092092629978553290700663189204;
2701 0 : wkronrod->ptr.p_double[2] = 0.104790010322250183839876322541518;
2702 0 : wkronrod->ptr.p_double[3] = 0.140653259715525918745189590510238;
2703 0 : wkronrod->ptr.p_double[4] = 0.169004726639267902826583426598550;
2704 0 : wkronrod->ptr.p_double[5] = 0.190350578064785409913256402421014;
2705 0 : wkronrod->ptr.p_double[6] = 0.204432940075298892414161999234649;
2706 0 : wkronrod->ptr.p_double[7] = 0.209482141084727828012999174891714;
2707 : }
2708 0 : if( n==21 )
2709 : {
2710 0 : ng = 5;
2711 0 : wgauss->ptr.p_double[0] = 0.066671344308688137593568809893332;
2712 0 : wgauss->ptr.p_double[1] = 0.149451349150580593145776339657697;
2713 0 : wgauss->ptr.p_double[2] = 0.219086362515982043995534934228163;
2714 0 : wgauss->ptr.p_double[3] = 0.269266719309996355091226921569469;
2715 0 : wgauss->ptr.p_double[4] = 0.295524224714752870173892994651338;
2716 0 : x->ptr.p_double[0] = 0.995657163025808080735527280689003;
2717 0 : x->ptr.p_double[1] = 0.973906528517171720077964012084452;
2718 0 : x->ptr.p_double[2] = 0.930157491355708226001207180059508;
2719 0 : x->ptr.p_double[3] = 0.865063366688984510732096688423493;
2720 0 : x->ptr.p_double[4] = 0.780817726586416897063717578345042;
2721 0 : x->ptr.p_double[5] = 0.679409568299024406234327365114874;
2722 0 : x->ptr.p_double[6] = 0.562757134668604683339000099272694;
2723 0 : x->ptr.p_double[7] = 0.433395394129247190799265943165784;
2724 0 : x->ptr.p_double[8] = 0.294392862701460198131126603103866;
2725 0 : x->ptr.p_double[9] = 0.148874338981631210884826001129720;
2726 0 : x->ptr.p_double[10] = 0.000000000000000000000000000000000;
2727 0 : wkronrod->ptr.p_double[0] = 0.011694638867371874278064396062192;
2728 0 : wkronrod->ptr.p_double[1] = 0.032558162307964727478818972459390;
2729 0 : wkronrod->ptr.p_double[2] = 0.054755896574351996031381300244580;
2730 0 : wkronrod->ptr.p_double[3] = 0.075039674810919952767043140916190;
2731 0 : wkronrod->ptr.p_double[4] = 0.093125454583697605535065465083366;
2732 0 : wkronrod->ptr.p_double[5] = 0.109387158802297641899210590325805;
2733 0 : wkronrod->ptr.p_double[6] = 0.123491976262065851077958109831074;
2734 0 : wkronrod->ptr.p_double[7] = 0.134709217311473325928054001771707;
2735 0 : wkronrod->ptr.p_double[8] = 0.142775938577060080797094273138717;
2736 0 : wkronrod->ptr.p_double[9] = 0.147739104901338491374841515972068;
2737 0 : wkronrod->ptr.p_double[10] = 0.149445554002916905664936468389821;
2738 : }
2739 0 : if( n==31 )
2740 : {
2741 0 : ng = 8;
2742 0 : wgauss->ptr.p_double[0] = 0.030753241996117268354628393577204;
2743 0 : wgauss->ptr.p_double[1] = 0.070366047488108124709267416450667;
2744 0 : wgauss->ptr.p_double[2] = 0.107159220467171935011869546685869;
2745 0 : wgauss->ptr.p_double[3] = 0.139570677926154314447804794511028;
2746 0 : wgauss->ptr.p_double[4] = 0.166269205816993933553200860481209;
2747 0 : wgauss->ptr.p_double[5] = 0.186161000015562211026800561866423;
2748 0 : wgauss->ptr.p_double[6] = 0.198431485327111576456118326443839;
2749 0 : wgauss->ptr.p_double[7] = 0.202578241925561272880620199967519;
2750 0 : x->ptr.p_double[0] = 0.998002298693397060285172840152271;
2751 0 : x->ptr.p_double[1] = 0.987992518020485428489565718586613;
2752 0 : x->ptr.p_double[2] = 0.967739075679139134257347978784337;
2753 0 : x->ptr.p_double[3] = 0.937273392400705904307758947710209;
2754 0 : x->ptr.p_double[4] = 0.897264532344081900882509656454496;
2755 0 : x->ptr.p_double[5] = 0.848206583410427216200648320774217;
2756 0 : x->ptr.p_double[6] = 0.790418501442465932967649294817947;
2757 0 : x->ptr.p_double[7] = 0.724417731360170047416186054613938;
2758 0 : x->ptr.p_double[8] = 0.650996741297416970533735895313275;
2759 0 : x->ptr.p_double[9] = 0.570972172608538847537226737253911;
2760 0 : x->ptr.p_double[10] = 0.485081863640239680693655740232351;
2761 0 : x->ptr.p_double[11] = 0.394151347077563369897207370981045;
2762 0 : x->ptr.p_double[12] = 0.299180007153168812166780024266389;
2763 0 : x->ptr.p_double[13] = 0.201194093997434522300628303394596;
2764 0 : x->ptr.p_double[14] = 0.101142066918717499027074231447392;
2765 0 : x->ptr.p_double[15] = 0.000000000000000000000000000000000;
2766 0 : wkronrod->ptr.p_double[0] = 0.005377479872923348987792051430128;
2767 0 : wkronrod->ptr.p_double[1] = 0.015007947329316122538374763075807;
2768 0 : wkronrod->ptr.p_double[2] = 0.025460847326715320186874001019653;
2769 0 : wkronrod->ptr.p_double[3] = 0.035346360791375846222037948478360;
2770 0 : wkronrod->ptr.p_double[4] = 0.044589751324764876608227299373280;
2771 0 : wkronrod->ptr.p_double[5] = 0.053481524690928087265343147239430;
2772 0 : wkronrod->ptr.p_double[6] = 0.062009567800670640285139230960803;
2773 0 : wkronrod->ptr.p_double[7] = 0.069854121318728258709520077099147;
2774 0 : wkronrod->ptr.p_double[8] = 0.076849680757720378894432777482659;
2775 0 : wkronrod->ptr.p_double[9] = 0.083080502823133021038289247286104;
2776 0 : wkronrod->ptr.p_double[10] = 0.088564443056211770647275443693774;
2777 0 : wkronrod->ptr.p_double[11] = 0.093126598170825321225486872747346;
2778 0 : wkronrod->ptr.p_double[12] = 0.096642726983623678505179907627589;
2779 0 : wkronrod->ptr.p_double[13] = 0.099173598721791959332393173484603;
2780 0 : wkronrod->ptr.p_double[14] = 0.100769845523875595044946662617570;
2781 0 : wkronrod->ptr.p_double[15] = 0.101330007014791549017374792767493;
2782 : }
2783 0 : if( n==41 )
2784 : {
2785 0 : ng = 10;
2786 0 : wgauss->ptr.p_double[0] = 0.017614007139152118311861962351853;
2787 0 : wgauss->ptr.p_double[1] = 0.040601429800386941331039952274932;
2788 0 : wgauss->ptr.p_double[2] = 0.062672048334109063569506535187042;
2789 0 : wgauss->ptr.p_double[3] = 0.083276741576704748724758143222046;
2790 0 : wgauss->ptr.p_double[4] = 0.101930119817240435036750135480350;
2791 0 : wgauss->ptr.p_double[5] = 0.118194531961518417312377377711382;
2792 0 : wgauss->ptr.p_double[6] = 0.131688638449176626898494499748163;
2793 0 : wgauss->ptr.p_double[7] = 0.142096109318382051329298325067165;
2794 0 : wgauss->ptr.p_double[8] = 0.149172986472603746787828737001969;
2795 0 : wgauss->ptr.p_double[9] = 0.152753387130725850698084331955098;
2796 0 : x->ptr.p_double[0] = 0.998859031588277663838315576545863;
2797 0 : x->ptr.p_double[1] = 0.993128599185094924786122388471320;
2798 0 : x->ptr.p_double[2] = 0.981507877450250259193342994720217;
2799 0 : x->ptr.p_double[3] = 0.963971927277913791267666131197277;
2800 0 : x->ptr.p_double[4] = 0.940822633831754753519982722212443;
2801 0 : x->ptr.p_double[5] = 0.912234428251325905867752441203298;
2802 0 : x->ptr.p_double[6] = 0.878276811252281976077442995113078;
2803 0 : x->ptr.p_double[7] = 0.839116971822218823394529061701521;
2804 0 : x->ptr.p_double[8] = 0.795041428837551198350638833272788;
2805 0 : x->ptr.p_double[9] = 0.746331906460150792614305070355642;
2806 0 : x->ptr.p_double[10] = 0.693237656334751384805490711845932;
2807 0 : x->ptr.p_double[11] = 0.636053680726515025452836696226286;
2808 0 : x->ptr.p_double[12] = 0.575140446819710315342946036586425;
2809 0 : x->ptr.p_double[13] = 0.510867001950827098004364050955251;
2810 0 : x->ptr.p_double[14] = 0.443593175238725103199992213492640;
2811 0 : x->ptr.p_double[15] = 0.373706088715419560672548177024927;
2812 0 : x->ptr.p_double[16] = 0.301627868114913004320555356858592;
2813 0 : x->ptr.p_double[17] = 0.227785851141645078080496195368575;
2814 0 : x->ptr.p_double[18] = 0.152605465240922675505220241022678;
2815 0 : x->ptr.p_double[19] = 0.076526521133497333754640409398838;
2816 0 : x->ptr.p_double[20] = 0.000000000000000000000000000000000;
2817 0 : wkronrod->ptr.p_double[0] = 0.003073583718520531501218293246031;
2818 0 : wkronrod->ptr.p_double[1] = 0.008600269855642942198661787950102;
2819 0 : wkronrod->ptr.p_double[2] = 0.014626169256971252983787960308868;
2820 0 : wkronrod->ptr.p_double[3] = 0.020388373461266523598010231432755;
2821 0 : wkronrod->ptr.p_double[4] = 0.025882133604951158834505067096153;
2822 0 : wkronrod->ptr.p_double[5] = 0.031287306777032798958543119323801;
2823 0 : wkronrod->ptr.p_double[6] = 0.036600169758200798030557240707211;
2824 0 : wkronrod->ptr.p_double[7] = 0.041668873327973686263788305936895;
2825 0 : wkronrod->ptr.p_double[8] = 0.046434821867497674720231880926108;
2826 0 : wkronrod->ptr.p_double[9] = 0.050944573923728691932707670050345;
2827 0 : wkronrod->ptr.p_double[10] = 0.055195105348285994744832372419777;
2828 0 : wkronrod->ptr.p_double[11] = 0.059111400880639572374967220648594;
2829 0 : wkronrod->ptr.p_double[12] = 0.062653237554781168025870122174255;
2830 0 : wkronrod->ptr.p_double[13] = 0.065834597133618422111563556969398;
2831 0 : wkronrod->ptr.p_double[14] = 0.068648672928521619345623411885368;
2832 0 : wkronrod->ptr.p_double[15] = 0.071054423553444068305790361723210;
2833 0 : wkronrod->ptr.p_double[16] = 0.073030690332786667495189417658913;
2834 0 : wkronrod->ptr.p_double[17] = 0.074582875400499188986581418362488;
2835 0 : wkronrod->ptr.p_double[18] = 0.075704497684556674659542775376617;
2836 0 : wkronrod->ptr.p_double[19] = 0.076377867672080736705502835038061;
2837 0 : wkronrod->ptr.p_double[20] = 0.076600711917999656445049901530102;
2838 : }
2839 0 : if( n==51 )
2840 : {
2841 0 : ng = 13;
2842 0 : wgauss->ptr.p_double[0] = 0.011393798501026287947902964113235;
2843 0 : wgauss->ptr.p_double[1] = 0.026354986615032137261901815295299;
2844 0 : wgauss->ptr.p_double[2] = 0.040939156701306312655623487711646;
2845 0 : wgauss->ptr.p_double[3] = 0.054904695975835191925936891540473;
2846 0 : wgauss->ptr.p_double[4] = 0.068038333812356917207187185656708;
2847 0 : wgauss->ptr.p_double[5] = 0.080140700335001018013234959669111;
2848 0 : wgauss->ptr.p_double[6] = 0.091028261982963649811497220702892;
2849 0 : wgauss->ptr.p_double[7] = 0.100535949067050644202206890392686;
2850 0 : wgauss->ptr.p_double[8] = 0.108519624474263653116093957050117;
2851 0 : wgauss->ptr.p_double[9] = 0.114858259145711648339325545869556;
2852 0 : wgauss->ptr.p_double[10] = 0.119455763535784772228178126512901;
2853 0 : wgauss->ptr.p_double[11] = 0.122242442990310041688959518945852;
2854 0 : wgauss->ptr.p_double[12] = 0.123176053726715451203902873079050;
2855 0 : x->ptr.p_double[0] = 0.999262104992609834193457486540341;
2856 0 : x->ptr.p_double[1] = 0.995556969790498097908784946893902;
2857 0 : x->ptr.p_double[2] = 0.988035794534077247637331014577406;
2858 0 : x->ptr.p_double[3] = 0.976663921459517511498315386479594;
2859 0 : x->ptr.p_double[4] = 0.961614986425842512418130033660167;
2860 0 : x->ptr.p_double[5] = 0.942974571228974339414011169658471;
2861 0 : x->ptr.p_double[6] = 0.920747115281701561746346084546331;
2862 0 : x->ptr.p_double[7] = 0.894991997878275368851042006782805;
2863 0 : x->ptr.p_double[8] = 0.865847065293275595448996969588340;
2864 0 : x->ptr.p_double[9] = 0.833442628760834001421021108693570;
2865 0 : x->ptr.p_double[10] = 0.797873797998500059410410904994307;
2866 0 : x->ptr.p_double[11] = 0.759259263037357630577282865204361;
2867 0 : x->ptr.p_double[12] = 0.717766406813084388186654079773298;
2868 0 : x->ptr.p_double[13] = 0.673566368473468364485120633247622;
2869 0 : x->ptr.p_double[14] = 0.626810099010317412788122681624518;
2870 0 : x->ptr.p_double[15] = 0.577662930241222967723689841612654;
2871 0 : x->ptr.p_double[16] = 0.526325284334719182599623778158010;
2872 0 : x->ptr.p_double[17] = 0.473002731445714960522182115009192;
2873 0 : x->ptr.p_double[18] = 0.417885382193037748851814394594572;
2874 0 : x->ptr.p_double[19] = 0.361172305809387837735821730127641;
2875 0 : x->ptr.p_double[20] = 0.303089538931107830167478909980339;
2876 0 : x->ptr.p_double[21] = 0.243866883720988432045190362797452;
2877 0 : x->ptr.p_double[22] = 0.183718939421048892015969888759528;
2878 0 : x->ptr.p_double[23] = 0.122864692610710396387359818808037;
2879 0 : x->ptr.p_double[24] = 0.061544483005685078886546392366797;
2880 0 : x->ptr.p_double[25] = 0.000000000000000000000000000000000;
2881 0 : wkronrod->ptr.p_double[0] = 0.001987383892330315926507851882843;
2882 0 : wkronrod->ptr.p_double[1] = 0.005561932135356713758040236901066;
2883 0 : wkronrod->ptr.p_double[2] = 0.009473973386174151607207710523655;
2884 0 : wkronrod->ptr.p_double[3] = 0.013236229195571674813656405846976;
2885 0 : wkronrod->ptr.p_double[4] = 0.016847817709128298231516667536336;
2886 0 : wkronrod->ptr.p_double[5] = 0.020435371145882835456568292235939;
2887 0 : wkronrod->ptr.p_double[6] = 0.024009945606953216220092489164881;
2888 0 : wkronrod->ptr.p_double[7] = 0.027475317587851737802948455517811;
2889 0 : wkronrod->ptr.p_double[8] = 0.030792300167387488891109020215229;
2890 0 : wkronrod->ptr.p_double[9] = 0.034002130274329337836748795229551;
2891 0 : wkronrod->ptr.p_double[10] = 0.037116271483415543560330625367620;
2892 0 : wkronrod->ptr.p_double[11] = 0.040083825504032382074839284467076;
2893 0 : wkronrod->ptr.p_double[12] = 0.042872845020170049476895792439495;
2894 0 : wkronrod->ptr.p_double[13] = 0.045502913049921788909870584752660;
2895 0 : wkronrod->ptr.p_double[14] = 0.047982537138836713906392255756915;
2896 0 : wkronrod->ptr.p_double[15] = 0.050277679080715671963325259433440;
2897 0 : wkronrod->ptr.p_double[16] = 0.052362885806407475864366712137873;
2898 0 : wkronrod->ptr.p_double[17] = 0.054251129888545490144543370459876;
2899 0 : wkronrod->ptr.p_double[18] = 0.055950811220412317308240686382747;
2900 0 : wkronrod->ptr.p_double[19] = 0.057437116361567832853582693939506;
2901 0 : wkronrod->ptr.p_double[20] = 0.058689680022394207961974175856788;
2902 0 : wkronrod->ptr.p_double[21] = 0.059720340324174059979099291932562;
2903 0 : wkronrod->ptr.p_double[22] = 0.060539455376045862945360267517565;
2904 0 : wkronrod->ptr.p_double[23] = 0.061128509717053048305859030416293;
2905 0 : wkronrod->ptr.p_double[24] = 0.061471189871425316661544131965264;
2906 0 : wkronrod->ptr.p_double[25] = 0.061580818067832935078759824240055;
2907 : }
2908 0 : if( n==61 )
2909 : {
2910 0 : ng = 15;
2911 0 : wgauss->ptr.p_double[0] = 0.007968192496166605615465883474674;
2912 0 : wgauss->ptr.p_double[1] = 0.018466468311090959142302131912047;
2913 0 : wgauss->ptr.p_double[2] = 0.028784707883323369349719179611292;
2914 0 : wgauss->ptr.p_double[3] = 0.038799192569627049596801936446348;
2915 0 : wgauss->ptr.p_double[4] = 0.048402672830594052902938140422808;
2916 0 : wgauss->ptr.p_double[5] = 0.057493156217619066481721689402056;
2917 0 : wgauss->ptr.p_double[6] = 0.065974229882180495128128515115962;
2918 0 : wgauss->ptr.p_double[7] = 0.073755974737705206268243850022191;
2919 0 : wgauss->ptr.p_double[8] = 0.080755895229420215354694938460530;
2920 0 : wgauss->ptr.p_double[9] = 0.086899787201082979802387530715126;
2921 0 : wgauss->ptr.p_double[10] = 0.092122522237786128717632707087619;
2922 0 : wgauss->ptr.p_double[11] = 0.096368737174644259639468626351810;
2923 0 : wgauss->ptr.p_double[12] = 0.099593420586795267062780282103569;
2924 0 : wgauss->ptr.p_double[13] = 0.101762389748405504596428952168554;
2925 0 : wgauss->ptr.p_double[14] = 0.102852652893558840341285636705415;
2926 0 : x->ptr.p_double[0] = 0.999484410050490637571325895705811;
2927 0 : x->ptr.p_double[1] = 0.996893484074649540271630050918695;
2928 0 : x->ptr.p_double[2] = 0.991630996870404594858628366109486;
2929 0 : x->ptr.p_double[3] = 0.983668123279747209970032581605663;
2930 0 : x->ptr.p_double[4] = 0.973116322501126268374693868423707;
2931 0 : x->ptr.p_double[5] = 0.960021864968307512216871025581798;
2932 0 : x->ptr.p_double[6] = 0.944374444748559979415831324037439;
2933 0 : x->ptr.p_double[7] = 0.926200047429274325879324277080474;
2934 0 : x->ptr.p_double[8] = 0.905573307699907798546522558925958;
2935 0 : x->ptr.p_double[9] = 0.882560535792052681543116462530226;
2936 0 : x->ptr.p_double[10] = 0.857205233546061098958658510658944;
2937 0 : x->ptr.p_double[11] = 0.829565762382768397442898119732502;
2938 0 : x->ptr.p_double[12] = 0.799727835821839083013668942322683;
2939 0 : x->ptr.p_double[13] = 0.767777432104826194917977340974503;
2940 0 : x->ptr.p_double[14] = 0.733790062453226804726171131369528;
2941 0 : x->ptr.p_double[15] = 0.697850494793315796932292388026640;
2942 0 : x->ptr.p_double[16] = 0.660061064126626961370053668149271;
2943 0 : x->ptr.p_double[17] = 0.620526182989242861140477556431189;
2944 0 : x->ptr.p_double[18] = 0.579345235826361691756024932172540;
2945 0 : x->ptr.p_double[19] = 0.536624148142019899264169793311073;
2946 0 : x->ptr.p_double[20] = 0.492480467861778574993693061207709;
2947 0 : x->ptr.p_double[21] = 0.447033769538089176780609900322854;
2948 0 : x->ptr.p_double[22] = 0.400401254830394392535476211542661;
2949 0 : x->ptr.p_double[23] = 0.352704725530878113471037207089374;
2950 0 : x->ptr.p_double[24] = 0.304073202273625077372677107199257;
2951 0 : x->ptr.p_double[25] = 0.254636926167889846439805129817805;
2952 0 : x->ptr.p_double[26] = 0.204525116682309891438957671002025;
2953 0 : x->ptr.p_double[27] = 0.153869913608583546963794672743256;
2954 0 : x->ptr.p_double[28] = 0.102806937966737030147096751318001;
2955 0 : x->ptr.p_double[29] = 0.051471842555317695833025213166723;
2956 0 : x->ptr.p_double[30] = 0.000000000000000000000000000000000;
2957 0 : wkronrod->ptr.p_double[0] = 0.001389013698677007624551591226760;
2958 0 : wkronrod->ptr.p_double[1] = 0.003890461127099884051267201844516;
2959 0 : wkronrod->ptr.p_double[2] = 0.006630703915931292173319826369750;
2960 0 : wkronrod->ptr.p_double[3] = 0.009273279659517763428441146892024;
2961 0 : wkronrod->ptr.p_double[4] = 0.011823015253496341742232898853251;
2962 0 : wkronrod->ptr.p_double[5] = 0.014369729507045804812451432443580;
2963 0 : wkronrod->ptr.p_double[6] = 0.016920889189053272627572289420322;
2964 0 : wkronrod->ptr.p_double[7] = 0.019414141193942381173408951050128;
2965 0 : wkronrod->ptr.p_double[8] = 0.021828035821609192297167485738339;
2966 0 : wkronrod->ptr.p_double[9] = 0.024191162078080601365686370725232;
2967 0 : wkronrod->ptr.p_double[10] = 0.026509954882333101610601709335075;
2968 0 : wkronrod->ptr.p_double[11] = 0.028754048765041292843978785354334;
2969 0 : wkronrod->ptr.p_double[12] = 0.030907257562387762472884252943092;
2970 0 : wkronrod->ptr.p_double[13] = 0.032981447057483726031814191016854;
2971 0 : wkronrod->ptr.p_double[14] = 0.034979338028060024137499670731468;
2972 0 : wkronrod->ptr.p_double[15] = 0.036882364651821229223911065617136;
2973 0 : wkronrod->ptr.p_double[16] = 0.038678945624727592950348651532281;
2974 0 : wkronrod->ptr.p_double[17] = 0.040374538951535959111995279752468;
2975 0 : wkronrod->ptr.p_double[18] = 0.041969810215164246147147541285970;
2976 0 : wkronrod->ptr.p_double[19] = 0.043452539701356069316831728117073;
2977 0 : wkronrod->ptr.p_double[20] = 0.044814800133162663192355551616723;
2978 0 : wkronrod->ptr.p_double[21] = 0.046059238271006988116271735559374;
2979 0 : wkronrod->ptr.p_double[22] = 0.047185546569299153945261478181099;
2980 0 : wkronrod->ptr.p_double[23] = 0.048185861757087129140779492298305;
2981 0 : wkronrod->ptr.p_double[24] = 0.049055434555029778887528165367238;
2982 0 : wkronrod->ptr.p_double[25] = 0.049795683427074206357811569379942;
2983 0 : wkronrod->ptr.p_double[26] = 0.050405921402782346840893085653585;
2984 0 : wkronrod->ptr.p_double[27] = 0.050881795898749606492297473049805;
2985 0 : wkronrod->ptr.p_double[28] = 0.051221547849258772170656282604944;
2986 0 : wkronrod->ptr.p_double[29] = 0.051426128537459025933862879215781;
2987 0 : wkronrod->ptr.p_double[30] = 0.051494729429451567558340433647099;
2988 : }
2989 :
2990 : /*
2991 : * copy nodes
2992 : */
2993 0 : for(i=n-1; i>=n/2; i--)
2994 : {
2995 0 : x->ptr.p_double[i] = -x->ptr.p_double[n-1-i];
2996 : }
2997 :
2998 : /*
2999 : * copy Kronrod weights
3000 : */
3001 0 : for(i=n-1; i>=n/2; i--)
3002 : {
3003 0 : wkronrod->ptr.p_double[i] = wkronrod->ptr.p_double[n-1-i];
3004 : }
3005 :
3006 : /*
3007 : * copy Gauss weights
3008 : */
3009 0 : for(i=ng-1; i>=0; i--)
3010 : {
3011 0 : wgauss->ptr.p_double[n-2-2*i] = wgauss->ptr.p_double[i];
3012 0 : wgauss->ptr.p_double[1+2*i] = wgauss->ptr.p_double[i];
3013 : }
3014 0 : for(i=0; i<=n/2; i++)
3015 : {
3016 0 : wgauss->ptr.p_double[2*i] = (double)(0);
3017 : }
3018 :
3019 : /*
3020 : * reorder
3021 : */
3022 0 : tagsort(x, n, &p1, &p2, _state);
3023 0 : for(i=0; i<=n-1; i++)
3024 : {
3025 0 : tmp = wkronrod->ptr.p_double[i];
3026 0 : wkronrod->ptr.p_double[i] = wkronrod->ptr.p_double[p2.ptr.p_int[i]];
3027 0 : wkronrod->ptr.p_double[p2.ptr.p_int[i]] = tmp;
3028 0 : tmp = wgauss->ptr.p_double[i];
3029 0 : wgauss->ptr.p_double[i] = wgauss->ptr.p_double[p2.ptr.p_int[i]];
3030 0 : wgauss->ptr.p_double[p2.ptr.p_int[i]] = tmp;
3031 : }
3032 0 : ae_frame_leave(_state);
3033 0 : }
3034 :
3035 :
3036 : #endif
3037 : #if defined(AE_COMPILE_AUTOGK) || !defined(AE_PARTIAL_BUILD)
3038 :
3039 :
3040 : /*************************************************************************
3041 : Integration of a smooth function F(x) on a finite interval [a,b].
3042 :
3043 : Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
3044 : is calculated with accuracy close to the machine precision.
3045 :
3046 : Algorithm works well only with smooth integrands. It may be used with
3047 : continuous non-smooth integrands, but with less performance.
3048 :
3049 : It should never be used with integrands which have integrable singularities
3050 : at lower or upper limits - algorithm may crash. Use AutoGKSingular in such
3051 : cases.
3052 :
3053 : INPUT PARAMETERS:
3054 : A, B - interval boundaries (A<B, A=B or A>B)
3055 :
3056 : OUTPUT PARAMETERS
3057 : State - structure which stores algorithm state
3058 :
3059 : SEE ALSO
3060 : AutoGKSmoothW, AutoGKSingular, AutoGKResults.
3061 :
3062 :
3063 : -- ALGLIB --
3064 : Copyright 06.05.2009 by Bochkanov Sergey
3065 : *************************************************************************/
3066 0 : void autogksmooth(double a,
3067 : double b,
3068 : autogkstate* state,
3069 : ae_state *_state)
3070 : {
3071 :
3072 0 : _autogkstate_clear(state);
3073 :
3074 0 : ae_assert(ae_isfinite(a, _state), "AutoGKSmooth: A is not finite!", _state);
3075 0 : ae_assert(ae_isfinite(b, _state), "AutoGKSmooth: B is not finite!", _state);
3076 0 : autogksmoothw(a, b, 0.0, state, _state);
3077 0 : }
3078 :
3079 :
3080 : /*************************************************************************
3081 : Integration of a smooth function F(x) on a finite interval [a,b].
3082 :
3083 : This subroutine is same as AutoGKSmooth(), but it guarantees that interval
3084 : [a,b] is partitioned into subintervals which have width at most XWidth.
3085 :
3086 : Subroutine can be used when integrating nearly-constant function with
3087 : narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth
3088 : subroutine can overlook them.
3089 :
3090 : INPUT PARAMETERS:
3091 : A, B - interval boundaries (A<B, A=B or A>B)
3092 :
3093 : OUTPUT PARAMETERS
3094 : State - structure which stores algorithm state
3095 :
3096 : SEE ALSO
3097 : AutoGKSmooth, AutoGKSingular, AutoGKResults.
3098 :
3099 :
3100 : -- ALGLIB --
3101 : Copyright 06.05.2009 by Bochkanov Sergey
3102 : *************************************************************************/
3103 0 : void autogksmoothw(double a,
3104 : double b,
3105 : double xwidth,
3106 : autogkstate* state,
3107 : ae_state *_state)
3108 : {
3109 :
3110 0 : _autogkstate_clear(state);
3111 :
3112 0 : ae_assert(ae_isfinite(a, _state), "AutoGKSmoothW: A is not finite!", _state);
3113 0 : ae_assert(ae_isfinite(b, _state), "AutoGKSmoothW: B is not finite!", _state);
3114 0 : ae_assert(ae_isfinite(xwidth, _state), "AutoGKSmoothW: XWidth is not finite!", _state);
3115 0 : state->wrappermode = 0;
3116 0 : state->a = a;
3117 0 : state->b = b;
3118 0 : state->xwidth = xwidth;
3119 0 : state->needf = ae_false;
3120 0 : ae_vector_set_length(&state->rstate.ra, 10+1, _state);
3121 0 : state->rstate.stage = -1;
3122 0 : }
3123 :
3124 :
3125 : /*************************************************************************
3126 : Integration on a finite interval [A,B].
3127 : Integrand have integrable singularities at A/B.
3128 :
3129 : F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known
3130 : alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates
3131 : from below can be used (but these estimates should be greater than -1 too).
3132 :
3133 : One of alpha/beta variables (or even both alpha/beta) may be equal to 0,
3134 : which means than function F(x) is non-singular at A/B. Anyway (singular at
3135 : bounds or not), function F(x) is supposed to be continuous on (A,B).
3136 :
3137 : Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result
3138 : is calculated with accuracy close to the machine precision.
3139 :
3140 : INPUT PARAMETERS:
3141 : A, B - interval boundaries (A<B, A=B or A>B)
3142 : Alpha - power-law coefficient of the F(x) at A,
3143 : Alpha>-1
3144 : Beta - power-law coefficient of the F(x) at B,
3145 : Beta>-1
3146 :
3147 : OUTPUT PARAMETERS
3148 : State - structure which stores algorithm state
3149 :
3150 : SEE ALSO
3151 : AutoGKSmooth, AutoGKSmoothW, AutoGKResults.
3152 :
3153 :
3154 : -- ALGLIB --
3155 : Copyright 06.05.2009 by Bochkanov Sergey
3156 : *************************************************************************/
3157 0 : void autogksingular(double a,
3158 : double b,
3159 : double alpha,
3160 : double beta,
3161 : autogkstate* state,
3162 : ae_state *_state)
3163 : {
3164 :
3165 0 : _autogkstate_clear(state);
3166 :
3167 0 : ae_assert(ae_isfinite(a, _state), "AutoGKSingular: A is not finite!", _state);
3168 0 : ae_assert(ae_isfinite(b, _state), "AutoGKSingular: B is not finite!", _state);
3169 0 : ae_assert(ae_isfinite(alpha, _state), "AutoGKSingular: Alpha is not finite!", _state);
3170 0 : ae_assert(ae_isfinite(beta, _state), "AutoGKSingular: Beta is not finite!", _state);
3171 0 : state->wrappermode = 1;
3172 0 : state->a = a;
3173 0 : state->b = b;
3174 0 : state->alpha = alpha;
3175 0 : state->beta = beta;
3176 0 : state->xwidth = 0.0;
3177 0 : state->needf = ae_false;
3178 0 : ae_vector_set_length(&state->rstate.ra, 10+1, _state);
3179 0 : state->rstate.stage = -1;
3180 0 : }
3181 :
3182 :
3183 : /*************************************************************************
3184 :
3185 : -- ALGLIB --
3186 : Copyright 07.05.2009 by Bochkanov Sergey
3187 : *************************************************************************/
3188 0 : ae_bool autogkiteration(autogkstate* state, ae_state *_state)
3189 : {
3190 : double s;
3191 : double tmp;
3192 : double eps;
3193 : double a;
3194 : double b;
3195 : double x;
3196 : double t;
3197 : double alpha;
3198 : double beta;
3199 : double v1;
3200 : double v2;
3201 : ae_bool result;
3202 :
3203 :
3204 :
3205 : /*
3206 : * Reverse communication preparations
3207 : * I know it looks ugly, but it works the same way
3208 : * anywhere from C++ to Python.
3209 : *
3210 : * This code initializes locals by:
3211 : * * random values determined during code
3212 : * generation - on first subroutine call
3213 : * * values from previous call - on subsequent calls
3214 : */
3215 0 : if( state->rstate.stage>=0 )
3216 : {
3217 0 : s = state->rstate.ra.ptr.p_double[0];
3218 0 : tmp = state->rstate.ra.ptr.p_double[1];
3219 0 : eps = state->rstate.ra.ptr.p_double[2];
3220 0 : a = state->rstate.ra.ptr.p_double[3];
3221 0 : b = state->rstate.ra.ptr.p_double[4];
3222 0 : x = state->rstate.ra.ptr.p_double[5];
3223 0 : t = state->rstate.ra.ptr.p_double[6];
3224 0 : alpha = state->rstate.ra.ptr.p_double[7];
3225 0 : beta = state->rstate.ra.ptr.p_double[8];
3226 0 : v1 = state->rstate.ra.ptr.p_double[9];
3227 0 : v2 = state->rstate.ra.ptr.p_double[10];
3228 : }
3229 : else
3230 : {
3231 0 : s = 359;
3232 0 : tmp = -58;
3233 0 : eps = -919;
3234 0 : a = -909;
3235 0 : b = 81;
3236 0 : x = 255;
3237 0 : t = 74;
3238 0 : alpha = -788;
3239 0 : beta = 809;
3240 0 : v1 = 205;
3241 0 : v2 = -838;
3242 : }
3243 0 : if( state->rstate.stage==0 )
3244 : {
3245 0 : goto lbl_0;
3246 : }
3247 0 : if( state->rstate.stage==1 )
3248 : {
3249 0 : goto lbl_1;
3250 : }
3251 0 : if( state->rstate.stage==2 )
3252 : {
3253 0 : goto lbl_2;
3254 : }
3255 :
3256 : /*
3257 : * Routine body
3258 : */
3259 0 : eps = (double)(0);
3260 0 : a = state->a;
3261 0 : b = state->b;
3262 0 : alpha = state->alpha;
3263 0 : beta = state->beta;
3264 0 : state->terminationtype = -1;
3265 0 : state->nfev = 0;
3266 0 : state->nintervals = 0;
3267 :
3268 : /*
3269 : * smooth function at a finite interval
3270 : */
3271 0 : if( state->wrappermode!=0 )
3272 : {
3273 0 : goto lbl_3;
3274 : }
3275 :
3276 : /*
3277 : * special case
3278 : */
3279 0 : if( ae_fp_eq(a,b) )
3280 : {
3281 0 : state->terminationtype = 1;
3282 0 : state->v = (double)(0);
3283 0 : result = ae_false;
3284 0 : return result;
3285 : }
3286 :
3287 : /*
3288 : * general case
3289 : */
3290 0 : autogk_autogkinternalprepare(a, b, eps, state->xwidth, &state->internalstate, _state);
3291 0 : lbl_5:
3292 0 : if( !autogk_autogkinternaliteration(&state->internalstate, _state) )
3293 : {
3294 0 : goto lbl_6;
3295 : }
3296 0 : x = state->internalstate.x;
3297 0 : state->x = x;
3298 0 : state->xminusa = x-a;
3299 0 : state->bminusx = b-x;
3300 0 : state->needf = ae_true;
3301 0 : state->rstate.stage = 0;
3302 0 : goto lbl_rcomm;
3303 0 : lbl_0:
3304 0 : state->needf = ae_false;
3305 0 : state->nfev = state->nfev+1;
3306 0 : state->internalstate.f = state->f;
3307 0 : goto lbl_5;
3308 0 : lbl_6:
3309 0 : state->v = state->internalstate.r;
3310 0 : state->terminationtype = state->internalstate.info;
3311 0 : state->nintervals = state->internalstate.heapused;
3312 0 : result = ae_false;
3313 0 : return result;
3314 0 : lbl_3:
3315 :
3316 : /*
3317 : * function with power-law singularities at the ends of a finite interval
3318 : */
3319 0 : if( state->wrappermode!=1 )
3320 : {
3321 0 : goto lbl_7;
3322 : }
3323 :
3324 : /*
3325 : * test coefficients
3326 : */
3327 0 : if( ae_fp_less_eq(alpha,(double)(-1))||ae_fp_less_eq(beta,(double)(-1)) )
3328 : {
3329 0 : state->terminationtype = -1;
3330 0 : state->v = (double)(0);
3331 0 : result = ae_false;
3332 0 : return result;
3333 : }
3334 :
3335 : /*
3336 : * special cases
3337 : */
3338 0 : if( ae_fp_eq(a,b) )
3339 : {
3340 0 : state->terminationtype = 1;
3341 0 : state->v = (double)(0);
3342 0 : result = ae_false;
3343 0 : return result;
3344 : }
3345 :
3346 : /*
3347 : * reduction to general form
3348 : */
3349 0 : if( ae_fp_less(a,b) )
3350 : {
3351 0 : s = (double)(1);
3352 : }
3353 : else
3354 : {
3355 0 : s = (double)(-1);
3356 0 : tmp = a;
3357 0 : a = b;
3358 0 : b = tmp;
3359 0 : tmp = alpha;
3360 0 : alpha = beta;
3361 0 : beta = tmp;
3362 : }
3363 0 : alpha = ae_minreal(alpha, (double)(0), _state);
3364 0 : beta = ae_minreal(beta, (double)(0), _state);
3365 :
3366 : /*
3367 : * first, integrate left half of [a,b]:
3368 : * integral(f(x)dx, a, (b+a)/2) =
3369 : * = 1/(1+alpha) * integral(t^(-alpha/(1+alpha))*f(a+t^(1/(1+alpha)))dt, 0, (0.5*(b-a))^(1+alpha))
3370 : */
3371 0 : autogk_autogkinternalprepare((double)(0), ae_pow(0.5*(b-a), 1+alpha, _state), eps, state->xwidth, &state->internalstate, _state);
3372 0 : lbl_9:
3373 0 : if( !autogk_autogkinternaliteration(&state->internalstate, _state) )
3374 : {
3375 0 : goto lbl_10;
3376 : }
3377 :
3378 : /*
3379 : * Fill State.X, State.XMinusA, State.BMinusX.
3380 : * Latter two are filled correctly even if B<A.
3381 : */
3382 0 : x = state->internalstate.x;
3383 0 : t = ae_pow(x, 1/(1+alpha), _state);
3384 0 : state->x = a+t;
3385 0 : if( ae_fp_greater(s,(double)(0)) )
3386 : {
3387 0 : state->xminusa = t;
3388 0 : state->bminusx = b-(a+t);
3389 : }
3390 : else
3391 : {
3392 0 : state->xminusa = a+t-b;
3393 0 : state->bminusx = -t;
3394 : }
3395 0 : state->needf = ae_true;
3396 0 : state->rstate.stage = 1;
3397 0 : goto lbl_rcomm;
3398 0 : lbl_1:
3399 0 : state->needf = ae_false;
3400 0 : if( ae_fp_neq(alpha,(double)(0)) )
3401 : {
3402 0 : state->internalstate.f = state->f*ae_pow(x, -alpha/(1+alpha), _state)/(1+alpha);
3403 : }
3404 : else
3405 : {
3406 0 : state->internalstate.f = state->f;
3407 : }
3408 0 : state->nfev = state->nfev+1;
3409 0 : goto lbl_9;
3410 0 : lbl_10:
3411 0 : v1 = state->internalstate.r;
3412 0 : state->nintervals = state->nintervals+state->internalstate.heapused;
3413 :
3414 : /*
3415 : * then, integrate right half of [a,b]:
3416 : * integral(f(x)dx, (b+a)/2, b) =
3417 : * = 1/(1+beta) * integral(t^(-beta/(1+beta))*f(b-t^(1/(1+beta)))dt, 0, (0.5*(b-a))^(1+beta))
3418 : */
3419 0 : autogk_autogkinternalprepare((double)(0), ae_pow(0.5*(b-a), 1+beta, _state), eps, state->xwidth, &state->internalstate, _state);
3420 0 : lbl_11:
3421 0 : if( !autogk_autogkinternaliteration(&state->internalstate, _state) )
3422 : {
3423 0 : goto lbl_12;
3424 : }
3425 :
3426 : /*
3427 : * Fill State.X, State.XMinusA, State.BMinusX.
3428 : * Latter two are filled correctly (X-A, B-X) even if B<A.
3429 : */
3430 0 : x = state->internalstate.x;
3431 0 : t = ae_pow(x, 1/(1+beta), _state);
3432 0 : state->x = b-t;
3433 0 : if( ae_fp_greater(s,(double)(0)) )
3434 : {
3435 0 : state->xminusa = b-t-a;
3436 0 : state->bminusx = t;
3437 : }
3438 : else
3439 : {
3440 0 : state->xminusa = -t;
3441 0 : state->bminusx = a-(b-t);
3442 : }
3443 0 : state->needf = ae_true;
3444 0 : state->rstate.stage = 2;
3445 0 : goto lbl_rcomm;
3446 0 : lbl_2:
3447 0 : state->needf = ae_false;
3448 0 : if( ae_fp_neq(beta,(double)(0)) )
3449 : {
3450 0 : state->internalstate.f = state->f*ae_pow(x, -beta/(1+beta), _state)/(1+beta);
3451 : }
3452 : else
3453 : {
3454 0 : state->internalstate.f = state->f;
3455 : }
3456 0 : state->nfev = state->nfev+1;
3457 0 : goto lbl_11;
3458 0 : lbl_12:
3459 0 : v2 = state->internalstate.r;
3460 0 : state->nintervals = state->nintervals+state->internalstate.heapused;
3461 :
3462 : /*
3463 : * final result
3464 : */
3465 0 : state->v = s*(v1+v2);
3466 0 : state->terminationtype = 1;
3467 0 : result = ae_false;
3468 0 : return result;
3469 0 : lbl_7:
3470 0 : result = ae_false;
3471 0 : return result;
3472 :
3473 : /*
3474 : * Saving state
3475 : */
3476 0 : lbl_rcomm:
3477 0 : result = ae_true;
3478 0 : state->rstate.ra.ptr.p_double[0] = s;
3479 0 : state->rstate.ra.ptr.p_double[1] = tmp;
3480 0 : state->rstate.ra.ptr.p_double[2] = eps;
3481 0 : state->rstate.ra.ptr.p_double[3] = a;
3482 0 : state->rstate.ra.ptr.p_double[4] = b;
3483 0 : state->rstate.ra.ptr.p_double[5] = x;
3484 0 : state->rstate.ra.ptr.p_double[6] = t;
3485 0 : state->rstate.ra.ptr.p_double[7] = alpha;
3486 0 : state->rstate.ra.ptr.p_double[8] = beta;
3487 0 : state->rstate.ra.ptr.p_double[9] = v1;
3488 0 : state->rstate.ra.ptr.p_double[10] = v2;
3489 0 : return result;
3490 : }
3491 :
3492 :
3493 : /*************************************************************************
3494 : Adaptive integration results
3495 :
3496 : Called after AutoGKIteration returned False.
3497 :
3498 : Input parameters:
3499 : State - algorithm state (used by AutoGKIteration).
3500 :
3501 : Output parameters:
3502 : V - integral(f(x)dx,a,b)
3503 : Rep - optimization report (see AutoGKReport description)
3504 :
3505 : -- ALGLIB --
3506 : Copyright 14.11.2007 by Bochkanov Sergey
3507 : *************************************************************************/
3508 0 : void autogkresults(autogkstate* state,
3509 : double* v,
3510 : autogkreport* rep,
3511 : ae_state *_state)
3512 : {
3513 :
3514 0 : *v = 0;
3515 0 : _autogkreport_clear(rep);
3516 :
3517 0 : *v = state->v;
3518 0 : rep->terminationtype = state->terminationtype;
3519 0 : rep->nfev = state->nfev;
3520 0 : rep->nintervals = state->nintervals;
3521 0 : }
3522 :
3523 :
3524 : /*************************************************************************
3525 : Internal AutoGK subroutine
3526 : eps<0 - error
3527 : eps=0 - automatic eps selection
3528 :
3529 : width<0 - error
3530 : width=0 - no width requirements
3531 : *************************************************************************/
3532 0 : static void autogk_autogkinternalprepare(double a,
3533 : double b,
3534 : double eps,
3535 : double xwidth,
3536 : autogkinternalstate* state,
3537 : ae_state *_state)
3538 : {
3539 :
3540 :
3541 :
3542 : /*
3543 : * Save settings
3544 : */
3545 0 : state->a = a;
3546 0 : state->b = b;
3547 0 : state->eps = eps;
3548 0 : state->xwidth = xwidth;
3549 :
3550 : /*
3551 : * Prepare RComm structure
3552 : */
3553 0 : ae_vector_set_length(&state->rstate.ia, 3+1, _state);
3554 0 : ae_vector_set_length(&state->rstate.ra, 8+1, _state);
3555 0 : state->rstate.stage = -1;
3556 0 : }
3557 :
3558 :
3559 : /*************************************************************************
3560 : Internal AutoGK subroutine
3561 : *************************************************************************/
3562 0 : static ae_bool autogk_autogkinternaliteration(autogkinternalstate* state,
3563 : ae_state *_state)
3564 : {
3565 : double c1;
3566 : double c2;
3567 : ae_int_t i;
3568 : ae_int_t j;
3569 : double intg;
3570 : double intk;
3571 : double inta;
3572 : double v;
3573 : double ta;
3574 : double tb;
3575 : ae_int_t ns;
3576 : double qeps;
3577 : ae_int_t info;
3578 : ae_bool result;
3579 :
3580 :
3581 :
3582 : /*
3583 : * Reverse communication preparations
3584 : * I know it looks ugly, but it works the same way
3585 : * anywhere from C++ to Python.
3586 : *
3587 : * This code initializes locals by:
3588 : * * random values determined during code
3589 : * generation - on first subroutine call
3590 : * * values from previous call - on subsequent calls
3591 : */
3592 0 : if( state->rstate.stage>=0 )
3593 : {
3594 0 : i = state->rstate.ia.ptr.p_int[0];
3595 0 : j = state->rstate.ia.ptr.p_int[1];
3596 0 : ns = state->rstate.ia.ptr.p_int[2];
3597 0 : info = state->rstate.ia.ptr.p_int[3];
3598 0 : c1 = state->rstate.ra.ptr.p_double[0];
3599 0 : c2 = state->rstate.ra.ptr.p_double[1];
3600 0 : intg = state->rstate.ra.ptr.p_double[2];
3601 0 : intk = state->rstate.ra.ptr.p_double[3];
3602 0 : inta = state->rstate.ra.ptr.p_double[4];
3603 0 : v = state->rstate.ra.ptr.p_double[5];
3604 0 : ta = state->rstate.ra.ptr.p_double[6];
3605 0 : tb = state->rstate.ra.ptr.p_double[7];
3606 0 : qeps = state->rstate.ra.ptr.p_double[8];
3607 : }
3608 : else
3609 : {
3610 0 : i = 939;
3611 0 : j = -526;
3612 0 : ns = 763;
3613 0 : info = -541;
3614 0 : c1 = -698;
3615 0 : c2 = -900;
3616 0 : intg = -318;
3617 0 : intk = -940;
3618 0 : inta = 1016;
3619 0 : v = -229;
3620 0 : ta = -536;
3621 0 : tb = 487;
3622 0 : qeps = -115;
3623 : }
3624 0 : if( state->rstate.stage==0 )
3625 : {
3626 0 : goto lbl_0;
3627 : }
3628 0 : if( state->rstate.stage==1 )
3629 : {
3630 0 : goto lbl_1;
3631 : }
3632 0 : if( state->rstate.stage==2 )
3633 : {
3634 0 : goto lbl_2;
3635 : }
3636 :
3637 : /*
3638 : * Routine body
3639 : */
3640 :
3641 : /*
3642 : * initialize quadratures.
3643 : * use 15-point Gauss-Kronrod formula.
3644 : */
3645 0 : state->n = 15;
3646 0 : gkqgenerategausslegendre(state->n, &info, &state->qn, &state->wk, &state->wg, _state);
3647 0 : if( info<0 )
3648 : {
3649 0 : state->info = -5;
3650 0 : state->r = (double)(0);
3651 0 : result = ae_false;
3652 0 : return result;
3653 : }
3654 0 : ae_vector_set_length(&state->wr, state->n, _state);
3655 0 : for(i=0; i<=state->n-1; i++)
3656 : {
3657 0 : if( i==0 )
3658 : {
3659 0 : state->wr.ptr.p_double[i] = 0.5*ae_fabs(state->qn.ptr.p_double[1]-state->qn.ptr.p_double[0], _state);
3660 0 : continue;
3661 : }
3662 0 : if( i==state->n-1 )
3663 : {
3664 0 : state->wr.ptr.p_double[state->n-1] = 0.5*ae_fabs(state->qn.ptr.p_double[state->n-1]-state->qn.ptr.p_double[state->n-2], _state);
3665 0 : continue;
3666 : }
3667 0 : state->wr.ptr.p_double[i] = 0.5*ae_fabs(state->qn.ptr.p_double[i-1]-state->qn.ptr.p_double[i+1], _state);
3668 : }
3669 :
3670 : /*
3671 : * special case
3672 : */
3673 0 : if( ae_fp_eq(state->a,state->b) )
3674 : {
3675 0 : state->info = 1;
3676 0 : state->r = (double)(0);
3677 0 : result = ae_false;
3678 0 : return result;
3679 : }
3680 :
3681 : /*
3682 : * test parameters
3683 : */
3684 0 : if( ae_fp_less(state->eps,(double)(0))||ae_fp_less(state->xwidth,(double)(0)) )
3685 : {
3686 0 : state->info = -1;
3687 0 : state->r = (double)(0);
3688 0 : result = ae_false;
3689 0 : return result;
3690 : }
3691 0 : state->info = 1;
3692 0 : if( ae_fp_eq(state->eps,(double)(0)) )
3693 : {
3694 0 : state->eps = 100000*ae_machineepsilon;
3695 : }
3696 :
3697 : /*
3698 : * First, prepare heap
3699 : * * column 0 - absolute error
3700 : * * column 1 - integral of a F(x) (calculated using Kronrod extension nodes)
3701 : * * column 2 - integral of a |F(x)| (calculated using modified rect. method)
3702 : * * column 3 - left boundary of a subinterval
3703 : * * column 4 - right boundary of a subinterval
3704 : */
3705 0 : if( ae_fp_neq(state->xwidth,(double)(0)) )
3706 : {
3707 0 : goto lbl_3;
3708 : }
3709 :
3710 : /*
3711 : * no maximum width requirements
3712 : * start from one big subinterval
3713 : */
3714 0 : state->heapwidth = 5;
3715 0 : state->heapsize = 1;
3716 0 : state->heapused = 1;
3717 0 : ae_matrix_set_length(&state->heap, state->heapsize, state->heapwidth, _state);
3718 0 : c1 = 0.5*(state->b-state->a);
3719 0 : c2 = 0.5*(state->b+state->a);
3720 0 : intg = (double)(0);
3721 0 : intk = (double)(0);
3722 0 : inta = (double)(0);
3723 0 : i = 0;
3724 0 : lbl_5:
3725 0 : if( i>state->n-1 )
3726 : {
3727 0 : goto lbl_7;
3728 : }
3729 :
3730 : /*
3731 : * obtain F
3732 : */
3733 0 : state->x = c1*state->qn.ptr.p_double[i]+c2;
3734 0 : state->rstate.stage = 0;
3735 0 : goto lbl_rcomm;
3736 0 : lbl_0:
3737 0 : v = state->f;
3738 :
3739 : /*
3740 : * Gauss-Kronrod formula
3741 : */
3742 0 : intk = intk+v*state->wk.ptr.p_double[i];
3743 0 : if( i%2==1 )
3744 : {
3745 0 : intg = intg+v*state->wg.ptr.p_double[i];
3746 : }
3747 :
3748 : /*
3749 : * Integral |F(x)|
3750 : * Use rectangles method
3751 : */
3752 0 : inta = inta+ae_fabs(v, _state)*state->wr.ptr.p_double[i];
3753 0 : i = i+1;
3754 0 : goto lbl_5;
3755 0 : lbl_7:
3756 0 : intk = intk*(state->b-state->a)*0.5;
3757 0 : intg = intg*(state->b-state->a)*0.5;
3758 0 : inta = inta*(state->b-state->a)*0.5;
3759 0 : state->heap.ptr.pp_double[0][0] = ae_fabs(intg-intk, _state);
3760 0 : state->heap.ptr.pp_double[0][1] = intk;
3761 0 : state->heap.ptr.pp_double[0][2] = inta;
3762 0 : state->heap.ptr.pp_double[0][3] = state->a;
3763 0 : state->heap.ptr.pp_double[0][4] = state->b;
3764 0 : state->sumerr = state->heap.ptr.pp_double[0][0];
3765 0 : state->sumabs = ae_fabs(inta, _state);
3766 0 : goto lbl_4;
3767 0 : lbl_3:
3768 :
3769 : /*
3770 : * maximum subinterval should be no more than XWidth.
3771 : * so we create Ceil((B-A)/XWidth)+1 small subintervals
3772 : */
3773 0 : ns = ae_iceil(ae_fabs(state->b-state->a, _state)/state->xwidth, _state)+1;
3774 0 : state->heapsize = ns;
3775 0 : state->heapused = ns;
3776 0 : state->heapwidth = 5;
3777 0 : ae_matrix_set_length(&state->heap, state->heapsize, state->heapwidth, _state);
3778 0 : state->sumerr = (double)(0);
3779 0 : state->sumabs = (double)(0);
3780 0 : j = 0;
3781 0 : lbl_8:
3782 0 : if( j>ns-1 )
3783 : {
3784 0 : goto lbl_10;
3785 : }
3786 0 : ta = state->a+j*(state->b-state->a)/ns;
3787 0 : tb = state->a+(j+1)*(state->b-state->a)/ns;
3788 0 : c1 = 0.5*(tb-ta);
3789 0 : c2 = 0.5*(tb+ta);
3790 0 : intg = (double)(0);
3791 0 : intk = (double)(0);
3792 0 : inta = (double)(0);
3793 0 : i = 0;
3794 0 : lbl_11:
3795 0 : if( i>state->n-1 )
3796 : {
3797 0 : goto lbl_13;
3798 : }
3799 :
3800 : /*
3801 : * obtain F
3802 : */
3803 0 : state->x = c1*state->qn.ptr.p_double[i]+c2;
3804 0 : state->rstate.stage = 1;
3805 0 : goto lbl_rcomm;
3806 0 : lbl_1:
3807 0 : v = state->f;
3808 :
3809 : /*
3810 : * Gauss-Kronrod formula
3811 : */
3812 0 : intk = intk+v*state->wk.ptr.p_double[i];
3813 0 : if( i%2==1 )
3814 : {
3815 0 : intg = intg+v*state->wg.ptr.p_double[i];
3816 : }
3817 :
3818 : /*
3819 : * Integral |F(x)|
3820 : * Use rectangles method
3821 : */
3822 0 : inta = inta+ae_fabs(v, _state)*state->wr.ptr.p_double[i];
3823 0 : i = i+1;
3824 0 : goto lbl_11;
3825 0 : lbl_13:
3826 0 : intk = intk*(tb-ta)*0.5;
3827 0 : intg = intg*(tb-ta)*0.5;
3828 0 : inta = inta*(tb-ta)*0.5;
3829 0 : state->heap.ptr.pp_double[j][0] = ae_fabs(intg-intk, _state);
3830 0 : state->heap.ptr.pp_double[j][1] = intk;
3831 0 : state->heap.ptr.pp_double[j][2] = inta;
3832 0 : state->heap.ptr.pp_double[j][3] = ta;
3833 0 : state->heap.ptr.pp_double[j][4] = tb;
3834 0 : state->sumerr = state->sumerr+state->heap.ptr.pp_double[j][0];
3835 0 : state->sumabs = state->sumabs+ae_fabs(inta, _state);
3836 0 : j = j+1;
3837 0 : goto lbl_8;
3838 0 : lbl_10:
3839 0 : lbl_4:
3840 :
3841 : /*
3842 : * method iterations
3843 : */
3844 0 : lbl_14:
3845 : if( ae_false )
3846 : {
3847 : goto lbl_15;
3848 : }
3849 :
3850 : /*
3851 : * additional memory if needed
3852 : */
3853 0 : if( state->heapused==state->heapsize )
3854 : {
3855 0 : autogk_mheapresize(&state->heap, &state->heapsize, 4*state->heapsize, state->heapwidth, _state);
3856 : }
3857 :
3858 : /*
3859 : * TODO: every 20 iterations recalculate errors/sums
3860 : */
3861 0 : if( ae_fp_less_eq(state->sumerr,state->eps*state->sumabs)||state->heapused>=autogk_maxsubintervals )
3862 : {
3863 0 : state->r = (double)(0);
3864 0 : for(j=0; j<=state->heapused-1; j++)
3865 : {
3866 0 : state->r = state->r+state->heap.ptr.pp_double[j][1];
3867 : }
3868 0 : result = ae_false;
3869 0 : return result;
3870 : }
3871 :
3872 : /*
3873 : * Exclude interval with maximum absolute error
3874 : */
3875 0 : autogk_mheappop(&state->heap, state->heapused, state->heapwidth, _state);
3876 0 : state->sumerr = state->sumerr-state->heap.ptr.pp_double[state->heapused-1][0];
3877 0 : state->sumabs = state->sumabs-state->heap.ptr.pp_double[state->heapused-1][2];
3878 :
3879 : /*
3880 : * Divide interval, create subintervals
3881 : */
3882 0 : ta = state->heap.ptr.pp_double[state->heapused-1][3];
3883 0 : tb = state->heap.ptr.pp_double[state->heapused-1][4];
3884 0 : state->heap.ptr.pp_double[state->heapused-1][3] = ta;
3885 0 : state->heap.ptr.pp_double[state->heapused-1][4] = 0.5*(ta+tb);
3886 0 : state->heap.ptr.pp_double[state->heapused][3] = 0.5*(ta+tb);
3887 0 : state->heap.ptr.pp_double[state->heapused][4] = tb;
3888 0 : j = state->heapused-1;
3889 0 : lbl_16:
3890 0 : if( j>state->heapused )
3891 : {
3892 0 : goto lbl_18;
3893 : }
3894 0 : c1 = 0.5*(state->heap.ptr.pp_double[j][4]-state->heap.ptr.pp_double[j][3]);
3895 0 : c2 = 0.5*(state->heap.ptr.pp_double[j][4]+state->heap.ptr.pp_double[j][3]);
3896 0 : intg = (double)(0);
3897 0 : intk = (double)(0);
3898 0 : inta = (double)(0);
3899 0 : i = 0;
3900 0 : lbl_19:
3901 0 : if( i>state->n-1 )
3902 : {
3903 0 : goto lbl_21;
3904 : }
3905 :
3906 : /*
3907 : * F(x)
3908 : */
3909 0 : state->x = c1*state->qn.ptr.p_double[i]+c2;
3910 0 : state->rstate.stage = 2;
3911 0 : goto lbl_rcomm;
3912 0 : lbl_2:
3913 0 : v = state->f;
3914 :
3915 : /*
3916 : * Gauss-Kronrod formula
3917 : */
3918 0 : intk = intk+v*state->wk.ptr.p_double[i];
3919 0 : if( i%2==1 )
3920 : {
3921 0 : intg = intg+v*state->wg.ptr.p_double[i];
3922 : }
3923 :
3924 : /*
3925 : * Integral |F(x)|
3926 : * Use rectangles method
3927 : */
3928 0 : inta = inta+ae_fabs(v, _state)*state->wr.ptr.p_double[i];
3929 0 : i = i+1;
3930 0 : goto lbl_19;
3931 0 : lbl_21:
3932 0 : intk = intk*(state->heap.ptr.pp_double[j][4]-state->heap.ptr.pp_double[j][3])*0.5;
3933 0 : intg = intg*(state->heap.ptr.pp_double[j][4]-state->heap.ptr.pp_double[j][3])*0.5;
3934 0 : inta = inta*(state->heap.ptr.pp_double[j][4]-state->heap.ptr.pp_double[j][3])*0.5;
3935 0 : state->heap.ptr.pp_double[j][0] = ae_fabs(intg-intk, _state);
3936 0 : state->heap.ptr.pp_double[j][1] = intk;
3937 0 : state->heap.ptr.pp_double[j][2] = inta;
3938 0 : state->sumerr = state->sumerr+state->heap.ptr.pp_double[j][0];
3939 0 : state->sumabs = state->sumabs+state->heap.ptr.pp_double[j][2];
3940 0 : j = j+1;
3941 0 : goto lbl_16;
3942 0 : lbl_18:
3943 0 : autogk_mheappush(&state->heap, state->heapused-1, state->heapwidth, _state);
3944 0 : autogk_mheappush(&state->heap, state->heapused, state->heapwidth, _state);
3945 0 : state->heapused = state->heapused+1;
3946 0 : goto lbl_14;
3947 : lbl_15:
3948 : result = ae_false;
3949 : return result;
3950 :
3951 : /*
3952 : * Saving state
3953 : */
3954 0 : lbl_rcomm:
3955 0 : result = ae_true;
3956 0 : state->rstate.ia.ptr.p_int[0] = i;
3957 0 : state->rstate.ia.ptr.p_int[1] = j;
3958 0 : state->rstate.ia.ptr.p_int[2] = ns;
3959 0 : state->rstate.ia.ptr.p_int[3] = info;
3960 0 : state->rstate.ra.ptr.p_double[0] = c1;
3961 0 : state->rstate.ra.ptr.p_double[1] = c2;
3962 0 : state->rstate.ra.ptr.p_double[2] = intg;
3963 0 : state->rstate.ra.ptr.p_double[3] = intk;
3964 0 : state->rstate.ra.ptr.p_double[4] = inta;
3965 0 : state->rstate.ra.ptr.p_double[5] = v;
3966 0 : state->rstate.ra.ptr.p_double[6] = ta;
3967 0 : state->rstate.ra.ptr.p_double[7] = tb;
3968 0 : state->rstate.ra.ptr.p_double[8] = qeps;
3969 0 : return result;
3970 : }
3971 :
3972 :
3973 0 : static void autogk_mheappop(/* Real */ ae_matrix* heap,
3974 : ae_int_t heapsize,
3975 : ae_int_t heapwidth,
3976 : ae_state *_state)
3977 : {
3978 : ae_int_t i;
3979 : ae_int_t p;
3980 : double t;
3981 : ae_int_t maxcp;
3982 :
3983 :
3984 0 : if( heapsize==1 )
3985 : {
3986 0 : return;
3987 : }
3988 0 : for(i=0; i<=heapwidth-1; i++)
3989 : {
3990 0 : t = heap->ptr.pp_double[heapsize-1][i];
3991 0 : heap->ptr.pp_double[heapsize-1][i] = heap->ptr.pp_double[0][i];
3992 0 : heap->ptr.pp_double[0][i] = t;
3993 : }
3994 0 : p = 0;
3995 0 : while(2*p+1<heapsize-1)
3996 : {
3997 0 : maxcp = 2*p+1;
3998 0 : if( 2*p+2<heapsize-1 )
3999 : {
4000 0 : if( ae_fp_greater(heap->ptr.pp_double[2*p+2][0],heap->ptr.pp_double[2*p+1][0]) )
4001 : {
4002 0 : maxcp = 2*p+2;
4003 : }
4004 : }
4005 0 : if( ae_fp_less(heap->ptr.pp_double[p][0],heap->ptr.pp_double[maxcp][0]) )
4006 : {
4007 0 : for(i=0; i<=heapwidth-1; i++)
4008 : {
4009 0 : t = heap->ptr.pp_double[p][i];
4010 0 : heap->ptr.pp_double[p][i] = heap->ptr.pp_double[maxcp][i];
4011 0 : heap->ptr.pp_double[maxcp][i] = t;
4012 : }
4013 0 : p = maxcp;
4014 : }
4015 : else
4016 : {
4017 0 : break;
4018 : }
4019 : }
4020 : }
4021 :
4022 :
4023 0 : static void autogk_mheappush(/* Real */ ae_matrix* heap,
4024 : ae_int_t heapsize,
4025 : ae_int_t heapwidth,
4026 : ae_state *_state)
4027 : {
4028 : ae_int_t i;
4029 : ae_int_t p;
4030 : double t;
4031 : ae_int_t parent;
4032 :
4033 :
4034 0 : if( heapsize==0 )
4035 : {
4036 0 : return;
4037 : }
4038 0 : p = heapsize;
4039 0 : while(p!=0)
4040 : {
4041 0 : parent = (p-1)/2;
4042 0 : if( ae_fp_greater(heap->ptr.pp_double[p][0],heap->ptr.pp_double[parent][0]) )
4043 : {
4044 0 : for(i=0; i<=heapwidth-1; i++)
4045 : {
4046 0 : t = heap->ptr.pp_double[p][i];
4047 0 : heap->ptr.pp_double[p][i] = heap->ptr.pp_double[parent][i];
4048 0 : heap->ptr.pp_double[parent][i] = t;
4049 : }
4050 0 : p = parent;
4051 : }
4052 : else
4053 : {
4054 0 : break;
4055 : }
4056 : }
4057 : }
4058 :
4059 :
4060 0 : static void autogk_mheapresize(/* Real */ ae_matrix* heap,
4061 : ae_int_t* heapsize,
4062 : ae_int_t newheapsize,
4063 : ae_int_t heapwidth,
4064 : ae_state *_state)
4065 : {
4066 : ae_frame _frame_block;
4067 : ae_matrix tmp;
4068 : ae_int_t i;
4069 :
4070 0 : ae_frame_make(_state, &_frame_block);
4071 0 : memset(&tmp, 0, sizeof(tmp));
4072 0 : ae_matrix_init(&tmp, 0, 0, DT_REAL, _state, ae_true);
4073 :
4074 0 : ae_matrix_set_length(&tmp, *heapsize, heapwidth, _state);
4075 0 : for(i=0; i<=*heapsize-1; i++)
4076 : {
4077 0 : ae_v_move(&tmp.ptr.pp_double[i][0], 1, &heap->ptr.pp_double[i][0], 1, ae_v_len(0,heapwidth-1));
4078 : }
4079 0 : ae_matrix_set_length(heap, newheapsize, heapwidth, _state);
4080 0 : for(i=0; i<=*heapsize-1; i++)
4081 : {
4082 0 : ae_v_move(&heap->ptr.pp_double[i][0], 1, &tmp.ptr.pp_double[i][0], 1, ae_v_len(0,heapwidth-1));
4083 : }
4084 0 : *heapsize = newheapsize;
4085 0 : ae_frame_leave(_state);
4086 0 : }
4087 :
4088 :
4089 0 : void _autogkreport_init(void* _p, ae_state *_state, ae_bool make_automatic)
4090 : {
4091 0 : autogkreport *p = (autogkreport*)_p;
4092 0 : ae_touch_ptr((void*)p);
4093 0 : }
4094 :
4095 :
4096 0 : void _autogkreport_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic)
4097 : {
4098 0 : autogkreport *dst = (autogkreport*)_dst;
4099 0 : autogkreport *src = (autogkreport*)_src;
4100 0 : dst->terminationtype = src->terminationtype;
4101 0 : dst->nfev = src->nfev;
4102 0 : dst->nintervals = src->nintervals;
4103 0 : }
4104 :
4105 :
4106 0 : void _autogkreport_clear(void* _p)
4107 : {
4108 0 : autogkreport *p = (autogkreport*)_p;
4109 0 : ae_touch_ptr((void*)p);
4110 0 : }
4111 :
4112 :
4113 0 : void _autogkreport_destroy(void* _p)
4114 : {
4115 0 : autogkreport *p = (autogkreport*)_p;
4116 0 : ae_touch_ptr((void*)p);
4117 0 : }
4118 :
4119 :
4120 0 : void _autogkinternalstate_init(void* _p, ae_state *_state, ae_bool make_automatic)
4121 : {
4122 0 : autogkinternalstate *p = (autogkinternalstate*)_p;
4123 0 : ae_touch_ptr((void*)p);
4124 0 : ae_matrix_init(&p->heap, 0, 0, DT_REAL, _state, make_automatic);
4125 0 : ae_vector_init(&p->qn, 0, DT_REAL, _state, make_automatic);
4126 0 : ae_vector_init(&p->wg, 0, DT_REAL, _state, make_automatic);
4127 0 : ae_vector_init(&p->wk, 0, DT_REAL, _state, make_automatic);
4128 0 : ae_vector_init(&p->wr, 0, DT_REAL, _state, make_automatic);
4129 0 : _rcommstate_init(&p->rstate, _state, make_automatic);
4130 0 : }
4131 :
4132 :
4133 0 : void _autogkinternalstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic)
4134 : {
4135 0 : autogkinternalstate *dst = (autogkinternalstate*)_dst;
4136 0 : autogkinternalstate *src = (autogkinternalstate*)_src;
4137 0 : dst->a = src->a;
4138 0 : dst->b = src->b;
4139 0 : dst->eps = src->eps;
4140 0 : dst->xwidth = src->xwidth;
4141 0 : dst->x = src->x;
4142 0 : dst->f = src->f;
4143 0 : dst->info = src->info;
4144 0 : dst->r = src->r;
4145 0 : ae_matrix_init_copy(&dst->heap, &src->heap, _state, make_automatic);
4146 0 : dst->heapsize = src->heapsize;
4147 0 : dst->heapwidth = src->heapwidth;
4148 0 : dst->heapused = src->heapused;
4149 0 : dst->sumerr = src->sumerr;
4150 0 : dst->sumabs = src->sumabs;
4151 0 : ae_vector_init_copy(&dst->qn, &src->qn, _state, make_automatic);
4152 0 : ae_vector_init_copy(&dst->wg, &src->wg, _state, make_automatic);
4153 0 : ae_vector_init_copy(&dst->wk, &src->wk, _state, make_automatic);
4154 0 : ae_vector_init_copy(&dst->wr, &src->wr, _state, make_automatic);
4155 0 : dst->n = src->n;
4156 0 : _rcommstate_init_copy(&dst->rstate, &src->rstate, _state, make_automatic);
4157 0 : }
4158 :
4159 :
4160 0 : void _autogkinternalstate_clear(void* _p)
4161 : {
4162 0 : autogkinternalstate *p = (autogkinternalstate*)_p;
4163 0 : ae_touch_ptr((void*)p);
4164 0 : ae_matrix_clear(&p->heap);
4165 0 : ae_vector_clear(&p->qn);
4166 0 : ae_vector_clear(&p->wg);
4167 0 : ae_vector_clear(&p->wk);
4168 0 : ae_vector_clear(&p->wr);
4169 0 : _rcommstate_clear(&p->rstate);
4170 0 : }
4171 :
4172 :
4173 0 : void _autogkinternalstate_destroy(void* _p)
4174 : {
4175 0 : autogkinternalstate *p = (autogkinternalstate*)_p;
4176 0 : ae_touch_ptr((void*)p);
4177 0 : ae_matrix_destroy(&p->heap);
4178 0 : ae_vector_destroy(&p->qn);
4179 0 : ae_vector_destroy(&p->wg);
4180 0 : ae_vector_destroy(&p->wk);
4181 0 : ae_vector_destroy(&p->wr);
4182 0 : _rcommstate_destroy(&p->rstate);
4183 0 : }
4184 :
4185 :
4186 0 : void _autogkstate_init(void* _p, ae_state *_state, ae_bool make_automatic)
4187 : {
4188 0 : autogkstate *p = (autogkstate*)_p;
4189 0 : ae_touch_ptr((void*)p);
4190 0 : _autogkinternalstate_init(&p->internalstate, _state, make_automatic);
4191 0 : _rcommstate_init(&p->rstate, _state, make_automatic);
4192 0 : }
4193 :
4194 :
4195 0 : void _autogkstate_init_copy(void* _dst, void* _src, ae_state *_state, ae_bool make_automatic)
4196 : {
4197 0 : autogkstate *dst = (autogkstate*)_dst;
4198 0 : autogkstate *src = (autogkstate*)_src;
4199 0 : dst->a = src->a;
4200 0 : dst->b = src->b;
4201 0 : dst->alpha = src->alpha;
4202 0 : dst->beta = src->beta;
4203 0 : dst->xwidth = src->xwidth;
4204 0 : dst->x = src->x;
4205 0 : dst->xminusa = src->xminusa;
4206 0 : dst->bminusx = src->bminusx;
4207 0 : dst->needf = src->needf;
4208 0 : dst->f = src->f;
4209 0 : dst->wrappermode = src->wrappermode;
4210 0 : _autogkinternalstate_init_copy(&dst->internalstate, &src->internalstate, _state, make_automatic);
4211 0 : _rcommstate_init_copy(&dst->rstate, &src->rstate, _state, make_automatic);
4212 0 : dst->v = src->v;
4213 0 : dst->terminationtype = src->terminationtype;
4214 0 : dst->nfev = src->nfev;
4215 0 : dst->nintervals = src->nintervals;
4216 0 : }
4217 :
4218 :
4219 0 : void _autogkstate_clear(void* _p)
4220 : {
4221 0 : autogkstate *p = (autogkstate*)_p;
4222 0 : ae_touch_ptr((void*)p);
4223 0 : _autogkinternalstate_clear(&p->internalstate);
4224 0 : _rcommstate_clear(&p->rstate);
4225 0 : }
4226 :
4227 :
4228 0 : void _autogkstate_destroy(void* _p)
4229 : {
4230 0 : autogkstate *p = (autogkstate*)_p;
4231 0 : ae_touch_ptr((void*)p);
4232 0 : _autogkinternalstate_destroy(&p->internalstate);
4233 0 : _rcommstate_destroy(&p->rstate);
4234 0 : }
4235 :
4236 :
4237 : #endif
4238 :
4239 : }
4240 :
|
