(*
   Funktionen, auf die Jaochim Mohr in seinen Programmen
   nicht mehr verzichten möchte.
   (c) Joachim Mohr, Rottenburg am Neckar

   Die Unit darf frei in nichtkommerziellen Programmen verwendet
   werden, wenn der CopyRight-Vermerk nicht entfernt wird.
   Kritik, Anregungen, Verbesserungsvorschläge bitte an

   Joachim Mohr Homepage
   www.joachimmohr.de
*)

unit mathMohr;

interface

FUNCTION TermToReal(s: string;x:extended): extended; overload;
Function TermToReal(s: string)       : extended; overload;
FUNCTION ReellZuBruch(const q:extended):string;
FUNCTION ReellZuGemZahl(const q:extended):string;
function AlsWurzeloderPi(x:extended):string;

implementation

uses controls,dialogs,sysutils;

var  AnzahlDerPrimzahlenBisher: integer;
     Primzahlen: array of integer;

//Allgemeine Funktionen
procedure showmessage0(s:string);
begin
  showmessage('Fehler:'+s);
end;

function copyab(const s:string; const i:integer):string;  //Rest von s ab i. em Zeichen
  begin result:=copy(s,i,length(s)-i+1) end;

function ohneLeerzeichenf(const s:string):string;
   var n:integer;
 begin
   result:=s;
   repeat
     n:=pos(' ',result);
     if n>0 then result:=copy(result,1,n-1)+copyab(result,n+1)
   until n=0;
 end;

//----------- im Parser vorkommende Funktionen -----------------

function hoch(const x:extended; n:integer):extended; //x^n  n El Z
  begin
    if n<0 then result:=1/hoch(x,-n) else //x^n=1/x^(-n) für n<0
      if n=0 then  result:=1 else
        result:=x*hoch(x,n-1)//rekursiv
  end;


function ArcTan2(Y, X: Extended): Extended;
asm
        FLD     Y
        FLD     X
        FPATAN
        FWAIT
end;

function ArcCos(X: Extended): Extended;
begin
  Result := ArcTan2(Sqrt(1 - X*X), X);
end;

function ArcSin(X: Extended): Extended;
begin
  Result := ArcTan2(X, Sqrt(1 - X*X))
end;


function asn(x:extended):extended; //arcsin
begin
  result:=0;
  if abs(x-1)<1E-15 then result:=pi/2 else     //asn(1)=Pi/2
    if abs(x+1)<1E-15 then result:=-Pi/2 else  //asn(-1)=-Pi/2
    if x*x<1 then result:=arcsin(x) //siehe unit tttmath
     else showmessage0('asn'+FloatToStr(x))
end;

function acs(x:extended):extended;   //Arccos
begin
 result:=0;
 if abs(x-1)<1E-15 then result:=0 else
   if abs(x+1)<1E-15 then result:=Pi else
     if x*x<1 then result:=arcCos(x) //siehe unit tttmath
     else showmessage0('asn'+floatToStr(x))
end;

function tan(x:extended):extended;
begin
 result:=0;
 if cos(x)<>0 then tan:=sin(x)/cos(x)
    else showmessage0('tan'+FloatToStr(x))
end;

function wurzel(x:extended):extended;
begin
  result:=0;
  if x>=-1E-15 //nicht negativ bis auf Rundungsfehler
    then wurzel:=sqrt(x)
      else showmessage0('sqr'+floatToStr(x))
end;



//----------------- Es folgt ein Parser ---------------------

FUNCTION TermToReal(s: string;x:extended): extended; overload;//rekursiv!
                                  // s beliebiger Term ohne wissensch Not
                                  // d.h. 1E-17 nicht erlaubt
    var u,v,p,q: string;
    function e(c:char): boolean; //falls möglich wird s zerlegt in s=u+c+v
            //z.B. s=(5+4)*(3+2) Bei c='+' wird u=(5+4) und v=(3+2). Dann e=true
            //    Bei c='*' wird wegen der Klammern e=false
          var i,k:integer;
      begin k:=0; // k zählt die Klammern
            i:=length(s)+1;
            repeat  dec(i);
              if s[i]=')' then inc(k);
              if s[i]='(' then k:=pred(k)
            until (i=1) or ((k=0) and (s[i]=c));
            if i=1 then Begin
              if k<>0 then
                showmessage0(s+'|'+floatToStr(k)+' Klammern zu viel!');
              result:=false
            End else Begin
              result:=true;
              u:=copy(s,1,i-1);
              v:=copyab(s,i+1)
            End
      end;
 begin
   result:=0;
   s:=OhneLeerzeichenf(s);
   u:=copy(s,1,3);
   v:=copyab(s,4);
   p:=copy(s,1,2);
   q:=copyab(s,3);
   if s='' then Begin result:=0; exit End;
   if s[1]='-' then s:='0'+s; //zB. s='-7/3x+14' -> s='0-7/3x+14'
 try {*****}
 {Punkt- vor Strichrechnung und vor Potenzen}
 if e('+') then result:=TermToReal(u,x)+TermToReal(v,x)        else
 if e('-') then result:=TermToReal(u,x)-TermToReal(v,x)        else
 if e('*') then result:=TermToReal(u,x)*TermToReal(v,x)        else
 if e('·') then result:=TermToReal(u,x)*TermToReal(v,x)        else //'·'=#183
 if e('/') or e(':') then Begin
   try result:=TermToReal(u,x)/TermToReal(v,x);
   except showmessage0(u+'/'+v) end
 End else
  if e('^') then result:=hoch(TermToReal(u,x),round(TermToReal(v,x)))  else
  if p='lg' then Begin
    try result:=ln(TermToReal(q,x))/ln(10)
    Except showmessage0('Lg'+q) End;
  End else
  if p='ln' then Begin
    try result:=ln(TermToReal(q,x))
    except showmessage0('ln'+q) End
  End else
  if p='lb' then Begin
    try result:=ln(TermToReal(q,x))/ln(2);
    except showmessage0('lb'+q) End
  End else
  if u='sin' then result:=sin(TermToReal(v,x))             else //sin im Bogenmaß
  if u='cos' then result:=cos(TermToReal(v,x))             else
  if u='si_' then result:=sin(Pi/180*TermToReal(v,x))      else //sin im Gradmaß
  if u='atn' then Begin
    try result:=arctan(TermToReal(v,x));
    except showmessage0('atn'+v) End
  End else
  if u='asn' then result:=asn(TermToReal(v,x)) {Fehler s.o.} else
  if u='acs' then result:=acs(TermToReal(v,x)) {Fehler bei asn} else
  if u='at_' then Begin  //arTan im Gradmaß
    try result:=arctan(TermToReal(v,x))*180/Pi;
    except showmessage0('at_'+v) End
  End else
  if u='as_' then Begin //arcsin im Gradmaß
    try result:=asn(TermToReal(v,x))*180/Pi;
    except showmessage0('as_'+v) End
  End else
  if u='ac_' then Begin //arccos im Gradmaß
    try result:=acs(TermToReal(v,x))*180/Pi;
    except showmessage0('ac_'+v) End
  End      else
  if u='co_' then Begin  //arccos im Gradmaß
    try result:=cos(Pi/180*TermToReal(v,x));
    except showmessage0('co_'+v) End
  End      else
  if u='tan' then Begin
    try result:=tan(TermToReal(v,x));
    except showmessage0('tan'+v) End
  End             else
  if u='ta_' then Begin
    try result:=tan(Pi/180*TermToReal(v,x));
    except showmessage0('ta_'+v) End
  End else
  if u='abs' then result:=abs(TermToReal(v,x))             else
  if u='exp' then result:=exp(TermToReal(v,x))             else
  if u='wur' then result:=wurzel(TermToReal(v,x))          else //Dort Fehler
  if u='int' then result:=int(TermToReal(v,x))             else
  if u='rou' then result:=int(TermToReal(v,x)+0.5)         else
  if s[1]='(' then Begin
     while (length(s)>1) and (s[length(s)]<>')') do s:=copy(s,1,length(s)-1);
    //z.B. 12*(250+50) DM
    s:=copy(s,2,length(s)-2);
    result:=TermToReal(s,x)
  End else if s='Pi' then result:=Pi else
    if s='x' then result:=x else //Variable x
    BEgin
      try result:=strToFloat(s); except showmessage0('Syntaxf. in '+s); result:=0; exit end;
    ENd
  except {******} showmessage0('Fehler in '+s) end;
 end;

Function TermToReal(s: string): extended; overload;
  const dummy=0;
begin
  result:=termToReal(s,dummy);
end;

//------------------------------ Es folgt Bruc darstellung  ------------------------
function char_l(const s:string; const k:integer):char; //falls k zu gross ' '
begin if (k<=0) or (k>length(s)) then result:=' ' else result:=s[k] end;

function char_last(const s:string):char;
begin result:=char_l(s,length(s)) end;

procedure kup(var s:string); overload;
begin s:=copy(s,1,length(s)-1) end;

FUNCTION ggtInt(a,b:longint):longint;
  begin if b=0 then result:=a
        else result:=ggtInt(b,a mod b); //...ist göttlich
  end;

FUNCTION ggtReal(a,b:Extended):Extended;
  begin if (a<maxlongint) and (b<maxlongint) then
           result:=ggTInt(round(a),round(b))
        else Begin
          if abs(b)<0.5 then result:=a
          else result:=ggtReal(b,a-b*int(a/b));
        end;
end;

PROCEDURE kuerzeReal(var a,b:Extended);
    var t:Extended;
  begin t:=ggtReal(a,b);
        a:=a/t;
        b:=b/t
  end;

function ErmittleBruch(const q:extended;var a,b:extended; eps_Genauigkeit:Extended):boolean;
       //=true, falls erfolgreich
    const MaxIteration=50;
    var q0,Qabs,GanzZahligerAnteilVonQ,NachkommaVonQ,am,amm,bmm,bm:extended;
       {am=a index -1 amm index -2}
        zaehl:integer;
begin //Methode: Kettenbruchentwicklung
  result:=true;
  qabs:=abs(q); //Es wird mit positiven Werten gerechnet
  q0:=qabs;
  a:=1;
  b:=0;
  am:=0; //altes a
  bm:=1; //altes b
  zaehl:=0;
  try //damit nach exit Finally ausgeführt wird
    repeat
      inc(zaehl);
      if zaehl>MaxIteration then exit; //falls eps_Genauigkeit zu knapp
      GanzZahligerAnteilVonQ:=int(q0+eps_Genauigkeit); //wg. Rundungsfeler muss z.B.  4.99.. zu 5.00.. werden
      NachkommaVonQ:=q0-GanzZahligerAnteilVonQ;        //=fraq(qabs)
      amm:=am; //uraltes a
      bmm:=bm; //uraltes b
      am:=a;   //altes   a
      bm:=b;   //altes   a
      try // möglich overflow (kommt nur vor, falls eps_Gen. zu klein)
        a:=GanzZahligerAnteilVonQ*am+amm; //Anfang: ak=GanzZahligerAnteilVonQ
        b:=GanzZahligerAnteilVonQ*bm+bmm; //        bk=1
        kuerzeReal(a,b);
      Except exit End;
      if NachkommaVonq<eps_Genauigkeit then exit;//endlicher Kettenbruch
      q0:=1/NachkommaVonq; //Jetzt ja NachkommaVonq>=Genauigkeit
    until abs(qabs-a/b)<eps_Genauigkeit;
  Finally
    if abs(qabs-a/b)>=eps_Genauigkeit then result:=false;
    if q<0 then a:=-a;
  End;
end;



function ReellZuBruch_(const q:extended;const alsGemZahl:boolean):string;
  const Hochkomma=#39;
  var a,b,g,r:extended;
      Erf:boolean;
  function ohneE(const s:string):string; //z.b. 1,4379...E18
      var p,n:integer;
          a:string;
    begin
       p:=pos('E',s);
       a:=copy(s,1,p-1);
       n:=StrToInt(copyab(s,p+1));
       result:=char_l(a,1)+copyab(a,3);
       while length(result)<n do result:=result+'0';
    end;
begin
  erf:=ErmittleBruch(q,a,b,1E-15);  //1/2+1/7 wird mit Fehler>1E-16 berechnet
  //Teste auch 9,999 999 999  und Pi/4
  if not erf or (abs(a)*abs(b)>1E14) then  Begin
    result:=trim(FloatToStrF(q,ffFixed,18,10));
    if pos('E',result)>0 then result:=OhneE(result) else
      while (pos(Decimalseparator,result)>0) and (char_last(result)='0')
        and (pos(Decimalseparator,result)<length(result)-1) do kup(result);
  End else if b=1 then result:=FloatToStr(a) else Begin
    if alsGemZahl and (abs(q)>1) then BEgin
      g:=int(abs(q));     //g ganzzahliger Teil von q
      r:=abs(a)-g*b;     //da a=g*b+r  q=a/b=g+r/b
      result:=FloatToStr(g)+Hochkomma+FloatToStr(r)+'/'+FloatToStr(b);
      if q<0 then result:='-'+result;
    ENd else result:=FloatToStr(a)+'/'+FloatToStr(b);
  End;
end;

function ReellZuBruch(const q:extended):string;
begin result:=ReellZuBruch_(q,false) end;

function ReellZuGemZahl(const q:extended):string;
begin result:=ReellZuBruch_(q,true) end;

function RealIsIntegerToStr(x:extended):string;
begin
  result:=trim(floatToStrf(x,ffGeneral,4,18));
end;

function VielfachesVonPi(x:extended):string;
   var z,n:extended;
begin
  result:=''; //Misserfolg
  if ErmittleBruch(x/Pi,z,n,1E-15) then Begin
     if abs(n)<1000000 then BEgin
       if n=1 then result:=RealIsIntegerToStr(z)+'*Pi' else
         result:=RealIsIntegerToStr(z)+'/'+RealIsIntegerToStr(n)+'*Pi'
     ENd;
  End;
end;

procedure findenaechstePrimzahl;
      var groesstePZ:longint;
    function hatKeineTeiler:boolean; //Versuch, ob groesstePrimzahl wirklich PZ ist, d.h.
              //nicht durch bisherige PZ geteilt werden kann
              //Prüfung:primzahl[j] Teiler von GroesstePrimzahl?  bis
              //        primzahl[j]>Wurzel(groesstePrimzahl) <beide hoch zwei ... s.u.

        var j:integer;
      begin
        result:=true; //falls bis zur Wurzel alle Primzahlen durchprobiert
        for j:=1 to AnzahlDerPrimzahlenBisher do Begin //ist 3,5,7 ... Teiler
          if groesstePZ<primzahlen[j]*primzahlen[j] then exit;
          if abs(groesstePZ-int(GroesstePZ/primzahlen[j])*primzahlen[j])<1/2 then
            Begin result:=false; exit End; //Teiler gefunden
        End;
      end;
    begin
        groesstePZ:=primzahlen[AnzahlDerPrimzahlenBisher]; //Größte bisher gefundene PZ
        repeat //In Prozeur, die diese aufruft, muß Beginkontr stehen!
          inc(groesstePZ,2); //Eine noch groesserPZ wird gesucht, evl. PZ-Zwilling
        until hatkeineTeiler;
        inc(AnzahlDerPrimzahlenBisher); //Jetzt gibt's eine PZ mehr
        if length(primzahlen)<=AnzahlDerPrimzahlenBisher+1 then
          setlength(primzahlen,AnzahlDerPrimzahlenBisher+500); //dyn Array
        primzahlen[AnzahlDerPrimzahlenBisher]:=groesstePZ; //Neue größte PZ
    end;

function FrageAbbruch:boolean; //=true, wenn Abgebrochen werden soll
begin
  if MessageDlg('Suche nach Wurzelausdruck kann lange dauern. Suche jetzt abbrechen?',
                 mtConfirmation, [mbYes, mbNo], 0) = mrYes then
                   result:=true else result:=false;
end;

function IsInteger(const x:extended;eps_Genauigkeit:extended):boolean;
   //eps_Genauigkeit>=1E-19
   var x0:Extended;
begin
  x0:=abs(x);
  while x0>1 do Begin
    x0:=x0/2;
    eps_Genauigkeit:=eps_Genauigkeit*2;
  End; // 1000000,001 oder 1,000000001 werden gleich behandelt
  result:=frac(abs(x)+eps_Genauigkeit)<eps_Genauigkeit*2;
end;

function AlsWurzeloderPi(x:extended):string;
         //es wird versucht x=a/b*wur(c)
         //-> x^2=z/n Nenner klein genug, dann x=1/n*wur(z*n)
        var c,xq,z,n,qu:extended;
            NrPZ:integer;
    function mitWurzel(x,z,n,c:extended):string;
        var wu:string;
      begin
        kuerzeReal(z,n);
        if abs(c-1)<1E-15 then wu:='' else wu:='*sqrt('+RealIsIntegerToStr(c)+')';
        if n=1 then Begin
          if z=1 then result:=copyab(wu,2) else result:=RealIsIntegerToStr(z)+wu
        End else result:=RealIsIntegerToStr(z)+'/'+
                             RealIsIntegerToStr(n)+wu;
        if x<0 then result:='-'+result;
      end;
  begin
      result:='';
      xq:=x*x;
      if abs(xq)<1E-15 then Begin
        result:=VielfachesVonPi(x);
        exit;
      End;
      if ErmittleBruch(xq,z,n,1E-15) then if abs(n)<1000000000 then Begin
           //x^2=z/n => x=1/n*wur(z*n);
         //z.B. 123/124*wur(2/19) hoch 2=  15 129*146 072
        c:=z*n; //x=1/n*wur(c)
      //Teilweise radizieren!
        z:=1;
        if c>1000000 then Begin result:=VielfachesVonPi(x); exit End;
        NrPZ:=0;
        repeat
          if NrPZ<AnzahlDerPrimzahlenBisher then findenaechstePrimzahl;
          if AnzahlDerPrimzahlenBisher mod 500 =0 then if FrageAbbruch then Begin
            VielfachesVonPi(x);
            exit;
          End;
          qu:=primzahlen[NrPz]*primzahlen[NrPz];
          while IsInteger(c/qu,1E-15) do Begin
            z:=z*primzahlen[NrPz];
            c:=c/qu;
          End;
          inc(NrPz);
        until (NrPz>1000) or (qu>c);
        result:=mitWurzel(x,z,n,c);
      End else result:=VielfachesVonPi(x); //
  end;


begin
  //initialisierung wird nur ebraucht für AlsWurzeloderPi
  setlength(Primzahlen,3);
  Primzahlen[0]:=2;
  primzahlen[1]:=3;
  AnzahlDerPrimzahlenBisher:=1; //Zählung ab 0
end.