unit matheunit;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils{$IFNDEF win64}{$IFNDEF win32}, gmp{$ENDIF}{$ENDIF}, 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
    chefs,knechte: array of tKnownValues;
    kvs:           array of tKnownValue;
    function       finde(nam: string): longint; inline;
  public
    constructor create(chef: tKnownValues);
    constructor createFromOriginal(original: tKnownValues);
    destructor destroy; override;
    procedure add(val: tKnownValue; ignoriereWennBereitsVorhanden: boolean = false); inline; overload;
    procedure add(nam: string; val: extended; ignoriereWennBereitsVorhanden: boolean = false); overload;
    function rem(nam: string): boolean;
    function extract(nam: string; out val: extended): boolean; inline;
    procedure dump(prot: tProtokollant; prefix: string); overload;
    procedure dump(prefix: string; out s: string); overload;
    procedure neuerKnecht(knecht: tKnownValues);
    procedure neuerChef(chef: tKnownValues; ignoriereWennBereitsVorhanden: boolean = false);
    procedure entferneKnecht(knecht: tKnownValues);
    procedure entferneChef(chef: tKnownValues);
    procedure clear;
  end;

function plus(a,b: tExtPoint): tExtPoint;
function durch(a: tExtPoint; b: extended): tExtPoint;
function myFrac(x: extended): extended;
{$IFNDEF win64}
{$IFNDEF win32}
function mpfToStr(f: mpf_t): string;
function mpfMyRoot(rad: mpf_t; wzlExp: int64): extended;
{$ENDIF}
{$ENDIF}
function signSqr(x: extended): extended; inline;
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 intToFixpoint(x,stellen: longint): string;
function floatToFixpoint(x: extended; vk,nk: longint): string;
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 exprToBool(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): boolean;
function formelnAuswerten(s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): string;
function knownValue(nam: string; val: extended): tKnownValue;
function berechneEinheitsZelle(invarianz,modulus: tInt64Point): tInt64Point;
function ggT(a,b: int64): int64;
function intRoot(n: int64): int64;
function ermittleAnstieg(dat: string; xSpalte,ySpalte: longestOrdinal): extended;
function ermittleMittelwert(dat: string; werteSpalte,gewichteSpalte: longestOrdinal): extended;
procedure sortiereNachWert(var maxima: tLongintArray; werte: tExtendedArray; logarithmischesPivot: boolean);
procedure sortiereNachDominanz(var maxima: tLongintArray; werte: tExtendedArray; nullteSpektrometerordnung: extended);
function RotationsMatrix(alpha: extended): t2x2Extended;
function einheitsVektor(alpha: extended): tExtPoint;

implementation

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

constructor tKnownValues.create(chef: tKnownValues);
begin
  inherited create;
  fillchar(kvs,sizeOf(kvs),#0);
  fillchar(knechte,sizeOf(knechte),#0);
  add('π',pi);
  fillchar(chefs,sizeOf(chefs),#0);
  if assigned(chef) then
    neuerChef(chef);
  setLength(knechte,0);
end;

constructor tKnownValues.createFromOriginal(original: tKnownValues);
var
  i: longint;
begin
  inherited create;
  fillchar(kvs,sizeOf(kvs),#0);
  fillchar(knechte,sizeOf(knechte),#0);
  for i:=0 to length(original.kvs)-1 do
    add(original.kvs[i]);
  for i:=0 to length(original.chefs)-1 do
    neuerChef(original.chefs[i]);
end;

destructor tKnownValues.destroy;
var
  i,j: longint;
begin
  for i:=0 to length(kvs)-1 do
    for j:=0 to length(knechte)-1 do
      knechte[j].add(kvs[i],true);
  setLength(kvs,0);
  for i:=length(knechte)-1 downto 0 do begin
    for j:=0 to length(chefs)-1 do
      knechte[i].neuerChef(chefs[j],true);
    knechte[i].entferneChef(self);
  end;
  setLength(knechte,0);
  for i:=length(chefs)-1 downto 0 do
    entferneChef(chefs[i]);
  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; ignoriereWennBereitsVorhanden: boolean = false);
begin
  add(val.name,val.value,ignoriereWennBereitsVorhanden);
end;

procedure tKnownValues.add(nam: string; val: extended; ignoriereWennBereitsVorhanden: boolean = false);
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
  else
    if ignoriereWennBereitsVorhanden then
      exit;
  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;
  i:=0;
  while (not result) and (i<length(chefs)) do begin
    result:=chefs[i].extract(nam,val);
    inc(i);
  end;
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));
  for i:=0 to length(chefs)-1 do
    chefs[i].dump(prot,prefix+'chef['+intToStr(i)+']: ');
end;

procedure tKnownValues.dump(prefix: string; out s: string);
var
  i: longint;
  t: string;
begin
  s:='';
  for i:=0 to length(kvs)-1 do
    s:=s+#10+prefix+kvs[i].name+' = '+floattostr(kvs[i].value);
  for i:=0 to length(chefs)-1 do begin
    chefs[i].dump(prefix+'chefs['+intToStr(i)+'] ',t);
    s:=s+#10+t;
  end;
  delete(s,1,1);
end;

procedure tKnownValues.neuerKnecht(knecht: tKnownValues);
var
  i: longint;
begin
  for i:=0 to length(knechte)-1 do
    if knechte[i]=knecht then
      fehler('tKnownValues soll einen Knecht hinzufügen, den es schon hat!');
  setLength(knechte,length(knechte)+1);
  knechte[length(knechte)-1]:=knecht;
end;

procedure tKnownValues.neuerChef(chef: tKnownValues; ignoriereWennBereitsVorhanden: boolean = false);
var
  i: longint;
begin
  for i:=0 to length(chefs)-1 do
    if chefs[i]=chef then begin
      if ignoriereWennBereitsVorhanden then
        exit;
      fehler('tKnownValues soll einen Chef hinzufügen, den es schon hat!');
    end;
  setLength(chefs,length(chefs)+1);
  chefs[length(chefs)-1]:=chef;
  chefs[length(chefs)-1].neuerKnecht(self);
end;

procedure tKnownValues.entferneKnecht(knecht: tKnownValues);
var
  i: longint;
begin
  i:=0;
  while i<length(knechte) do begin
    if knechte[i]=knecht then begin
      while i<length(knechte)-1 do begin
        knechte[i]:=knechte[i+1];
        inc(i);
      end;
      setLength(knechte,length(knechte)-1);
      exit;
    end;
    inc(i);
  end;
  fehler('tKnownValues soll einen Knecht entfernen, den es nicht hat!');
end;

procedure tKnownValues.entferneChef(chef: tKnownValues);
var
  i: longint;
begin
  i:=0;
  while i<length(chefs) do begin
    if chefs[i]=chef then begin
      while i<length(chefs)-1 do begin
        chefs[i]:=chefs[i+1];
        inc(i);
      end;
      setLength(chefs,length(chefs)-1);
      chef.entferneKnecht(self);
      exit;
    end;
    inc(i);
  end;
  fehler('tKnownValues soll einen Chef entfernen, den es nicht hat!');
end;

procedure tKnownValues.clear;
begin
  setLength(kvs,0);
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;

{$IFNDEF win64}
{$IFNDEF win32}
function mpfToStr(f: mpf_t): string;
var
  ex:  mp_exp_t;
  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 mpfMyRoot(rad: mpf_t; wzlExp: int64): extended;
var
  ex: mp_exp_t;
begin
  result:=power(mpf_get_d_2exp(ex,rad),1/wzlExp);
  result:=result*power(2,ex/wzlExp);
end;
{$ENDIF}
{$ENDIF}

function signSqr(x: extended): extended;
begin
 result:=sign(x)*sqr(x);
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 intToFixpoint(x,stellen: longint): string;
var
  neg: boolean;
begin
  neg:=x<0;
  if neg then
    x:=-x;
  result:=intToStr(x);
  while length(result)+byte(neg)<stellen do
    result:='0'+result;
  if neg then
    result:='-'+result;
end;

function floatToFixpoint(x: extended; vk,nk: longint): string;
begin
  result:=intToFixpoint(round(x*power(10,nk)),vk+nk);
  result:=leftStr(result,length(result)-nk)+'.'+rightStr(result,nk);
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..10] of string = ('exp','sin','cos','tan','sqr','sqrt','ln','round','floor','ceil','abs');
  fkt2: array[0..2] of string = ('min','max','mod');
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);
        7: val1:=round(val1);
        8: val1:=floor(val1);
        9: val1:=ceil(val1);
        10: val1:=abs(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);
        2:
          begin
            if val2<0 then
              val2:=-val2;
            assert(val2>0,'Der Modul darf nicht 0 sein.');
            while val1<0 do
              val1:=val1+val2;
            while val1>=val2 do
              val1:=val1-val2;
          end;
      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 assigned(cbgv) then
    while pos('[',s)>0 do begin
      i:=pos('[',s)-1;      // i wird Index des letzten Zeichens vor "Werte[bla].blubb" o. ä.
      while (i>0) and not (s[i] in ['+','-','*','/','^']) do
        dec(i);
      j:=pos(']',s)+1;      // j wird Index des ersten Zeichens nach "Werte[bla].blubb" o. ä.
      while (j<=length(s)) and not (s[j] in ['+','-','*','/','^']) do
        inc(j);

      s:=copy(s,1,i)+                         // vorher
         floattostr(cbgv(copy(s,i+1,j-i-1)))+ // zwischen
         copy(s,j,length(s)-j+1);             // nachher
    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 (assigned(kvs) and kvs.extract(s,result)) then begin
    if st then result:=1
    else if assigned(cbgv) then result:=cbgv(s)
    else begin
      result:=nan;
      fehler(''''+s+''' kann ich nicht auflösen!');
    end;
  end;

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

function simpleExprToBool(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): boolean;
var
  i: longint;
  t: string;
const
  binOps: array[0..12] of string =
    ('<=','>=','<>','≤','≥','=','≠','<','>',' in ','∈',' notIn ','∉');
begin
  s:=trim(s);
  while (leftStr(s,1)='(') and (rightStr(s,1)=')') do
    s:=trim(copy(s,2,length(s)-2));

  if uppercase(s)='TRUE' then begin
    result:=true;
    exit;
  end;
  if uppercase(s)='FALSE' then begin
    result:=false;
    exit;
  end;

  for i:=0 to length(binOps)-1 do
    if pos(binOps[i],s)>0 then begin
      t:=erstesArgument(s,binOps[i]);
      case binOps[i] of
        '≤','<=': result:=exprToFloat(st,t,kvs,cbgv)<=exprToFloat(st,s,kvs,cbgv);
        '≥','>=': result:=exprToFloat(st,t,kvs,cbgv)>=exprToFloat(st,s,kvs,cbgv);
        '=':
          try
            result:=exprToFloat(st,t,kvs,cbgv)=exprToFloat(st,s,kvs,cbgv);
          except
            result:=t=s;
          end;
        '≠','<>':
          try
            result:=exprToFloat(st,t,kvs,cbgv)<>exprToFloat(st,s,kvs,cbgv);
          except
            result:=t<>s;
          end;
        '<': result:=exprToFloat(st,t,kvs,cbgv)<exprToFloat(st,s,kvs,cbgv);
        '>': result:=exprToFloat(st,t,kvs,cbgv)>exprToFloat(st,s,kvs,cbgv);
        ' in ','∈': begin
          result:=false;
          while (s<>'') and not result do
            result:=erstesArgument(s)=t;
        end;
        ' notIn ','∉': begin
          result:=true;
          while (s<>'') and result do
            result:=erstesArgument(s)<>t;
        end;
        else
          fehler('Operator '''+binOps[i]+''' ist nicht implementiert!');
      end{of case};
      exit;
    end;
  fehler(''''+s+''' ist kein gültiger logischer Ausdruck!');
end;

function exprToBool(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): boolean;
var
  i,j,anz:       longint;
  operatoren,
  ausdruecke:    tIntPointArray;
  werte:         tBooleanArray;
  klammernZahl:  tLongintArray;
const
  logOpsNamen: array[0..6] of string =
    ('not ',' and ','&&',' or ','||',' xor ',' equiv ');
  logOpsAktionen: array[0..6,boolean,boolean] of boolean =
    (((true,false),(true,false)),  // not (1. Argument ignoriert!)
     ((false,false),(false,true)), // and
     ((false,false),(false,true)), // &&
     ((false,true),(true,true)),   // or
     ((false,true),(true,true)),   // ||
     ((false,true),(true,false)),  // xor
     ((true,false),(false,true))); // equiv

begin
  i:=1;
  setLength(operatoren,0);
  while i<=length(s) do begin
    for j:=0 to length(logOpsNamen)-1 do
      if copy(s,i,length(logOpsNamen[j]))=logOpsNamen[j] then begin
        setLength(operatoren,length(operatoren)+1);
        operatoren[length(operatoren)-1,'x']:=j;
        operatoren[length(operatoren)-1,'y']:=i;
        i:=i+length(logOpsNamen[j]);
        break;
      end;
    inc(i);
  end;

  setLength(ausdruecke,length(operatoren)+1);
  setLength(klammernZahl,length(operatoren)+1);

  for i:=0 to length(ausdruecke)-1 do begin
    if i=0 then
      ausdruecke[i,'x']:=1
    else
      ausdruecke[i,'x']:=operatoren[i-1,'y']+length(logOpsNamen[operatoren[i-1,'x']]);
    if i=length(operatoren) then
      ausdruecke[i,'y']:=length(s)
    else
      ausdruecke[i,'y']:=operatoren[i,'y']-1;

    anz:=0;
    for j:=ausdruecke[i,'x'] to ausdruecke[i,'y'] do
      if s[j]='(' then
        inc(anz)
      else if s[j]=')' then
        dec(anz);

    klammernZahl[i]:=anz;

    while (anz<0) and (ausdruecke[i,'y']>=ausdruecke[i,'x']) do begin // zu viele schließende Klammern (am rechten Rand)
      if s[ausdruecke[i,'y']]=')' then
        inc(anz);
      dec(ausdruecke[i,'y']);
    end;
    while (anz>0) and (ausdruecke[i,'y']>=ausdruecke[i,'x']) do begin // zu viele öffnende Klammern (am linken Rand)
      if s[ausdruecke[i,'x']]='(' then
        dec(anz);
      inc(ausdruecke[i,'x']);
    end;
    if (ausdruecke[i,'y']<ausdruecke[i,'x']) and ((i=length(operatoren)) or (operatoren[i,'x']<>0)) then // das linke Argument von "not" darf leer sein
      fehler('Klammerfehler in '''+s+''' ('+intToStr(ausdruecke[i,'y'])+'<'+intToStr(ausdruecke[i,'x'])+')!');
  end;

  anz:=0;
  for i:=0 to length(klammernZahl)-1 do begin
    anz:=anz+klammernZahl[i];
    if anz<0 then
      fehler('zu viele oder zu zeitige schließende Klammern in '''+s+'''!');
  end;
  if anz<>0 then
    fehler('Klammern in '''+s+''' sind nicht ausgeglichen!');

  setLength(werte,length(ausdruecke));
  for i:=0 to length(werte)-1 do
    if (i=length(werte)-1) or (operatoren[i,'x']<>0) then // kein linkes Argument von "not"
      werte[i]:=simpleExprToBool(st,copy(s,ausdruecke[i,'x'],ausdruecke[i,'y']-ausdruecke[i,'x']+1),kvs,cbgv)
    else
      werte[i]:=false;
  setLength(ausdruecke,0);
  i:=length(werte)-2;
  while i>=0 do
    if ((klammernZahl[i]>0) or (i=0)) and (klammernZahl[i+1]<=0) then begin
      result:=true;
      werte[i]:=logOpsAktionen[operatoren[i,'x'],werte[i],werte[i+1]];
      klammernZahl[i]:=klammernZahl[i]+klammernZahl[i+1];
      for j:=i+1 to length(operatoren)-1 do
        operatoren[j-1]:=operatoren[j];
      for j:=i+2 to length(werte)-1 do begin
        werte[j-1]:=werte[j];
        klammernZahl[j-1]:=klammernZahl[j];
      end;
      setLength(operatoren,length(operatoren)-1);
      setLength(werte,length(werte)-1);
      setLength(klammernZahl,length(klammernZahl)-1);
      if i>length(werte)-2 then
        dec(i);
    end
    else
      dec(i);
  if length(werte)<>1 then
    fehler('Kann kein Paar von öffnenden und schließenden Klammern in '''+s+''' mehr finden - das sollte eigentlich nicht passieren können!');

  result:=werte[0];
end;

function formelnAuswerten(s: string; kvs: tKnownValues; cbgv: tCallBackGetValue): string;
var
  i,start,mitte,
    vk,nk:       longint;
  t:             string;
  ok:            boolean;
  wert:          extended;
begin
  ok:=false;
  i:=1;
  mitte:=0;
  start:=0;
  while i<=length(s) do begin
    if copy(s,i,2)='$[' then begin
      ok:=true;
      i:=i+2;
      start:=i;
      continue;
    end;
    if ok then
      case s[i] of
        ';':
          mitte:=i;
        ']': begin
          if mitte<=start then
            mitte:=i;
          t:=trim(copy(s,start,mitte-start));
          try
            wert:=exprToFloat(false,t,kvs,cbgv);
          except
            ok:=false;
            inc(i);
            continue;
          end;
          if mitte<>i then begin
            t:=trim(copy(s,mitte+1,i-mitte-1));
            if pos('.',t)=0 then
              t:=t+'.0';
            try
              vk:=strToInt(erstesArgument(t,'.'));
              nk:=strToInt(t);
            except
              ok:=false;
              inc(i);
              continue;
            end;
            if nk=0 then
              t:=intToFixpoint(round(wert),vk)
            else
              t:=floatToFixpoint(wert,vk,nk);
          end
          else // keine Formatierung => generisch als float umwandeln
            t:=floattostr(wert);

          s:=
            leftStr(s,start-3) + t + rightStr(s,length(s)-i);
          i:=start-3 + length(t);
          continue;
        end;
        '$': ok:=false;
      end{of case};
    inc(i);
  end;
  result:=s;
end;

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

function berechneEinheitsZelle(invarianz,modulus: tInt64Point): tInt64Point;
var
  c:    char;
  fak:  tInt64Point;
  verh: tExtPoint;
begin
  for c:='x' to 'y' do begin
    invarianz[c]:=invarianz[c] mod modulus[c];
    result[c]:=ggT(invarianz[c],modulus[c]);
    invarianz[c]:=invarianz[c] div result[c];
    modulus[c]:=modulus[c] div result[c];
  end;

  // x*y = p (konstant)
  // für x>y wollen wir x/y maximieren (lang und schmal)
  // => x^2 ist zu maximieren, falls x>y

  fak['x']:=ggT(modulus['y']*invarianz['x'],modulus['x']);
  fak['y']:=ggT(modulus['x']*invarianz['y'],modulus['y']);

  for c:='x' to 'y' do begin
    verh[c]:=result[c]*fak[c]/result[char(ord('y')-byte(c='y'))];
    if verh[c]<1 then
      verh[c]:=1/verh[c];
  end;

  if verh['x']<verh['y'] then
    result['y']:=result['y']*fak['y']
  else
    result['x']:=result['x']*fak['x'];
end;

function ggT(a,b: int64): int64;
var
  tmp: int64;
begin  // a <= b
  if a>b then begin
    tmp:=abs(a);
    a:=abs(b);
    b:=tmp;
  end
  else begin
    a:=abs(a);
    b:=abs(b);
  end;
  while a>0 do begin
    tmp:=b mod a;
    if tmp>a then
      b:=tmp
    else begin
      b:=a;
      a:=tmp;
    end;
  end;
  result:=b;
end;

function intRoot(n: int64): int64;
begin
  if n<=0 then begin
    result:=0;
    exit;
  end;
  result:=floor(sqrt(n));
  while n mod result > 0 do
    dec(result);
end;

function ermittleAnstieg(dat: string; xSpalte,ySpalte: longestOrdinal): extended;
var
  f:                       textFile;
  s:                       string;
  xSum,xQdrSum,ySum,xySum: extended;
  wert:                    tExtPoint;
  i,j,anz:                 longestOrdinal;
  reihenfolge:             string[2];
  positionen:              array[0..2] of longestOrdinal;
begin
  if xSpalte=ySpalte then begin
    result:=1;
    exit;
  end;

  positionen[0]:=-1;
  positionen[1]:=min(xSpalte,ySpalte);
  positionen[2]:=max(xSpalte,ySpalte);

  if xSpalte<ySpalte then
    reihenfolge:='xy'
  else
    reihenfolge:='yx';

  xSum:=0;
  xQdrSum:=0;
  ySum:=0;
  xySum:=0;
  anz:=0;

  assignFile(f,dat);
  reset(f);
  while not eof(f) do begin
    readln(f,s);
    for j:=1 to 2 do begin
      for i:=positionen[j-1]+1 to positionen[j]-1 do
        erstesArgument(s);
      wert[reihenfolge[j]]:=
        strToFloat(erstesArgument(s));
    end;
    xSum:=xSum+wert['x'];
    ySum:=ySum+wert['y'];
    xySum:=xySum+wert['x']*wert['y'];
    xQdrSum:=xQdrSum+sqr(wert['x']);
    inc(anz);
  end;
  closeFile(f);
  if anz=0 then
    fehler('Keine Werte in Anstiegs-Datei'''+dat+'''!');
  if anz*xQdrSum=sqr(xSum) then
    fehler('Anstieg aus '''+dat+''' nicht bestimmbar!');

  result:=(anz*xySum-xSum*ySum)/(anz*xQdrSum-sqr(xSum));
end;

function ermittleMittelwert(dat: string; werteSpalte,gewichteSpalte: longestOrdinal): extended;
var
  f:           textFile;
  s:           string;
  wSum,gSum:   extended;
  wert:        tExtPoint;
  i,j:         longestOrdinal;
  reihenfolge: string[2];
  positionen:  array[0..2] of longestOrdinal;
begin
  if werteSpalte=gewichteSpalte then
    fehler('Werte- und Gewichtespalte dürfen nicht übereinstimmen!');

  wSum:=0;
  gSum:=0;
  assignFile(f,dat);
  reset(f);
  if gewichteSpalte < 0 then begin
    while not eof(f) do begin
      readln(f,s);
      for i:=1 to werteSpalte-1 do
        erstesArgument(s);
      wSum:=wSum+strToFloat(erstesArgument(s));
      gSum:=gSum+1;
    end;
  end
  else begin
    positionen[0]:=-1;
    positionen[1]:=min(werteSpalte,gewichteSpalte);
    positionen[2]:=max(werteSpalte,gewichteSpalte);
    if werteSpalte<gewichteSpalte then
      reihenfolge:='xy'
    else
      reihenfolge:='yx';
    while not eof(f) do begin
      readln(f,s);
      for j:=1 to 2 do begin
        for i:=positionen[j-1]+1 to positionen[j]-1 do
          erstesArgument(s);
        wert[reihenfolge[j]]:=
          strToFloat(erstesArgument(s));
      end;
      wSum:=wSum+wert['x']*wert['y'];
      gSum:=gSum+wert['y'];
    end;
  end;
  closeFile(f);
  if gSum=0 then
    fehler('Keine Werte in Mittelwert-Datei'''+dat+'''!');
  result:=wSum/gSum;
end;

procedure sortiereNachWert(var maxima: tLongintArray; werte: tExtendedArray; logarithmischesPivot: boolean);
var
  mins,maxs:                tExtendedArray;
  pivot,wert,wertLi,wertRe: extended;
  vons,biss:                tLongintArray;
  i,li,re,cnt,tmp:          int64;
begin
  setLength(vons,1);
  vons[0]:=0;
  setLength(biss,1);
  biss[0]:=length(maxima)-1;
  setLength(mins,1);
  setLength(maxs,1);
  mins[0]:=0;
  maxs[0]:=0;
  for i:=vons[0] to biss[0] do begin
    wert:=werte[maxima[i]];
    if (i=vons[0]) or (wert>maxs[0]) then
      maxs[0]:=wert;
    if (i=vons[0]) or (wert<mins[0]) then
      mins[0]:=wert;
  end;
  cnt:=1;

  while cnt>0 do begin
    li:=vons[cnt-1];
    re:=biss[cnt-1];

    if (li>=re) or (mins[cnt-1]=maxs[cnt-1]) then begin
      dec(cnt);
      continue;
    end;

    if cnt>=length(vons) then begin
      setLength(vons,cnt+100);
      setLength(biss,cnt+100);
      setLength(mins,cnt+100);
      setLength(maxs,cnt+100);
    end;

    if logarithmischesPivot and (maxs[cnt-1]*mins[cnt-1]>0) then begin
      pivot:=sqrt(maxs[cnt-1]*mins[cnt-1]);
      if (pivot<=mins[cnt-1]) or (pivot>=maxs[cnt-1]) then
        pivot:=(maxs[cnt-1]+mins[cnt-1])/2;
    end
    else
      pivot:=(maxs[cnt-1]+mins[cnt-1])/2;
    if pivot>=maxs[cnt-1] then
      pivot:=mins[cnt-1];

    mins[cnt]:=mins[cnt-1];
    biss[cnt]:=biss[cnt-1];
    maxs[cnt]:=mins[cnt];
    mins[cnt-1]:=maxs[cnt-1];

    while li<=re do begin
      wertLi:=werte[maxima[li]];
      if wertLi>pivot then begin
        if wertLi<mins[cnt-1] then
          mins[cnt-1]:=wertLi;
        inc(li);
        continue;
      end;
      wertRe:=werte[maxima[re]];
      if wertRe<=pivot then begin
        if wertRe>maxs[cnt] then
          maxs[cnt]:=wertRe;
        dec(re);
        continue;
      end;
      if wertLi>maxs[cnt] then
        maxs[cnt]:=wertLi;
      if wertRe<mins[cnt-1] then
        mins[cnt-1]:=wertRe;
      tmp:=maxima[re];
      maxima[re]:=maxima[li];
      maxima[li]:=tmp;
      inc(li);
      dec(re);
    end;

    vons[cnt]:=li;
    biss[cnt-1]:=re;

    inc(cnt);
  end;

  for i:=1 to length(maxima)-1 do
    if werte[maxima[i-1]]<werte[maxima[i]] then
      fehler('Interner Fehler: Werte sind nicht sortiert, sollten sie aber sein!');
end;

procedure sortiereNachDominanz(var maxima: tLongintArray; werte: tExtendedArray; nullteSpektrometerordnung: extended);
var
  i,j,dom,schranke: int64;
  dominanzen:       tExtendedArray;
  tmp:              tLongintArray;
begin
  // zuerst sortieren wir die Maxima nach Höhe
  sortiereNachWert(maxima,werte,true);
  // dann berechnen wir die Dominanzen aller Maxima
  setLength(dominanzen,length(maxima));
  for i:=0 to length(dominanzen)-1 do begin
    schranke:=length(werte);
    for j:=0 to i-1 do
      schranke:=min(schranke,abs(maxima[i]-maxima[j]));
    dom:=1;
    while (dom<schranke) and
      ((maxima[i]+dom>=length(werte)) or (werte[maxima[i]+dom]<=werte[maxima[i]])) and
      ((maxima[i]-dom<0) or (werte[maxima[i]-dom]<=werte[maxima[i]])) do
      inc(dom);
    dominanzen[i]:=dom;
    if (maxima[i]+dom < length(werte)) and (werte[maxima[i]]<werte[maxima[i]+dom]) then // der rechte Rand hat getriggert
      dominanzen[i]:=dom - 1 + 1 / (1 - (werte[maxima[i]+dom]-werte[maxima[i]])/(werte[maxima[i]+dom-1]-werte[maxima[i]]));
    if (maxima[i]-dom >= 0) and (werte[maxima[i]]<werte[maxima[i]-dom]) then // der linke Rand hat getriggert
      dominanzen[i]:=min(dominanzen[i],dom - 1 + 1 / (1 - (werte[maxima[i]-dom]-werte[maxima[i]])/(werte[maxima[i]-dom+1]-werte[maxima[i]])));
    if not isNaN(nullteSpektrometerordnung) then
      dominanzen[i]:=dominanzen[i] / sqr(i - nullteSpektrometerordnung);
  end;

  // und sortieren nach Dominanz
  setLength(tmp,length(maxima));
  for i:=0 to length(tmp)-1 do begin
    tmp[i]:=maxima[i];
    maxima[i]:=i;
  end;
  sortiereNachWert(maxima,dominanzen,false);

  for i:=0 to length(tmp)-1 do
    maxima[i]:=tmp[maxima[i]];
  setLength(tmp,0);
  setLength(dominanzen,0);
end;

function RotationsMatrix(alpha: extended): t2x2Extended;
begin
  result:=_2x2Extended(
    cos(alpha),-sin(alpha),
    sin(alpha), cos(alpha)
  )
end;

function einheitsVektor(alpha: extended): tExtPoint;
begin
  result:=extPoint(
    cos(alpha),sin(alpha)
  );
end;

end.