Unit Mathe; //02.08.2011 by Heizkoerper

Interface

Const ErgebnisFehler=1E2200;

Function  SgnPas(x:Extended):Extended;
//Signum (Vorzeichen)
Function  FloorPas(x:Extended):Extended;
//abrunden	
Function  CeilPas(x:Extended):Extended;
//aufrunden
Function  ModPas(x,y:Extended):Extended;
//Modulos
Function  OddPas(x:Extended):Extended;
//ungerade=1,sonst 0 
Function  CubPas(x:Extended):Extended;
//x^3
Function  CurPas(x:Extended):Extended;
//dritte Wurzel
Function  LnTwoPas(x:Extended):Extended;
//Logarithmus zur Basis 2
Function  LogPas(x:Extended):Extended;
//Logarithmus zur Basis 10
Function  LnBasisPas(Zahl,Basis:Extended):Extended;
//Logarithmus zur Basis Basic         
Function  TwoPas(x:Extended):Extended;
//2^x
Function  TenPas(x:Extended):Extended;
//10^x
Function  RundenPas(x,Nk:Extended):Extended;
//runden mit Nk Nachkommastellen
Function  PrimPas(x:Extended):Extended;
//Primzahl=1, sonst 0
Function  AnzPrimPas(x:Extended):Extended;
//Anzahl der Primzahlen bis x
Function  ntePrimPas(x:Extended):Extended;
//die x-te Primzahl
Function  FakPas(x:Extended):Extended;
//Fakultaet mit Ganzzahl
Function  PotPas(x,y:Extended):Extended;
//x^y
Function  BernoulliPas(x:Extended):Extended;
//x.te Bernoullische-Zahl
Function  EulerPas(x:Extended):Extended;
//x.te Eulersche-Zahl
Function  SiPas(x:Extended):Extended;
//Integral-Sinus
Function  CiPas(x:Extended):Extended;
//Integral-Cosinus
Function  TanPas(x:Extended):Extended;
//Tangens
Function  TaniPas(x:Extended):Extended;
//Integral-Tangens
Function  CotPas(x:Extended):Extended;
//Cotangens
Function  CotiPas(x:Extended):Extended;
//Integral-Cotangens
Function  FSPas(x:Extended):Extended;
//Fresnel-Sin-Integral
Function  FCPas(x:Extended):Extended;
//Fresnel-Cos-Integral
Function  ASinPas(x:Extended):Extended;
//ArkusSinus
Function  ACosPas(x:Extended):Extended;
//ArcusCosinus
Function  ACotPas(x:Extended):Extended;
//ArcusCotangens 
Function  SinhPas(x:Extended):Extended;
//SinusHyperbelicus
Function  SihPas(x:Extended):Extended;
//Integral-SinusHyperbelicus 
Function  CoshPas(x:Extended):Extended;
//CosinusHyperbelicus
Function  CihPas(x:Extended):Extended;
//Integral-CosinusHyperbelicus
Function  TanhPas(x:Extended):Extended;
//TangensHyperbelicus
Function  TanhiPas(x:Extended):Extended;
//Integral-TangensHyperbelicus
Function  CothPas(x:Extended):Extended;
//CotangensHyperbelicus
Function  CothiPas(x:Extended):Extended;
//Integral-CotangensHyperbelicus 
Function  ASinhPas(x:Extended):Extended;
//AreaSinusHyperbelicus
Function  ACoshPas(x:Extended):Extended;
//AreaCosinusHyperbelicus
Function  ATanhPas(x:Extended):Extended;
//AreaTangensHyperbelicus
Function  ACothPas(x:Extended):Extended;
//AreaCotangensHyperbelicus 
Function  AnnuitaetPas(Zinsen,Laufzeit:Extended):Extended;
//Annuitaet in Prozent 
Function  LaufzeitPas(Zinsen,Annuitaet:Extended):Extended;
//Laufzeit in Jahren
Function  ZinsenPas(Laufzeit,Annuitaet:Extended):Extended;
//Zinsen in Prozent
Function  GammaPas(x:Extended):Extended;
//Fakultaet für reale Zahlen (Gamma(x)=Fak(x-1)
Function  VariationenPas(x,y:Extended):Extended;
//Variationen von y Elementen zur x-ten Klasse ohne Wiederholung 
Function  KombinationenPas(x,y:Extended):Extended;
//Kombinationen von y Elementen zur x-ten Klasse ohne Wiederholung 
Function  ggTPas(x,y:Extended):Extended;
//groesster gemeinsamer Teiler
Function  kgVPas(x,y:Extended):Extended;
//kleinstes gemeinsames Vielfaches 
Function  MandelPas(Bildx,Bildy:Extended;Tiefe,Schranke:Word):Extended;
//Farbtiefe der Mandelbrotmenge 
Function  EiPas(x:Extended):Extended;
//ExponentialIntegral
Function  LiPas(x:Extended):Extended;
//IntegralLogarithmus
Function  PhiPas(x:Extended):Extended;
//Gaussches Fehlerintegral
Function  LegendrePas(x:Extended;n,m:Byte):Extended;
//n,m-tes Legendre-Polynom
Function  HermitePas(x:Extended;n,m:Byte):Extended;
//n,m-tes Hermitesches Polynom 
Function  BesselPas(x:Extended;n:Byte):Extended;
//n-tes Bessel-Polynom
Function  TschebyschewPas(x:Extended;n:Byte):Extended;
//n-tes Tschebyschew-Polynom
Function  FibPas(x:Extended):Extended;
//Fibonacci-Zahlen
Function  BinomialPas(n,k:Extended):Extended;
//Binomialkoeffizienten x ueber y
Function  HofPas(k:Extended):Extended;
//Hofstadter-Funktion
Function  ZetaPas(x:Extended):Extended;
//Zeta-Funktion (bis 4 nur angenaehert)
Function  GammaBerechnenPas:Extended;
//Gammazahl
Function  GoldbachPas(x:Extended):Extended;
//Goldbach-Funktion

Implementation

Const EulerKonstante=0.577215664901532860606;

Function SgnPas(x:Extended):Extended;
Begin
 If x>0 Then SgnPas:=1 Else SgnPas:=-1;
 If x=0 Then SgnPas:=0
End;

Function FloorPas(x:Extended):Extended;
Begin
 If Frac(x)<0 Then FloorPas:=Int(x)-1 Else FloorPas:=Int(x)
End;

Function CeilPas(x:Extended):Extended;
Begin
 If Frac(x)>0 Then CeilPas:=Int(x)+1 Else CeilPas:=Int(x)
End;

Function ModPas(x,y:Extended):Extended;
 Begin
  if y<>0 Then ModPas:=x-Int(x/y)*y Else ModPas:=ErgebnisFehler
 End;

Function OddPas(x:Extended):Extended;
Begin OddPas:=ModPas(x,2) End;

Function CubPas(x:Extended):Extended;
Begin CubPas:=x*Sqr(x) End;

Function CurPas(x:Extended):Extended;
Begin
 If x<>0 Then CurPas:=SgnPas(x)*Exp(1/3*Ln(Abs(x))) Else CurPas:=0
End;

Function LnTwoPas(x:Extended):Extended;
Begin If x>0 Then LnTwoPas:=Ln(x)/Ln(2) Else LnTwoPas:=ErgebnisFehler End;

Function LogPas(x:Extended):Extended;
Begin If x>0 Then LogPas:=Ln(x)/Ln(10) Else LogPas:=ErgebnisFehler End;

Function LnBasisPas(Zahl,Basis:Extended):Extended;
Begin
 If (Zahl>0) And (Basis>1) Then LnBasisPas:=Ln(Zahl)/Ln(Basis)
                           Else LnBasisPas:=ErgebnisFehler
End;

Function TwoPas(x:Extended):Extended;
Begin TwoPas:=Exp(x*Ln(2)) End;

Function TenPas(x:Extended):Extended;
Begin TenPas:=Exp(x*Ln(10)) End;

Function RundenPas(x,Nk:Extended):Extended;
Begin
 Nk:=Round(Nk);
 If (Nk>=0) And (Nk<15) Then RundenPas:=Round(x*TenPas(Nk))/TenPas(Nk)
                        Else RundenPas:=ErgebnisFehler
 End;

Function PrimPas(x:Extended):Extended;
Var Nenner,Wurzel:Extended;
Begin
 PrimPas:=0;x:=Round(x);
 If (x<2) Or (x>1e16) Then Exit;
 If (x=2) Or (x=3) Then Begin PrimPas:=1;Exit End;
 If OddPas(x)=0 Then Exit;
 Nenner:=3;Wurzel:=Sqrt(x);
 Repeat
  If ModPas(x,Nenner)=0 Then Exit;
  Nenner:=Nenner+2
 Until Nenner>Wurzel;
 PrimPas:=1
End;

Function AnzPrimPas(x:Extended):Extended;
Var xi :Extended;
    Anz:LongInt;
Begin
 If Abs(x)>1E5 Then Begin AnzPrimPas:=ErgebnisFehler;Exit End;
 Anz:=0;
 If x>=2 Then
  Begin
   xi:=2;
   While xi<=x Do
    Begin
     If PrimPas(xi)=1 Then Inc(Anz);
     xi:=xi+1
    End
  End;
 AnzPrimPas:=Anz
End;

Function ntePrimPas(x:Extended):Extended;
Var xi,Anz:LongInt;
Begin
 If (Abs(x)>1E5) Or (x<1) Then Begin ntePrimPas:=ErgebnisFehler;Exit End;
 Anz:=0;x:=Int(x);xi:=1;
 While Anz<x Do
  Begin
   Inc(xi);
   If PrimPas(xi)=1 Then Inc(Anz)
  End;
 ntePrimPas:=xi
End;

Function FakPas(x:Extended):Extended;
Var i,Ergebnis,vz:Extended;
Begin
 If x=0 Then FakPas:=1
        Else
  If Abs(x)>=1000 Then FakPas:=ErgebnisFehler
                  Else
   Begin
    Ergebnis:=1;i:=1;vz:=SgnPas(x);x:=Int(Abs(x));
    If ModPas(x,2)=0 Then vz:=1;
    While i<=x Do Begin Ergebnis:=Ergebnis*i;i:=i+1 End;
    FakPas:=Ergebnis*vz
   End
End;

Function PotPas(x,y:Extended):Extended;
Var vz:Extended;
Begin
 If y=0 Then Begin PotPas:=1;Exit End;
 If x=0 Then
  Begin
   If y>0 Then PotPas:=0 Else PotPas:=ErgebnisFehler;
   Exit
  End;
 If y*Ln(Abs(x))>5000 Then Begin PotPas:=ErgebnisFehler;Exit End;
 vz:=SgnPas(x);
 If x<0 Then
  If Frac(y)=0 Then
   If OddPas(Abs(y))=0 Then vz:=1 Else
        Else
   If (OddPas(1/Abs(y))=0) Or (Frac(1/Abs(y))<>0) Then
    Begin PotPas:=ErgebnisFehler;Exit End;
 PotPas:=vz*Exp(y*Ln(Abs(x)))
End;

Function BernoulliPas(x:Extended):Extended;
Const Epsilon=1E-18;
Var   Summe,Delta,Teiler:Extended;
      Zahl              :LongInt;
      i,k               :Word;
Begin
 If (x<1) Or (x>500) Then Begin BernoulliPas:=ErgebnisFehler;Exit End;
 k:=Round(x);
 Case k Of
  1:Begin BernoulliPas:=1/6;Exit End;
  2:Begin BernoulliPas:=1/30;Exit End;
  3:Begin BernoulliPas:=1/42;Exit End;
  4:Begin BernoulliPas:=1/30;Exit End;
  5:Begin BernoulliPas:=5/66;Exit End
    Else
     Begin
      Zahl:=2;Summe:=1;
      Repeat
       Teiler:=Zahl;
       For i:=2 To 2*k Do Teiler:=Teiler*Zahl;
       Delta:=1/Teiler;Summe:=Summe+Delta;Inc(Zahl)
      Until Delta<Epsilon
     End
 End; 
 BernoulliPas:=Summe/PotPas(Pi,2*k)/PotPas(2,2*k-1)*FakPas(2*k)
End;

Function EulerPas(x:Extended):Extended;
Const Epsilon=1E-18;
Var   Summe,Delta,Teiler:Extended;
      Zahl              :LongInt;
      i,k               :Word;
      Vz                :ShortInt;
Begin
 If (x<1) Or (x>300) Then Begin EulerPas:=ErgebnisFehler;Exit End;
 k:=Round(x);
 Case k Of
  1:Begin EulerPas:=1;Exit End;
  2:Begin EulerPas:=5;Exit End;
  3:Begin EulerPas:=61;Exit End;
  4:Begin EulerPas:=1385;Exit End;
  5:Begin EulerPas:=50521;Exit End
    Else
     Begin
      Zahl:=3;Summe:=1;Vz:=-1;
      Repeat
       Teiler:=Zahl;
       For i:=2 To 2*k+1 Do Teiler:=Teiler*Zahl;
       Delta:=1/Teiler;
       Summe:=Summe+Vz*Delta;Vz:=-Vz;Inc(Zahl,2)
      Until Delta<Epsilon
     End
 End; 
 EulerPas:=Summe/PotPas(Pi,2*k+1)*PotPas(2,2*k+2)*FakPas(2*k)
End;

Function SiPas(x:Extended):Extended;
Var Ins,n,Fak,Pot,Delta:Extended;
    Vz                 :-1..1;
Begin
 If Abs(x)>100 Then SiPas:=ErgebnisFehler
               Else
  Begin
   Ins:=x;n:=1;Fak:=1;Pot:=x;x:=Sqr(x);Vz:=-1;
   Repeat
    Pot:=Pot*x;Delta:=1;n:=n+Delta;Fak:=Fak*n;n:=n+Delta;Fak:=Fak*n;
    Delta:=Pot;Delta:=Delta/Fak;Delta:=Delta/n;
    Ins:=Ins+vz*Delta;vz:=-vz;Delta:=Abs(Delta)
   Until Delta<1e-15;
   SiPas:=Ins
  End
End;

Function CiPas(x:Extended):Extended;
Var Ic,n,Fak,Pot,Delta:Extended;
    Vz                :-1..1;
Begin
 If (x>0) And (x<=100) Then
  Begin
   Ic:=Ln(x);n:=0;Fak:=1;Pot:=1;x:=Sqr(x);Vz:=-1;
   Repeat
    Pot:=Pot*x;Delta:=1;n:=n+Delta;Fak:=Fak*n;n:=n+Delta;Fak:=Fak*n;
    Delta:=Pot;Delta:=Delta/Fak;Delta:=Delta/n;Ic:=Ic+vz*Delta;vz:=-vz
   Until Delta<1e-15;
   CiPas:=Ic+EulerKonstante
  End
                       Else CiPas:=ErgebnisFehler
End;

Function TanPas(x:Extended):Extended;
Begin
 If Cos(x)<>0 Then TanPas:=Sin(x)/Cos(x) Else TanPas:=ErgebnisFehler
End;

Function TaniPas(x:Extended):Extended;
Const Epsilon=1E-12;
Var   n            :Byte;
      Delta,Zwei,
      Pot,Fak,Summe:Extended;
Begin
 If (x<-Pi/2+0.1) Or (x>Pi/2-0.1) Then Begin TaniPas:=ErgebnisFehler;Exit End;
 n:=2;Summe:=x;Zwei:=4;Pot:=x;Fak:=2;
 Repeat
  Zwei:=Zwei*4;Pot:=Pot*x*x;Fak:=Fak*(2*n-1)*(2*n);
  Delta:=Zwei*(Zwei-1)*BernoulliPas(n)/Fak/(2*n-1)*Pot;
  Summe:=Summe+Delta;Inc(n)
 Until Delta<Epsilon;
 TaniPas:=Summe
End;

Function CotPas(x:Extended):Extended;
Begin
 If Sin(x)<>0 Then CotPas:=Cos(x)/Sin(x) Else CotPas:=ErgebnisFehler
End;

Function CotiPas(x:Extended):Extended;
Const Epsilon=1E-12;
Var   n            :Byte;
      Delta,Zwei,
      Pot,Fak,Summe:Extended;
Begin
 If x=0 Then Begin CotiPas:=ErgebnisFehler;Exit End;
 If (x<-2.5) Or (x>2.5) Then Begin CotiPas:=ErgebnisFehler;Exit End;
 n:=2;Summe:=-1/x;Zwei:=4;Pot:=x;Fak:=2;
 Repeat
  Zwei:=Zwei*4;Pot:=Pot*x*x;Fak:=Fak*(2*n-1)*(2*n);
  Delta:=-Zwei*BernoulliPas(n)/Fak/(2*n-1)*Pot;
  Summe:=Summe+Delta;Inc(n);
 Until Delta<Epsilon;
 CotiPas:=Summe
End;

Function Integral(von,bis:Extended;Art:Boolean):Extended;
 Function Formel(x:Extended;Art:Boolean):Extended;
 Begin
  If Art Then Formel:=Sin(Pi/2*Sqr(x)) Else Formel:=Cos(Pi/2*Sqr(x))
 End;
 Function dy(Y1,Y2:Extended):Extended;
 Begin
  If (Abs(Y1)>1000) And (Abs(Y2)>1000) Then dy:=Abs(Y1/Y2-1)
                                       Else dy:=Abs(Y1-Y2)
 End;
Const GroesseIntegral=30;
      Schranke       =1E-6;
      Durchlaeufe    =2000000;
Var   h,s,n1:Extended;
      k,l,r :0..GroesseIntegral;
      i,n   :LongInt;
      Feld  :Array[0..GroesseIntegral,0..GroesseIntegral] Of Extended;
Begin
 k:=0;h:=bis-von;Integral:=ErgebnisFehler;
 s:=Formel(von,Art)+Formel(bis,Art);
 Feld[k,0]:=0.5*h*s;n:=1;
 Repeat
  Inc(k);h:=h*0.5;s:=0;
  If k=GroesseIntegral Then Exit;
  For i:=1 To n Do s:=s+Formel(von+(2*i-1)*h,Art);
  Feld[k,0]:=0.5*Feld[k-1,0]+h*s;l:=1;n:=n shl 1;n1:=4;
  If n>Durchlaeufe Then Exit;
  For r:=k-1 DownTo 0 Do
   Begin Feld[r,l]:=(n1*Feld[r+1,l-1]-Feld[r,l-1])/(n1-1);Inc(l);n1:=n1*4 End
 Until (dy(Feld[0,l-1],Feld[0,l-2])<Schranke) And (n>16);
 Integral:=Feld[0,l-1]
End;

Function FSPas(x:Extended):Extended;
Begin
 If Abs(x)>100 Then FSPas:=ErgebnisFehler Else FSPas:=Integral(0,x,False)
End;

Function FCPas(x:Extended):Extended;
Begin
 If Abs(x)>100 Then FCPas:=ErgebnisFehler Else FCPas:=Integral(0,x,True)
End;

Function ASinPas(x:Extended):Extended;
Begin
 If Abs(x)>1 Then ASinPas:=ErgebnisFehler
             Else
  If Abs(x)=1 Then ASinPas:=SgnPas(x)*Pi/2 Else ASinPas:=ArcTan(x/Sqrt(1-Sqr(x)))
 End;

Function ACosPas(x:Extended):Extended;
Begin
 If Abs(x)>1 Then ACosPas:=ErgebnisFehler Else ACosPas:=Pi/2-ASinPas(x)
 End;

Function ACotPas(x:Extended):Extended;
Begin ACotPas:=Pi/2-ArcTan(x) End;

Function SinhPas(x:Extended):Extended;
Begin SinhPas:=(Exp(x)-Exp(-x))/2 End;

Function SihPas(x:Extended):Extended;
Var Ish,n,Fak,Pot,Delta:Extended;
Begin
 If Abs(x)>100 Then SihPas:=ErgebnisFehler
               Else
  Begin
   Ish:=x;n:=1;Fak:=1;Pot:=x;x:=Sqr(x);
   Repeat
    Pot:=Pot*x;Delta:=1;n:=n+Delta;Fak:=Fak*n;n:=n+Delta;Fak:=Fak*n;
    Delta:=Pot;Delta:=Delta/Fak;Delta:=Delta/n;Ish:=Ish+Delta;
    Delta:=Abs(Delta)
   Until Delta<1e-15;
   SihPas:=Ish
  End
End;

Function CoshPas(x:Extended):Extended;
Begin CoshPas:=(Exp(x)+Exp(-x))/2 End;

Function CihPas(x:Extended):Extended;
Var Ich,n,Fak,Pot,Delta:Extended;
Begin
 If (x>0) And (x<=100) Then
  Begin
   Ich:=Ln(x);n:=0;Fak:=1;Pot:=1;x:=Sqr(x);
   Repeat
    Pot:=Pot*x;Delta:=1;n:=n+Delta;Fak:=Fak*n;n:=n+Delta;Fak:=Fak*n;
    Delta:=Pot;Delta:=Delta/Fak;Delta:=Delta/n;Ich:=Ich+Delta
   Until Delta<1e-15;
   CihPas:=Ich+EulerKonstante
  End
                       Else CihPas:=ErgebnisFehler
End;

Function TanhPas(x:Extended):Extended;
Begin TanhPas:=(Exp(x)-Exp(-x))/(Exp(x)+Exp(-x)) End;

Function TanhiPas(x:Extended):Extended;
Const Epsilon=1E-12;
Var   n            :Byte;
      Delta,Zwei,
      Pot,Fak,Summe:Extended;
      Vz           :-1..1;
Begin
 If (x<-Pi/2+0.1) Or (x>Pi/2-0.1) Then Begin TanhiPas:=ErgebnisFehler;Exit End;
 n:=2;Summe:=x;Zwei:=4;Pot:=x;Fak:=2;Vz:=-1;
 Repeat
  Zwei:=Zwei*4;Pot:=Pot*x*x;Fak:=Fak*(2*n-1)*(2*n);
  Delta:=Zwei*(Zwei-1)*BernoulliPas(n)/Fak/(2*n-1)*Pot;
  Summe:=Summe+Vz*Delta;Inc(n);Vz:=-Vz
 Until Delta<Epsilon;
 TanhiPas:=Summe
End;

Function CothPas(x:Extended):Extended;
Begin
 If x<>0 Then CothPas:=(Exp(x)+Exp(-x))/(Exp(x)-Exp(-x)) Else CothPas:=ErgebnisFehler
 End;

Function CothiPas(x:Extended):Extended;
Const Epsilon=1E-12;
Var   n            :Byte;
      Delta,Zwei,
      Pot,Fak,Summe:Extended;
      Vz           :-1..1;
Begin
 If (x<-3) Or (x>3) Or (x=0) Then Begin CothiPas:=ErgebnisFehler;Exit End;
 n:=2;Summe:=-1/x+x/3;Zwei:=4;Pot:=x;Fak:=2;Vz:=-1;
 Repeat
  Zwei:=Zwei*4;Pot:=Pot*x*x;Fak:=Fak*(2*n-1)*(2*n);
  Delta:=Zwei*BernoulliPas(n)/Fak/(2*n-1)*Pot;
  Summe:=Summe+Vz*Delta;Inc(n);Vz:=-Vz
 Until Delta<Epsilon;
 CothiPas:=Summe
End;

Function ASinhPas(x:Extended):Extended;
Begin ASinhPas:=Ln(x+Sqrt(Sqr(x)+1)) End;

Function ACoshPas(x:Extended):Extended;
Begin
 If x>=1 Then ACoshPas:=ASinhPas(Sqrt(Sqr(x)-1))
         Else ACoshPas:=ErgebnisFehler
End;

Function ATanhPas(x:Extended):Extended;
Begin
 If Abs(x)>=1 Then ATanhPas:=ErgebnisFehler
              Else ATanhPas:=ASinhPas(x/Sqrt(1-Sqr(x)))
End;

Function ACothPas(x:Extended):Extended;
Begin
 If Abs(x)>1 Then ACothPas:=ATanhPas(1/x)
             Else ACothPas:=ErgebnisFehler
End;

Function AnnuitaetPas(Zinsen,Laufzeit:Extended):Extended;
Begin
 If (Zinsen<0) Or (Laufzeit<=0) Then AnnuitaetPas:=ErgebnisFehler
                                Else
  If Zinsen=0 Then AnnuitaetPas:=1/Laufzeit
              Else AnnuitaetPas:=Zinsen/(1-Exp(-Laufzeit*Ln(1+Zinsen)))
End;

Function LaufzeitPas(Zinsen,Annuitaet:Extended):Extended;
Begin
 If (Zinsen<0) Or (Annuitaet<=0) Or (Zinsen>=Annuitaet) Then
  LaufzeitPas:=ErgebnisFehler
                                                        Else
  If Zinsen=0 Then LaufzeitPas:=1/Annuitaet
              Else LaufzeitPas:=Ln(1-(Zinsen/Annuitaet))/Ln(1/(1+Zinsen))
End;
Function ZinsenPas(Laufzeit,Annuitaet:Extended):Extended;
Var zAlt,zNeu:Extended;
    Zaehler  :Word;
Begin
 If (Laufzeit<=0) Or (Annuitaet<1/Laufzeit) Then ZinsenPas:=ErgebnisFehler
                                            Else
  If Annuitaet=1/Laufzeit Then ZinsenPas:=0
                          Else
   Begin
    zNeu:=0.5;Zaehler:=0;
    Repeat
     zAlt:=zNeu;
     zNeu:=Annuitaet*(1-Exp(-Laufzeit*Ln(1+zAlt)));
     Inc(Zaehler);
     If Zaehler>1000 Then Begin ZinsenPas:=ErgebnisFehler;Exit End
    Until Abs(zNeu-zAlt)<1e-8;
    ZinsenPas:=zNeu
   End
End;

Function GammaPas(x:Extended):Extended;
Label 1,2,3;
Var   g,g1,zAlt,zNeu,n:Extended;
Begin
 If (x=0) Or (Abs(x)>=1000) Then GammaPas:=ErgebnisFehler
                            Else
  Begin
   g:=1;g1:=1;zNeu:=1;n:=1;
 1:If x<=1.5 Then Goto 2;
   x:=x-1;g1:=g1*x;Goto 1;
 2:If x>=0.5 Then Goto 3;
   g1:=g1/x;x:=x+1;
   If x<>0 Then Goto 2;
   GammaPas:=ErgebnisFehler;Exit;
 3:zAlt:=zNeu;zNeu:=n/(x+n-1);g:=g*zNeu;n:=n+1;
   If Abs(zNeu-zAlt)>1e-8 Then Goto 3;
   GammaPas:=g1*g*Exp((x-1)*Ln(n))
  End
End;

Function VariationenPas(x,y:Extended):Extended;
Var h:Extended;
Begin
 If (x<1) Or (y<1) Or (x>=1000) Or (y>=1000) Or (x=y) Then
  Begin VariationenPas:=ErgebnisFehler;Exit End;
 If x<y Then Begin h:=x;x:=y;y:=h End;
 VariationenPas:=FakPas(x)/FakPas(x-y)
End;

Function KombinationenPas(x,y:Extended):Extended;
Var h:Extended;
Begin
 If (x<1) Or (y<1) Or (x>=1000) Or (y>=1000) Or (x=y) Then
  Begin KombinationenPas:=ErgebnisFehler;Exit End;
 If x<y Then Begin h:=x;x:=y;y:=h End;
 KombinationenPas:=VariationenPas(x,y)/FakPas(y)
End;

Function ggTPas(x,y:Extended):Extended;
Label 1,2;
Var   h:Extended;
Begin
  x:=Int(x+0.5);y:=Int(y+0.5);
  If (x<1) Or (y<1) Or (x>1e9) Or (y>1e9) Then
   Begin ggTPas:=ErgebnisFehler;Exit End;
  If y<x Then Begin h:=x;x:=y;y:=h End;
1: h:=ModPas(y,x);
   If h=0 Then Goto 2;
   y:=x;x:=h;Goto 1;
2: ggTPas:=x
End;

Function kgVPas(x,y:Extended):Extended;
Var h:Extended;
Begin
 x:=Int(x+0.5);y:=Int(y+0.5);
 If (x<1) Or (y<1) Or (x>1e9) Or (y>1e9) Then
  Begin kgVPas:=ErgebnisFehler;Exit End;
 If y<x Then Begin h:=x;x:=y;y:=h End;
 kgVPas:=x*y/ggTPas(x,y)
End;

Function MandelPas(Bildx,Bildy:Extended;Tiefe,Schranke:Word):Extended;
Var Farbe        :Word;
    Test         :Boolean;
    u,t,x,x2,y,y2:Extended;
Begin
 Test:=False;
 If Bildx>-0.75 Then
  Begin
   u:=Bildx-0.25;t:=u*u+Bildy*Bildy;
   If (2*t+u)*(2*t+u)<=t Then Test:=True
  End
                 Else
  If (Bildx+1)*(Bildx+1)+Bildy*Bildy<=0.0625 Then Test:=True;
 If Not(Test) Then
  Begin
   x:=0;y:=0;x2:=0;y2:=0;Farbe:=0;
   While (x2+y2<Schranke) And (Farbe<Tiefe) Do
    Begin y:=2*x*y+Bildy;x:=x2-y2+Bildx;Inc(Farbe);x2:=Sqr(x);y2:=y*y End;
   MandelPas:=-Farbe
  End
              Else MandelPas:=-Tiefe
End;
Function EiPas(x:Extended):Extended;
Var n,FEi,Pot,Fak,Delta:Extended;
Begin
 If (x=0) Or (Abs(x)>50) Then Begin EiPas:=ErgebnisFehler;Exit End;
 n:=1;FEi:=EulerKonstante+Ln(Abs(x));Pot:=1;Fak:=1;
 Repeat
  Pot:=Pot*x;Fak:=Fak*n;Delta:=Pot/Fak/n;FEi:=FEi+Delta;n:=n+1
 Until Abs(Delta)<1e-15;
 EiPas:=FEi
End;

Function LiPas(x:Extended):Extended;
Begin
 If (Abs(x)=1) Or (Abs(x)>50) Then Begin LiPas:=ErgebnisFehler;Exit End;
 If x=0 Then LiPas:=0 Else LiPas:=EiPas(Ln(Abs(x)))
End;

Function PhiPas(x:Extended):Extended;
Var n,FPhi,Pot,Fak,Delta:Extended;
    vz                  :-1..1;
Begin
 If Abs(x)>25 Then Begin PhiPas:=ErgebnisFehler;Exit End;
 n:=1;FPhi:=x;Pot:=x;Fak:=1;Vz:=1;
 Repeat
  Vz:=-Vz;Pot:=Pot*Sqr(x);Fak:=Fak*n;Delta:=Pot/Fak/(n+n+1);
  FPhi:=FPhi+Vz*Delta;n:=n+1
 Until Abs(Delta)<=1e-15;
 PhiPas:=2/Sqrt(Pi)*FPhi
End;

Function LegendrePas(x:Extended;n,m:Byte):Extended;
Begin
 If n<m Then LegendrePas:=0
        Else
  If n=m+1 Then LegendrePas:=FakPas(2*m+1)/(PotPas(2,m)*FakPas(m))*x
           Else
   If n=m Then LegendrePas:=FakPas(2*m)/(PotPas(2,m)*FakPas(m))
          Else LegendrePas:=((2*n-1)*x*LegendrePas(x,n-1,m)-(n-1+m)*LegendrePas(x,n-2,m))/(n-m)
End;

Function HermitePas(x:Extended;n,m:Byte):Extended;
 Function Hermite(x:Extended;n:Byte):Extended;
 Begin
  If n=0 Then Hermite:=1
         Else If n=1 Then Hermite:=2*x
                     Else Hermite:=2*x*Hermite(x,n-1)-2*(n-1)*Hermite(x,n-2)
 End;
Begin
 If m=0 Then HermitePas:=Hermite(x,n)
        Else If n<m Then HermitePas:=0
                    Else If m=1 Then HermitePas:=2*n*Hermite(x,n-1)
                                Else HermitePas:=2*n*HermitePas(x,n-1,m-1)
End;

Function BesselPas(x:Extended;n:Byte):Extended;
Var k                    :Byte;
    b,Delta,Fak1,Fak2,Pot:Extended;
    vz                   :-1..1;
Begin
 If n>=0 Then
  Begin
   b:=0;vz:=1;Fak1:=1;Fak2:=1;Pot:=1;
   For k:=1 To n Do Begin Fak2:=Fak2*k;Pot:=Pot*0.5*x End;
   k:=0;
   Repeat
    Delta:=vz*Pot/(Fak1*Fak2);b:=b+Delta;Inc(k);vz:=-vz;
    Fak1:=Fak1*k;Fak2:=Fak2*(n+k);Pot:=Pot*Sqr(x)/4
   Until Abs(Delta)<1e-15;
   BesselPas:=b
  End
         Else BesselPas:=ErgebnisFehler
End;

Function TschebyschewPas(x:Extended;n:Byte):Extended;
Var z:Extended;
Begin
 If Abs(x)<1 Then TschebyschewPas:=ErgebnisFehler
              Else
  Begin
   z:=Sqrt(Sqr(x)-1);
   TschebyschewPas:=(PotPas(x+z,n)+PotPas(x-z,n))/PotPas(2,n)
  End
End;

Function FibPas(x:Extended):Extended;
Begin
 If (x<0) Or (x>20) Then FibPas:=ErgebnisFehler
                    Else
  Begin
   x:=Round(x);
   IF x<2 Then FibPas:=x
          Else FibPas:=FibPas(x-1)+FibPas(x-2)
  End
End;

Function BinomialPas(n,k:Extended):Extended;
Var i:Byte;
    B:Extended;
Begin
 BinomialPas:=ErgebnisFehler;
 n:=Round(n);k:=Round(k);
 If (n<1) Or (n>100) Then Exit;
 If (k<0) Or (k>n) Then Exit;
 B:=1;
 For i:=Round(n) DownTo Round(n)-Round(k)+1 Do B:=B*i;
 BinomialPas:=B/FakPas(k)
End;

Function HofPas(k:Extended):Extended;
Begin
 If (Round(k)=1) Or (Round(k)=2) Then HofPas:=1
                                 Else
  HofPas:=HofPas(k-HofPas(k-1))+HofPas(k-HofPas(k-2))
End;

Function ZetaPas(x:Extended):Extended;
Const A        =1.312251738456346;
      B        =0.851136628120492;
      Ausgleich=0.007578105045971;
      n        =1.39;
      Schranke =1E-18;
Var   Delta,S:Extended;
      i      :LongInt;
Begin
 If x<=1 Then Begin ZetaPas:=ErgebnisFehler;Exit End;
 If x>50 Then Begin ZetaPas:=1;Exit End;
 If x>=4 Then
  Begin
   i:=2;S:=1;
   Repeat
    Delta:=1/Exp(Ln(i)*x);S:=S+Delta;Inc(i)
   Until Delta<Schranke;
   ZetaPas:=S
  End
         Else ZetaPas:=A/(PotPas(x,n)-1)+B+Ausgleich
End;

Function GammaBerechnenPas:Extended;
Const Schranke=1E-18;
Var   Delta,C:Extended;
      i      :LongInt;
Begin
 i:=3;C:=(Sqr(Pi)/6-1)/(i-1);
 Repeat
  Delta:=(ZetaPas(i)-1)/i;C:=C+Delta;Inc(i)
 Until Delta<Schranke;
 GammaBerechnenPas:=1-C
End;

Function GoldbachPas(x:Extended):Extended;
Var Prim,Gold:LongInt;
Begin
 If Abs(x)>1E6 Then Begin GoldbachPas:=ErgebnisFehler;Exit End;
 If x<4 Then Begin GoldbachPas:=0;Exit End;
 Gold:=0;
 For Prim:=2 To Trunc(x) Div 2 Do
  If (PrimPas(Prim)<>0) And (PrimPas(x-Prim)<>0) Then Inc(Gold);
 GoldbachPas:=Gold
End;

Initialization

Randomize

End.