Plugins
Loading...
Searching...
No Matches
ln_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
34#include "lib_disc/common/groups_util.h"
35#include "lib_disc/common/function_group.h"
36#include "lib_grid/tools/subset_group.h"
37#include "bridge/bridge.h"
38#include "lib_disc/domain.h"
39#include "lib_disc/domain_traits.h"
40#include "lib_algebra/cpu_algebra_types.h"
41#include "lib_disc/quadrature/quadrature.h"
42#include "lib_disc/local_finite_element/local_finite_element_provider.h"
43#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
44#include "lib_disc/spatial_disc/user_data/user_data.h"
45#include "common/util/provider.h"
46#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
47#include "lib_disc/spatial_disc/disc_util/hfv1_geom.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
56#include "bindings/lua/lua_user_data.h"
57#endif
58
59namespace ug{
60namespace d3f{
61
62
63/*functions copied from: ln_gamma_function.c*/
64
66// double Ln_Gamma_Function( double x ) //
67// //
68// Description: //
69// This function calculates the natural log of Gamma(x) for positive real //
70// x. //
71// Assuming that Gamma_Function_Max_Arg() = 171, then //
72// If 0 < x <= 171, then ln(gamma(x)) is calculated by taking the natural //
73// log of the result from Gamma_Function(x). If x > 171, then //
74// ln(gamma(x)) is calculated using the asymptotic expansion //
75// ln(gamma(x)) ~ ln(2sqrt(2pi)) - x + (x - 1/2) ln x + //
76// Sum B[2j] / [ 2j * (2j-1) * x^(2j-1) ], summed over //
77// j from 1 to 3 and where B[2j] is the 2j-th Bernoulli number. //
78// //
79// Arguments: //
80// double x Argument of the ln Gamma function. The argument x must be //
81// positive. //
82// //
83// Return Values: //
84// ln(Gamma(x)) where x > 0. //
85// //
86// Example: //
87// double x, g; //
88// //
89// g = Ln_Gamma_Function( x ); //
91
92template <typename TData, int dim, typename TRet>
93double FractalField<TData,dim,TRet>::Ln_Gamma_Function(double x)
94{
95 // For a positive argument, 0 < x <= Gamma_Function_Max_Arg() //
96 // then return log Gamma(x). //
97 if (x <= Gamma_Function_Max_Arg())
98 return log(Gamma_Function(x));
99 // otherwise return result from asymptotic expansion of ln Gamma(x). //
100 return (double) xLnGamma_Asymptotic_Expansion( (long double) x );
101}
102
103
105// long double xLn_Gamma_Function( long double x ) //
106// //
107// Description: //
108// This function calculates the natural log of Gamma(x) for positive real //
109// x. //
110// Assuming that Gamma_Function_Max_Arg() = 171, then //
111// If 0 < x <= 171, then ln(gamma(x)) is calculated by taking the natural //
112// log of the result from Gamma_Function(x). If x > 171, then //
113// ln(gamma(x)) is calculated using the asymptotic expansion //
114// ln(gamma(x)) ~ ln(2sqrt(2pi)) - x + (x - 1/2) ln x + //
115// Sum B[2j] / [ 2j * (2j-1) * x^(2j-1) ], summed over //
116// j from 1 to 3 and where B[2j] is the 2j-th Bernoulli number. //
117// //
118// Arguments: //
119// long double x Argument of the ln Gamma function. The argument x must //
120// be positive. //
121// //
122// Return Values: //
123// ln(Gamma(x)) where x > 0. //
124// //
125// Example: //
126// double x; //
127// long double g; //
128// //
129// g = xLn_Gamma_Function( x ); //
131template <typename TData, int dim, typename TRet>
132long double FractalField<TData,dim,TRet>::xLn_Gamma_Function(long double x)
133{
134 // For a positive argument, 0 < x <= Gamma_Function_Max_Arg() //
135 // then return log Gamma(x). //
136 if (x <= Gamma_Function_Max_Arg())
137 return logl(xGamma_Function(x));
138 // otherwise return result from asymptotic expansion of ln Gamma(x). //
139 return xLnGamma_Asymptotic_Expansion( x );
140}
141
142
144// static long double xLnGamma_Asymptotic_Expansion( long double x ) //
145// //
146// Description: //
147// This function estimates log(gamma(x)) by evaluating the asymptotic //
148// expression: //
149// ln(Gamma(x)) ~ ln(2sqrt(2pi)) - x + (x - 1/2) ln x + //
150// Sum B[2j] / [ 2j * (2j-1) * x^(2j-1) ], summed over //
151// j from 1 to 3 and where B[2j] is the 2j-th Bernoulli number. //
152// //
153// Arguments: //
154// long double x Argument of the ln Gamma function. The argument x must //
155// be positive. //
156// //
157// Return Values: //
158// ln(Gamma(x)) where x > Gamma_Function_Max_Arg() //
159// //
160// Example: //
161// double x; //
162// long double g; //
163// //
164// g = xlnGamma_Asymptotic_Expansion( x ); //
166
167
168
169template <typename TData, int dim, typename TRet>
170long double FractalField<TData,dim,TRet>::xLnGamma_Asymptotic_Expansion( long double x ) {
171 const int m = 3;
172 long double term[3];
173 long double sum = 0.0L;
174 long double xx = x * x;
175 long double xj = x;
176 long double lngamma;
177 int i;
178
179 lngamma = log_sqrt_2pi - xj + (xj - 0.5L) * logl(xj);
180 for (i = 0; i < m; i++)
181 {
182 term[i] = B[i] / xj; xj *= xx;
183 }
184 for (i = m - 1; i >= 0; i--)
185 sum += term[i];
186 return lngamma + sum;
187}
188
189
190
191}
192}
ug::d3f::FractalField::Ln_Gamma_Function
double Ln_Gamma_Function(double x)
Definition ln_gamma_function.h:93
ug::d3f::FractalField::xLn_Gamma_Function
long double xLn_Gamma_Function(long double x)
Definition ln_gamma_function.h:132
ug::d3f::FractalField::xLnGamma_Asymptotic_Expansion
long double xLnGamma_Asymptotic_Expansion(long double x)
Definition ln_gamma_function.h:170
cpu_algebra_types.h
domain_traits.h
function_group.h
fv1_geom.h
groups_util.h
hfv1_geom.h
user_data.h
local_finite_element_provider.h
lua_user_data.h
subset_group.h