floor Implementierung



  • Hallo,

    hab mich gerade gefragt wie floor in C implementiert ist?
    Denn das naheliegendste, nämlich return (double)((int)d); ist ja eigentlich Schwachsinn weil ja der int Zahlenbereich bei weitem kleiner als der double Zahlenbereich ist.


  • Mod

    Bei den oftmals recht hoch optimierten Mathebibliotheken werden das allgemein irgendwelche fiesen Bithacks sein und kein allgemeiner, portabler Code. Außerdem habe ich gehört, dass es Open Source Software geben soll (unvorstellbar!), die aus der Implementierung ihrer Standardbibliothek kein Geheimnis macht (Blasphemie!). Leider ist so etwas aus dem im ersten Satz genannten Grund in der Regel recht unleserlich.



  • floor Implementierung schrieb:

    Hallo,

    hab mich gerade gefragt wie floor in C implementiert ist?
    Denn das naheliegendste, nämlich return (double)((int)d); ist ja eigentlich Schwachsinn weil ja der int Zahlenbereich bei weitem kleiner als der double Zahlenbereich ist.

    Vermutlich mit Hilfe der modf()-Funktion, die eine Double-Zahl in ihren ganzzahligen und den Nachkommteil trennen kann.



  • Also wird - ohne jetzt ins Detail zu gehen - die Mantisse ab einer bestimmten Stelle mit 0ern überschrieben, sofern der Exponent in einem speziellen Bereich liegt?

    Prinzipiell also nach folgendem Schema:

    2.345*10^20 => überschrieben nicht nötig weil Exponent so groß dass die double Zahl eigentlich eine Ganzzahl ist.

    2.345*10^0 => Mantisse 2345... ab zweiter Stelle mit 0ern überschreiben=>2.0



  • x86 hat dafür eine frndint-Instruktion. Üblicherweise (ich kann das nur für die glibc sicher sagen, aber es sollte mich wundern, wenn das anderswo groß anders gehandhabt würde) wird es darauf hinauslaufen, in der FPU das Runden-in-Richtung-minus-unendlich-Bit zu setzen und frndint hinterherzuschieben. Auf vielen anderen Architekturen sieht es ähnlich aus.

    Da das aber nicht ist, was du hören wolltest - es gibt in der glibc eine Fallback-Implementation für verschiedene ieee754-Floats, ursprünglich anscheinend von Sun entwickelt. Für 64-Bit-Floats (double) liest sich das wie folgt:

    /* @(#)s_floor.c 5.1 93/09/24 */
    /*
     * ====================================================
     * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     *
     * Developed at SunPro, a Sun Microsystems, Inc. business.
     * Permission to use, copy, modify, and distribute this
     * software is freely granted, provided that this notice
     * is preserved.
     * ====================================================
     */
    
    #if defined(LIBM_SCCS) && !defined(lint)
    static char rcsid[] = "$NetBSD: s_floor.c,v 1.8 1995/05/10 20:47:20 jtc Exp $";
    #endif
    
    /*
     * floor(x)
     * Return x rounded toward -inf to integral value
     * Method:
     *        Bit twiddling.
     * Exception:
     *        Inexact flag raised if x not equal to floor(x).
     */
    
    #include "math.h"
    #include "math_private.h"
    
    #ifdef __STDC__
    static const double huge = 1.0e300;
    #else
    static double huge = 1.0e300;
    #endif
    
    #ifdef __STDC__
            double __floor(double x)
    #else
            double __floor(x)
            double x;
    #endif
    {
            int32_t i0,i1,j0;
            u_int32_t i,j;
            EXTRACT_WORDS(i0,i1,x);
            j0 = ((i0>>20)&0x7ff)-0x3ff;
            if(j0<20) {
                if(j0<0) {         /* raise inexact if x != 0 */
                    if(huge+x>0.0) {/* return 0*sign(x) if |x|<1 */
                        if(i0>=0) {i0=i1=0;}
                        else if(((i0&0x7fffffff)|i1)!=0)
                            { i0=0xbff00000;i1=0;}
                    }
                } else {
                    i = (0x000fffff)>>j0;
                    if(((i0&i)|i1)==0) return x; /* x is integral */
                    if(huge+x>0.0) {        /* raise inexact flag */
                        if(i0<0) i0 += (0x00100000)>>j0;
                        i0 &= (~i); i1=0;
                    }
                }
            } else if (j0>51) {
                if(j0==0x400) return x+x;        /* inf or NaN */
                else return x;                /* x is integral */
            } else {
                i = ((u_int32_t)(0xffffffff))>>(j0-20);
                if((i1&i)==0) return x;        /* x is integral */
                if(huge+x>0.0) {                 /* raise inexact flag */
                    if(i0<0) {
                        if(j0==20) i0+=1;
                        else {
                            j = i1+(1<<(52-j0));
                            if(j<i1) i0 +=1 ;         /* got a carry */
                            i1=j;
                        }
                    }
                    i1 &= (~i);
                }
            }
            INSERT_WORDS(x,i0,i1);
            return x;
    }
    weak_alias (__floor, floor)
    #ifdef NO_LONG_DOUBLE
    strong_alias (__floor, __floorl)
    weak_alias (__floor, floorl)
    #endif
    

    Dabei sind EXTRACT_WORDS und INSERT_WORDS Makros, die x in zwei Maschinenworte aufteilen bzw. daraus zusammensetzen. i0 ist hier das niedrige Wort, i1 das höhere, und in Makros verpackt ist das aus Endian-Gründen.

    j0 wird der Exponent sein, und abgesehen von Sonderbehandlungen für unendlich und NaN sowie eine Unterscheidung zwischen +0 und -0 läuft es im Wesentlichen in der Tat darauf hinaus, die überflüssigen Stellen der Mantisse abzusäbeln.



  • Danke dir für Code + Erklärung!



  • Ist ja schlimmer als erwartet 😮

    Leute, die solchen Code schreiben, gehören entlassen.



  • 314159265358979 schrieb:

    Ist ja schlimmer als erwartet 😮

    Leute, die solchen Code schreiben, gehören entlassen.

    Dass man hexadezimal rechnet lass ich mir ja noch einreden, denn teilweise sieht man damit gewisse Dinge wirklich besser als dezimal.

    ABER: Ich werde Leute nie verstehen, die so blöde Variablennamen hernehmen.
    Wieso kann man nicht einfach sprechende Variablen hernehmen???

    Auf der Arbeit musste ich mal ein altes C Programm umbauen. Die meiste Zeit habe ich damit verbracht, zu erraten, was s0, s1, s2, i0, i1, i2 denn nun bedeuten soll. Und wenn mans dann mal weiß muss man es sich aufschreiben weil ich mir die Bedeutung nicht merken kann von so doofen Namen.
    So ... das war jetzt offtopic aber das musste mal raus 😃



  • 314159265358979 schrieb:

    Ist ja schlimmer als erwartet 😮

    Leute, die solchen Code schreiben, gehören entlassen.

    Ach bitte...

    Derart grundlegende Mathefunktionen dürfen so schlimm aussehen, wenn es der Performance zu Gute kommt.



  • Naja, es geht hier um einen Sachverhalt, der extrem low-level ist. Zum einen ist eine sinnvolle Benennung in dem Sinne, wie wir sie üblicherweise verstehen (nämlich eine, die die Rolle der Variable im abstrakten Sachverhalt, den wir gerade modellieren, erklärt) ziemlich schwierig, weil es wenig Abstraktion zu erklären gibt - j0 hätte man exponent nennen können, i0 und i1 vielleicht hiword und loword, i gerade noch mask und bei j wird's dann richtig schwierig.

    Ich denke, man muss das so sehen: Wenn man in IEEE-754 zuhause ist, erkennt man einen Ausdruck der Form

    x = ((y >> kleine_zahl) & 0x[137f]f+) - 0x[137f]f+;
    

    sofort als Extraktion des Exponenten. i0 und i1 als Bezeichner für Teile des selben Wertes sind so schlimm für mein Verständnis eh nicht, und der Rest ist elendes Bitmaskengefrickel, bei dem kein Variablenname hilft. Ich gehe ferner davon aus, dass der Code von jemandem geschrieben wurde, der auf dem Gradienten zwischen Mathematiker und Programmierer irgendwo in der Nähe von Don Knuth steht - und Mathematiker sind Variablennamen gewohnt, die genau einen Buchstaben haben. Wenn im Code "hiword" steht, parst der das zunächst als "h * i * w * o * r * d".



  • Der Code hätte auch besser dokumentiert werden können. Beispielsweise wäre eine "Grafik" des Speicherlayouts von doubles am Anfang ziemlich nett.

    Und nicht nur die Variablennamen sind schrecklich, auch die Klammernsetzung und die fehlenden Abstände. Einrückungen und Leerzeilen wurden auch gespart.



  • 314159265358979 schrieb:

    Der Code hätte auch besser dokumentiert werden können. Beispielsweise wäre eine "Grafik" des Speicherlayouts von doubles am Anfang ziemlich nett.

    Junge, wenn du in so einer Funktion herumschreiben willst, darf davon ausgegangen werden, dass du IEEE-754 auswendig kennst.



  • 314159265358979 schrieb:

    Leute, die solchen Code schreiben, gehören entlassen.

    Ich würd ihn dann für die harten Probleme einstellen.



  • seldon poste bitte mal noch die ein oder andere Funktion aus dieser Bibliothek.



  • DevilMayCry schrieb:

    seldon poste bitte mal noch die ein oder andere Funktion aus dieser Bibliothek.

    http://sourceware.org/git/?p=glibc.git;a=tree;f=sysdeps/ieee754;hb=HEAD

    dbl-64 ist für double, flt-32 für float und die ldbl für diverse long double Varianten.

    Implementierungen optimiert für x86-64 sind hier http://sourceware.org/git/?p=glibc.git;a=tree;f=sysdeps/x86_64/fpu;h=7ed9aacca047946dc03db5fd79078c79623f6c1e;hb=HEAD



  • Du kannst dir die ganze Bibliothek herunterholen. Liegt unter ftp://ftp.gnu.org/gnu/glibc/ - mit langer Versionshistorie; die derzeit neueste ist 2.14.1.

    Die in C geschriebenen ieee-754-Funktionen liegen darin unter sysdeps/ieee754.



  • Man darf bei der GNU Clib auch nicht vergessen dass der Code auch von einem 20 Jahre alten GCC gefressen werden können soll.

    Das hat so einige Auswirkungen auf den Code.


Log in to reply