www.pudn.com > 1.rar > Mathutil.cpp
///////////////////////////////////////////////////////////////////////////////
// //
// Dateiname: MATHUTIL.CPP //
// //
// Autor: Andreas Jäger, Friedrich-Schiller-Universität Jena //
// //
// System: WIN_RWPM.EXE und WINLRWPM.EXE //
// //
// Beschreibung: Mathematische Routinen //
// //
// Hinweise: //
// //
///////////////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include
#include "mathutil.h"
long Fixl(long double x)
{
#pragma warning ( disable : 4244 )
return x;
#pragma warning ( default : 4244 )
}
long double Reml(long double x, long double y)
{
return x-y*Fixl(x/y);
}
long Fix(double x)
{
#pragma warning ( disable : 4244 )
return x;
#pragma warning ( default : 4244 )
}
double Rem(double x, double y)
{
return x-y*Fix(x/y);
}
bool Rational(long double x, long maxtest, long double TOL, long& zaehler, long& nenner)
{
#pragma warning ( disable : 4244 )
long y = x;
#pragma warning ( default : 4244 )
if (fabsl(y - x) < TOL)
{
zaehler = y;
nenner = 1L;
return true;
}
for (long p=2; p<=maxtest; p++)
{
#pragma warning ( disable : 4244 )
y = x*p;
#pragma warning ( default : 4244 )
if (fabsl(y - p*x) < TOL)
{
zaehler = y;
nenner = p;
return true;
}
}
return false;
}
long ggT(long x, long y)
{
for (long r = x%y; r!=0; x=y, y=r, r=x%y)
{
}
return abs(y);
}
long kgV(long x, long y)
{
long z = x;
z /= ggT(x,y);
return z*y;
}
long Fak(long x)
{
long s = 1;
for (int i=2;i <= x;i++) s *= i;
return s;
}
// Diese Version nur für gutmütige n, m verwenden
long BinKoeff(long n, long m)
{
long p = 1;
long i;
for (i=n; i>=n-m+1; i--)
{
p *= i;
}
for (i=2; i<=m; i++)
{
p /= i;
}
return p;
}
long double BinKoeffEx(long double ny, long k)
{
long double t;
long i;
t = 1;
for (i=1; i<=k; i++)
{
t *= (ny-i+1)/i;
}
return t;
}
#define Fib_Rho (0.5*(1+sqrtl(5)))
long double Fibonacci(int n)
{
return (powl(Fib_Rho, n) - powl(1-Fib_Rho, n)) / (2.0 * Fib_Rho - 1);
}
// Berechnet eine nichtkrumme Zahl im Intervall x1<=x<=x2
// z.B. 1 für x1=0.6, x2=3.5
long double NichtKrumm(long double x1, long double x2)
{
if (x1*x2<=0) return 0.0L;
long double dist = fabsl(x1 - x2);
if (dist <= 1E-300) return x1; // Sicherheit vor Überläufen
#pragma warning ( disable : 4244 )
int l = floorl(log10l(dist));
#pragma warning ( default : 4244 )
long double d1 = ((l<0) ? (1.0/powl(10.0, -l)) : powl(10.0, l)); // "Maßeinheit"
while (floorl(min(x1,x2)/(10.0*d1)) 0) and (Rem(k,2)=0) then t := t + 2*bk;
end;
if (k = 0) then t := t + bk;
{ Normalizing condition, J0(x) + 2*J2(x) + 2*J4(x) + ... = 1}
J := bn/t;
{Restore results for x = 0; J0(0) = 1, Jn(0) = 0 for n > 0}
J := (1-z)*J + z*Integer(n=0);
Besseln:=J;
end;
*/
long double FGammal(long double x)
{
long double MaxGammaNumber = 1700;
long double inf = 233;
long double sig,fact,y,y1,Res,z,xnum,xden,gam,ysq,sum;
int i,n;
long double eps=1e-20L;
//{Constants for Aproximation of Gamma(x)}
long double p [] = {0,-1.71618513886549492533811e+0, 2.47656508055759199108314e+1,
-3.79804256470945635097577e+2, 6.29331155312818442661052e+2,
8.66966202790413211295064e+2, -3.14512729688483675254357e+4,
-3.61444134186911729807069e+4, 6.64561438202405440627855e+4};
long double q[] = {0,-3.08402300119738975254353e+1, 3.15350626979604161529144e+2,
-1.01515636749021914166146e+3, -3.10777167157231109440444e+3,
2.25381184209801510330112e+4, 4.75584627752788110767815e+3,
-1.34659959864969306392456e+5, -1.15132259675553483497211e+5};
long double c[] = {0, -1.910444077728e-03 , 8.4171387781295e-04,
-5.952379913043012e-04 , 7.93650793500350248e-04,
-2.777777777777681622553e-03 , 8.333333333333333331554247e-02,
5.7083835261e-03};
sig = 1;
fact = 1;
y = x;
//{ Negative x }
if (y <= 0)
{
y = -x;
if (Reml(Fixl(y),2) != 0)
sig = -1;
Res = Reml(y,1);
if (Res != 0)
fact = -PI/sinl(PI*Res);
else
{
//{Pole at negative integers}
fact = inf;
}
y = y + 1;
}
if (y < eps)
gam=1/y; //{ Small x }
if ((y > eps) && (y <= 12))
{
//{ eps < x <= 12 }
y1 = y;
if (y < 1)
{
n = 0;
z = y;
y = y + 1;
}
else
{
n = Fixl(y) - 1;
y = y - n;
z = y - 1;
}
//{Rational approximation for 1.0 < x < 2.0}
xnum = 0;
xden = 1;
for (i = 1; i<=8; i++)
{
xnum = (xnum + p[i]) * z;
xden = xden * z + q[i];
}
gam = xnum / xden + 1;
if (y > y1)
{
//{ Adjust result for 0.0 < x < 1.0}
gam = gam / y1;
}
else
{
if (y < y1)
{
//{ Adjust result for 2.0 < x < 12.0}
for (i=1; i<=n; i++)
{
gam = gam * y;
y = y + 1;
}
}
}
}
//{ Different rational approximation for 12 < x < xbig}
if ((y > 12) && (y <= MaxGammaNumber))
{
ysq = y * y;
sum = c[7];
for (i=1; i<=6; i++)
{
sum = sum/ysq + c[i];
}
sum = sum/y - y + logl(2*PI)/2 + (y-0.5L)*logl(y);
gam = expl(sum); //{Large x}
}
//{Final adjustments }
if (y > MaxGammaNumber)
gam = inf;
gam = sig*gam;
if (fact != 1 )
gam = fact / gam;
if (x==(int)x)
gam = (int)gam;
return gam;
}
double FGamma(double x)
{
double MaxGammaNumber = 1700;
double inf = 233;
double sig,fact,y,y1,Res,z,xnum,xden,gam,ysq,sum;
int i,n;
double eps=1e-20L;
//{Constants for Aproximation of Gamma(x)}
double p [] = {0,-1.71618513886549492533811e+0, 2.47656508055759199108314e+1,
-3.79804256470945635097577e+2, 6.29331155312818442661052e+2,
8.66966202790413211295064e+2, -3.14512729688483675254357e+4,
-3.61444134186911729807069e+4, 6.64561438202405440627855e+4};
double q[] = {0,-3.08402300119738975254353e+1, 3.15350626979604161529144e+2,
-1.01515636749021914166146e+3, -3.10777167157231109440444e+3,
2.25381184209801510330112e+4, 4.75584627752788110767815e+3,
-1.34659959864969306392456e+5, -1.15132259675553483497211e+5};
double c[] = {0, -1.910444077728e-03 , 8.4171387781295e-04,
-5.952379913043012e-04 , 7.93650793500350248e-04,
-2.777777777777681622553e-03 , 8.333333333333333331554247e-02,
5.7083835261e-03};
sig = 1;
fact = 1;
y = x;
//{ Negative x }
if (y <= 0)
{
y = -x;
if (Rem(Fix(y),2) != 0)
sig = -1;
Res = Rem(y,1);
if (Res != 0)
fact = -PI_/sin(PI_*Res);
else
{
//{Pole at negative integers}
fact = inf;
}
y = y + 1;
}
if (y < eps)
gam=1/y; //{ Small x }
if ((y > eps) && (y <= 12))
{
//{ eps < x <= 12 }
y1 = y;
if (y < 1)
{
n = 0;
z = y;
y = y + 1;
}
else
{
n = Fix(y) - 1;
y = y - n;
z = y - 1;
}
//{Rational approximation for 1.0 < x < 2.0}
xnum = 0;
xden = 1;
for (i = 1; i<=8; i++)
{
xnum = (xnum + p[i]) * z;
xden = xden * z + q[i];
}
gam = xnum / xden + 1;
if (y > y1)
{
//{ Adjust result for 0.0 < x < 1.0}
gam = gam / y1;
}
else
{
if (y < y1)
{
//{ Adjust result for 2.0 < x < 12.0}
for (i=1; i<=n; i++)
{
gam = gam * y;
y = y + 1;
}
}
}
}
//{ Different rational approximation for 12 < x < xbig}
if ((y > 12) && (y <= MaxGammaNumber))
{
ysq = y * y;
sum = c[7];
for (i=1; i<=6; i++)
{
sum = sum/ysq + c[i];
}
sum = sum/y - y + log(2*PI)/2 + (y-0.5)*log(y);
gam = exp(sum); //{Large x}
}
//{Final adjustments }
if (y > MaxGammaNumber)
gam = inf;
gam = sig*gam;
if (fact != 1 )
gam = fact / gam;
if (x==(int)x)
gam = (int)gam;
return gam;
}
/*
long double Bessela(long double nu, long double x);
{
long double xz,KMIN,z,t,q1,p1,r,mu,l,k,s,y,J,Ja,da,Jp,dp,digits;
int oka,okp;
xz = (x==0) ? 1 : 0;
x = x + xz;
//{ Asymptotic series }
KMIN = 12;//{ Try the series for at least KMIN/2 terms.}
z = (-1)/Powl(8*x,2.0L);
t = 0*x;
q1 = t;
t = t + 1;
p1 = t;
r = abs(t);
mu = 2*nu;
k = 0;
l = -1;
oka = Integer(x > 0.9*nu);
while ((FAbs(t) > eps*Max(FAbs(p1),FAbs(q1))) and Boolean(oka)) do begin
s := t;
k := k+1;
l := l+2;
t := (mu-l)*(mu+l)/k*t;
q1 := q1 + t;
k := k+1;
l := l+2;
t := (mu-l)*(mu+l)/k*t*z;
p1 := p1+t;
r := Max(r,FAbs(t));
if k >= 2*KMIN then oka := oka * Integer(FAbs(t) <= 5*FAbs(s));
end;
q1 := q1/(8*x);
y := x - (2*nu+1)*PI/4;
Ja := FSqrt((2)/(PI*x)) * (p1*FCos(y) - q1*FSin(y));
da := Max(oka*FLg(Max(FAbs(p1),FAbs(q1))/(eps*FAbs(r))),0);
oka := oka*Integer(da > 0);
{ If necessary, use power series expansion.
Requires the gamma function from another M file}
okp := Integer(da < 12);
if Boolean(okp) then begin
t := okp * FPower(x/2,nu) / Gamma(nu+1);
s := t;
r := abs(t);
z := -Sqr(x/2);
k := 0;
while ((abs(t) > eps*FAbs(s)) and Boolean(okp)) do begin
k := k + 1;
t := z*t/(k*(nu+k));
s := s + t;
r := Max(r,abs(t));
okp := okp * Integer((eps*FAbs(r) <= FAbs(s)));
end;
Jp := s;
dp := Max(okp*FLg((1-okp+FAbs(s))/(1-okp+eps*FAbs(r))),0);
okp := okp*Integer(dp > 0)*Integer(dp>da);
oka := oka*Integer(da>=dp);
digits := oka*da + okp*dp;
J := oka*Ja + okp*Jp + (0)/digits;end
else begin
J := Ja;
digits := da;
end;
{ Restore results for x = 0; J0(0) = 1, Jnu(0) = 0 for nu > 0.}
Bessela := (1-xz)*J + xz*Integer(nu=0);
end;
long double Bessel(long double n, long double x);
{
int nn = n;
if (nn == n)
return Besseln(n, x);
else
return Bessela(n, x);
}
*/
long double ArSinhl(long double x)
{
if (x < 0) return -logl(-x+sqrtl(sqrl(-x)+1));
return logl(x+sqrtl(sqrl(x)+1));
}
long double ArCoshl(long double x)
{
return logl(x+sqrtl(sqrl(x)-1));
}
long double ArTanhl(long double x)
{
return 0.5L*logl((1+x)/(1-x));
}
long double ArCothl(long double x)
{
return 0.5L*logl((x+1)/(x-1));
}
double ArSinh(double x)
{
if (x < 0) return -log(-x+sqrt(sqr(-x)+1));
return log(x+sqrt(sqr(x)+1));
}
double ArCosh(double x)
{
return log(x+sqrt(sqr(x)-1));
}
double ArTanh(double x)
{
return 0.5*logl((1+x)/(1-x));
}
double ArCoth(double x)
{
return 0.5*log((x+1)/(x-1));
}