Nächstgrößere, nächstkleinere Fließkommazahl

[Ekkehard]

Hallo,
ich suche eine FPC-Funktion, die der C-Funktion nextafter entspricht.
Diese Funktion liefert zu einer gegebenen Fließkommazahl die nächstgrößere bzw. nächstkleinere zurück.
Im Prinzip läuft es darauf hinaus, in der Mantisse das unterste Bit zu manipulieren, muss aber dabei natürlich die Überläufe im Exponenten berücksichtigen.

Der Sinn der Sache ist folgender.
Will man eine Ebene lückenlos mit Rechtecken pflastern, so dürfen keine Punkte der Ebene in mehr als einem Rechteck liegen.
Hierzu definiert man das Rechteck so, dass die jeweils rechte (und untere) Seite nicht mehr zum Rechteck gehören.
Beispiel (nur die X-Achse betrachtet):
Rechteck 1: Links: 10.0, Rechts: 16.0
Rechteck 2: Links: 16.0, Rechts: 30.0

Ein Punkt mit der x-Koordinate 16.0 liegt nach dieser Definition in Rechteck 2.
Das lässt sich einfach programmieren:

Code: Alles auswählen

  if (x >= r.Left) and (x < r.Right) then Result := True;

Jetzt aber die Frage: Welches x ist das Größte (also am weitesten rechts liegende), welches gerade noch im Rechteck liegt?
Die Antwort darauf liefert eben genau die gesuchte Funktion:
d := nextafter(16.0,10.0);

Durch geduldiges Spielen im Debugger kommt man auf das Ergebnis: 15,9999999999999980
Und eine kleinere Differenz ist nicht mehr möglich.
Nur hätte ich das gerne berechnet und nicht ausprobiert.

Vielen Dank!

Hier mein Bastelprogramm:
Man braucht ein Form, einen Button und ein Memo

Code: Alles auswählen

procedure TForm1.Button1Click(Sender: TObject);
var
  d1 : Double;  // Die zu manipulierende Double
  b : array[0..SizeOf(d1)-1] of Byte absolute d1; // Die Bits des Double
  d0 : Double; // Die Ausgangs double
  i : Integer;
  s : String;
begin
  s := '';
  d0 := 16.0; // Geht nur mit sehr speziellen Zahlen, bspw. 10.0 oder 20.0!!
  memo1.Lines.Add(Format('%16.16f',[d0])); // 16,0000000000000000
  d1 := d0; // Die Ausgangszahl der manipulierenden Double zuweisen
  // durch die Deklaration "absolute" liegen die Werte auch im array b
  for i := 0 to High(b) do
    s := s + Format(' %2.2x',[b[i]]); // Hex zahlen zusammenstellen
  memo1.Lines.Add(Format('%1.16f als HexDump %s',[d1,s])); // ausgeben
  // Jetzt die Bits manipulieren
  for i := 0 to 5 do
    b[i] := $FF;
  b[6] := b[6] - 1;
  s := '';
  for i := 0 to High(b) do
    s := s + Format(' %2.2x',[b[i]]);
  memo1.Lines.Add(Format('%1.16f als HexDump %s',[d1,s])); // ausgeben
  memo1.Lines.Add(Format('Differenz: %1.6g',[d0-d1])); // -1,77636E-15
end;     

Die Ausgabe des Programms sieht dann so aus:
16,0000000000000000
16,0000000000000000 als HexDump 00 00 00 00 00 00 30 40
15,9999999999999980 als HexDump FF FF FF FF FF FF 2F 40
Differenz: 1,77636E-15

[Mathias]

Code: Alles auswählen

uses
  ctypes;

function nextafter(x, y: cdouble): cdouble;cdecl;external 'm';  

[Ekkehard]

Besten Dank für die Antwort

external 'm';

DLL (Windows) wird nicht gefunden?

Wo soll die sein?

[Mathias]

Versuche mal dies:

Code: Alles auswählen

function nextafter(x,y:cdouble):cdouble;cdecl;external 'msvcrt.dll';   

[af0815]

Naja, man kommt da schnell auf ähnliche Probleme wie hier diskutiert viewtopic.php?f=10&t=15869

Dieses nextafter ist ja genaugenommen vom Zahlentyp abhängig, weil es geht ja darum den entsprechend kleinsten vorigen Wert der möglich ist zu ermitteln. Und das ist halt bei ganzzahligen und den Float Typen unterschiedlich, je nach Genauigkeit.

[Sieben]

Sicher, mathematisch gesehen ist die Fragestellung sinnlos, da zwischen zwei beliebigen Fließkommazahlen, seien sie noch so dicht bei einander, immer unendlich viele weitere Zahlen liegen. Aber dazu hat der Fragesteller ja auch einen Kontext geliefert, wenn wir auch nicht wissen, wozu diese 'Pflasterung' letztlich dienen soll. Vielleicht ist die Genauigkeit ja sogar schon zu hoch?

[af0815]

Er hat schon ein Beispiel geliefert

Der Sinn der Sache ist folgender.
Will man eine Ebene lückenlos mit Rechtecken pflastern, so dürfen keine Punkte der Ebene in mehr als einem Rechteck liegen.
Hierzu definiert man das Rechteck so, dass die jeweils rechte (und untere) Seite nicht mehr zum Rechteck gehören.

Wenn das ganze auf einen Canvas liegt, so ist es klar - ein Bildpunkt weniger Wenn es aber keinen begrenzeten Kontext bzw. bekannten gibt, so müsst es wirklich in der Auflösung der Zahl der nächste kleinstmögliche Schritt zurück sein. Und das ist einmal abhängig von dem Zahlentyp (4,8,10 Byte). Vielleicht gibt es sowas in der unit math. Die ja schon immer einige Überraschungen hat (was da alles drinnen ist). https://www.freepascal.org/docs-html/rt ... dex-5.html

Hier wird auch über sowas diskutiert https://stackoverflow.com/questions/290 ... -to-and-to
Interessanter Ansatz: Die Float als integer betrachten und einfach die Mantisse um eins manipulieren.
Zitat:

By design, for IEEE754 data types, you can simply treat the value as an integer and increment the value. Or decrement it if the value is negative.

[Zvoni]

Selbst zusammenbauen?
https://codebrowser.dev/glibc/glibc/mat ... ter.c.html

[Ekkehard]

Versuche mal dies:

Code: Alles auswählen

function nextafter(x,y:cdouble):cdouble;cdecl;external 'msvcrt.dll';   

Danke, dass hat (fasz auf Anhieb) geklappt.
Die Funktion in der DLL heißt "_nextafter", aber das macht ja nichts, sie funktioniert trotzdem!
Gibt es irgendwo eine Import-Unit für diese DLL? Dann könnte ich die Arbeit sparen, diese dynamisch zu laden.
Ich mag das nicht, wenn ein Programm wegen einer fehlenden/inkompaiblen DLL ohne wirkliche Fehlermeldung einfach beim Start den Dienst quittiert.

[Zvoni]

Ich würds versuchen es selbst zu bauen gemäss dem C-Source-Code, um sich diese externe Abhängigkeit zu sparen.

Das einzige was ich nicht auf Anhieb verstanden habe ist die Funktionsweise von math_opt_barrier und math_force_eval
Vielleicht hat Freepascal (Math-Unit) ja schon was passendes

[Ekkehard]

Sicher, mathematisch gesehen ist die Fragestellung sinnlos, da zwischen zwei beliebigen Fließkommazahlen, seien sie noch so dicht bei einander, immer unendlich viele weitere Zahlen liegen. Aber dazu hat der Fragesteller ja auch einen Kontext geliefert, wenn wir auch nicht wissen, wozu diese 'Pflasterung' letztlich dienen soll. Vielleicht ist die Genauigkeit ja sogar schon zu hoch?

Ich will das Problem noch einmal genauer anhand dessen beschreiben, was ich da tue.
Derzeit baue ich an einer größeren Komponente (eigentlich nur eine Klasse, da nicht von TComponent abgeleitet) die die Funktion eines Quadratur-Suchbaumes übernimmt.
Für mich war es sehr überraschend, dass es so etwas nicht schon als OpenSource gibt.
Der Einsatz ist immer dort sinnvoll, wo viele Elemente auf einer Ebene verstreut liegen und irgendwie ausgewertet werden müssen.
Landkarten sind da ein gutes Beispiel. Nehmen wir an, wir hätten eine Liste aller Haltestellen in Deutschland und würden jetzt alle Haltestellen im Abstand von 500m um einen Standort wissen wollen, dann ist es wenig sinnvoll, alle Haltestellen zu durchforsten um in 99,9999% der Fälle "Außerhalb!" zu erfahren.
Deshalb teilt man die Welt in Rechtecke ein, so dass immer vier Rechtecke zu einem größeren zusammengefasst werden, bzw andersherum ein Rechteck von vier Rechtecken vollständig und nicht überlappend bedeckt wird.
So kann man sehr schnell weite Teile des Suchraums ausschließen, weil man immr nur das nächst kleinere Rechteck mit den Koordinaten wählen muss, bis man am Ziel ist. So bleiben dann statt der Millionen von Koordinaten-Vergleichen, nur einige Dutzend übrig.
Jetzt ist bspw die Erdorberfläche mit den Koordinaten (-180; 90; 180; -90) beschrieben und es gibt auch kein Problem mit den Unterteilungen. Denn man kann 360 problemlos 26 mal durch 2 teilen ohne dabei an eine Begrenzung der Double-Darstellung zu gelangen (ein Meter hat eine Länge von etwa 0,000009 Grad, 26 Teilungen = 26 Bit).
Das sieht aber vielleicht anders aus, wenn jemand die Galaxie kartieren will, oder einen Kristall auf atomarer Ebene. Für diese Fälle möchte ich nicht in Zufälligkeiten geraten sondern klare Fehlermeldungen haben.
All das lässt sich ohne die Fragestellung beantworten, denn dazu gibt es die Vergleichsfunktionen, die auf "kleiner als" prüfen können.

Jetzt bin ich ziemlich am Ende meiner Arbeit und wollte langsam Richtung Veröffentlichung schreiten, da ist mir eine Kleinigkeit aufgefallen, die leider doch ärgerlich ist.
Wie am Anfang des Threads ausgeführt, enthält ein Rechteck nicht den rechten und unteren Rand seiner Definition. Das bedeutet, dass ein Rechteck mit den Koordinaten (0; 0; 0; 0) strenggenommen so punktförmig ist, dass kein Punkt hinein passt. Denn ein Element mit den Koordinaten (0; 0) kann da nicht reinpassen, da der rechte und untere Rand (0; 0) nicht Teil davon ist.
Um es digital auszudrücken (Das Beispiel Canvas war ja genannt) ich suche ein Rechteck, welches genau die eine Double-Koordinate (0;0) aufnehmen kann.
Die Werte dafür lauten: (0; 0; 4,94065645841247E-324; 4,94065645841247E-324)!
In dieses Rechteck passt nur der Punkt 0;0, weil jeder Punkt, der weiter rechts oder drunter liegt schon außerhalb ist (links und oben sowieso).

Ein Rechteck, welches nur einen Punkt umfasst, man könnte denken, das ist nun wirklich kein Problem.
Aber das sieht allgemeiner formuliert anders aus.
Beispiel: Wenn ich jetzt eine Funktion habe, die ein Rechteck berechnen soll, in dem alle Elemente liegen, so kann ich für den rechten und den unteren Rand eben nicht die Koordinaten des am weitesten rechts bzw weitesten unten liegenden Punktes nehmen, denn der ist ja dann per Definition nicht Teil des Rechtecks.
Das umfassende Rechteck muss also ein infinitisimalen Teil (genauer ein Bit Unterschied in der Double Darstellung) größer sein. Und genau diese Größe möchte ich bestimmen.

Jetzt könnte man meinen, dass dieser Ansatz nur in Chaos führen kann. Aber dem ist nicht so. Wenn beispielsweise ein Element am rechten Rand eines der oben genannten Rechtecke des Baumes liegt (größte Zahl die gerade noch kleiner ist als der Wert des rechten Randes) und dieses Element jetzt in die zu umfassende Menge gelangt, so wird das etwas vergrößerte Rechteck genau den Wert erhalten, den der rechte Rand des Rechtecks hat, in dem es ursprünglich eingefügt war.

Die oben angegebene Funktion "nextafter" liefert das genau und so muss ich keinerlei Bitgebastel machen, welches dann doch an irgendeiner Stelle klemmt.

Also Vielen Dank für den Input!

[af0815]

Die oben angegebene Funktion "nextafter" liefert das genau und so muss ich keinerlei Bitgebastel machen, welches dann doch an irgendeiner Stelle klemmt.

Es ist immer ein Bitgebastel, allerdings ausserhalb deiner Reichweite

Denn nur über das Bitgebastel ist es möglich das Problem technisch zu lösen. Alleine die Recherche dazu war schon den Aufwand wert, es zu verstehen Danke für die Problemstellung.

[wp_xyz]

Du brauchst in der Fragestellung mit den Karten sicher keine Mikrometer-Auflösung und darunter. Wäre es nicht einfacher und allgemeiner mit Integern zu arbeiten? Einfach die Float-Koordinaten mit dem entsprechenden Faktor multiplizieren und die Dezimalstellen weglassen?

Die Verwendung von Float-Zahlen hat auch den Nachteil, dass, wenn du so einen Grenz-Punkt hast und nur eine Kleinigkeit mit den Koordinaten rechnest, durch die Rundungsfehler der Punkt plötzlich in einen anderen Quadranten landen kann.

[Mathias]

Versuche mal dies:

Code: Alles auswählen

function nextafter(x,y:cdouble):cdouble;cdecl;external 'msvcrt.dll';   

Danke, dass hat (fasz auf Anhieb) geklappt.
Die Funktion in der DLL heißt "_nextafter", aber das macht ja nichts, sie funktioniert trotzdem!
Gibt es irgendwo eine Import-Unit für diese DLL? Dann könnte ich die Arbeit sparen, diese dynamisch zu laden.
Ich mag das nicht, wenn ein Programm wegen einer fehlenden/inkompaiblen DLL ohne wirkliche Fehlermeldung einfach beim Start den Dienst quittiert.

Jetzt geht es auf beiden Plattformen.

Code: Alles auswählen

 function nextafter(x, y: cdouble): cdouble; cdecl; external {$ifdef linux} 'm' {$else} 'msvcrt.dll' Name '_nextafter' {$endif};    

Das die DLL anders heisst kann ich noch verstehen, aber das Unterline geht mir gegen jede Logik.
Auf diese Idee muss man erst mal kommen.

Auch da bin ich kürzlich über dieses komische Unterline gestossen.
https://learn.microsoft.com/de-de/cpp/c ... w=msvc-170

[Warf]

Die LibC variante wird zwar vermutlich deutlich effizienter implementiert sein, man kann das ganze aber auch recht trivial selbst in Pascal bauen:

Code: Alles auswählen

function nextAfter(const AValue: Single): Single;
var
  Sign, Exp, Mantissa: LongWord;
begin
  // Single precision float: 8 bit exponent 23 bit mantissa
  Sign := LongWord(aValue) and (1 shl 31);
  Exp := LongWord(aValue) and (255 shl 23);
  Mantissa := LongWord(aValue) and ((1 shl 23) - 1);
  if Mantissa + 1 < 1 shl 23 then
    PLongWord(@Result)^ := Sign Or Exp Or (Mantissa + 1)
  else if Exp + 1 < 256 then
    PLongWord(@Result)^ := Sign Or (Exp+1)
  else if Sign = 0 then
    Result := Infinity
  else
    Result := NegInfinity;
end;

function nextAfter(const AValue: Double): Double;
var
  Sign, Exp, Mantissa: QWord;
begin
  // Double precision float: 11 bit exponent 52 bit mantissa
  Sign := QWord(aValue) and (1 shl 63);
  Exp := QWord(aValue) and (2047 shl 52);
  Mantissa := QWord(aValue) and ((1 shl 52) - 1);
  if Mantissa + 1 < 1 shl 52 then
    PQWord(@Result)^ := Sign Or Exp Or (Mantissa + 1)
  else if Exp + 1 < 2048 then
    PQWord(@Result)^ := Sign Or (Exp+1)
  else if Sign = 0 then
    Result := Infinity
  else
    Result := NegInfinity;
end;