matheunit.pas: berechneEinheitsZelle und ggT neu

1 files changed, 59 insertions, 0 deletions

diff --git a/matheunit.pas b/matheunit.pas
index 80c4be7..d186e27 100644
--- a/matheunit.pas
+++ b/matheunit.pas
@@ -50,6 +50,8 @@ function exprToFloat(st: boolean; s: string; kvs: tKnownValues; cbgv: tCallBackG
 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: tIntPoint): tIntPoint;
+function ggT(a,b: int64): int64;
 
 implementation
 
@@ -782,5 +784,62 @@ begin
   result.value:=val;
 end;
 
+function berechneEinheitsZelle(invarianz,modulus: tIntPoint): tIntPoint;
+var
+  c:    char;
+  fak:  tIntPoint;
+  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;
+
 end.