Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.ALGLIB-2.5.0/ALGLIB-2.5.0/trigintegrals.cs @ 5192

Last change on this file since 5192 was 3839, checked in by mkommend, 14 years ago

implemented first version of LR (ticket #1012)
File size: 22.3 KB
Line 
1/*************************************************************************
2Cephes Math Library Release 2.8:  June, 2000
3Copyright by Stephen L. Moshier
4
5Contributors:
6    * Sergey Bochkanov (ALGLIB project). Translation from C to
7      pseudocode.
8
9See subroutines comments for additional copyrights.
10
11>>> SOURCE LICENSE >>>
12This program is free software; you can redistribute it and/or modify
13it under the terms of the GNU General Public License as published by
14the Free Software Foundation (www.fsf.org); either version 2 of the
15License, or (at your option) any later version.
16
17This program is distributed in the hope that it will be useful,
18but WITHOUT ANY WARRANTY; without even the implied warranty of
19MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20GNU General Public License for more details.
21
22A copy of the GNU General Public License is available at
23http://www.fsf.org/licensing/licenses
24
25>>> END OF LICENSE >>>
26*************************************************************************/
27
28using System;
29
30namespace alglib
31{
32    public class trigintegrals
33    {
34        /*************************************************************************
35        Sine and cosine integrals
36
37        Evaluates the integrals
38
39                                 x
40                                 -
41                                |  cos t - 1
42          Ci(x) = eul + ln x +  |  --------- dt,
43                                |      t
44                               -
45                                0
46                    x
47                    -
48                   |  sin t
49          Si(x) =  |  ----- dt
50                   |    t
51                  -
52                   0
53
54        where eul = 0.57721566490153286061 is Euler's constant.
55        The integrals are approximated by rational functions.
56        For x > 8 auxiliary functions f(x) and g(x) are employed
57        such that
58
59        Ci(x) = f(x) sin(x) - g(x) cos(x)
60        Si(x) = pi/2 - f(x) cos(x) - g(x) sin(x)
61
62
63        ACCURACY:
64           Test interval = [0,50].
65        Absolute error, except relative when > 1:
66        arithmetic   function   # trials      peak         rms
67           IEEE        Si        30000       4.4e-16     7.3e-17
68           IEEE        Ci        30000       6.9e-16     5.1e-17
69
70        Cephes Math Library Release 2.1:  January, 1989
71        Copyright 1984, 1987, 1989 by Stephen L. Moshier
72        *************************************************************************/
73        public static void sinecosineintegrals(double x,
74            ref double si,
75            ref double ci)
76        {
77            double z = 0;
78            double c = 0;
79            double s = 0;
80            double f = 0;
81            double g = 0;
82            int sg = 0;
83            double sn = 0;
84            double sd = 0;
85            double cn = 0;
86            double cd = 0;
87            double fn = 0;
88            double fd = 0;
89            double gn = 0;
90            double gd = 0;
91
92            if( (double)(x)<(double)(0) )
93            {
94                sg = -1;
95                x = -x;
96            }
97            else
98            {
99                sg = 0;
100            }
101            if( (double)(x)==(double)(0) )
102            {
103                si = 0;
104                ci = -AP.Math.MaxRealNumber;
105                return;
106            }
107            if( (double)(x)>(double)(1.0E9) )
108            {
109                si = 1.570796326794896619-Math.Cos(x)/x;
110                ci = Math.Sin(x)/x;
111                return;
112            }
113            if( (double)(x)<=(double)(4) )
114            {
115                z = x*x;
116                sn = -8.39167827910303881427E-11;
117                sn = sn*z+4.62591714427012837309E-8;
118                sn = sn*z-9.75759303843632795789E-6;
119                sn = sn*z+9.76945438170435310816E-4;
120                sn = sn*z-4.13470316229406538752E-2;
121                sn = sn*z+1.00000000000000000302E0;
122                sd = 2.03269266195951942049E-12;
123                sd = sd*z+1.27997891179943299903E-9;
124                sd = sd*z+4.41827842801218905784E-7;
125                sd = sd*z+9.96412122043875552487E-5;
126                sd = sd*z+1.42085239326149893930E-2;
127                sd = sd*z+9.99999999999999996984E-1;
128                s = x*sn/sd;
129                cn = 2.02524002389102268789E-11;
130                cn = cn*z-1.35249504915790756375E-8;
131                cn = cn*z+3.59325051419993077021E-6;
132                cn = cn*z-4.74007206873407909465E-4;
133                cn = cn*z+2.89159652607555242092E-2;
134                cn = cn*z-1.00000000000000000080E0;
135                cd = 4.07746040061880559506E-12;
136                cd = cd*z+3.06780997581887812692E-9;
137                cd = cd*z+1.23210355685883423679E-6;
138                cd = cd*z+3.17442024775032769882E-4;
139                cd = cd*z+5.10028056236446052392E-2;
140                cd = cd*z+4.00000000000000000080E0;
141                c = z*cn/cd;
142                if( sg!=0 )
143                {
144                    s = -s;
145                }
146                si = s;
147                ci = 0.57721566490153286061+Math.Log(x)+c;
148                return;
149            }
150            s = Math.Sin(x);
151            c = Math.Cos(x);
152            z = 1.0/(x*x);
153            if( (double)(x)<(double)(8) )
154            {
155                fn = 4.23612862892216586994E0;
156                fn = fn*z+5.45937717161812843388E0;
157                fn = fn*z+1.62083287701538329132E0;
158                fn = fn*z+1.67006611831323023771E-1;
159                fn = fn*z+6.81020132472518137426E-3;
160                fn = fn*z+1.08936580650328664411E-4;
161                fn = fn*z+5.48900223421373614008E-7;
162                fd = 1.00000000000000000000E0;
163                fd = fd*z+8.16496634205391016773E0;
164                fd = fd*z+7.30828822505564552187E0;
165                fd = fd*z+1.86792257950184183883E0;
166                fd = fd*z+1.78792052963149907262E-1;
167                fd = fd*z+7.01710668322789753610E-3;
168                fd = fd*z+1.10034357153915731354E-4;
169                fd = fd*z+5.48900252756255700982E-7;
170                f = fn/(x*fd);
171                gn = 8.71001698973114191777E-2;
172                gn = gn*z+6.11379109952219284151E-1;
173                gn = gn*z+3.97180296392337498885E-1;
174                gn = gn*z+7.48527737628469092119E-2;
175                gn = gn*z+5.38868681462177273157E-3;
176                gn = gn*z+1.61999794598934024525E-4;
177                gn = gn*z+1.97963874140963632189E-6;
178                gn = gn*z+7.82579040744090311069E-9;
179                gd = 1.00000000000000000000E0;
180                gd = gd*z+1.64402202413355338886E0;
181                gd = gd*z+6.66296701268987968381E-1;
182                gd = gd*z+9.88771761277688796203E-2;
183                gd = gd*z+6.22396345441768420760E-3;
184                gd = gd*z+1.73221081474177119497E-4;
185                gd = gd*z+2.02659182086343991969E-6;
186                gd = gd*z+7.82579218933534490868E-9;
187                g = z*gn/gd;
188            }
189            else
190            {
191                fn = 4.55880873470465315206E-1;
192                fn = fn*z+7.13715274100146711374E-1;
193                fn = fn*z+1.60300158222319456320E-1;
194                fn = fn*z+1.16064229408124407915E-2;
195                fn = fn*z+3.49556442447859055605E-4;
196                fn = fn*z+4.86215430826454749482E-6;
197                fn = fn*z+3.20092790091004902806E-8;
198                fn = fn*z+9.41779576128512936592E-11;
199                fn = fn*z+9.70507110881952024631E-14;
200                fd = 1.00000000000000000000E0;
201                fd = fd*z+9.17463611873684053703E-1;
202                fd = fd*z+1.78685545332074536321E-1;
203                fd = fd*z+1.22253594771971293032E-2;
204                fd = fd*z+3.58696481881851580297E-4;
205                fd = fd*z+4.92435064317881464393E-6;
206                fd = fd*z+3.21956939101046018377E-8;
207                fd = fd*z+9.43720590350276732376E-11;
208                fd = fd*z+9.70507110881952025725E-14;
209                f = fn/(x*fd);
210                gn = 6.97359953443276214934E-1;
211                gn = gn*z+3.30410979305632063225E-1;
212                gn = gn*z+3.84878767649974295920E-2;
213                gn = gn*z+1.71718239052347903558E-3;
214                gn = gn*z+3.48941165502279436777E-5;
215                gn = gn*z+3.47131167084116673800E-7;
216                gn = gn*z+1.70404452782044526189E-9;
217                gn = gn*z+3.85945925430276600453E-12;
218                gn = gn*z+3.14040098946363334640E-15;
219                gd = 1.00000000000000000000E0;
220                gd = gd*z+1.68548898811011640017E0;
221                gd = gd*z+4.87852258695304967486E-1;
222                gd = gd*z+4.67913194259625806320E-2;
223                gd = gd*z+1.90284426674399523638E-3;
224                gd = gd*z+3.68475504442561108162E-5;
225                gd = gd*z+3.57043223443740838771E-7;
226                gd = gd*z+1.72693748966316146736E-9;
227                gd = gd*z+3.87830166023954706752E-12;
228                gd = gd*z+3.14040098946363335242E-15;
229                g = z*gn/gd;
230            }
231            si = 1.570796326794896619-f*c-g*s;
232            if( sg!=0 )
233            {
234                si = -si;
235            }
236            ci = f*s-g*c;
237        }
238
239
240        /*************************************************************************
241        Hyperbolic sine and cosine integrals
242
243        Approximates the integrals
244
245                                   x
246                                   -
247                                  | |   cosh t - 1
248          Chi(x) = eul + ln x +   |    -----------  dt,
249                                | |          t
250                                 -
251                                 0
252
253                      x
254                      -
255                     | |  sinh t
256          Shi(x) =   |    ------  dt
257                   | |       t
258                    -
259                    0
260
261        where eul = 0.57721566490153286061 is Euler's constant.
262        The integrals are evaluated by power series for x < 8
263        and by Chebyshev expansions for x between 8 and 88.
264        For large x, both functions approach exp(x)/2x.
265        Arguments greater than 88 in magnitude return MAXNUM.
266
267
268        ACCURACY:
269
270        Test interval 0 to 88.
271                             Relative error:
272        arithmetic   function  # trials      peak         rms
273           IEEE         Shi      30000       6.9e-16     1.6e-16
274               Absolute error, except relative when |Chi| > 1:
275           IEEE         Chi      30000       8.4e-16     1.4e-16
276
277        Cephes Math Library Release 2.8:  June, 2000
278        Copyright 1984, 1987, 2000 by Stephen L. Moshier
279        *************************************************************************/
280        public static void hyperbolicsinecosineintegrals(double x,
281            ref double shi,
282            ref double chi)
283        {
284            double k = 0;
285            double z = 0;
286            double c = 0;
287            double s = 0;
288            double a = 0;
289            int sg = 0;
290            double b0 = 0;
291            double b1 = 0;
292            double b2 = 0;
293
294            if( (double)(x)<(double)(0) )
295            {
296                sg = -1;
297                x = -x;
298            }
299            else
300            {
301                sg = 0;
302            }
303            if( (double)(x)==(double)(0) )
304            {
305                shi = 0;
306                chi = -AP.Math.MaxRealNumber;
307                return;
308            }
309            if( (double)(x)<(double)(8.0) )
310            {
311                z = x*x;
312                a = 1.0;
313                s = 1.0;
314                c = 0.0;
315                k = 2.0;
316                do
317                {
318                    a = a*z/k;
319                    c = c+a/k;
320                    k = k+1.0;
321                    a = a/k;
322                    s = s+a/k;
323                    k = k+1.0;
324                }
325                while( (double)(Math.Abs(a/s))>=(double)(AP.Math.MachineEpsilon) );
326                s = s*x;
327            }
328            else
329            {
330                if( (double)(x)<(double)(18.0) )
331                {
332                    a = (576.0/x-52.0)/10.0;
333                    k = Math.Exp(x)/x;
334                    b0 = 1.83889230173399459482E-17;
335                    b1 = 0.0;
336                    chebiterationshichi(a, -9.55485532279655569575E-17, ref b0, ref b1, ref b2);
337                    chebiterationshichi(a, 2.04326105980879882648E-16, ref b0, ref b1, ref b2);
338                    chebiterationshichi(a, 1.09896949074905343022E-15, ref b0, ref b1, ref b2);
339                    chebiterationshichi(a, -1.31313534344092599234E-14, ref b0, ref b1, ref b2);
340                    chebiterationshichi(a, 5.93976226264314278932E-14, ref b0, ref b1, ref b2);
341                    chebiterationshichi(a, -3.47197010497749154755E-14, ref b0, ref b1, ref b2);
342                    chebiterationshichi(a, -1.40059764613117131000E-12, ref b0, ref b1, ref b2);
343                    chebiterationshichi(a, 9.49044626224223543299E-12, ref b0, ref b1, ref b2);
344                    chebiterationshichi(a, -1.61596181145435454033E-11, ref b0, ref b1, ref b2);
345                    chebiterationshichi(a, -1.77899784436430310321E-10, ref b0, ref b1, ref b2);
346                    chebiterationshichi(a, 1.35455469767246947469E-9, ref b0, ref b1, ref b2);
347                    chebiterationshichi(a, -1.03257121792819495123E-9, ref b0, ref b1, ref b2);
348                    chebiterationshichi(a, -3.56699611114982536845E-8, ref b0, ref b1, ref b2);
349                    chebiterationshichi(a, 1.44818877384267342057E-7, ref b0, ref b1, ref b2);
350                    chebiterationshichi(a, 7.82018215184051295296E-7, ref b0, ref b1, ref b2);
351                    chebiterationshichi(a, -5.39919118403805073710E-6, ref b0, ref b1, ref b2);
352                    chebiterationshichi(a, -3.12458202168959833422E-5, ref b0, ref b1, ref b2);
353                    chebiterationshichi(a, 8.90136741950727517826E-5, ref b0, ref b1, ref b2);
354                    chebiterationshichi(a, 2.02558474743846862168E-3, ref b0, ref b1, ref b2);
355                    chebiterationshichi(a, 2.96064440855633256972E-2, ref b0, ref b1, ref b2);
356                    chebiterationshichi(a, 1.11847751047257036625E0, ref b0, ref b1, ref b2);
357                    s = k*0.5*(b0-b2);
358                    b0 = -8.12435385225864036372E-18;
359                    b1 = 0.0;
360                    chebiterationshichi(a, 2.17586413290339214377E-17, ref b0, ref b1, ref b2);
361                    chebiterationshichi(a, 5.22624394924072204667E-17, ref b0, ref b1, ref b2);
362                    chebiterationshichi(a, -9.48812110591690559363E-16, ref b0, ref b1, ref b2);
363                    chebiterationshichi(a, 5.35546311647465209166E-15, ref b0, ref b1, ref b2);
364                    chebiterationshichi(a, -1.21009970113732918701E-14, ref b0, ref b1, ref b2);
365                    chebiterationshichi(a, -6.00865178553447437951E-14, ref b0, ref b1, ref b2);
366                    chebiterationshichi(a, 7.16339649156028587775E-13, ref b0, ref b1, ref b2);
367                    chebiterationshichi(a, -2.93496072607599856104E-12, ref b0, ref b1, ref b2);
368                    chebiterationshichi(a, -1.40359438136491256904E-12, ref b0, ref b1, ref b2);
369                    chebiterationshichi(a, 8.76302288609054966081E-11, ref b0, ref b1, ref b2);
370                    chebiterationshichi(a, -4.40092476213282340617E-10, ref b0, ref b1, ref b2);
371                    chebiterationshichi(a, -1.87992075640569295479E-10, ref b0, ref b1, ref b2);
372                    chebiterationshichi(a, 1.31458150989474594064E-8, ref b0, ref b1, ref b2);
373                    chebiterationshichi(a, -4.75513930924765465590E-8, ref b0, ref b1, ref b2);
374                    chebiterationshichi(a, -2.21775018801848880741E-7, ref b0, ref b1, ref b2);
375                    chebiterationshichi(a, 1.94635531373272490962E-6, ref b0, ref b1, ref b2);
376                    chebiterationshichi(a, 4.33505889257316408893E-6, ref b0, ref b1, ref b2);
377                    chebiterationshichi(a, -6.13387001076494349496E-5, ref b0, ref b1, ref b2);
378                    chebiterationshichi(a, -3.13085477492997465138E-4, ref b0, ref b1, ref b2);
379                    chebiterationshichi(a, 4.97164789823116062801E-4, ref b0, ref b1, ref b2);
380                    chebiterationshichi(a, 2.64347496031374526641E-2, ref b0, ref b1, ref b2);
381                    chebiterationshichi(a, 1.11446150876699213025E0, ref b0, ref b1, ref b2);
382                    c = k*0.5*(b0-b2);
383                }
384                else
385                {
386                    if( (double)(x)<=(double)(88.0) )
387                    {
388                        a = (6336.0/x-212.0)/70.0;
389                        k = Math.Exp(x)/x;
390                        b0 = -1.05311574154850938805E-17;
391                        b1 = 0.0;
392                        chebiterationshichi(a, 2.62446095596355225821E-17, ref b0, ref b1, ref b2);
393                        chebiterationshichi(a, 8.82090135625368160657E-17, ref b0, ref b1, ref b2);
394                        chebiterationshichi(a, -3.38459811878103047136E-16, ref b0, ref b1, ref b2);
395                        chebiterationshichi(a, -8.30608026366935789136E-16, ref b0, ref b1, ref b2);
396                        chebiterationshichi(a, 3.93397875437050071776E-15, ref b0, ref b1, ref b2);
397                        chebiterationshichi(a, 1.01765565969729044505E-14, ref b0, ref b1, ref b2);
398                        chebiterationshichi(a, -4.21128170307640802703E-14, ref b0, ref b1, ref b2);
399                        chebiterationshichi(a, -1.60818204519802480035E-13, ref b0, ref b1, ref b2);
400                        chebiterationshichi(a, 3.34714954175994481761E-13, ref b0, ref b1, ref b2);
401                        chebiterationshichi(a, 2.72600352129153073807E-12, ref b0, ref b1, ref b2);
402                        chebiterationshichi(a, 1.66894954752839083608E-12, ref b0, ref b1, ref b2);
403                        chebiterationshichi(a, -3.49278141024730899554E-11, ref b0, ref b1, ref b2);
404                        chebiterationshichi(a, -1.58580661666482709598E-10, ref b0, ref b1, ref b2);
405                        chebiterationshichi(a, -1.79289437183355633342E-10, ref b0, ref b1, ref b2);
406                        chebiterationshichi(a, 1.76281629144264523277E-9, ref b0, ref b1, ref b2);
407                        chebiterationshichi(a, 1.69050228879421288846E-8, ref b0, ref b1, ref b2);
408                        chebiterationshichi(a, 1.25391771228487041649E-7, ref b0, ref b1, ref b2);
409                        chebiterationshichi(a, 1.16229947068677338732E-6, ref b0, ref b1, ref b2);
410                        chebiterationshichi(a, 1.61038260117376323993E-5, ref b0, ref b1, ref b2);
411                        chebiterationshichi(a, 3.49810375601053973070E-4, ref b0, ref b1, ref b2);
412                        chebiterationshichi(a, 1.28478065259647610779E-2, ref b0, ref b1, ref b2);
413                        chebiterationshichi(a, 1.03665722588798326712E0, ref b0, ref b1, ref b2);
414                        s = k*0.5*(b0-b2);
415                        b0 = 8.06913408255155572081E-18;
416                        b1 = 0.0;
417                        chebiterationshichi(a, -2.08074168180148170312E-17, ref b0, ref b1, ref b2);
418                        chebiterationshichi(a, -5.98111329658272336816E-17, ref b0, ref b1, ref b2);
419                        chebiterationshichi(a, 2.68533951085945765591E-16, ref b0, ref b1, ref b2);
420                        chebiterationshichi(a, 4.52313941698904694774E-16, ref b0, ref b1, ref b2);
421                        chebiterationshichi(a, -3.10734917335299464535E-15, ref b0, ref b1, ref b2);
422                        chebiterationshichi(a, -4.42823207332531972288E-15, ref b0, ref b1, ref b2);
423                        chebiterationshichi(a, 3.49639695410806959872E-14, ref b0, ref b1, ref b2);
424                        chebiterationshichi(a, 6.63406731718911586609E-14, ref b0, ref b1, ref b2);
425                        chebiterationshichi(a, -3.71902448093119218395E-13, ref b0, ref b1, ref b2);
426                        chebiterationshichi(a, -1.27135418132338309016E-12, ref b0, ref b1, ref b2);
427                        chebiterationshichi(a, 2.74851141935315395333E-12, ref b0, ref b1, ref b2);
428                        chebiterationshichi(a, 2.33781843985453438400E-11, ref b0, ref b1, ref b2);
429                        chebiterationshichi(a, 2.71436006377612442764E-11, ref b0, ref b1, ref b2);
430                        chebiterationshichi(a, -2.56600180000355990529E-10, ref b0, ref b1, ref b2);
431                        chebiterationshichi(a, -1.61021375163803438552E-9, ref b0, ref b1, ref b2);
432                        chebiterationshichi(a, -4.72543064876271773512E-9, ref b0, ref b1, ref b2);
433                        chebiterationshichi(a, -3.00095178028681682282E-9, ref b0, ref b1, ref b2);
434                        chebiterationshichi(a, 7.79387474390914922337E-8, ref b0, ref b1, ref b2);
435                        chebiterationshichi(a, 1.06942765566401507066E-6, ref b0, ref b1, ref b2);
436                        chebiterationshichi(a, 1.59503164802313196374E-5, ref b0, ref b1, ref b2);
437                        chebiterationshichi(a, 3.49592575153777996871E-4, ref b0, ref b1, ref b2);
438                        chebiterationshichi(a, 1.28475387530065247392E-2, ref b0, ref b1, ref b2);
439                        chebiterationshichi(a, 1.03665693917934275131E0, ref b0, ref b1, ref b2);
440                        c = k*0.5*(b0-b2);
441                    }
442                    else
443                    {
444                        if( sg!=0 )
445                        {
446                            shi = -AP.Math.MaxRealNumber;
447                        }
448                        else
449                        {
450                            shi = AP.Math.MaxRealNumber;
451                        }
452                        chi = AP.Math.MaxRealNumber;
453                        return;
454                    }
455                }
456            }
457            if( sg!=0 )
458            {
459                s = -s;
460            }
461            shi = s;
462            chi = 0.57721566490153286061+Math.Log(x)+c;
463        }
464
465
466        private static void chebiterationshichi(double x,
467            double c,
468            ref double b0,
469            ref double b1,
470            ref double b2)
471        {
472            b2 = b1;
473            b1 = b0;
474            b0 = x*b1-b2+c;
475        }
476    }
477}
Note: See TracBrowser for help on using the repository browser.