unit matheunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, Gmp, Math, lowlevelunit, protokollunit;

type
  tCallBackGetValue = function(name: string): extended of object;
  tExprToFloat = function(syntaxtest: boolean; name: string): extended of object;
  tKnownValue = record
    name:  string;
    value: extended;
  end;
  tKnownValues = class
  private
    kvs: array of tKnownValue;
    function finde(nam: string): longint; inline;
  public
    constructor create;
    destructor destroy; override;
    procedure add(val: tKnownValue); inline; overload;
    procedure add(nam: string; val: extended); overload;
    function rem(nam: string): boolean;
    function extract(nam: string; out val: extended): boolean; inline;
    procedure dump(prot: tProtokollant; prefix: string);
  end;

function plus(a,b: tExtPoint): tExtPoint;
function durch(a: tExtPoint; b: extended): tExtPoint;
function myFrac(x: extended): extended;
function mpfToStr(f: mpf_t): string;
function signSqr(x: extended): extended; inline;
function mpfMyRoot(rad: mpf_t; wzlExp: int64): extended;
function myTimeToStr(t: extended): string;
function cmpStr(s1,s2: string): longint;
function mitte(s1,s2: string): string;
function myFloatToStr(x: extended): string;
function myStrToFloat(s: string): extended;
function istGanzZahl(s: string): boolean;
procedure copyArray(i: tExtPointArray; out o: tExtPointArray); overload;
procedure copyArray(i: tLongintArray; out o: tLongintArray); overload;
procedure copyArray(i: tExtendedArray; out o: tExtendedArray); overload;
function nullfunktion(x: extended): extended;
function exprtofloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): extended;
function knownValue(nam: string; val: extended): tKnownValue;
procedure fft(var res,ims: tExtendedArray; inv: boolean); inline;

implementation

// tKnownValueArray ************************************************************

constructor tKnownValues.create;
begin
  inherited create;
  fillchar(kvs,sizeof(kvs),#0);
  add('π',pi);
end;

destructor tKnownValues.destroy;
begin
  setlength(kvs,0);
  inherited destroy;
end;

function tKnownValues.finde(nam: string): longint;
var
  i: longint;
begin
  result:=-1;
  for i:=0 to length(kvs)-1 do
    if kvs[i].name=nam then begin
      result:=i;
      exit;
    end;
end;

procedure tKnownValues.add(val: tKnownValue);
begin
  add(val.name,val.value);
end;

procedure tKnownValues.add(nam: string; val: extended);
var
  i: longint;
begin
  i:=finde(nam);
  if i<0 then begin
    i:=length(kvs);
    setlength(kvs,length(kvs)+1);
    kvs[i].name:=nam;
  end;
  kvs[i].value:=val;
end;

function tKnownValues.rem(nam: string): boolean;
var
  i: longint;
begin
  i:=finde(nam);
  result:=i>=0;
  if result then begin
    while i<length(kvs)-1 do begin
      kvs[i].name:=kvs[i+1].name;
      kvs[i].value:=kvs[i+1].value;
      inc(i);
    end;
    setlength(kvs,length(kvs)-1);
  end;
end;

function tKnownValues.extract(nam: string; out val: extended): boolean;
var
  i: longint;
begin
  i:=finde(nam);
  result:=i>=0;
  if result then
    val:=kvs[i].value;
end;

procedure tKnownValues.dump(prot: tProtokollant; prefix: string);
var
  i: longint;
begin
  for i:=0 to length(kvs)-1 do
    prot.schreibe(prefix+kvs[i].name+' = '+floattostr(kvs[i].value));
end;

// allgemeine Funktionen *******************************************************

function plus(a,b: tExtPoint): tExtPoint;
var
  c: char;
begin
  for c:='x' to 'y' do
    result[c]:=a[c]+b[c];
end;

function durch(a: tExtPoint; b: extended): tExtPoint;
var
  c: char;
begin
  for c:='x' to 'y' do
    result[c]:=a[c]/b;
end;

function myFrac(x: extended): extended;
begin
  result:=frac(x);
  while result<0 do
    result:=result+1;
end;

function mpfToStr(f: mpf_t): string;
var
  ex:  int64;
  off: byte;
begin
 result:=mpf_get_str(nil,ex,10,0,f);
 off:=1+byte(pos('-',result)=1);
 if result='' then
   result:='0'
 else if ex=1 then
   result:=copy(result,1,off)+','+copy(result,off+1,length(result)-off)
 else
   result:=copy(result,1,off)+','+copy(result,off+1,length(result)-off)+' * 10^'+inttostr(ex-1);
end;

function signSqr(x: extended): extended;
begin
 result:=sign(x)*sqr(x);
end;

function mpfMyRoot(rad: mpf_t; wzlExp: int64): extended;
var
  ex: int64;
begin
  result:=power(mpf_get_d_2exp(ex,rad),1/wzlExp);
  result:=result*power(2,ex/wzlExp);
end;

function myTimeToStr(t: extended): string;
var
  tim: int64;
begin
 tim:=floor(t*24*60*60);
 result:=inttostr(tim mod 10)+'s';
 tim:=tim div 10;
 if tim=0 then exit;
 result:=inttostr(tim mod 6)+result;
 tim:=tim div 6;
 if tim=0 then exit;
 result:=inttostr(tim mod 10)+'min '+result;
 tim:=tim div 10;
 if tim=0 then exit;
 result:=inttostr(tim mod 6)+result;
 tim:=tim div 6;
 if tim=0 then exit;
 result:=inttostr(tim mod 24)+'h '+result;
 tim:=tim div 24;
 if tim=0 then exit;
 result:=' '+result;
 if (tim mod 7)<>1 then
   result:='e'+result;
 result:=inttostr(tim mod 7)+'Tag'+result;
 tim:=tim div 7;
 if tim=0 then exit;
 result:=' '+result;
 if tim<>1 then
   result:='n'+result;
 result:=inttostr(tim)+'Woche'+result;
end;

function cmpStr(s1,s2: string): longint;
var
  i: longint;
begin
 for i:=1 to min(length(s1),length(s2)) do
   if s1[i]<>s2[i] then begin
     result:=2*byte(s1[i]>s2[i])-1;
     exit;
   end;
 if length(s1)<>length(s2) then begin
   result:=2*byte(length(s1)>length(s2))-1;
   exit;
 end;
 result:=0;
end;

function mitte(s1,s2: string): string;
var
  i:    longint;
  w,nw: word;
begin
 setlength(result,max(length(s1),length(s2)));
 w:=0;
 for i:=length(result) downto 1 do begin // result:= "s1+s2";
   if i<=length(s1) then
     w:=w+byte(s1[i]);
   if i<=length(s2) then
     w:=w+byte(s2[i]);
   result[i]:=char(w and $ff);
   w:=w shr 8;
 end;
 result:=char(w)+result;
 w:=0;
 for i:=1 to length(result) do begin
   nw:=byte(odd(byte(result[i])+w));
   result[i]:=char((byte(result[i])+w) div 2);
   w:=nw shl 8;
 end;
 if w<>0 then
   result:=result+char(w div 2);
 if result[1]<>#0 then begin
   writeln('Fehler bei der Mittenfindeung!');
   halt;
 end;
 delete(result,1,1);
end;

function myFloatToStr(x: extended): string;
var
  i,e: longint;
begin
 e:=0;
 if x<0 then begin
   result:='-';
   x:=-x;
 end
 else
   result:='';
 if x=0 then begin
   result:='0';
   exit;
 end;
 while x<1 do begin
   dec(e);
   x:=x*10;
 end;
 while x>=10 do begin
   inc(e);
   x:=x/10;
 end;
 result:=result+char(ord('0')+floor(x))+'.';
 for i:=0 to 20 do begin
   x:=(x-floor(x))*10;
   result:=result+char(ord('0')+floor(x));
 end;
 if e<>0 then
   result:=result+'E'+inttostr(e);
end;

function myStrToFloat(s: string): extended;
var
  i,e: longint;
  neg: boolean;
begin
 if pos('E',s)>0 then begin
   e:=strtoint(rightStr(s,length(s)-pos('E',s)));
   delete(s,pos('E',s),length(s));
 end
 else e:=0;
 if pos('.',s)=0 then begin
   result:=strtoint(s)*power(10,e);
   exit;
 end;
 neg:=leftStr(s,1)='-';
 if neg then
   delete(s,1,1);
 if pos('.',s)=2 then begin
   result:=0;
   for i:=length(s) downto 3 do
     result:=result/10 + ord(s[i])-ord('0');
   result:=result/10 + ord(s[1])-ord('0');
 end
 else result:=strtofloat(s);
 result:=result*power(10,e);
 if neg then
   result:=-result;
end;

function istGanzZahl(s: string): boolean;
var
  i:              longint;
begin
  result:=length(s)>0;
  result:=result and (length(s) > byte(s[1] in ['-','+']));
  if result then
    for i:=byte(s[1] in ['-','+'])+1 to length(s) do
      result:=result and (s[i] in ['0'..'9']);
end;

procedure copyArray(i: tExtPointArray; out o: tExtPointArray);
var
  j: longint;
  c: char;
begin
  setlength(o,length(i));
  for j:=0 to length(o)-1 do
    for c:='x' to 'y' do
      o[j,c]:=i[j,c];
end;

procedure copyArray(i: tLongintArray; out o: tLongintArray);
var
  j: longint;
begin
  setlength(o,length(i));
  for j:=0 to length(o)-1 do
    o[j]:=i[j];
end;

procedure copyArray(i: tExtendedArray; out o: tExtendedArray);
var
  j: longint;
begin
  setlength(o,length(i));
  for j:=0 to length(o)-1 do
    o[j]:=i[j];
end;

function Nullfunktion(x: extended): extended;
begin
  result:=0*x;
end;

function exprtofloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): extended;
var i,j,k,l,m:   longint;
    inv,neg,cbv: boolean;
    val1,val2:   extended;
const
    fkt1: array[0..6] of string = ('exp','sin','cos','tan','sqr','sqrt','ln');
    fkt2: array[0..1] of string = ('min','max');
begin
  s:=trimAll(s);

  for i:=0 to length(fkt1)-1 do
    while fktpos(fkt1[i],s)>0 do begin
      j:=fktpos(fkt1[i],s)+length(fkt1[i]);
      while (j<=length(s)) and (s[j]<>'(') do
        inc(j);
      m:=j+1;                // erstes Zeichen innerhalb der Klammern
      k:=1;
      while (j<length(s)) and (k>0) do begin
        inc(j);
        case s[j] of
          '(': inc(k);
          ')': dec(k);
        end;
      end;
      k:=fktpos(fkt1[i],s); // erstes Zeichen des Funktionsnamens
      val1:=exprtofloat(st,copy(s,m,j-m),kvs,cbgv);
      case i of
        0: val1:=exp(val1);
        1: val1:=sin(val1);
        2: val1:=cos(val1);
        3: val1:=tan(val1);
        4: val1:=sqr(val1);
        5: val1:=sqrt(val1);
        6: val1:=ln(val1);
      end{of case};
      s:=copy(s,1,k-1) + floattostr(val1) + copy(s,j+1,length(s));
    end;

  for i:=0 to length(fkt2)-1 do
    while fktpos(fkt2[i],s)>0 do begin
      j:=fktpos(fkt2[i],s)+length(fkt2[i]);
      while (j<=length(s)) and (s[j]<>'(') do
        inc(j);
      m:=j+1;                // erstes Zeichen innerhalb der Klammern
      k:=1;
      l:=-1;
      while (j<length(s)) and (k>0) do begin
        if (k=1) and (s[j] in [',',';']) then
          l:=j;              // das Komma/Semikolon
        inc(j);
        case s[j] of
          '(': inc(k);
          ')': dec(k);
        end;
      end;
      k:=fktpos(fkt1[i],s); // erstes Zeichen des Funktionsnamens
      val1:=exprtofloat(st,copy(s,m,l-m),kvs,cbgv);
      val2:=exprtofloat(st,copy(s,l+1,j-l-1),kvs,cbgv);
      case i of
        0: val1:=min(val1,val2);
        1: val1:=max(val1,val2);
      end{of case};
      s:=copy(s,1,k-1) + floattostr(val1) + copy(s,j+1,length(s));
    end;

  while pos('(',s)>0 do begin
    i:=pos('(',s);
    j:=1;
    while j>0 do begin
      inc(i);
      case s[i] of
        '(': inc(j);
        ')': dec(j);
      end;
    end;
    s:=copy(s,1,pos('(',s)-1)+
       floattostr(exprtofloat(st,copy(s,pos('(',s)+1,i-pos('(',s)-1),kvs,cbgv))+
       copy(s,i+1,length(s)-i);
  end;
  if (binOpPos('+',s)>0) or (binOpPos('-',s)>0) then begin
    result:=0;
    inv:=false;
    repeat
      i:=min(binOpPos('+',s),binOpPos('-',s));
      if i=0 then i:=max(binOpPos('+',s),binOpPos('-',s));
      if i=0 then i:=length(s)+1;
      if inv then result:=result-exprtofloat(st,copy(s,1,i-1),kvs,cbgv)
             else result:=result+exprtofloat(st,copy(s,1,i-1),kvs,cbgv);
      inv:=s[i-byte(i>length(s))]='-';
      delete(s,1,i);
    until s='';
    exit;
  end;
  if (binOpPos('*',s)>0) or (binOpPos('/',s)>0) then begin
    result:=1;
    inv:=false;
    repeat
      i:=min(binOpPos('*',s),binOpPos('/',s));
      if i=0 then i:=max(binOpPos('*',s),binOpPos('/',s));
      if i=0 then i:=length(s)+1;
      if inv then result:=result/exprtofloat(st,copy(s,1,i-1),kvs,cbgv)
             else result:=result*exprtofloat(st,copy(s,1,i-1),kvs,cbgv);
      inv:=s[i-byte(i>length(s))]='/';
      delete(s,1,i);
    until s='';
    exit;
  end;
  if binOpPos('^',s)>0 then begin
    i:=binOpPos('^',s);
    result:=power(exprtofloat(st,copy(s,1,i-1),kvs,cbgv),
                  exprtofloat(st,copy(s,i+1,length(s)-i),kvs,cbgv));
    exit
  end;
  neg:=startetMit('-',s);
  cbv:=false;
  for i:=1 to length(s) do
    cbv:=cbv or not (s[i] in ['-','0'..'9','.',',','e','E']);
  if not cbv then result:=strtofloat(s)
  else if not kvs.extract(s,result) then begin
    if st then result:=1
    else if assigned(cbgv) then result:=cbgv(s)
    else begin
      result:=nan;
      gibAus(''''+s+''' kann ich nicht auflösen!',3);
      halt(1);
    end;
  end;
(*   if s='np' then result:=params.np
  else if s='maxw' then result:=params.maxW
  else if s='minw' then result:=params.minW
  else if s='beta' then result:=params.beta
  else if s='xstart' then result:=params.xstart
  else if s='xstop' then result:=params.xstop
  else if s='tstart' then result:=params.tstart
  else if s='tstop' then result:=params.tstop
  else if st then result:=1
  else if assigned(cbgv) then result:=cbgv(s)
  else result:=nan; *)

  result:=result*(2*byte(not neg)-1);
end;

function knownValue(nam: string; val: extended): tKnownValue;
begin
  result.name:=nam;
  result.value:=val;
end;

procedure fft(var res,ims: tExtendedArray; inv: boolean);
var
  haL,i,j,k,n,dist,wstep,absch,wnum: longint;
  umsortierung:                      tLongintArray;
  wRe,wIm:                           tExtendedArray;
  t1,t2:                             extended;
begin
  haL:=length(res) div 2;
  n:=round(ln(length(res))/ln(2));

  if round(power(2,n))<>length(res) then
    raise exception.create(inttostr(length(res))+' ist keine Zweierpotenz!');

  setlength(umsortierung,2*haL);
  for j:=0 to 2*haL-1 do begin
    k:=0;
    for i:=0 to n-1 do             // Bits umkehren
      k:=2*k + byte(odd(j shr i));
    umsortierung[j]:=k;
  end;
  setlength(wRe,haL);
  setlength(wIm,haL);
  for i:=0 to haL-1 do begin
    wRe[i]:=cos(pi*i/haL);
    wIm[i]:=sin(pi*i/haL);
  end;

  for j:=0 to 2*haL-1 do begin
    if umsortierung[j]>j then begin
      t1:=res[j];
      res[j]:=res[umsortierung[j]];
      res[umsortierung[j]]:=t1;
      t1:=ims[j];
      ims[j]:=ims[umsortierung[j]] * (1-2*byte(inv));
      ims[umsortierung[j]]:=t1 * (1-2*byte(inv));
    end;
    if (umsortierung[j]=j) and inv then
      ims[j]:=-ims[j];
  end;

  dist:=1;
  wstep:=haL;

  while dist<2*haL do begin
    absch:=0;
    while absch<2*haL do begin
      wnum:=0;
      for j:=0 to dist-1 do begin
        // x_links:  [absch+j]
        // x_rechts: [absch+j+dist]
        t1:=wRe[wnum]*res[absch+j+dist] - wIm[wnum]*ims[absch+j+dist];
        t2:=wRe[wnum]*ims[absch+j+dist] + wIm[wnum]*res[absch+j+dist];

        res[absch+j+dist]:=res[absch+j]-t1;
        ims[absch+j+dist]:=ims[absch+j]-t2;
        res[absch+j]:=res[absch+j]+t1;
        ims[absch+j]:=ims[absch+j]+t2;
        wnum:=wnum+wstep;
      end;
      absch:=absch+2*dist;
    end;
    wstep:=wstep div 2;
    dist:=dist*2;
  end;

  t1:=1/sqrt(length(res));
  for i:=0 to 2*haL-1 do begin
    res[i]:=res[i]*t1;
    ims[i]:=ims[i]*t1;
  end;
end;

end.