Gyors logaritmus fixpontos aritmetikával

Fórumok

Sziasztok!

Egész számoknak keresem a logaritmusát és mivel kernelben működne a program ezért nem használhatok lebegőpontos számokat. Van pl. egy 1 és 100e9 közötti tartomány, ezt szeretném logaritmussal leképezni egy kettő hatvány méretű területre. Teszem azt 10 bitem van, legnagyobb értékem az 1024, ez lenne a 100e9 és a hozzá közel eső számok. Az ehhez tartozó logaritmus a 1.025 alapú, mert ha x^1024 = 100e9 akkor x ~ 1.025. Mivel x bármilyen értéket felvehet 1 és 100e9 között én az ehhez legközelebbi logaritmus értéket keresem.

A jelenlegi megoldásom úgy néz ki, hogy előre kigenerálom az értékeket egy táblázatba, az első ezer számra ez így néz ki:
x : log_1.025(x)
1   : 0
2   : 28
3   : 44
4   : 56
....
912-934: 276
935-957: 277
958-981: 278
982-1000: 279

És ebben a táblázatban max 10 lépésben csinálok bináris keresést x-re és a kapott index lesz a logaritmus.

A kérdésem: lehetne-e ezt elegánsabban? Tudnátok-e mondani valami olyan módszert, amivel esetleg megspórolom a táblázatot meg a bináris keresést és képes lennék "egy lépésben" mondani x-hez egy log_1.025(x)-et? Szerencsére minden ismert fordítási időben x-et leszámítva, tehát megvan mi a legkisebb és legnagyobb x, megvan hogy hány bitre kell leképezni és milyen alapú logaritmussal.

Előre is köszönöm az ötleteket.

Hozzászólások

Szerkesztve: 2021. 07. 20., k – 10:37

Első megközelítésben a log2, az a bitek száma, amennyin ábrázolod a számodat binárisan, azaz a legbaloldalibb 1-es pozíciója. (A legbaloldalibb 1-es pozíciójára egyes processzorokon van 1 lépéses művelet is, de nem bonyolult ciklussal sem megcsinálni. Sok bites ábrázolás esetén ezt a keresést is meg lehet csinálni felezgetősre, hogy ne kelljen pl 64-szer iterálni 64 bites számokra.)

Ezután normalizálod a számod ábrázolását (az eredeti szám osztása elméletileg, de megvalósításban csak shiftelés, az osztás, az logaritmusban nézve kivonás ugye) mondjuk fix pontos ábrázolásban az 1-2 tartományra, és ezen belül szerintem már táblázatot kell használni, jobb ötletem nincsen. Az első bited mindig 1, az ez után jövőkből szintén bitműveletekkel lehet csinálni egy kereső-indexként használható számot (shift+maszkolás és művelettel) a táblázathoz.

A táblázatban a logaritmus értékeket fix pontos ábrázolásban tárolod, és a logaritmus értéknek a kettedes pont utáni része lesz az, ami a táblázatból jön.

 

(Mivel a log2(x) függvény értékeit keressük, ennek a függvénynek pedig az 1-2 szakaszon nem túlságosan változik a meredeksége, ezért a táblázatos lekérdezése az értéknek kellően pontos tud lenni. Lehet becsülni, hogy mennyit fogsz tévedni, és a kívánt pontosság szerinti bit számot lehet alkalmazni a táblázat címzésére. Emiatt, hogy ezen a szakaszon nem változik nagyon a meredekség, emiatt lehet kispórolni a bináris keresést szerintem.)

(fix pontos ábrázolás: sima int-ben van a számod, amibe oda képzelsz egy kettedes pontot valahova, mindjuk a 8. bitre. Ekkor minden számot úgy kell érteni, x/256 van a változóban, nem x. De ezt csak te tudod, a program sima int-ekkel számol.)

Ami kijött a táblázatból ahhoz hozzáadod az eredeti becslést (megint csak shift művelettel a helyére tolva a fix pontos ábrázolásban), ami a logartimus kettedes pont előtti része ugye (a logaritmusban hozzáadás, az ugye az eredetiben szorzás) és kész vagy.

Egy nagyságrendileg 1024 elemű táblázat elfér RAM-ban ha nem mikrovezérlőről beszélünk, én így csinálnám.

Köszönöm, hogy összefoglaltad elsőre nem esett le de most már értem. Ez tök jól néz ki, igazából memória van, ez egyszer lesz csak benne globálisan a memóriában akkor sincs baj ha monduk pár megabájtra meghízik (sebesség kevésbé fontos nekem most, mint a pontosság). És asch-nak is köszönöm a részletes leírást, alapvetően jónak tűnik ez az algoritmus, talán optimalizáltabb is mint amire szükségem van :-)

A 100e9 az 1e9-et jelent, vagy 1e11-et? Az előbbi 30 biten, az utóbbi 37 biten fér el, tehát 20 illetve 27 bites jobbraléptetés jutna eszembe először.

Ha van 1024 lehetőséged a logaritmusra, akkor ennyiszer 5 byteban elfér a hozzá tartozó 0...100e9 döntés határ érték.
Ebben 10 lépéses bináris keresésből megvagy, hogy melyik értéknek van éppen alatta, azaz mely logaritmusérték kell neked.

A log(1+x) sorfejtese egyszeru, es ha megvan a tablazat egy y ertekre, akkor a log(y+x)=log(y)+log(1+x/y) -t erdemes szamolni. Ha ketto hatvany az y, akkor az osztas is egyszeru. 

log(1+x) = x  -x*x/2+  x*x*x/3 ... -(-x)^n /n 

Szerencsére ez nem mikrokontroller, rendes PC csak nem lehet lebegőpontos számokat használni. De igen, alapvetően ezért is írtam ide, mert reméltem hogy vannak erre civilizált megoldások viszont a libfixmath-ben én nem látok csak egy log2 függvényt ami alapból is van a kernelben.

Ez jo, mert pl. egy 8 bites Arduinonak nem okoz gondot a 32 bites float. Lehet, hogy nem lesz olyan gyors, mint a hardware-es, de a fordito osszerakja minden problema nelkul. Persze azt te tudod, mennyire szamit neked a sebesseg. (logaritmus mondjuk pont nincs benne alapbol, hatvany, szogfuggvenyek es hasonlok vannak)

Amugy ugy emlekszem, az AGC is Taylor-sorral szamolta ezeket, szoval az eleg lehet neked is.

When you tear out a man's tongue, you are not proving him a liar, you're only telling the world that you fear what he might say. -George R.R. Martin

Nem Taylor sor, van gyorsabb is:

float _fast_log2(const float val) {
    float mp = 0.346607f; 
    float result = (float)*((int*)&val); 
    result *= 1.0/(1<<23); 
    result = result - 127;   
    float tmp = result - floorf(result); 
    tmp = (tmp - tmp*tmp) * mp; 
    return tmp + result;
}

Másik logaritmus alapra áttérés: az eredményt egyszerűen szorozd meg float-os konstanssal.

Ez mintha a float definiciojat hasznalna ki, mint a quake gyors 1/sqrt-je. A kitevot kinyeri az utolso nehany bitbol, es az egesz resz logaritmusanak hasznalja. A maradek reszre pedig kvadratikusan kozelit. A tmp-tmp*tmp masodfoiug epp az ln(1+tmp) Taylor sora. Ez van leosztva ln(2)/2-vel, hogy a log_2(1+tmp)-t kozelitse.

A kernelben miért ne lehetne lebegőpontos számokat használni?

Nem titkos, egy korábban publikált algoritmusunkat próbálom könnyen használható formára hozni hogy bárki akit érdekel ki tudja próbálni egy laptopon vagy több fizikai/virtuális géppel is. Egyik komponense egy linux Qdisc, ez elég egyszerű, a másik komponens egy csomag jelölő, ami a kérdésemben részletezett módon fogja és a sebességet leképezi értékekre amit berak a csomagok fejlécébe. Ennek van egy DPDK-ban működő verziója, ott lehet lebegőpontos számokat is használni, de sajnos sok más felesleges dolog is van benne és belsős céges kód.

Először AF_XDP-vel kezdtem el a mostanit, hasonló a DPDK-hoz. Viszont ez macerás volt, a kód 95%a boilerplate, ráadásul vagy támogatja a kernel/hálókártya a zero-copy-t vagy nem így szűkül a kompatibilis környezetek tere. A linux tc eBPF része pedig jó ideje stabil, meg szerintem elegáns is itt csomagokat markolni anélkül hogy járatnád a csomagot kernel és userspace között. Ami miatt ez jó még, hogy vannak mapek, amiket pl. bpftool-al akár bash command lineból fel tudsz tölteni. Egy ilyen map-et lát a userspace és a kernel spaceben futó eBPF progi is így ha én hirtelen váltanék a csomag-markolási stratégia között, akkor kényelmesen beállítom másra az 1024 értéket, lockolás, kommunikáció stb. kódolásával nem kell foglalkozni meg aki esetleg letölti a kódom annak lesz esélye a boilerplate között látni valamit a működésből.

Elvileg van még nfqueue, meg AF_PACKET de ezeknek szintén vannak performancia a problémák, ami ha valaki laptopon namespacek között futtatja a kódot akkor nem feltétlen gond, de ha fizikai gépek között ott már lehetnek bajok vele. Ráadásul mindkettő elég régi API, és felesleges oda-vissza másolást csinálnak.

Bármi egyéb ötletre nyitott vagyok, de beüzemelési komplexitásnak egy sudo apt install, make, meg 3-4 parancs amit el tudok képzelni. Efölött olyanokat riaszthat el akik nem értenek a linuxhoz de magát a kódot kipróbálnák. Sajnos más alternatív userpsace stackek beüzemelése nem egyszerű, már egy DPDK-é sem feltétlen az.

Szerkesztve: 2021. 07. 23., p – 05:31

Csak mert előjött: Linux kernel alatt nem igazán megengedett a lebegőpontos aritmerika használata, a leggyakoribb ok, amit felhoznak, az az, hogy az FPU kontextusának mentése költséges feladat, nem akarjuk minden rendszerhívásnál megtenni, bőven elég task váltásakor. Meg lehet csinálni (emlékeim szerint valamelyik digital rendering driver csinál ilyet), de okkal ritka.

A sokadik driverem megírását magam mögött tudva (disclaimer: kb. 2009 óta csinálom, ma ebből élek, volt minden, ha van érdeklődés rá, talán blogolhatnék róla) az a tapasztalatom, hogy a kernelen belül kábé minden utility függvényre van jó (értsd: használható, általában használt, "elég gyors") implementáció, nem érdemes sajátot behozni. Ha kétségeid vannak a log2 függvénnyel, mérd ki. Integer aritmetikára abszolút igaz.

Egyébként bármi kétség esetén érdemes még egy jó cross referencet megnézni, pl elixir.bootlin.com . Látni fogod, hogyan csinálják mások. Enélkül szerintem ma már nem nagyon lehet könnyen tanulni.

 

Szerk: https://elixir.bootlin.com/linux/latest/A/ident/double , tényleg GPUnál, de valami BPF sample is van.

Értem, köszönöm az infót, nem tudtam hogy ezért nem preferálják a lebegőpontos számokat.

kábé minden utility függvényre van jó [..] implementáció, nem érdemes sajátot behozni

Nem szeretnék semmiképp sajátot behozni, a log2-vel sincs bajom csak nem elég pontos ahhoz amire nekem kell, ezért indítottam ezt a threadet, hogy mások milyen trükkökkel csinálnak logaritmust egész számokkal (nem feltétlen a kernelben, hanem általában). És jött is szerencsére sok tanács sőt az is kiderült nekem, hogy lehet a log2-t felhasználni mégis a saját problémámhoz. Sajnos ami BPF/XDP sample-k vannak azok minden esetben a user-space kódok, BPF restriktívebb és nem engedi a lebegőpontos aritmetikát. Ettől függetlenül levlistákon volt szó róla, hogy lehetne de a mai napig nem csinálta meg senki.

ha van érdeklődés rá, talán blogolhatnék róla

Azért nagyon hálás lennék! Akkor is ha épp nem tolonganak az emberek, meg fogja találni a közönségét itt (engem mindenképpen).

Mondjuk igazad van abban, a shift a logaritmusra jól használható. A [0..1024] értékekkészletért a végén osztunk: (longlong >> 17)/745
Az osztást kerülni akartam, de a konstanssal való osztás valójában trükkösen kiváltható. Lásd 64 bites ARM esetén a GCC által fordított kódot:

mylog_prec:
        mov     x1, 31405
        lsr     x0, x0, 17
        movk    x1, 0xdae3, lsl 16
        movk    x1, 0x818b, lsl 32
        movk    x1, 0xafef, lsl 48
        umulh   x0, x0, x1
        ubfx    x0, x0, 9, 32
        ret

Tehát végülis ez lesz a jó megoldás.

További optimalizáció, amely főként 32 bites rendszeren sokat segít. A fenti kódrészt 32 bitre lefordítva látszik a rondaság:

mylog_prec:
        lsrs    r0, r0, #17
        movw    r2, #745
        push    {r3, lr}
        orr     r0, r0, r1, lsl #15
        movs    r3, #0
        lsrs    r1, r1, #17
        bl      __aeabi_uldivmod(PLT)    ; hoppá
        pop     {r3, pc}

Ezen segíthetünk. Kasztoljuk át osztás előtt integerre az eredményt, hiszen a maximuma bőven belefér ebbe is.
(unsigned int)(longlong >> 17)/745

mylog_prec:
        lsrs    r0, r0, #17
        movw    r3, #57443
        orr     r0, r0, r1, lsl #15
        movt    r3, 11259
        umull   r3, r0, r3, r0
        lsrs    r0, r0, #7
        bx      lr

Egyébként a 64 bites architektúránál is nyerünk ezzel. Vesd össze az előző ARM64 assembly kóddal:

mylog_prec:
        lsr     x0, x0, 17
        mov     w1, 57443
        movk    w1, 0x2bfb, lsl 16
        umull   x0, w0, w1
        lsr     x0, x0, 39
        ret