Plugins
Loading...
Searching...
No Matches
gamma_function.h
Go to the documentation of this file.
1/*
2 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
3 * SPDX-License-Identifier: EUPL-1.2
4 * SPDX-FileContributor: Sabine Stichel
5 * SPDX-FileContributor: Goethe Universität Frankfurt
6 * SPDX-FileType: SOURCE
7 *
8 * This file is part of d3f++.
9 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
10 */
11
12
13/*
14 * Original by J.Musuuza
15 */
16
17// extern headers
18#include <iostream>
19#include <sstream>
20#include <fstream>
21#include <string>
22#include <cmath>
23#include <vector>
24
25
26#define _USE_MATH_DEFINES
27#include <math.h>
28#include <float.h> //for DBL_MAX
29#include <iomanip>
30#include <cstdlib>
31#include <time.h>
32#include "common/common.h"
33
37#include "bridge/bridge.h"
38#include "lib_disc/domain.h"
48#include "common/parser/rapidxml/rapidxml_print.hpp"
49#include "common/parser/rapidxml/rapidxml.hpp"
50#include "../plugins/experimental/biogas/save_load.h"
51
52
53//using namespace std;
54//using namespace rapidxml;
55#ifdef UG_FOR_LUA
57#endif
58
59namespace ug{
60namespace d3f{
61
62// Internally Defined Constants //
63template <typename TData, int dim, typename TRet>
65{
66 // Bernoulli numbers B(2),B(4),B(6),...,B(20). Only B(2),...,B(6) currently used.//
67 B[0] = 1.0L / (long double)(6 * 2 * 1);
68 B[1] = -1.0L / (long double)(30 * 4 * 3);
69 B[2] = 1.0L / (long double)(42 * 6 * 5);
70 B[3] = -1.0L / (long double)(30 * 8 * 7);
71 B[4] = 5.0L / (long double)(66 * 10 * 9);
72 B[5] = -691.0L / (long double)(2730 * 12 * 11);
73 B[6] = 7.0L / (long double)(6 * 14 * 13);
74 B[7] = -3617.0L / (long double)(510 * 16 * 15);
75 B[8] = 43867.0L / (long double)(796 * 18 * 17);
76 B[9] = -174611.0L / (long double)(330 * 20 * 19);
77
78 n = sizeof(B) / sizeof(long double);
79
80 a[0] = +1.14400529453851095667309e+4L;
81 a[1] = -3.23988020152318335053598e+4L;
82 a[2] = +3.50514523505571666566083e+4L;
83 a[3] = -1.81641309541260702610647e+4L;
84 a[4] = +4.63232990536666818409138e+3L;
85 a[5] = -5.36976777703356780555748e+2L;
86 a[6] = +2.28754473395181007645155e+1L;
87 a[7] = -2.17925748738865115560082e-1L;
88 a[8] = +1.08314836272589368860689e-4L;
89
90 e = 2.71828182845904523536028747L;
91 log_sqrt_2pi = 9.18938533204672741780329736e-1L;
92 g = 9.65657815377331589457187L;
93 exp_g_o_sqrt_2pi = +6.23316569877722552586386e+3L;
94 max_double_arg = 171.0;
95 max_long_double_arg = 1755.5L;
96
97 factorials[0] = 1.000000000000000000000e+0L; // 0!
98 factorials[1] = 1.000000000000000000000e+0L; // 1!
99 factorials[2] = 2.000000000000000000000e+0L; // 2!
100 factorials[3] = 6.000000000000000000000e+0L; // 3!
101 factorials[4] = 2.400000000000000000000e+1L; // 4!
102 factorials[5] = 1.200000000000000000000e+2L; // 5!
103 factorials[6] =7.200000000000000000000e+2L; // 6!
104 factorials[7] =5.040000000000000000000e+3L; // 7!
105 factorials[8] =4.032000000000000000000e+4L; // 8!
106 factorials[9] =3.628800000000000000000e+5L; // 9!
107 factorials[10] =3.628800000000000000000e+6L; // 10!
108 factorials[11] = 3.991680000000000000000e+7L; // 11!
109 factorials[12] = 4.790016000000000000000e+8L; // 12!
110 factorials[13] = 6.227020800000000000000e+9L; // 13!
111 factorials[14] = 8.717829120000000000000e+10L; // 14!
112 factorials[15] = 1.307674368000000000000e+12L; // 15!
113 factorials[16] = 2.092278988800000000000e+13L; // 16!
114 factorials[17] = 3.556874280960000000000e+14L; // 17!
115 factorials[18] = 6.402373705728000000000e+15L; // 18!
116 factorials[19] = 1.216451004088320000000e+17L; // 19!
117 factorials[20] = 2.432902008176640000000e+18L; // 20!
118 factorials[21] = 5.109094217170944000000e+19L; // 21!
119 factorials[22] = 1.124000727777607680000e+21L; // 22!
120 factorials[23] = 2.585201673888497664000e+22L; // 23!
121 factorials[24] = 6.204484017332394393600e+23L; // 24!
122 factorials[25] = 1.551121004333098598400e+25L; // 25!
123 factorials[26] = 4.032914611266056355840e+26L; // 26!
124 factorials[27] = 1.088886945041835216077e+28L; // 27!
125 factorials[28] = 3.048883446117138605015e+29L; // 28!
126 factorials[29] = 8.841761993739701954544e+30L; // 29!
127 factorials[30] = 2.652528598121910586363e+32L; // 30!
128 factorials[31] = 8.222838654177922817726e+33L; // 31!
129 factorials[32] = 2.631308369336935301672e+35L; // 32!
130 factorials[33] = 8.683317618811886495518e+36L; // 33!
131 factorials[34] = 2.952327990396041408476e+38L; // 34!
132 factorials[35] = 1.033314796638614492967e+40L; // 35!
133 factorials[36] = 3.719933267899012174680e+41L; // 36!
134 factorials[37] = 1.376375309122634504632e+43L; // 37!
135 factorials[38] = 5.230226174666011117600e+44L; // 38!
136 factorials[39] = 2.039788208119744335864e+46L; // 39!
137 factorials[40] = 8.159152832478977343456e+47L; // 40!
138 factorials[41] = 3.345252661316380710817e+49L; // 41!
139 factorials[42] = 1.405006117752879898543e+51L; // 42!
140 factorials[43] = 6.041526306337383563736e+52L; // 43!
141 factorials[44] = 2.658271574788448768044e+54L; // 44!
142 factorials[45] = 1.196222208654801945620e+56L; // 45!
143 factorials[46] = 5.502622159812088949850e+57L; // 46!
144 factorials[47] = 2.586232415111681806430e+59L; // 47!
145 factorials[48] = 1.241391559253607267086e+61L; // 48!
146 factorials[49] = 6.082818640342675608723e+62L; // 49!
147 factorials[50] = 3.041409320171337804361e+64L; // 50!
148 factorials[51] = 1.551118753287382280224e+66L; // 51!
149 factorials[52] = 8.065817517094387857166e+67L; // 52!
150 factorials[53] = 4.274883284060025564298e+69L; // 53!
151 factorials[54] = 2.308436973392413804721e+71L; // 54!
152 factorials[55] = 1.269640335365827592597e+73L; // 55!
153 factorials[56] = 7.109985878048634518540e+74L; // 56!
154 factorials[57] = 4.052691950487721675568e+76L; // 57!
155 factorials[58] = 2.350561331282878571829e+78L; // 58!
156 factorials[59] = 1.386831185456898357379e+80L; // 59!
157 factorials[60] = 8.320987112741390144276e+81L; // 60!
158 factorials[61] = 5.075802138772247988009e+83L; // 61!
159 factorials[62] = 3.146997326038793752565e+85L; // 62!
160 factorials[63] = 1.982608315404440064116e+87L; // 63!
161 factorials[64] = 1.268869321858841641034e+89L; // 64!
162 factorials[65] = 8.247650592082470666723e+90L; // 65!
163 factorials[66] = 5.443449390774430640037e+92L; // 66!
164 factorials[67] = 3.647111091818868528825e+94L; // 67!
165 factorials[68] = 2.480035542436830599601e+96L; // 68!
166 factorials[69] = 1.711224524281413113725e+98L; // 69!
167 factorials[70] = 1.197857166996989179607e+100L; // 70!
168 factorials[71] = 8.504785885678623175212e+101L; // 71!
169 factorials[72] = 6.123445837688608686152e+103L; // 72!
170 factorials[73] = 4.470115461512684340891e+105L; // 73!
171 factorials[74] = 3.307885441519386412260e+107L; // 74!
172 factorials[75] = 2.480914081139539809195e+109L; // 75!
173 factorials[76] = 1.885494701666050254988e+111L; // 76!
174 factorials[77] = 1.451830920282858696341e+113L; // 77!
175 factorials[78] = 1.132428117820629783146e+115L; // 78!
176 factorials[79] = 8.946182130782975286851e+116L; // 79!
177 factorials[80] = 7.156945704626380229481e+118L; // 80!
178 factorials[81] = 5.797126020747367985880e+120L; // 81!
179 factorials[82] = 4.753643337012841748421e+122L; // 82!
180 factorials[83] = 3.945523969720658651190e+124L; // 83!
181 factorials[84] = 3.314240134565353266999e+126L; // 84!
182 factorials[85] = 2.817104114380550276949e+128L; // 85!
183 factorials[86] = 2.422709538367273238177e+130L; // 86!
184 factorials[87] = 2.107757298379527717214e+132L; // 87!
185 factorials[88] = 1.854826422573984391148e+134L; // 88!
186 factorials[89] = 1.650795516090846108122e+136L; // 89!
187 factorials[90] = 1.485715964481761497310e+138L; // 90!
188 factorials[91] = 1.352001527678402962552e+140L; // 91!
189 factorials[92] = 1.243841405464130725548e+142L; // 92!
190 factorials[93] = 1.156772507081641574759e+144L; // 93!
191 factorials[94] = 1.087366156656743080274e+146L; // 94!
192 factorials[95] = 1.032997848823905926260e+148L; // 95!
193 factorials[96] = 9.916779348709496892096e+149L; // 96!
194 factorials[97] = 9.619275968248211985333e+151L; // 97!
195 factorials[98] = 9.426890448883247745626e+153L; // 98!
196 factorials[99] = 9.332621544394415268170e+155L; // 99!
197 factorials[100] = 9.332621544394415268170e+157L; // 100!
198 factorials[101] = 9.425947759838359420852e+159L; // 101!
199 factorials[102] = 9.614466715035126609269e+161L; // 102!
200 factorials[103] = 9.902900716486180407547e+163L; // 103!
201 factorials[104] = 1.029901674514562762385e+166L; // 104!
202 factorials[105] = 1.081396758240290900504e+168L; // 105!
203 factorials[106] = 1.146280563734708354534e+170L; // 106!
204 factorials[107] = 1.226520203196137939352e+172L; // 107!
205 factorials[108] = 1.324641819451828974500e+174L; // 108!
206 factorials[109] = 1.443859583202493582205e+176L; // 109!
207 factorials[110] = 1.588245541522742940425e+178L; // 110!
208 factorials[111] = 1.762952551090244663872e+180L; // 111!
209 factorials[112] = 1.974506857221074023537e+182L; // 112!
210 factorials[113] = 2.231192748659813646597e+184L; // 113!
211 factorials[114] = 2.543559733472187557120e+186L; // 114!
212 factorials[115] = 2.925093693493015690688e+188L; // 115!
213 factorials[116] = 3.393108684451898201198e+190L; // 116!
214 factorials[117] = 3.969937160808720895402e+192L; // 117!
215 factorials[118] = 4.684525849754290656574e+194L; // 118!
216 factorials[119] = 5.574585761207605881323e+196L; // 119!
217 factorials[120] = 6.689502913449127057588e+198L; // 120!
218 factorials[121] = 8.094298525273443739682e+200L; // 121!
219 factorials[122] = 9.875044200833601362412e+202L; // 122!
220 factorials[123] = 1.214630436702532967577e+205L; // 123!
221 factorials[124] = 1.506141741511140879795e+207L; // 124!
222 factorials[125] = 1.882677176888926099744e+209L; // 125!
223 factorials[126] = 2.372173242880046885677e+211L; // 126!
224 factorials[127] = 3.012660018457659544810e+213L; // 127!
225 factorials[128] = 3.856204823625804217357e+215L; // 128!
226 factorials[129] = 4.974504222477287440390e+217L; // 129!
227 factorials[130] = 6.466855489220473672507e+219L; // 130!
228 factorials[131] = 8.471580690878820510985e+221L; // 131!
229 factorials[132] = 1.118248651196004307450e+224L; // 132!
230 factorials[133] = 1.487270706090685728908e+226L; // 133!
231 factorials[134] = 1.992942746161518876737e+228L; // 134!
232 factorials[135] = 2.690472707318050483595e+230L; // 135!
233 factorials[136] = 3.659042881952548657690e+232L; // 136!
234 factorials[137] = 5.012888748274991661035e+234L; // 137!
235 factorials[138] = 6.917786472619488492228e+236L; // 138!
236 factorials[139] = 9.615723196941089004197e+238L; // 139!
237 factorials[140] = 1.346201247571752460588e+241L; // 140!
238 factorials[141] = 1.898143759076170969429e+243L; // 141!
239 factorials[142] = 2.695364137888162776589e+245L; // 142!
240 factorials[143] = 3.854370717180072770522e+247L; // 143!
241 factorials[144] = 5.550293832739304789551e+249L; // 144!
242 factorials[145] = 8.047926057471991944849e+251L; // 145!
243 factorials[146] = 1.174997204390910823948e+254L; // 146!
244 factorials[147] = 1.727245890454638911203e+256L; // 147!
245 factorials[148] = 2.556323917872865588581e+258L; // 148!
246 factorials[149] = 3.808922637630569726986e+260L; // 149!
247 factorials[150] = 5.713383956445854590479e+262L; // 150!
248 factorials[151] = 8.627209774233240431623e+264L; // 151!
249 factorials[152] = 1.311335885683452545607e+267L; // 152!
250 factorials[153] = 2.006343905095682394778e+269L; // 153!
251 factorials[154] = 3.089769613847350887959e+271L; // 154!
252 factorials[155] = 4.789142901463393876336e+273L; // 155!
253 factorials[156] = 7.471062926282894447084e+275L; // 156!
254 factorials[157] = 1.172956879426414428192e+278L; // 157!
255 factorials[158] = 1.853271869493734796544e+280L; // 158!
256 factorials[159] = 2.946702272495038326504e+282L; // 159!
257 factorials[160] = 4.714723635992061322407e+284L; // 160!
258 factorials[161] = 7.590705053947218729075e+286L; // 161!
259 factorials[162] = 1.229694218739449434110e+289L; // 162!
260 factorials[163] = 2.004401576545302577600e+291L; // 163!
261 factorials[164] = 3.287218585534296227263e+293L; // 164!
262 factorials[165] = 5.423910666131588774984e+295L; // 165!
263 factorials[166] = 9.003691705778437366474e+297L; // 166!
264 factorials[167] = 1.503616514864999040201e+300L; // 167!
265 factorials[168] = 2.526075744973198387538e+302L; // 168!
266 factorials[169] = 4.269068009004705274939e+304L; // 169!
267 factorials[170] = 7.257415615307998967397e+306L; // 170!
268
269 N = sizeof(factorials) / sizeof(long double);
270}
271
272
273/*functions copied from: gamma_function.c*/
274
275
276
278// double Gamma_Function( double x ) //
279// //
280// Description: //
281// This function uses Lanczos' expression to calculate Gamma(x) for real //
282// x, where -(max_double_arg - 1) < x < max_double_arg. //
283// Note the Gamma function is meromorphic in the complex plane and has //
284// poles at the nonpositive integers. //
285// Tests for x a positive integer or a half positive integer give a //
286// maximum absolute relative error of about 1.9e-16. //
287// //
288// If x > max_double_arg, then one should either use xGamma_Function(x) //
289// or calculate lnGamma(x). //
290// Note that for x < 0, ln (Gamma(x)) may be a complex number. //
291// //
292// Arguments: //
293// double x Argument of the Gamma function. //
294// //
295// Return Values: //
296// If x is positive and is less than max_double_arg then Gamma(x) is //
297// returned and if x > max_double_arg then DBL_MAX is returned. If x is //
298// a nonpositive integer i.e. x is a pole, then DBL_MAX is returned //
299// ( note that Gamma(x) changes sign on each side of the pole). If x is //
300// nonpositive nonintegral, then if Gamma(x) > DBL_MAX, then DBL_MAX is //
301// returned and if Gamma(x) < -DBL_MAX, then -DBL_MAX is returned. //
302// //
303// Example: //
304// double x, g; //
305// //
306// g = Gamma_Function( x ); //
308template <typename TData, int dim, typename TRet>
310{
311 long double g;
312
313 if ( x > max_double_arg )
314 return DBL_MAX;
315 g = xGamma_Function( (long double) x);
316 if (fabsl(g) < DBL_MAX)
317 return (double) g;
318 return (g < 0.0L) ? -DBL_MAX : DBL_MAX;
319
320}
321
322
324// long double xGamma_Function( long double x ) //
325// //
326// Description: //
327// This function uses Lanczos' expression to calculate Gamma(x) for real //
328// x, where -(max_long_double_arg - 1) < x < max_long_double_arg. //
329// Note the Gamma function is meromorphic in the complex plane and has //
330// poles at the nonpositive integers. //
331// Tests for x a positive integer or a half positive integer give a //
332// maximum absolute relative error of about 3.5e-16. //
333// //
334// If x > max_long_double_arg, then one should use lnGamma(x). //
335// Note that for x < 0, ln (Gamma(x)) may be a complex number. //
336// //
337// Arguments: //
338// long double x Argument of the Gamma function. //
339// //
340// Return Values: //
341// If x is positive and is less than max_long_double_arg, then Gamma(x) //
342// is returned and if x > max_long_double_arg, then LDBL_MAX is returned. //
343// If x is a nonpositive integer i.e. x is a pole, then LDBL_MAX is //
344// returned ( note that Gamma(x) changes sign on each side of the pole). //
345// If x is nonpositive nonintegral, then if x > -(max_long_double_arg + 1)//
346// then Gamma(x) is returned otherwise 0.0 is returned. //
347// //
348// Example: //
349// long double x, g; //
350// //
351// g = xGamma_Function( x ); //
353template <typename TData, int dim, typename TRet>
355{
356 long double sin_x;
357 long double rg;
358 long int ix;
359
360 // For a positive argument (x > 0) //
361 // if x <= max_long_double_arg return Gamma(x) //
362 // otherwise return LDBL_MAX. //
363
364 if (x > 0.0L)
365 {
366 if (x <= max_long_double_arg)
367 return xGamma(x);
368 else
369 return LDBL_MAX;
370 }
371
372 // For a nonpositive argument (x <= 0) //
373 // if x is a pole return LDBL_MAX //
374
375 if ( x > -(long double)LONG_MAX)
376 {
377 ix = (long int) x;
378 if ( x == (long double) ix)
379 return LDBL_MAX;
380 }
381 sin_x = sinl(M_PI * x);
382 if ( sin_x == 0.0L )
383 return LDBL_MAX;
384
385 // if x is not a pole and x < -(max_long_double_arg - 1) //
386 // then return 0.0L //
387
388 if ( x < -(max_long_double_arg - 1.0L) )
389 return 0.0L;
390
391 // if x is not a pole and x >= -(max_long_double - 1) //
392 // then return Gamma(x) //
393
394 rg = xGamma(1.0L - x) * sin_x / M_PI;
395 if ( rg != 0.0L )
396 return (1.0L / rg);
397 return LDBL_MAX;
398}
399
400
402// static long double xGamma( long double x ) //
403// //
404// Description: //
405// This function uses Lanczos' expression to calculate Gamma(x) for real //
406// x, where 0 < x <= 900. For 900 < x < 1755.5, the duplication formula //
407// is used. //
408// The major source of relative error is in the use of the c library //
409// function powl(). The results have a relative error of about 10^-16. //
410// except near x = 0. //
411// //
412// If x > 1755.5, then one should calculate lnGamma(x). //
413// //
414// Arguments: //
415// long double x Argument of the Gamma function. //
416// //
417// Return Values: //
418// If x is positive and is less than 1755.5 then Gamma(x) is returned and //
419// if x > 1755.5 then LDBL_MAX is returned. //
420// //
421// Example: //
422// long double x; //
423// long double g; //
424// //
425// g = xGamma_Function( x ); //
427template <typename TData, int dim, typename TRet>
428long double FractalField<TData,dim,TRet>::xGamma(long double x)
429{
430
431 long double xx = (x < 1.0L) ? x + 1.0L : x;
432 long double temp;
433 int const n = sizeof(a) / sizeof(long double);
434 int i;
435
436 if (x > 1755.5L) return LDBL_MAX;
437
438 if (x > 900.0L) return Duplication_Formula(x);
439
440 temp = 0.0L;
441 for (i = n-1; i >= 0; i--) {
442 temp += ( a[i] / (xx + (long double) i) );
443 }
444 temp += 1.0L;
445 temp *= ( powl((g + xx - 0.5L) / e, xx - 0.5L) / exp_g_o_sqrt_2pi );
446 return (x < 1.0L) ? temp / x : temp;
447}
448
449
451// static long double Duplication_Formula(long double two_x) //
452// //
453// Description: //
454// This function returns the Gamma(two_x) using the duplication formula //
455// Gamma(2x) = (2^(2x-1) / sqrt(pi)) Gamma(x) Gamma(x+1/2). //
456// //
457// Arguments: //
458// none //
459// //
460// Return Values: //
461// Gamma(two_x) //
462// //
463// Example: //
464// long double two_x, g; //
465// //
466// g = Duplication_Formula(two_x); //
468template <typename TData, int dim, typename TRet>
470{
471 long double x = 0.5L * two_x;
472 long double g;
473 //double two_n = 1.0;
474 int n = (int) two_x - 1;
475
476 g = powl(2.0L, two_x - 1.0L - (long double) n);
477 g = ldexpl(g,n);
478 g /= sqrt(M_PI);
479 g *= xGamma_Function(x);
480 g *= xGamma_Function(x + 0.5L);
481
482 return g;
483}
484
485
487// double Gamma_Function_Max_Arg( void ) //
488// //
489// Description: //
490// This function returns the maximum argument of Gamma_Function for which //
491// a number < DBL_MAX is returned, for arguments greater than 1. //
492// //
493// Arguments: //
494// none //
495// //
496// Return Values: //
497// max_double_arg (171.0). //
498// //
499// Example: //
500// double x; //
501// //
502// x = Gamma_Function_Max_Arg(); //
504template <typename TData, int dim, typename TRet>
506{
507 return max_double_arg;
508}
509
510
512// long double xGamma_Function_Max_Arg( void ) //
513// //
514// Description: //
515// This function returns the maximum argument of Gamma_Function for which //
516// a number < LDBL_MAX is returned, for arguments greater than 1. //
517// //
518// Arguments: //
519// none //
520// //
521// Return Values: //
522// max_long_double_arg (1755.5). //
523// //
524// Example: //
525// long double x; //
526// //
527// x = xGamma_Function_Max_Arg(); //
529template <typename TData, int dim, typename TRet>
531{
532 return max_long_double_arg;
533}
534
535
536}
537}
void InitConstants()
Definition gamma_function.h:64
long double xGamma(long double x)
Definition gamma_function.h:428
long double xGamma_Function_Max_Arg()
Definition gamma_function.h:530
long double Duplication_Formula(long double two_x)
Definition gamma_function.h:469
long double xGamma_Function(long double x)
Definition gamma_function.h:354
double Gamma_Function(double x)
Definition gamma_function.h:309
double Gamma_Function_Max_Arg()
Definition gamma_function.h:505