Changeset 7294 for trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.4.0/ALGLIB-3.4.0/integration.cs
- Timestamp:
- 01/09/12 11:42:08 (13 years ago)
- Location:
- trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.4.0
- Files:
-
- 1 edited
- 1 copied
- 1 moved
-
. (moved) (moved from trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.1.0)
-
ALGLIB-3.4.0 (copied) (copied from trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.1.0/ALGLIB-3.1.0)
-
ALGLIB-3.4.0/integration.cs (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/3.4.0/ALGLIB-3.4.0/integration.cs
r4977 r7294 21 21 using System; 22 22 23 public partial class alglib24 {25 26 27 /*************************************************************************28 Integration report:29 * TerminationType = completetion code:30 * -5 non-convergence of Gauss-Kronrod nodes31 calculation subroutine.32 * -1 incorrect parameters were specified33 * 1 OK34 * Rep.NFEV countains number of function calculations35 * Rep.NIntervals contains number of intervals [a,b]36 was partitioned into.37 *************************************************************************/38 public class autogkreport39 {40 //41 // Public declarations42 //43 public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } }44 public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } }45 public int nintervals { get { return _innerobj.nintervals; } set { _innerobj.nintervals = value; } }46 47 public autogkreport()48 {49 _innerobj = new autogk.autogkreport();50 }51 52 //53 // Although some of declarations below are public, you should not use them54 // They are intended for internal use only55 //56 private autogk.autogkreport _innerobj;57 public autogk.autogkreport innerobj { get { return _innerobj; } }58 public autogkreport(autogk.autogkreport obj)59 {60 _innerobj = obj;61 }62 }63 64 65 /*************************************************************************66 This structure stores state of the integration algorithm.67 68 Although this class has public fields, they are not intended for external69 use. You should use ALGLIB functions to work with this class:70 * autogksmooth()/AutoGKSmoothW()/... to create objects71 * autogkintegrate() to begin integration72 * autogkresults() to get results73 *************************************************************************/74 public class autogkstate75 {76 //77 // Public declarations78 //79 public bool needf { get { return _innerobj.needf; } set { _innerobj.needf = value; } }80 public double x { get { return _innerobj.x; } set { _innerobj.x = value; } }81 public double xminusa { get { return _innerobj.xminusa; } set { _innerobj.xminusa = value; } }82 public double bminusx { get { return _innerobj.bminusx; } set { _innerobj.bminusx = value; } }83 public double f { get { return _innerobj.f; } set { _innerobj.f = value; } }84 85 public autogkstate()86 {87 _innerobj = new autogk.autogkstate();88 }89 90 //91 // Although some of declarations below are public, you should not use them92 // They are intended for internal use only93 //94 private autogk.autogkstate _innerobj;95 public autogk.autogkstate innerobj { get { return _innerobj; } }96 public autogkstate(autogk.autogkstate obj)97 {98 _innerobj = obj;99 }100 }101 102 /*************************************************************************103 Integration of a smooth function F(x) on a finite interval [a,b].104 105 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result106 is calculated with accuracy close to the machine precision.107 108 Algorithm works well only with smooth integrands. It may be used with109 continuous non-smooth integrands, but with less performance.110 111 It should never be used with integrands which have integrable singularities112 at lower or upper limits - algorithm may crash. Use AutoGKSingular in such113 cases.114 115 INPUT PARAMETERS:116 A, B - interval boundaries (A<B, A=B or A>B)117 118 OUTPUT PARAMETERS119 State - structure which stores algorithm state120 121 SEE ALSO122 AutoGKSmoothW, AutoGKSingular, AutoGKResults.123 124 125 -- ALGLIB --126 Copyright 06.05.2009 by Bochkanov Sergey127 *************************************************************************/128 public static void autogksmooth(double a, double b, out autogkstate state)129 {130 state = new autogkstate();131 autogk.autogksmooth(a, b, state.innerobj);132 return;133 }134 135 /*************************************************************************136 Integration of a smooth function F(x) on a finite interval [a,b].137 138 This subroutine is same as AutoGKSmooth(), but it guarantees that interval139 [a,b] is partitioned into subintervals which have width at most XWidth.140 141 Subroutine can be used when integrating nearly-constant function with142 narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth143 subroutine can overlook them.144 145 INPUT PARAMETERS:146 A, B - interval boundaries (A<B, A=B or A>B)147 148 OUTPUT PARAMETERS149 State - structure which stores algorithm state150 151 SEE ALSO152 AutoGKSmooth, AutoGKSingular, AutoGKResults.153 154 155 -- ALGLIB --156 Copyright 06.05.2009 by Bochkanov Sergey157 *************************************************************************/158 public static void autogksmoothw(double a, double b, double xwidth, out autogkstate state)159 {160 state = new autogkstate();161 autogk.autogksmoothw(a, b, xwidth, state.innerobj);162 return;163 }164 165 /*************************************************************************166 Integration on a finite interval [A,B].167 Integrand have integrable singularities at A/B.168 169 F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known170 alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates171 from below can be used (but these estimates should be greater than -1 too).172 173 One of alpha/beta variables (or even both alpha/beta) may be equal to 0,174 which means than function F(x) is non-singular at A/B. Anyway (singular at175 bounds or not), function F(x) is supposed to be continuous on (A,B).176 177 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result178 is calculated with accuracy close to the machine precision.179 180 INPUT PARAMETERS:181 A, B - interval boundaries (A<B, A=B or A>B)182 Alpha - power-law coefficient of the F(x) at A,183 Alpha>-1184 Beta - power-law coefficient of the F(x) at B,185 Beta>-1186 187 OUTPUT PARAMETERS188 State - structure which stores algorithm state189 190 SEE ALSO191 AutoGKSmooth, AutoGKSmoothW, AutoGKResults.192 193 194 -- ALGLIB --195 Copyright 06.05.2009 by Bochkanov Sergey196 *************************************************************************/197 public static void autogksingular(double a, double b, double alpha, double beta, out autogkstate state)198 {199 state = new autogkstate();200 autogk.autogksingular(a, b, alpha, beta, state.innerobj);201 return;202 }203 204 /*************************************************************************205 This function provides reverse communication interface206 Reverse communication interface is not documented or recommended to use.207 See below for functions which provide better documented API208 *************************************************************************/209 public static bool autogkiteration(autogkstate state)210 {211 212 bool result = autogk.autogkiteration(state.innerobj);213 return result;214 }215 216 217 /*************************************************************************218 This function is used to launcn iterations of ODE solver219 220 It accepts following parameters:221 diff - callback which calculates dy/dx for given y and x222 obj - optional object which is passed to diff; can be NULL223 224 225 -- ALGLIB --226 Copyright 07.05.2009 by Bochkanov Sergey227 228 *************************************************************************/229 public static void autogkintegrate(autogkstate state, integrator1_func func, object obj)230 {231 if( func==null )232 throw new alglibexception("ALGLIB: error in 'autogkintegrate()' (func is null)");233 while( alglib.autogkiteration(state) )234 {235 if( state.needf )236 {237 func(state.innerobj.x, state.innerobj.xminusa, state.innerobj.bminusx, ref state.innerobj.f, obj);238 continue;239 }240 throw new alglibexception("ALGLIB: unexpected error in 'autogksolve'");241 }242 }243 244 /*************************************************************************245 Adaptive integration results246 247 Called after AutoGKIteration returned False.248 249 Input parameters:250 State - algorithm state (used by AutoGKIteration).251 252 Output parameters:253 V - integral(f(x)dx,a,b)254 Rep - optimization report (see AutoGKReport description)255 256 -- ALGLIB --257 Copyright 14.11.2007 by Bochkanov Sergey258 *************************************************************************/259 public static void autogkresults(autogkstate state, out double v, out autogkreport rep)260 {261 v = 0;262 rep = new autogkreport();263 autogk.autogkresults(state.innerobj, ref v, rep.innerobj);264 return;265 }266 267 }268 23 public partial class alglib 269 24 { … … 764 519 public partial class alglib 765 520 { 766 public class autogk 521 522 523 /************************************************************************* 524 Integration report: 525 * TerminationType = completetion code: 526 * -5 non-convergence of Gauss-Kronrod nodes 527 calculation subroutine. 528 * -1 incorrect parameters were specified 529 * 1 OK 530 * Rep.NFEV countains number of function calculations 531 * Rep.NIntervals contains number of intervals [a,b] 532 was partitioned into. 533 *************************************************************************/ 534 public class autogkreport 767 535 { 768 /************************************************************************* 769 Integration report: 770 * TerminationType = completetion code: 771 * -5 non-convergence of Gauss-Kronrod nodes 772 calculation subroutine. 773 * -1 incorrect parameters were specified 774 * 1 OK 775 * Rep.NFEV countains number of function calculations 776 * Rep.NIntervals contains number of intervals [a,b] 777 was partitioned into. 778 *************************************************************************/ 779 public class autogkreport 536 // 537 // Public declarations 538 // 539 public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } } 540 public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } } 541 public int nintervals { get { return _innerobj.nintervals; } set { _innerobj.nintervals = value; } } 542 543 public autogkreport() 780 544 { 781 public int terminationtype; 782 public int nfev; 783 public int nintervals; 784 }; 785 786 787 public class autogkinternalstate 545 _innerobj = new autogk.autogkreport(); 546 } 547 548 // 549 // Although some of declarations below are public, you should not use them 550 // They are intended for internal use only 551 // 552 private autogk.autogkreport _innerobj; 553 public autogk.autogkreport innerobj { get { return _innerobj; } } 554 public autogkreport(autogk.autogkreport obj) 788 555 { 789 public double a; 790 public double b; 791 public double eps; 792 public double xwidth; 793 public double x; 794 public double f; 795 public int info; 796 public double r; 797 public double[,] heap; 798 public int heapsize; 799 public int heapwidth; 800 public int heapused; 801 public double sumerr; 802 public double sumabs; 803 public double[] qn; 804 public double[] wg; 805 public double[] wk; 806 public double[] wr; 807 public int n; 808 public rcommstate rstate; 809 public autogkinternalstate() 810 { 811 heap = new double[0,0]; 812 qn = new double[0]; 813 wg = new double[0]; 814 wk = new double[0]; 815 wr = new double[0]; 816 rstate = new rcommstate(); 817 } 818 }; 819 820 821 /************************************************************************* 822 This structure stores state of the integration algorithm. 823 824 Although this class has public fields, they are not intended for external 825 use. You should use ALGLIB functions to work with this class: 826 * autogksmooth()/AutoGKSmoothW()/... to create objects 827 * autogkintegrate() to begin integration 828 * autogkresults() to get results 829 *************************************************************************/ 830 public class autogkstate 556 _innerobj = obj; 557 } 558 } 559 560 561 /************************************************************************* 562 This structure stores state of the integration algorithm. 563 564 Although this class has public fields, they are not intended for external 565 use. You should use ALGLIB functions to work with this class: 566 * autogksmooth()/AutoGKSmoothW()/... to create objects 567 * autogkintegrate() to begin integration 568 * autogkresults() to get results 569 *************************************************************************/ 570 public class autogkstate 571 { 572 // 573 // Public declarations 574 // 575 public bool needf { get { return _innerobj.needf; } set { _innerobj.needf = value; } } 576 public double x { get { return _innerobj.x; } set { _innerobj.x = value; } } 577 public double xminusa { get { return _innerobj.xminusa; } set { _innerobj.xminusa = value; } } 578 public double bminusx { get { return _innerobj.bminusx; } set { _innerobj.bminusx = value; } } 579 public double f { get { return _innerobj.f; } set { _innerobj.f = value; } } 580 581 public autogkstate() 831 582 { 832 public double a; 833 public double b; 834 public double alpha; 835 public double beta; 836 public double xwidth; 837 public double x; 838 public double xminusa; 839 public double bminusx; 840 public bool needf; 841 public double f; 842 public int wrappermode; 843 public autogkinternalstate internalstate; 844 public rcommstate rstate; 845 public double v; 846 public int terminationtype; 847 public int nfev; 848 public int nintervals; 849 public autogkstate() 850 { 851 internalstate = new autogkinternalstate(); 852 rstate = new rcommstate(); 853 } 854 }; 855 856 857 858 859 /************************************************************************* 860 Integration of a smooth function F(x) on a finite interval [a,b]. 861 862 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result 863 is calculated with accuracy close to the machine precision. 864 865 Algorithm works well only with smooth integrands. It may be used with 866 continuous non-smooth integrands, but with less performance. 867 868 It should never be used with integrands which have integrable singularities 869 at lower or upper limits - algorithm may crash. Use AutoGKSingular in such 870 cases. 871 872 INPUT PARAMETERS: 873 A, B - interval boundaries (A<B, A=B or A>B) 874 875 OUTPUT PARAMETERS 876 State - structure which stores algorithm state 877 878 SEE ALSO 879 AutoGKSmoothW, AutoGKSingular, AutoGKResults. 880 881 882 -- ALGLIB -- 883 Copyright 06.05.2009 by Bochkanov Sergey 884 *************************************************************************/ 885 public static void autogksmooth(double a, 886 double b, 887 autogkstate state) 583 _innerobj = new autogk.autogkstate(); 584 } 585 586 // 587 // Although some of declarations below are public, you should not use them 588 // They are intended for internal use only 589 // 590 private autogk.autogkstate _innerobj; 591 public autogk.autogkstate innerobj { get { return _innerobj; } } 592 public autogkstate(autogk.autogkstate obj) 888 593 { 889 ap.assert(math.isfinite(a), "AutoGKSmooth: A is not finite!"); 890 ap.assert(math.isfinite(b), "AutoGKSmooth: B is not finite!"); 891 autogksmoothw(a, b, 0.0, state); 594 _innerobj = obj; 892 595 } 893 894 895 /************************************************************************* 896 Integration of a smooth function F(x) on a finite interval [a,b]. 897 898 This subroutine is same as AutoGKSmooth(), but it guarantees that interval 899 [a,b] is partitioned into subintervals which have width at most XWidth. 900 901 Subroutine can be used when integrating nearly-constant function with 902 narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth 903 subroutine can overlook them. 904 905 INPUT PARAMETERS: 906 A, B - interval boundaries (A<B, A=B or A>B) 907 908 OUTPUT PARAMETERS 909 State - structure which stores algorithm state 910 911 SEE ALSO 912 AutoGKSmooth, AutoGKSingular, AutoGKResults. 913 914 915 -- ALGLIB -- 916 Copyright 06.05.2009 by Bochkanov Sergey 917 *************************************************************************/ 918 public static void autogksmoothw(double a, 919 double b, 920 double xwidth, 921 autogkstate state) 596 } 597 598 /************************************************************************* 599 Integration of a smooth function F(x) on a finite interval [a,b]. 600 601 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result 602 is calculated with accuracy close to the machine precision. 603 604 Algorithm works well only with smooth integrands. It may be used with 605 continuous non-smooth integrands, but with less performance. 606 607 It should never be used with integrands which have integrable singularities 608 at lower or upper limits - algorithm may crash. Use AutoGKSingular in such 609 cases. 610 611 INPUT PARAMETERS: 612 A, B - interval boundaries (A<B, A=B or A>B) 613 614 OUTPUT PARAMETERS 615 State - structure which stores algorithm state 616 617 SEE ALSO 618 AutoGKSmoothW, AutoGKSingular, AutoGKResults. 619 620 621 -- ALGLIB -- 622 Copyright 06.05.2009 by Bochkanov Sergey 623 *************************************************************************/ 624 public static void autogksmooth(double a, double b, out autogkstate state) 625 { 626 state = new autogkstate(); 627 autogk.autogksmooth(a, b, state.innerobj); 628 return; 629 } 630 631 /************************************************************************* 632 Integration of a smooth function F(x) on a finite interval [a,b]. 633 634 This subroutine is same as AutoGKSmooth(), but it guarantees that interval 635 [a,b] is partitioned into subintervals which have width at most XWidth. 636 637 Subroutine can be used when integrating nearly-constant function with 638 narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth 639 subroutine can overlook them. 640 641 INPUT PARAMETERS: 642 A, B - interval boundaries (A<B, A=B or A>B) 643 644 OUTPUT PARAMETERS 645 State - structure which stores algorithm state 646 647 SEE ALSO 648 AutoGKSmooth, AutoGKSingular, AutoGKResults. 649 650 651 -- ALGLIB -- 652 Copyright 06.05.2009 by Bochkanov Sergey 653 *************************************************************************/ 654 public static void autogksmoothw(double a, double b, double xwidth, out autogkstate state) 655 { 656 state = new autogkstate(); 657 autogk.autogksmoothw(a, b, xwidth, state.innerobj); 658 return; 659 } 660 661 /************************************************************************* 662 Integration on a finite interval [A,B]. 663 Integrand have integrable singularities at A/B. 664 665 F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known 666 alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates 667 from below can be used (but these estimates should be greater than -1 too). 668 669 One of alpha/beta variables (or even both alpha/beta) may be equal to 0, 670 which means than function F(x) is non-singular at A/B. Anyway (singular at 671 bounds or not), function F(x) is supposed to be continuous on (A,B). 672 673 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result 674 is calculated with accuracy close to the machine precision. 675 676 INPUT PARAMETERS: 677 A, B - interval boundaries (A<B, A=B or A>B) 678 Alpha - power-law coefficient of the F(x) at A, 679 Alpha>-1 680 Beta - power-law coefficient of the F(x) at B, 681 Beta>-1 682 683 OUTPUT PARAMETERS 684 State - structure which stores algorithm state 685 686 SEE ALSO 687 AutoGKSmooth, AutoGKSmoothW, AutoGKResults. 688 689 690 -- ALGLIB -- 691 Copyright 06.05.2009 by Bochkanov Sergey 692 *************************************************************************/ 693 public static void autogksingular(double a, double b, double alpha, double beta, out autogkstate state) 694 { 695 state = new autogkstate(); 696 autogk.autogksingular(a, b, alpha, beta, state.innerobj); 697 return; 698 } 699 700 /************************************************************************* 701 This function provides reverse communication interface 702 Reverse communication interface is not documented or recommended to use. 703 See below for functions which provide better documented API 704 *************************************************************************/ 705 public static bool autogkiteration(autogkstate state) 706 { 707 708 bool result = autogk.autogkiteration(state.innerobj); 709 return result; 710 } 711 712 713 /************************************************************************* 714 This function is used to launcn iterations of ODE solver 715 716 It accepts following parameters: 717 diff - callback which calculates dy/dx for given y and x 718 obj - optional object which is passed to diff; can be NULL 719 720 721 -- ALGLIB -- 722 Copyright 07.05.2009 by Bochkanov Sergey 723 724 *************************************************************************/ 725 public static void autogkintegrate(autogkstate state, integrator1_func func, object obj) 726 { 727 if( func==null ) 728 throw new alglibexception("ALGLIB: error in 'autogkintegrate()' (func is null)"); 729 while( alglib.autogkiteration(state) ) 922 730 { 923 ap.assert(math.isfinite(a), "AutoGKSmoothW: A is not finite!"); 924 ap.assert(math.isfinite(b), "AutoGKSmoothW: B is not finite!"); 925 ap.assert(math.isfinite(xwidth), "AutoGKSmoothW: XWidth is not finite!"); 926 state.wrappermode = 0; 927 state.a = a; 928 state.b = b; 929 state.xwidth = xwidth; 930 state.needf = false; 931 state.rstate.ra = new double[10+1]; 932 state.rstate.stage = -1; 731 if( state.needf ) 732 { 733 func(state.innerobj.x, state.innerobj.xminusa, state.innerobj.bminusx, ref state.innerobj.f, obj); 734 continue; 735 } 736 throw new alglibexception("ALGLIB: unexpected error in 'autogksolve'"); 933 737 } 934 935 936 /*************************************************************************937 Integration on a finite interval [A,B].938 Integrand have integrable singularities at A/B.939 940 F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known941 alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates942 from below can be used (but these estimates should be greater than -1 too).943 944 One of alpha/beta variables (or even both alpha/beta) may be equal to 0,945 which means than function F(x) is non-singular at A/B. Anyway (singular at946 bounds or not), function F(x) is supposed to be continuous on (A,B).947 948 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result949 is calculated with accuracy close to the machine precision.950 951 INPUT PARAMETERS:952 A, B - interval boundaries (A<B, A=B or A>B)953 Alpha - power-law coefficient of the F(x) at A,954 Alpha>-1955 Beta - power-law coefficient of the F(x) at B,956 Beta>-1957 958 OUTPUT PARAMETERS959 State - structure which stores algorithm state960 961 SEE ALSO962 AutoGKSmooth, AutoGKSmoothW, AutoGKResults.963 964 965 -- ALGLIB --966 Copyright 06.05.2009 by Bochkanov Sergey967 *************************************************************************/968 public static void autogksingular(double a,969 double b,970 double alpha,971 double beta,972 autogkstate state)973 {974 ap.assert(math.isfinite(a), "AutoGKSingular: A is not finite!");975 ap.assert(math.isfinite(b), "AutoGKSingular: B is not finite!");976 ap.assert(math.isfinite(alpha), "AutoGKSingular: Alpha is not finite!");977 ap.assert(math.isfinite(beta), "AutoGKSingular: Beta is not finite!");978 state.wrappermode = 1;979 state.a = a;980 state.b = b;981 state.alpha = alpha;982 state.beta = beta;983 state.xwidth = 0.0;984 state.needf = false;985 state.rstate.ra = new double[10+1];986 state.rstate.stage = -1;987 }988 989 990 /*************************************************************************991 992 -- ALGLIB --993 Copyright 07.05.2009 by Bochkanov Sergey994 *************************************************************************/995 public static bool autogkiteration(autogkstate state)996 {997 bool result = new bool();998 double s = 0;999 double tmp = 0;1000 double eps = 0;1001 double a = 0;1002 double b = 0;1003 double x = 0;1004 double t = 0;1005 double alpha = 0;1006 double beta = 0;1007 double v1 = 0;1008 double v2 = 0;1009 1010 1011 //1012 // Reverse communication preparations1013 // I know it looks ugly, but it works the same way1014 // anywhere from C++ to Python.1015 //1016 // This code initializes locals by:1017 // * random values determined during code1018 // generation - on first subroutine call1019 // * values from previous call - on subsequent calls1020 //1021 if( state.rstate.stage>=0 )1022 {1023 s = state.rstate.ra[0];1024 tmp = state.rstate.ra[1];1025 eps = state.rstate.ra[2];1026 a = state.rstate.ra[3];1027 b = state.rstate.ra[4];1028 x = state.rstate.ra[5];1029 t = state.rstate.ra[6];1030 alpha = state.rstate.ra[7];1031 beta = state.rstate.ra[8];1032 v1 = state.rstate.ra[9];1033 v2 = state.rstate.ra[10];1034 }1035 else1036 {1037 s = -983;1038 tmp = -989;1039 eps = -834;1040 a = 900;1041 b = -287;1042 x = 364;1043 t = 214;1044 alpha = -338;1045 beta = -686;1046 v1 = 912;1047 v2 = 585;1048 }1049 if( state.rstate.stage==0 )1050 {1051 goto lbl_0;1052 }1053 if( state.rstate.stage==1 )1054 {1055 goto lbl_1;1056 }1057 if( state.rstate.stage==2 )1058 {1059 goto lbl_2;1060 }1061 1062 //1063 // Routine body1064 //1065 eps = 0;1066 a = state.a;1067 b = state.b;1068 alpha = state.alpha;1069 beta = state.beta;1070 state.terminationtype = -1;1071 state.nfev = 0;1072 state.nintervals = 0;1073 1074 //1075 // smooth function at a finite interval1076 //1077 if( state.wrappermode!=0 )1078 {1079 goto lbl_3;1080 }1081 1082 //1083 // special case1084 //1085 if( (double)(a)==(double)(b) )1086 {1087 state.terminationtype = 1;1088 state.v = 0;1089 result = false;1090 return result;1091 }1092 1093 //1094 // general case1095 //1096 autogkinternalprepare(a, b, eps, state.xwidth, state.internalstate);1097 lbl_5:1098 if( !autogkinternaliteration(state.internalstate) )1099 {1100 goto lbl_6;1101 }1102 x = state.internalstate.x;1103 state.x = x;1104 state.xminusa = x-a;1105 state.bminusx = b-x;1106 state.needf = true;1107 state.rstate.stage = 0;1108 goto lbl_rcomm;1109 lbl_0:1110 state.needf = false;1111 state.nfev = state.nfev+1;1112 state.internalstate.f = state.f;1113 goto lbl_5;1114 lbl_6:1115 state.v = state.internalstate.r;1116 state.terminationtype = state.internalstate.info;1117 state.nintervals = state.internalstate.heapused;1118 result = false;1119 return result;1120 lbl_3:1121 1122 //1123 // function with power-law singularities at the ends of a finite interval1124 //1125 if( state.wrappermode!=1 )1126 {1127 goto lbl_7;1128 }1129 1130 //1131 // test coefficients1132 //1133 if( (double)(alpha)<=(double)(-1) | (double)(beta)<=(double)(-1) )1134 {1135 state.terminationtype = -1;1136 state.v = 0;1137 result = false;1138 return result;1139 }1140 1141 //1142 // special cases1143 //1144 if( (double)(a)==(double)(b) )1145 {1146 state.terminationtype = 1;1147 state.v = 0;1148 result = false;1149 return result;1150 }1151 1152 //1153 // reduction to general form1154 //1155 if( (double)(a)<(double)(b) )1156 {1157 s = 1;1158 }1159 else1160 {1161 s = -1;1162 tmp = a;1163 a = b;1164 b = tmp;1165 tmp = alpha;1166 alpha = beta;1167 beta = tmp;1168 }1169 alpha = Math.Min(alpha, 0);1170 beta = Math.Min(beta, 0);1171 1172 //1173 // first, integrate left half of [a,b]:1174 // integral(f(x)dx, a, (b+a)/2) =1175 // = 1/(1+alpha) * integral(t^(-alpha/(1+alpha))*f(a+t^(1/(1+alpha)))dt, 0, (0.5*(b-a))^(1+alpha))1176 //1177 autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+alpha), eps, state.xwidth, state.internalstate);1178 lbl_9:1179 if( !autogkinternaliteration(state.internalstate) )1180 {1181 goto lbl_10;1182 }1183 1184 //1185 // Fill State.X, State.XMinusA, State.BMinusX.1186 // Latter two are filled correctly even if B<A.1187 //1188 x = state.internalstate.x;1189 t = Math.Pow(x, 1/(1+alpha));1190 state.x = a+t;1191 if( (double)(s)>(double)(0) )1192 {1193 state.xminusa = t;1194 state.bminusx = b-(a+t);1195 }1196 else1197 {1198 state.xminusa = a+t-b;1199 state.bminusx = -t;1200 }1201 state.needf = true;1202 state.rstate.stage = 1;1203 goto lbl_rcomm;1204 lbl_1:1205 state.needf = false;1206 if( (double)(alpha)!=(double)(0) )1207 {1208 state.internalstate.f = state.f*Math.Pow(x, -(alpha/(1+alpha)))/(1+alpha);1209 }1210 else1211 {1212 state.internalstate.f = state.f;1213 }1214 state.nfev = state.nfev+1;1215 goto lbl_9;1216 lbl_10:1217 v1 = state.internalstate.r;1218 state.nintervals = state.nintervals+state.internalstate.heapused;1219 1220 //1221 // then, integrate right half of [a,b]:1222 // integral(f(x)dx, (b+a)/2, b) =1223 // = 1/(1+beta) * integral(t^(-beta/(1+beta))*f(b-t^(1/(1+beta)))dt, 0, (0.5*(b-a))^(1+beta))1224 //1225 autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+beta), eps, state.xwidth, state.internalstate);1226 lbl_11:1227 if( !autogkinternaliteration(state.internalstate) )1228 {1229 goto lbl_12;1230 }1231 1232 //1233 // Fill State.X, State.XMinusA, State.BMinusX.1234 // Latter two are filled correctly (X-A, B-X) even if B<A.1235 //1236 x = state.internalstate.x;1237 t = Math.Pow(x, 1/(1+beta));1238 state.x = b-t;1239 if( (double)(s)>(double)(0) )1240 {1241 state.xminusa = b-t-a;1242 state.bminusx = t;1243 }1244 else1245 {1246 state.xminusa = -t;1247 state.bminusx = a-(b-t);1248 }1249 state.needf = true;1250 state.rstate.stage = 2;1251 goto lbl_rcomm;1252 lbl_2:1253 state.needf = false;1254 if( (double)(beta)!=(double)(0) )1255 {1256 state.internalstate.f = state.f*Math.Pow(x, -(beta/(1+beta)))/(1+beta);1257 }1258 else1259 {1260 state.internalstate.f = state.f;1261 }1262 state.nfev = state.nfev+1;1263 goto lbl_11;1264 lbl_12:1265 v2 = state.internalstate.r;1266 state.nintervals = state.nintervals+state.internalstate.heapused;1267 1268 //1269 // final result1270 //1271 state.v = s*(v1+v2);1272 state.terminationtype = 1;1273 result = false;1274 return result;1275 lbl_7:1276 result = false;1277 return result;1278 1279 //1280 // Saving state1281 //1282 lbl_rcomm:1283 result = true;1284 state.rstate.ra[0] = s;1285 state.rstate.ra[1] = tmp;1286 state.rstate.ra[2] = eps;1287 state.rstate.ra[3] = a;1288 state.rstate.ra[4] = b;1289 state.rstate.ra[5] = x;1290 state.rstate.ra[6] = t;1291 state.rstate.ra[7] = alpha;1292 state.rstate.ra[8] = beta;1293 state.rstate.ra[9] = v1;1294 state.rstate.ra[10] = v2;1295 return result;1296 }1297 1298 1299 /*************************************************************************1300 Adaptive integration results1301 1302 Called after AutoGKIteration returned False.1303 1304 Input parameters:1305 State - algorithm state (used by AutoGKIteration).1306 1307 Output parameters:1308 V - integral(f(x)dx,a,b)1309 Rep - optimization report (see AutoGKReport description)1310 1311 -- ALGLIB --1312 Copyright 14.11.2007 by Bochkanov Sergey1313 *************************************************************************/1314 public static void autogkresults(autogkstate state,1315 ref double v,1316 autogkreport rep)1317 {1318 v = 0;1319 1320 v = state.v;1321 rep.terminationtype = state.terminationtype;1322 rep.nfev = state.nfev;1323 rep.nintervals = state.nintervals;1324 }1325 1326 1327 /*************************************************************************1328 Internal AutoGK subroutine1329 eps<0 - error1330 eps=0 - automatic eps selection1331 1332 width<0 - error1333 width=0 - no width requirements1334 *************************************************************************/1335 private static void autogkinternalprepare(double a,1336 double b,1337 double eps,1338 double xwidth,1339 autogkinternalstate state)1340 {1341 1342 //1343 // Save settings1344 //1345 state.a = a;1346 state.b = b;1347 state.eps = eps;1348 state.xwidth = xwidth;1349 1350 //1351 // Prepare RComm structure1352 //1353 state.rstate.ia = new int[3+1];1354 state.rstate.ra = new double[8+1];1355 state.rstate.stage = -1;1356 }1357 1358 1359 /*************************************************************************1360 Internal AutoGK subroutine1361 *************************************************************************/1362 private static bool autogkinternaliteration(autogkinternalstate state)1363 {1364 bool result = new bool();1365 double c1 = 0;1366 double c2 = 0;1367 int i = 0;1368 int j = 0;1369 double intg = 0;1370 double intk = 0;1371 double inta = 0;1372 double v = 0;1373 double ta = 0;1374 double tb = 0;1375 int ns = 0;1376 double qeps = 0;1377 int info = 0;1378 1379 1380 //1381 // Reverse communication preparations1382 // I know it looks ugly, but it works the same way1383 // anywhere from C++ to Python.1384 //1385 // This code initializes locals by:1386 // * random values determined during code1387 // generation - on first subroutine call1388 // * values from previous call - on subsequent calls1389 //1390 if( state.rstate.stage>=0 )1391 {1392 i = state.rstate.ia[0];1393 j = state.rstate.ia[1];1394 ns = state.rstate.ia[2];1395 info = state.rstate.ia[3];1396 c1 = state.rstate.ra[0];1397 c2 = state.rstate.ra[1];1398 intg = state.rstate.ra[2];1399 intk = state.rstate.ra[3];1400 inta = state.rstate.ra[4];1401 v = state.rstate.ra[5];1402 ta = state.rstate.ra[6];1403 tb = state.rstate.ra[7];1404 qeps = state.rstate.ra[8];1405 }1406 else1407 {1408 i = 497;1409 j = -271;1410 ns = -581;1411 info = 745;1412 c1 = -533;1413 c2 = -77;1414 intg = 678;1415 intk = -293;1416 inta = 316;1417 v = 647;1418 ta = -756;1419 tb = 830;1420 qeps = -871;1421 }1422 if( state.rstate.stage==0 )1423 {1424 goto lbl_0;1425 }1426 if( state.rstate.stage==1 )1427 {1428 goto lbl_1;1429 }1430 if( state.rstate.stage==2 )1431 {1432 goto lbl_2;1433 }1434 1435 //1436 // Routine body1437 //1438 1439 //1440 // initialize quadratures.1441 // use 15-point Gauss-Kronrod formula.1442 //1443 state.n = 15;1444 gkq.gkqgenerategausslegendre(state.n, ref info, ref state.qn, ref state.wk, ref state.wg);1445 if( info<0 )1446 {1447 state.info = -5;1448 state.r = 0;1449 result = false;1450 return result;1451 }1452 state.wr = new double[state.n];1453 for(i=0; i<=state.n-1; i++)1454 {1455 if( i==0 )1456 {1457 state.wr[i] = 0.5*Math.Abs(state.qn[1]-state.qn[0]);1458 continue;1459 }1460 if( i==state.n-1 )1461 {1462 state.wr[state.n-1] = 0.5*Math.Abs(state.qn[state.n-1]-state.qn[state.n-2]);1463 continue;1464 }1465 state.wr[i] = 0.5*Math.Abs(state.qn[i-1]-state.qn[i+1]);1466 }1467 1468 //1469 // special case1470 //1471 if( (double)(state.a)==(double)(state.b) )1472 {1473 state.info = 1;1474 state.r = 0;1475 result = false;1476 return result;1477 }1478 1479 //1480 // test parameters1481 //1482 if( (double)(state.eps)<(double)(0) | (double)(state.xwidth)<(double)(0) )1483 {1484 state.info = -1;1485 state.r = 0;1486 result = false;1487 return result;1488 }1489 state.info = 1;1490 if( (double)(state.eps)==(double)(0) )1491 {1492 state.eps = 1000*math.machineepsilon;1493 }1494 1495 //1496 // First, prepare heap1497 // * column 0 - absolute error1498 // * column 1 - integral of a F(x) (calculated using Kronrod extension nodes)1499 // * column 2 - integral of a |F(x)| (calculated using modified rect. method)1500 // * column 3 - left boundary of a subinterval1501 // * column 4 - right boundary of a subinterval1502 //1503 if( (double)(state.xwidth)!=(double)(0) )1504 {1505 goto lbl_3;1506 }1507 1508 //1509 // no maximum width requirements1510 // start from one big subinterval1511 //1512 state.heapwidth = 5;1513 state.heapsize = 1;1514 state.heapused = 1;1515 state.heap = new double[state.heapsize, state.heapwidth];1516 c1 = 0.5*(state.b-state.a);1517 c2 = 0.5*(state.b+state.a);1518 intg = 0;1519 intk = 0;1520 inta = 0;1521 i = 0;1522 lbl_5:1523 if( i>state.n-1 )1524 {1525 goto lbl_7;1526 }1527 1528 //1529 // obtain F1530 //1531 state.x = c1*state.qn[i]+c2;1532 state.rstate.stage = 0;1533 goto lbl_rcomm;1534 lbl_0:1535 v = state.f;1536 1537 //1538 // Gauss-Kronrod formula1539 //1540 intk = intk+v*state.wk[i];1541 if( i%2==1 )1542 {1543 intg = intg+v*state.wg[i];1544 }1545 1546 //1547 // Integral |F(x)|1548 // Use rectangles method1549 //1550 inta = inta+Math.Abs(v)*state.wr[i];1551 i = i+1;1552 goto lbl_5;1553 lbl_7:1554 intk = intk*(state.b-state.a)*0.5;1555 intg = intg*(state.b-state.a)*0.5;1556 inta = inta*(state.b-state.a)*0.5;1557 state.heap[0,0] = Math.Abs(intg-intk);1558 state.heap[0,1] = intk;1559 state.heap[0,2] = inta;1560 state.heap[0,3] = state.a;1561 state.heap[0,4] = state.b;1562 state.sumerr = state.heap[0,0];1563 state.sumabs = Math.Abs(inta);1564 goto lbl_4;1565 lbl_3:1566 1567 //1568 // maximum subinterval should be no more than XWidth.1569 // so we create Ceil((B-A)/XWidth)+1 small subintervals1570 //1571 ns = (int)Math.Ceiling(Math.Abs(state.b-state.a)/state.xwidth)+1;1572 state.heapsize = ns;1573 state.heapused = ns;1574 state.heapwidth = 5;1575 state.heap = new double[state.heapsize, state.heapwidth];1576 state.sumerr = 0;1577 state.sumabs = 0;1578 j = 0;1579 lbl_8:1580 if( j>ns-1 )1581 {1582 goto lbl_10;1583 }1584 ta = state.a+j*(state.b-state.a)/ns;1585 tb = state.a+(j+1)*(state.b-state.a)/ns;1586 c1 = 0.5*(tb-ta);1587 c2 = 0.5*(tb+ta);1588 intg = 0;1589 intk = 0;1590 inta = 0;1591 i = 0;1592 lbl_11:1593 if( i>state.n-1 )1594 {1595 goto lbl_13;1596 }1597 1598 //1599 // obtain F1600 //1601 state.x = c1*state.qn[i]+c2;1602 state.rstate.stage = 1;1603 goto lbl_rcomm;1604 lbl_1:1605 v = state.f;1606 1607 //1608 // Gauss-Kronrod formula1609 //1610 intk = intk+v*state.wk[i];1611 if( i%2==1 )1612 {1613 intg = intg+v*state.wg[i];1614 }1615 1616 //1617 // Integral |F(x)|1618 // Use rectangles method1619 //1620 inta = inta+Math.Abs(v)*state.wr[i];1621 i = i+1;1622 goto lbl_11;1623 lbl_13:1624 intk = intk*(tb-ta)*0.5;1625 intg = intg*(tb-ta)*0.5;1626 inta = inta*(tb-ta)*0.5;1627 state.heap[j,0] = Math.Abs(intg-intk);1628 state.heap[j,1] = intk;1629 state.heap[j,2] = inta;1630 state.heap[j,3] = ta;1631 state.heap[j,4] = tb;1632 state.sumerr = state.sumerr+state.heap[j,0];1633 state.sumabs = state.sumabs+Math.Abs(inta);1634 j = j+1;1635 goto lbl_8;1636 lbl_10:1637 lbl_4:1638 1639 //1640 // method iterations1641 //1642 lbl_14:1643 if( false )1644 {1645 goto lbl_15;1646 }1647 1648 //1649 // additional memory if needed1650 //1651 if( state.heapused==state.heapsize )1652 {1653 mheapresize(ref state.heap, ref state.heapsize, 4*state.heapsize, state.heapwidth);1654 }1655 1656 //1657 // TODO: every 20 iterations recalculate errors/sums1658 // TODO: one more criterion to prevent infinite loops with too strict Eps1659 //1660 if( (double)(state.sumerr)<=(double)(state.eps*state.sumabs) )1661 {1662 state.r = 0;1663 for(j=0; j<=state.heapused-1; j++)1664 {1665 state.r = state.r+state.heap[j,1];1666 }1667 result = false;1668 return result;1669 }1670 1671 //1672 // Exclude interval with maximum absolute error1673 //1674 mheappop(ref state.heap, state.heapused, state.heapwidth);1675 state.sumerr = state.sumerr-state.heap[state.heapused-1,0];1676 state.sumabs = state.sumabs-state.heap[state.heapused-1,2];1677 1678 //1679 // Divide interval, create subintervals1680 //1681 ta = state.heap[state.heapused-1,3];1682 tb = state.heap[state.heapused-1,4];1683 state.heap[state.heapused-1,3] = ta;1684 state.heap[state.heapused-1,4] = 0.5*(ta+tb);1685 state.heap[state.heapused,3] = 0.5*(ta+tb);1686 state.heap[state.heapused,4] = tb;1687 j = state.heapused-1;1688 lbl_16:1689 if( j>state.heapused )1690 {1691 goto lbl_18;1692 }1693 c1 = 0.5*(state.heap[j,4]-state.heap[j,3]);1694 c2 = 0.5*(state.heap[j,4]+state.heap[j,3]);1695 intg = 0;1696 intk = 0;1697 inta = 0;1698 i = 0;1699 lbl_19:1700 if( i>state.n-1 )1701 {1702 goto lbl_21;1703 }1704 1705 //1706 // F(x)1707 //1708 state.x = c1*state.qn[i]+c2;1709 state.rstate.stage = 2;1710 goto lbl_rcomm;1711 lbl_2:1712 v = state.f;1713 1714 //1715 // Gauss-Kronrod formula1716 //1717 intk = intk+v*state.wk[i];1718 if( i%2==1 )1719 {1720 intg = intg+v*state.wg[i];1721 }1722 1723 //1724 // Integral |F(x)|1725 // Use rectangles method1726 //1727 inta = inta+Math.Abs(v)*state.wr[i];1728 i = i+1;1729 goto lbl_19;1730 lbl_21:1731 intk = intk*(state.heap[j,4]-state.heap[j,3])*0.5;1732 intg = intg*(state.heap[j,4]-state.heap[j,3])*0.5;1733 inta = inta*(state.heap[j,4]-state.heap[j,3])*0.5;1734 state.heap[j,0] = Math.Abs(intg-intk);1735 state.heap[j,1] = intk;1736 state.heap[j,2] = inta;1737 state.sumerr = state.sumerr+state.heap[j,0];1738 state.sumabs = state.sumabs+state.heap[j,2];1739 j = j+1;1740 goto lbl_16;1741 lbl_18:1742 mheappush(ref state.heap, state.heapused-1, state.heapwidth);1743 mheappush(ref state.heap, state.heapused, state.heapwidth);1744 state.heapused = state.heapused+1;1745 goto lbl_14;1746 lbl_15:1747 result = false;1748 return result;1749 1750 //1751 // Saving state1752 //1753 lbl_rcomm:1754 result = true;1755 state.rstate.ia[0] = i;1756 state.rstate.ia[1] = j;1757 state.rstate.ia[2] = ns;1758 state.rstate.ia[3] = info;1759 state.rstate.ra[0] = c1;1760 state.rstate.ra[1] = c2;1761 state.rstate.ra[2] = intg;1762 state.rstate.ra[3] = intk;1763 state.rstate.ra[4] = inta;1764 state.rstate.ra[5] = v;1765 state.rstate.ra[6] = ta;1766 state.rstate.ra[7] = tb;1767 state.rstate.ra[8] = qeps;1768 return result;1769 }1770 1771 1772 private static void mheappop(ref double[,] heap,1773 int heapsize,1774 int heapwidth)1775 {1776 int i = 0;1777 int p = 0;1778 double t = 0;1779 int maxcp = 0;1780 1781 if( heapsize==1 )1782 {1783 return;1784 }1785 for(i=0; i<=heapwidth-1; i++)1786 {1787 t = heap[heapsize-1,i];1788 heap[heapsize-1,i] = heap[0,i];1789 heap[0,i] = t;1790 }1791 p = 0;1792 while( 2*p+1<heapsize-1 )1793 {1794 maxcp = 2*p+1;1795 if( 2*p+2<heapsize-1 )1796 {1797 if( (double)(heap[2*p+2,0])>(double)(heap[2*p+1,0]) )1798 {1799 maxcp = 2*p+2;1800 }1801 }1802 if( (double)(heap[p,0])<(double)(heap[maxcp,0]) )1803 {1804 for(i=0; i<=heapwidth-1; i++)1805 {1806 t = heap[p,i];1807 heap[p,i] = heap[maxcp,i];1808 heap[maxcp,i] = t;1809 }1810 p = maxcp;1811 }1812 else1813 {1814 break;1815 }1816 }1817 }1818 1819 1820 private static void mheappush(ref double[,] heap,1821 int heapsize,1822 int heapwidth)1823 {1824 int i = 0;1825 int p = 0;1826 double t = 0;1827 int parent = 0;1828 1829 if( heapsize==0 )1830 {1831 return;1832 }1833 p = heapsize;1834 while( p!=0 )1835 {1836 parent = (p-1)/2;1837 if( (double)(heap[p,0])>(double)(heap[parent,0]) )1838 {1839 for(i=0; i<=heapwidth-1; i++)1840 {1841 t = heap[p,i];1842 heap[p,i] = heap[parent,i];1843 heap[parent,i] = t;1844 }1845 p = parent;1846 }1847 else1848 {1849 break;1850 }1851 }1852 }1853 1854 1855 private static void mheapresize(ref double[,] heap,1856 ref int heapsize,1857 int newheapsize,1858 int heapwidth)1859 {1860 double[,] tmp = new double[0,0];1861 int i = 0;1862 int i_ = 0;1863 1864 tmp = new double[heapsize, heapwidth];1865 for(i=0; i<=heapsize-1; i++)1866 {1867 for(i_=0; i_<=heapwidth-1;i_++)1868 {1869 tmp[i,i_] = heap[i,i_];1870 }1871 }1872 heap = new double[newheapsize, heapwidth];1873 for(i=0; i<=heapsize-1; i++)1874 {1875 for(i_=0; i_<=heapwidth-1;i_++)1876 {1877 heap[i,i_] = tmp[i,i_];1878 }1879 }1880 heapsize = newheapsize;1881 }1882 1883 1884 738 } 739 740 /************************************************************************* 741 Adaptive integration results 742 743 Called after AutoGKIteration returned False. 744 745 Input parameters: 746 State - algorithm state (used by AutoGKIteration). 747 748 Output parameters: 749 V - integral(f(x)dx,a,b) 750 Rep - optimization report (see AutoGKReport description) 751 752 -- ALGLIB -- 753 Copyright 14.11.2007 by Bochkanov Sergey 754 *************************************************************************/ 755 public static void autogkresults(autogkstate state, out double v, out autogkreport rep) 756 { 757 v = 0; 758 rep = new autogkreport(); 759 autogk.autogkresults(state.innerobj, ref v, rep.innerobj); 760 return; 761 } 762 763 } 764 public partial class alglib 765 { 1885 766 public class gq 1886 767 { … … 3546 2427 3547 2428 } 2429 public class autogk 2430 { 2431 /************************************************************************* 2432 Integration report: 2433 * TerminationType = completetion code: 2434 * -5 non-convergence of Gauss-Kronrod nodes 2435 calculation subroutine. 2436 * -1 incorrect parameters were specified 2437 * 1 OK 2438 * Rep.NFEV countains number of function calculations 2439 * Rep.NIntervals contains number of intervals [a,b] 2440 was partitioned into. 2441 *************************************************************************/ 2442 public class autogkreport 2443 { 2444 public int terminationtype; 2445 public int nfev; 2446 public int nintervals; 2447 }; 2448 2449 2450 public class autogkinternalstate 2451 { 2452 public double a; 2453 public double b; 2454 public double eps; 2455 public double xwidth; 2456 public double x; 2457 public double f; 2458 public int info; 2459 public double r; 2460 public double[,] heap; 2461 public int heapsize; 2462 public int heapwidth; 2463 public int heapused; 2464 public double sumerr; 2465 public double sumabs; 2466 public double[] qn; 2467 public double[] wg; 2468 public double[] wk; 2469 public double[] wr; 2470 public int n; 2471 public rcommstate rstate; 2472 public autogkinternalstate() 2473 { 2474 heap = new double[0,0]; 2475 qn = new double[0]; 2476 wg = new double[0]; 2477 wk = new double[0]; 2478 wr = new double[0]; 2479 rstate = new rcommstate(); 2480 } 2481 }; 2482 2483 2484 /************************************************************************* 2485 This structure stores state of the integration algorithm. 2486 2487 Although this class has public fields, they are not intended for external 2488 use. You should use ALGLIB functions to work with this class: 2489 * autogksmooth()/AutoGKSmoothW()/... to create objects 2490 * autogkintegrate() to begin integration 2491 * autogkresults() to get results 2492 *************************************************************************/ 2493 public class autogkstate 2494 { 2495 public double a; 2496 public double b; 2497 public double alpha; 2498 public double beta; 2499 public double xwidth; 2500 public double x; 2501 public double xminusa; 2502 public double bminusx; 2503 public bool needf; 2504 public double f; 2505 public int wrappermode; 2506 public autogkinternalstate internalstate; 2507 public rcommstate rstate; 2508 public double v; 2509 public int terminationtype; 2510 public int nfev; 2511 public int nintervals; 2512 public autogkstate() 2513 { 2514 internalstate = new autogkinternalstate(); 2515 rstate = new rcommstate(); 2516 } 2517 }; 2518 2519 2520 2521 2522 public const int maxsubintervals = 10000; 2523 2524 2525 /************************************************************************* 2526 Integration of a smooth function F(x) on a finite interval [a,b]. 2527 2528 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result 2529 is calculated with accuracy close to the machine precision. 2530 2531 Algorithm works well only with smooth integrands. It may be used with 2532 continuous non-smooth integrands, but with less performance. 2533 2534 It should never be used with integrands which have integrable singularities 2535 at lower or upper limits - algorithm may crash. Use AutoGKSingular in such 2536 cases. 2537 2538 INPUT PARAMETERS: 2539 A, B - interval boundaries (A<B, A=B or A>B) 2540 2541 OUTPUT PARAMETERS 2542 State - structure which stores algorithm state 2543 2544 SEE ALSO 2545 AutoGKSmoothW, AutoGKSingular, AutoGKResults. 2546 2547 2548 -- ALGLIB -- 2549 Copyright 06.05.2009 by Bochkanov Sergey 2550 *************************************************************************/ 2551 public static void autogksmooth(double a, 2552 double b, 2553 autogkstate state) 2554 { 2555 ap.assert(math.isfinite(a), "AutoGKSmooth: A is not finite!"); 2556 ap.assert(math.isfinite(b), "AutoGKSmooth: B is not finite!"); 2557 autogksmoothw(a, b, 0.0, state); 2558 } 2559 2560 2561 /************************************************************************* 2562 Integration of a smooth function F(x) on a finite interval [a,b]. 2563 2564 This subroutine is same as AutoGKSmooth(), but it guarantees that interval 2565 [a,b] is partitioned into subintervals which have width at most XWidth. 2566 2567 Subroutine can be used when integrating nearly-constant function with 2568 narrow "bumps" (about XWidth wide). If "bumps" are too narrow, AutoGKSmooth 2569 subroutine can overlook them. 2570 2571 INPUT PARAMETERS: 2572 A, B - interval boundaries (A<B, A=B or A>B) 2573 2574 OUTPUT PARAMETERS 2575 State - structure which stores algorithm state 2576 2577 SEE ALSO 2578 AutoGKSmooth, AutoGKSingular, AutoGKResults. 2579 2580 2581 -- ALGLIB -- 2582 Copyright 06.05.2009 by Bochkanov Sergey 2583 *************************************************************************/ 2584 public static void autogksmoothw(double a, 2585 double b, 2586 double xwidth, 2587 autogkstate state) 2588 { 2589 ap.assert(math.isfinite(a), "AutoGKSmoothW: A is not finite!"); 2590 ap.assert(math.isfinite(b), "AutoGKSmoothW: B is not finite!"); 2591 ap.assert(math.isfinite(xwidth), "AutoGKSmoothW: XWidth is not finite!"); 2592 state.wrappermode = 0; 2593 state.a = a; 2594 state.b = b; 2595 state.xwidth = xwidth; 2596 state.needf = false; 2597 state.rstate.ra = new double[10+1]; 2598 state.rstate.stage = -1; 2599 } 2600 2601 2602 /************************************************************************* 2603 Integration on a finite interval [A,B]. 2604 Integrand have integrable singularities at A/B. 2605 2606 F(X) must diverge as "(x-A)^alpha" at A, as "(B-x)^beta" at B, with known 2607 alpha/beta (alpha>-1, beta>-1). If alpha/beta are not known, estimates 2608 from below can be used (but these estimates should be greater than -1 too). 2609 2610 One of alpha/beta variables (or even both alpha/beta) may be equal to 0, 2611 which means than function F(x) is non-singular at A/B. Anyway (singular at 2612 bounds or not), function F(x) is supposed to be continuous on (A,B). 2613 2614 Fast-convergent algorithm based on a Gauss-Kronrod formula is used. Result 2615 is calculated with accuracy close to the machine precision. 2616 2617 INPUT PARAMETERS: 2618 A, B - interval boundaries (A<B, A=B or A>B) 2619 Alpha - power-law coefficient of the F(x) at A, 2620 Alpha>-1 2621 Beta - power-law coefficient of the F(x) at B, 2622 Beta>-1 2623 2624 OUTPUT PARAMETERS 2625 State - structure which stores algorithm state 2626 2627 SEE ALSO 2628 AutoGKSmooth, AutoGKSmoothW, AutoGKResults. 2629 2630 2631 -- ALGLIB -- 2632 Copyright 06.05.2009 by Bochkanov Sergey 2633 *************************************************************************/ 2634 public static void autogksingular(double a, 2635 double b, 2636 double alpha, 2637 double beta, 2638 autogkstate state) 2639 { 2640 ap.assert(math.isfinite(a), "AutoGKSingular: A is not finite!"); 2641 ap.assert(math.isfinite(b), "AutoGKSingular: B is not finite!"); 2642 ap.assert(math.isfinite(alpha), "AutoGKSingular: Alpha is not finite!"); 2643 ap.assert(math.isfinite(beta), "AutoGKSingular: Beta is not finite!"); 2644 state.wrappermode = 1; 2645 state.a = a; 2646 state.b = b; 2647 state.alpha = alpha; 2648 state.beta = beta; 2649 state.xwidth = 0.0; 2650 state.needf = false; 2651 state.rstate.ra = new double[10+1]; 2652 state.rstate.stage = -1; 2653 } 2654 2655 2656 /************************************************************************* 2657 2658 -- ALGLIB -- 2659 Copyright 07.05.2009 by Bochkanov Sergey 2660 *************************************************************************/ 2661 public static bool autogkiteration(autogkstate state) 2662 { 2663 bool result = new bool(); 2664 double s = 0; 2665 double tmp = 0; 2666 double eps = 0; 2667 double a = 0; 2668 double b = 0; 2669 double x = 0; 2670 double t = 0; 2671 double alpha = 0; 2672 double beta = 0; 2673 double v1 = 0; 2674 double v2 = 0; 2675 2676 2677 // 2678 // Reverse communication preparations 2679 // I know it looks ugly, but it works the same way 2680 // anywhere from C++ to Python. 2681 // 2682 // This code initializes locals by: 2683 // * random values determined during code 2684 // generation - on first subroutine call 2685 // * values from previous call - on subsequent calls 2686 // 2687 if( state.rstate.stage>=0 ) 2688 { 2689 s = state.rstate.ra[0]; 2690 tmp = state.rstate.ra[1]; 2691 eps = state.rstate.ra[2]; 2692 a = state.rstate.ra[3]; 2693 b = state.rstate.ra[4]; 2694 x = state.rstate.ra[5]; 2695 t = state.rstate.ra[6]; 2696 alpha = state.rstate.ra[7]; 2697 beta = state.rstate.ra[8]; 2698 v1 = state.rstate.ra[9]; 2699 v2 = state.rstate.ra[10]; 2700 } 2701 else 2702 { 2703 s = -983; 2704 tmp = -989; 2705 eps = -834; 2706 a = 900; 2707 b = -287; 2708 x = 364; 2709 t = 214; 2710 alpha = -338; 2711 beta = -686; 2712 v1 = 912; 2713 v2 = 585; 2714 } 2715 if( state.rstate.stage==0 ) 2716 { 2717 goto lbl_0; 2718 } 2719 if( state.rstate.stage==1 ) 2720 { 2721 goto lbl_1; 2722 } 2723 if( state.rstate.stage==2 ) 2724 { 2725 goto lbl_2; 2726 } 2727 2728 // 2729 // Routine body 2730 // 2731 eps = 0; 2732 a = state.a; 2733 b = state.b; 2734 alpha = state.alpha; 2735 beta = state.beta; 2736 state.terminationtype = -1; 2737 state.nfev = 0; 2738 state.nintervals = 0; 2739 2740 // 2741 // smooth function at a finite interval 2742 // 2743 if( state.wrappermode!=0 ) 2744 { 2745 goto lbl_3; 2746 } 2747 2748 // 2749 // special case 2750 // 2751 if( (double)(a)==(double)(b) ) 2752 { 2753 state.terminationtype = 1; 2754 state.v = 0; 2755 result = false; 2756 return result; 2757 } 2758 2759 // 2760 // general case 2761 // 2762 autogkinternalprepare(a, b, eps, state.xwidth, state.internalstate); 2763 lbl_5: 2764 if( !autogkinternaliteration(state.internalstate) ) 2765 { 2766 goto lbl_6; 2767 } 2768 x = state.internalstate.x; 2769 state.x = x; 2770 state.xminusa = x-a; 2771 state.bminusx = b-x; 2772 state.needf = true; 2773 state.rstate.stage = 0; 2774 goto lbl_rcomm; 2775 lbl_0: 2776 state.needf = false; 2777 state.nfev = state.nfev+1; 2778 state.internalstate.f = state.f; 2779 goto lbl_5; 2780 lbl_6: 2781 state.v = state.internalstate.r; 2782 state.terminationtype = state.internalstate.info; 2783 state.nintervals = state.internalstate.heapused; 2784 result = false; 2785 return result; 2786 lbl_3: 2787 2788 // 2789 // function with power-law singularities at the ends of a finite interval 2790 // 2791 if( state.wrappermode!=1 ) 2792 { 2793 goto lbl_7; 2794 } 2795 2796 // 2797 // test coefficients 2798 // 2799 if( (double)(alpha)<=(double)(-1) | (double)(beta)<=(double)(-1) ) 2800 { 2801 state.terminationtype = -1; 2802 state.v = 0; 2803 result = false; 2804 return result; 2805 } 2806 2807 // 2808 // special cases 2809 // 2810 if( (double)(a)==(double)(b) ) 2811 { 2812 state.terminationtype = 1; 2813 state.v = 0; 2814 result = false; 2815 return result; 2816 } 2817 2818 // 2819 // reduction to general form 2820 // 2821 if( (double)(a)<(double)(b) ) 2822 { 2823 s = 1; 2824 } 2825 else 2826 { 2827 s = -1; 2828 tmp = a; 2829 a = b; 2830 b = tmp; 2831 tmp = alpha; 2832 alpha = beta; 2833 beta = tmp; 2834 } 2835 alpha = Math.Min(alpha, 0); 2836 beta = Math.Min(beta, 0); 2837 2838 // 2839 // first, integrate left half of [a,b]: 2840 // integral(f(x)dx, a, (b+a)/2) = 2841 // = 1/(1+alpha) * integral(t^(-alpha/(1+alpha))*f(a+t^(1/(1+alpha)))dt, 0, (0.5*(b-a))^(1+alpha)) 2842 // 2843 autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+alpha), eps, state.xwidth, state.internalstate); 2844 lbl_9: 2845 if( !autogkinternaliteration(state.internalstate) ) 2846 { 2847 goto lbl_10; 2848 } 2849 2850 // 2851 // Fill State.X, State.XMinusA, State.BMinusX. 2852 // Latter two are filled correctly even if B<A. 2853 // 2854 x = state.internalstate.x; 2855 t = Math.Pow(x, 1/(1+alpha)); 2856 state.x = a+t; 2857 if( (double)(s)>(double)(0) ) 2858 { 2859 state.xminusa = t; 2860 state.bminusx = b-(a+t); 2861 } 2862 else 2863 { 2864 state.xminusa = a+t-b; 2865 state.bminusx = -t; 2866 } 2867 state.needf = true; 2868 state.rstate.stage = 1; 2869 goto lbl_rcomm; 2870 lbl_1: 2871 state.needf = false; 2872 if( (double)(alpha)!=(double)(0) ) 2873 { 2874 state.internalstate.f = state.f*Math.Pow(x, -(alpha/(1+alpha)))/(1+alpha); 2875 } 2876 else 2877 { 2878 state.internalstate.f = state.f; 2879 } 2880 state.nfev = state.nfev+1; 2881 goto lbl_9; 2882 lbl_10: 2883 v1 = state.internalstate.r; 2884 state.nintervals = state.nintervals+state.internalstate.heapused; 2885 2886 // 2887 // then, integrate right half of [a,b]: 2888 // integral(f(x)dx, (b+a)/2, b) = 2889 // = 1/(1+beta) * integral(t^(-beta/(1+beta))*f(b-t^(1/(1+beta)))dt, 0, (0.5*(b-a))^(1+beta)) 2890 // 2891 autogkinternalprepare(0, Math.Pow(0.5*(b-a), 1+beta), eps, state.xwidth, state.internalstate); 2892 lbl_11: 2893 if( !autogkinternaliteration(state.internalstate) ) 2894 { 2895 goto lbl_12; 2896 } 2897 2898 // 2899 // Fill State.X, State.XMinusA, State.BMinusX. 2900 // Latter two are filled correctly (X-A, B-X) even if B<A. 2901 // 2902 x = state.internalstate.x; 2903 t = Math.Pow(x, 1/(1+beta)); 2904 state.x = b-t; 2905 if( (double)(s)>(double)(0) ) 2906 { 2907 state.xminusa = b-t-a; 2908 state.bminusx = t; 2909 } 2910 else 2911 { 2912 state.xminusa = -t; 2913 state.bminusx = a-(b-t); 2914 } 2915 state.needf = true; 2916 state.rstate.stage = 2; 2917 goto lbl_rcomm; 2918 lbl_2: 2919 state.needf = false; 2920 if( (double)(beta)!=(double)(0) ) 2921 { 2922 state.internalstate.f = state.f*Math.Pow(x, -(beta/(1+beta)))/(1+beta); 2923 } 2924 else 2925 { 2926 state.internalstate.f = state.f; 2927 } 2928 state.nfev = state.nfev+1; 2929 goto lbl_11; 2930 lbl_12: 2931 v2 = state.internalstate.r; 2932 state.nintervals = state.nintervals+state.internalstate.heapused; 2933 2934 // 2935 // final result 2936 // 2937 state.v = s*(v1+v2); 2938 state.terminationtype = 1; 2939 result = false; 2940 return result; 2941 lbl_7: 2942 result = false; 2943 return result; 2944 2945 // 2946 // Saving state 2947 // 2948 lbl_rcomm: 2949 result = true; 2950 state.rstate.ra[0] = s; 2951 state.rstate.ra[1] = tmp; 2952 state.rstate.ra[2] = eps; 2953 state.rstate.ra[3] = a; 2954 state.rstate.ra[4] = b; 2955 state.rstate.ra[5] = x; 2956 state.rstate.ra[6] = t; 2957 state.rstate.ra[7] = alpha; 2958 state.rstate.ra[8] = beta; 2959 state.rstate.ra[9] = v1; 2960 state.rstate.ra[10] = v2; 2961 return result; 2962 } 2963 2964 2965 /************************************************************************* 2966 Adaptive integration results 2967 2968 Called after AutoGKIteration returned False. 2969 2970 Input parameters: 2971 State - algorithm state (used by AutoGKIteration). 2972 2973 Output parameters: 2974 V - integral(f(x)dx,a,b) 2975 Rep - optimization report (see AutoGKReport description) 2976 2977 -- ALGLIB -- 2978 Copyright 14.11.2007 by Bochkanov Sergey 2979 *************************************************************************/ 2980 public static void autogkresults(autogkstate state, 2981 ref double v, 2982 autogkreport rep) 2983 { 2984 v = 0; 2985 2986 v = state.v; 2987 rep.terminationtype = state.terminationtype; 2988 rep.nfev = state.nfev; 2989 rep.nintervals = state.nintervals; 2990 } 2991 2992 2993 /************************************************************************* 2994 Internal AutoGK subroutine 2995 eps<0 - error 2996 eps=0 - automatic eps selection 2997 2998 width<0 - error 2999 width=0 - no width requirements 3000 *************************************************************************/ 3001 private static void autogkinternalprepare(double a, 3002 double b, 3003 double eps, 3004 double xwidth, 3005 autogkinternalstate state) 3006 { 3007 3008 // 3009 // Save settings 3010 // 3011 state.a = a; 3012 state.b = b; 3013 state.eps = eps; 3014 state.xwidth = xwidth; 3015 3016 // 3017 // Prepare RComm structure 3018 // 3019 state.rstate.ia = new int[3+1]; 3020 state.rstate.ra = new double[8+1]; 3021 state.rstate.stage = -1; 3022 } 3023 3024 3025 /************************************************************************* 3026 Internal AutoGK subroutine 3027 *************************************************************************/ 3028 private static bool autogkinternaliteration(autogkinternalstate state) 3029 { 3030 bool result = new bool(); 3031 double c1 = 0; 3032 double c2 = 0; 3033 int i = 0; 3034 int j = 0; 3035 double intg = 0; 3036 double intk = 0; 3037 double inta = 0; 3038 double v = 0; 3039 double ta = 0; 3040 double tb = 0; 3041 int ns = 0; 3042 double qeps = 0; 3043 int info = 0; 3044 3045 3046 // 3047 // Reverse communication preparations 3048 // I know it looks ugly, but it works the same way 3049 // anywhere from C++ to Python. 3050 // 3051 // This code initializes locals by: 3052 // * random values determined during code 3053 // generation - on first subroutine call 3054 // * values from previous call - on subsequent calls 3055 // 3056 if( state.rstate.stage>=0 ) 3057 { 3058 i = state.rstate.ia[0]; 3059 j = state.rstate.ia[1]; 3060 ns = state.rstate.ia[2]; 3061 info = state.rstate.ia[3]; 3062 c1 = state.rstate.ra[0]; 3063 c2 = state.rstate.ra[1]; 3064 intg = state.rstate.ra[2]; 3065 intk = state.rstate.ra[3]; 3066 inta = state.rstate.ra[4]; 3067 v = state.rstate.ra[5]; 3068 ta = state.rstate.ra[6]; 3069 tb = state.rstate.ra[7]; 3070 qeps = state.rstate.ra[8]; 3071 } 3072 else 3073 { 3074 i = 497; 3075 j = -271; 3076 ns = -581; 3077 info = 745; 3078 c1 = -533; 3079 c2 = -77; 3080 intg = 678; 3081 intk = -293; 3082 inta = 316; 3083 v = 647; 3084 ta = -756; 3085 tb = 830; 3086 qeps = -871; 3087 } 3088 if( state.rstate.stage==0 ) 3089 { 3090 goto lbl_0; 3091 } 3092 if( state.rstate.stage==1 ) 3093 { 3094 goto lbl_1; 3095 } 3096 if( state.rstate.stage==2 ) 3097 { 3098 goto lbl_2; 3099 } 3100 3101 // 3102 // Routine body 3103 // 3104 3105 // 3106 // initialize quadratures. 3107 // use 15-point Gauss-Kronrod formula. 3108 // 3109 state.n = 15; 3110 gkq.gkqgenerategausslegendre(state.n, ref info, ref state.qn, ref state.wk, ref state.wg); 3111 if( info<0 ) 3112 { 3113 state.info = -5; 3114 state.r = 0; 3115 result = false; 3116 return result; 3117 } 3118 state.wr = new double[state.n]; 3119 for(i=0; i<=state.n-1; i++) 3120 { 3121 if( i==0 ) 3122 { 3123 state.wr[i] = 0.5*Math.Abs(state.qn[1]-state.qn[0]); 3124 continue; 3125 } 3126 if( i==state.n-1 ) 3127 { 3128 state.wr[state.n-1] = 0.5*Math.Abs(state.qn[state.n-1]-state.qn[state.n-2]); 3129 continue; 3130 } 3131 state.wr[i] = 0.5*Math.Abs(state.qn[i-1]-state.qn[i+1]); 3132 } 3133 3134 // 3135 // special case 3136 // 3137 if( (double)(state.a)==(double)(state.b) ) 3138 { 3139 state.info = 1; 3140 state.r = 0; 3141 result = false; 3142 return result; 3143 } 3144 3145 // 3146 // test parameters 3147 // 3148 if( (double)(state.eps)<(double)(0) | (double)(state.xwidth)<(double)(0) ) 3149 { 3150 state.info = -1; 3151 state.r = 0; 3152 result = false; 3153 return result; 3154 } 3155 state.info = 1; 3156 if( (double)(state.eps)==(double)(0) ) 3157 { 3158 state.eps = 100000*math.machineepsilon; 3159 } 3160 3161 // 3162 // First, prepare heap 3163 // * column 0 - absolute error 3164 // * column 1 - integral of a F(x) (calculated using Kronrod extension nodes) 3165 // * column 2 - integral of a |F(x)| (calculated using modified rect. method) 3166 // * column 3 - left boundary of a subinterval 3167 // * column 4 - right boundary of a subinterval 3168 // 3169 if( (double)(state.xwidth)!=(double)(0) ) 3170 { 3171 goto lbl_3; 3172 } 3173 3174 // 3175 // no maximum width requirements 3176 // start from one big subinterval 3177 // 3178 state.heapwidth = 5; 3179 state.heapsize = 1; 3180 state.heapused = 1; 3181 state.heap = new double[state.heapsize, state.heapwidth]; 3182 c1 = 0.5*(state.b-state.a); 3183 c2 = 0.5*(state.b+state.a); 3184 intg = 0; 3185 intk = 0; 3186 inta = 0; 3187 i = 0; 3188 lbl_5: 3189 if( i>state.n-1 ) 3190 { 3191 goto lbl_7; 3192 } 3193 3194 // 3195 // obtain F 3196 // 3197 state.x = c1*state.qn[i]+c2; 3198 state.rstate.stage = 0; 3199 goto lbl_rcomm; 3200 lbl_0: 3201 v = state.f; 3202 3203 // 3204 // Gauss-Kronrod formula 3205 // 3206 intk = intk+v*state.wk[i]; 3207 if( i%2==1 ) 3208 { 3209 intg = intg+v*state.wg[i]; 3210 } 3211 3212 // 3213 // Integral |F(x)| 3214 // Use rectangles method 3215 // 3216 inta = inta+Math.Abs(v)*state.wr[i]; 3217 i = i+1; 3218 goto lbl_5; 3219 lbl_7: 3220 intk = intk*(state.b-state.a)*0.5; 3221 intg = intg*(state.b-state.a)*0.5; 3222 inta = inta*(state.b-state.a)*0.5; 3223 state.heap[0,0] = Math.Abs(intg-intk); 3224 state.heap[0,1] = intk; 3225 state.heap[0,2] = inta; 3226 state.heap[0,3] = state.a; 3227 state.heap[0,4] = state.b; 3228 state.sumerr = state.heap[0,0]; 3229 state.sumabs = Math.Abs(inta); 3230 goto lbl_4; 3231 lbl_3: 3232 3233 // 3234 // maximum subinterval should be no more than XWidth. 3235 // so we create Ceil((B-A)/XWidth)+1 small subintervals 3236 // 3237 ns = (int)Math.Ceiling(Math.Abs(state.b-state.a)/state.xwidth)+1; 3238 state.heapsize = ns; 3239 state.heapused = ns; 3240 state.heapwidth = 5; 3241 state.heap = new double[state.heapsize, state.heapwidth]; 3242 state.sumerr = 0; 3243 state.sumabs = 0; 3244 j = 0; 3245 lbl_8: 3246 if( j>ns-1 ) 3247 { 3248 goto lbl_10; 3249 } 3250 ta = state.a+j*(state.b-state.a)/ns; 3251 tb = state.a+(j+1)*(state.b-state.a)/ns; 3252 c1 = 0.5*(tb-ta); 3253 c2 = 0.5*(tb+ta); 3254 intg = 0; 3255 intk = 0; 3256 inta = 0; 3257 i = 0; 3258 lbl_11: 3259 if( i>state.n-1 ) 3260 { 3261 goto lbl_13; 3262 } 3263 3264 // 3265 // obtain F 3266 // 3267 state.x = c1*state.qn[i]+c2; 3268 state.rstate.stage = 1; 3269 goto lbl_rcomm; 3270 lbl_1: 3271 v = state.f; 3272 3273 // 3274 // Gauss-Kronrod formula 3275 // 3276 intk = intk+v*state.wk[i]; 3277 if( i%2==1 ) 3278 { 3279 intg = intg+v*state.wg[i]; 3280 } 3281 3282 // 3283 // Integral |F(x)| 3284 // Use rectangles method 3285 // 3286 inta = inta+Math.Abs(v)*state.wr[i]; 3287 i = i+1; 3288 goto lbl_11; 3289 lbl_13: 3290 intk = intk*(tb-ta)*0.5; 3291 intg = intg*(tb-ta)*0.5; 3292 inta = inta*(tb-ta)*0.5; 3293 state.heap[j,0] = Math.Abs(intg-intk); 3294 state.heap[j,1] = intk; 3295 state.heap[j,2] = inta; 3296 state.heap[j,3] = ta; 3297 state.heap[j,4] = tb; 3298 state.sumerr = state.sumerr+state.heap[j,0]; 3299 state.sumabs = state.sumabs+Math.Abs(inta); 3300 j = j+1; 3301 goto lbl_8; 3302 lbl_10: 3303 lbl_4: 3304 3305 // 3306 // method iterations 3307 // 3308 lbl_14: 3309 if( false ) 3310 { 3311 goto lbl_15; 3312 } 3313 3314 // 3315 // additional memory if needed 3316 // 3317 if( state.heapused==state.heapsize ) 3318 { 3319 mheapresize(ref state.heap, ref state.heapsize, 4*state.heapsize, state.heapwidth); 3320 } 3321 3322 // 3323 // TODO: every 20 iterations recalculate errors/sums 3324 // 3325 if( (double)(state.sumerr)<=(double)(state.eps*state.sumabs) | state.heapused>=maxsubintervals ) 3326 { 3327 state.r = 0; 3328 for(j=0; j<=state.heapused-1; j++) 3329 { 3330 state.r = state.r+state.heap[j,1]; 3331 } 3332 result = false; 3333 return result; 3334 } 3335 3336 // 3337 // Exclude interval with maximum absolute error 3338 // 3339 mheappop(ref state.heap, state.heapused, state.heapwidth); 3340 state.sumerr = state.sumerr-state.heap[state.heapused-1,0]; 3341 state.sumabs = state.sumabs-state.heap[state.heapused-1,2]; 3342 3343 // 3344 // Divide interval, create subintervals 3345 // 3346 ta = state.heap[state.heapused-1,3]; 3347 tb = state.heap[state.heapused-1,4]; 3348 state.heap[state.heapused-1,3] = ta; 3349 state.heap[state.heapused-1,4] = 0.5*(ta+tb); 3350 state.heap[state.heapused,3] = 0.5*(ta+tb); 3351 state.heap[state.heapused,4] = tb; 3352 j = state.heapused-1; 3353 lbl_16: 3354 if( j>state.heapused ) 3355 { 3356 goto lbl_18; 3357 } 3358 c1 = 0.5*(state.heap[j,4]-state.heap[j,3]); 3359 c2 = 0.5*(state.heap[j,4]+state.heap[j,3]); 3360 intg = 0; 3361 intk = 0; 3362 inta = 0; 3363 i = 0; 3364 lbl_19: 3365 if( i>state.n-1 ) 3366 { 3367 goto lbl_21; 3368 } 3369 3370 // 3371 // F(x) 3372 // 3373 state.x = c1*state.qn[i]+c2; 3374 state.rstate.stage = 2; 3375 goto lbl_rcomm; 3376 lbl_2: 3377 v = state.f; 3378 3379 // 3380 // Gauss-Kronrod formula 3381 // 3382 intk = intk+v*state.wk[i]; 3383 if( i%2==1 ) 3384 { 3385 intg = intg+v*state.wg[i]; 3386 } 3387 3388 // 3389 // Integral |F(x)| 3390 // Use rectangles method 3391 // 3392 inta = inta+Math.Abs(v)*state.wr[i]; 3393 i = i+1; 3394 goto lbl_19; 3395 lbl_21: 3396 intk = intk*(state.heap[j,4]-state.heap[j,3])*0.5; 3397 intg = intg*(state.heap[j,4]-state.heap[j,3])*0.5; 3398 inta = inta*(state.heap[j,4]-state.heap[j,3])*0.5; 3399 state.heap[j,0] = Math.Abs(intg-intk); 3400 state.heap[j,1] = intk; 3401 state.heap[j,2] = inta; 3402 state.sumerr = state.sumerr+state.heap[j,0]; 3403 state.sumabs = state.sumabs+state.heap[j,2]; 3404 j = j+1; 3405 goto lbl_16; 3406 lbl_18: 3407 mheappush(ref state.heap, state.heapused-1, state.heapwidth); 3408 mheappush(ref state.heap, state.heapused, state.heapwidth); 3409 state.heapused = state.heapused+1; 3410 goto lbl_14; 3411 lbl_15: 3412 result = false; 3413 return result; 3414 3415 // 3416 // Saving state 3417 // 3418 lbl_rcomm: 3419 result = true; 3420 state.rstate.ia[0] = i; 3421 state.rstate.ia[1] = j; 3422 state.rstate.ia[2] = ns; 3423 state.rstate.ia[3] = info; 3424 state.rstate.ra[0] = c1; 3425 state.rstate.ra[1] = c2; 3426 state.rstate.ra[2] = intg; 3427 state.rstate.ra[3] = intk; 3428 state.rstate.ra[4] = inta; 3429 state.rstate.ra[5] = v; 3430 state.rstate.ra[6] = ta; 3431 state.rstate.ra[7] = tb; 3432 state.rstate.ra[8] = qeps; 3433 return result; 3434 } 3435 3436 3437 private static void mheappop(ref double[,] heap, 3438 int heapsize, 3439 int heapwidth) 3440 { 3441 int i = 0; 3442 int p = 0; 3443 double t = 0; 3444 int maxcp = 0; 3445 3446 if( heapsize==1 ) 3447 { 3448 return; 3449 } 3450 for(i=0; i<=heapwidth-1; i++) 3451 { 3452 t = heap[heapsize-1,i]; 3453 heap[heapsize-1,i] = heap[0,i]; 3454 heap[0,i] = t; 3455 } 3456 p = 0; 3457 while( 2*p+1<heapsize-1 ) 3458 { 3459 maxcp = 2*p+1; 3460 if( 2*p+2<heapsize-1 ) 3461 { 3462 if( (double)(heap[2*p+2,0])>(double)(heap[2*p+1,0]) ) 3463 { 3464 maxcp = 2*p+2; 3465 } 3466 } 3467 if( (double)(heap[p,0])<(double)(heap[maxcp,0]) ) 3468 { 3469 for(i=0; i<=heapwidth-1; i++) 3470 { 3471 t = heap[p,i]; 3472 heap[p,i] = heap[maxcp,i]; 3473 heap[maxcp,i] = t; 3474 } 3475 p = maxcp; 3476 } 3477 else 3478 { 3479 break; 3480 } 3481 } 3482 } 3483 3484 3485 private static void mheappush(ref double[,] heap, 3486 int heapsize, 3487 int heapwidth) 3488 { 3489 int i = 0; 3490 int p = 0; 3491 double t = 0; 3492 int parent = 0; 3493 3494 if( heapsize==0 ) 3495 { 3496 return; 3497 } 3498 p = heapsize; 3499 while( p!=0 ) 3500 { 3501 parent = (p-1)/2; 3502 if( (double)(heap[p,0])>(double)(heap[parent,0]) ) 3503 { 3504 for(i=0; i<=heapwidth-1; i++) 3505 { 3506 t = heap[p,i]; 3507 heap[p,i] = heap[parent,i]; 3508 heap[parent,i] = t; 3509 } 3510 p = parent; 3511 } 3512 else 3513 { 3514 break; 3515 } 3516 } 3517 } 3518 3519 3520 private static void mheapresize(ref double[,] heap, 3521 ref int heapsize, 3522 int newheapsize, 3523 int heapwidth) 3524 { 3525 double[,] tmp = new double[0,0]; 3526 int i = 0; 3527 int i_ = 0; 3528 3529 tmp = new double[heapsize, heapwidth]; 3530 for(i=0; i<=heapsize-1; i++) 3531 { 3532 for(i_=0; i_<=heapwidth-1;i_++) 3533 { 3534 tmp[i,i_] = heap[i,i_]; 3535 } 3536 } 3537 heap = new double[newheapsize, heapwidth]; 3538 for(i=0; i<=heapsize-1; i++) 3539 { 3540 for(i_=0; i_<=heapwidth-1;i_++) 3541 { 3542 heap[i,i_] = tmp[i,i_]; 3543 } 3544 } 3545 heapsize = newheapsize; 3546 } 3547 3548 3549 } 3548 3550 } 3549 3551
Note: See TracChangeset
for help on using the changeset viewer.
