A programozó munkája a legkönnyebb...

... legalábbis a többi kolléga szerint, de azért egy nagyon egyszerű probléma megoldására készült nagyon egyszerű programot beidéznék (ami könnyen lehet, hogy még nem is fed le minden esetet): számoljuk ki a a^(p/q) értékét, ahol a valós, p és q egész.

/* mypow.c */

#include <math.h>
#include <stdio.h>

int lnko (int a, int b)
{
    if (a<0) a= -a;
    if (b<0) b= -b;
    if (b>a) { int tmp= a; a= b; b= tmp; }
    while (b!=0) {
        int tmp= a%b;
        a= b;
        b= tmp;
    }
    return a;
}

double mypow (double base, int p, int q)
{
    double retval;
    int d;
/* spec esetek:
   kitevő nevezője nulla: hiba
   kitevő számlálója nulla: a válasz 1.0
 */
    if (q==0) return NAN;
    if (p==0) return 1.0;
/* a továbbiakban p és q nem nulla */

    if (q<0) {
       p= -p;
       q= -q;
    }
/* a továbbiakban p nem nulla, q pozitív */

    d= lnko (p, q);
    if (d!=1) {
       p /= d;
       q /= d;
    }
/* a továbbiakban p nem nulla, q pozitív, relatív prímek */

    if (base<0 && p%2==0) {
        base= -base;
    }
/* a továbbiakban ha base<0 akkor p páratlan */

    if (base<0) { /* negatív alapnál q (is) legyen páratlan */
        if (q%2==0) return NAN;
        retval= - pow (-base, (double)p/(double)q);

    } else if (base==0.0) { /* ha az alap nulla, a kitevő legyen pozitív */
        if (p<0) return NAN;
        retval= 0.0;        /* p==0 eset korábban megvolt */

    } else {
        retval= pow (base, (double)p/(double)q);
    }
    return retval;
}

Hozzászólások

Szerkesztve: 2021. 04. 14., sze – 14:20

És erre miért van szükség? Ez úgy, ahogy van bugos.

if (base==0.0)

Lebegőpontos számra egyenlőségvizsgálat?

Akkora elégtelen, hogy kilóg a NEPTUN-ból...

"Share what you know. Learn what you don't."

Ha esetleg az lenne a helyzet, hogy a nulla pontosan tárolható egy ilyen lebegőpontos változóban, akkor a nullával való összehasonlítás értelmes művelet lehet.

Ettől még, amit a kolléga mond, az egy létező szabályból származik, durván ilyesmi: ha egy xn sorozattal közelítünk egy X értéket, akkor ne xn==X egyenlőségvizsgálattal ellenőrizzük a cél elérését, mert lehet, hogy sosem érünk célba, hanem |xn-X|<ε vizsgálattal, de ha X épp nem nulla, és programozósan akarjuk csinálni, akkor |xn-X|/|X|<ε még jobb lehet, ugyanis ezek a lebegőpontos számok nem egyenletesen oszlanak el (például van olyan, amihez egyet adva önmagát kapjuk).

Ez nem GCC-specifikus, ez a szabvány:https://en.wikipedia.org/wiki/Signed_zero#Comparisons

According to the IEEE 754 standard, negative zero and positive zero should compare as equal with the usual (numerical) comparison operators, like the == operators of C and Java.

    <meta http-equiv="CONTENT-TYPE" content="text/html; charset=windows-1250">
    <meta name="GENERATOR" content="LibreOffice 3.6  (Windows)">
  

Egy programozás tanár miért LibreOffice-szal készíti el a weboldalát? Ráadásul a karakterkódolás sem jó, mert a valódi charset az a windows-1250 de mivel a webszerver http headerben ezt mondja: "Content-Type: text/html; charset=UTF-8" ezért alapból helytelenül jelennek meg az ékezetes karakterek. (jah igen ezt te is írtad):(

> Mi van, ha p 0 és b is 0? Akkor 1-gyel tér vissza, pedig az értéke NaN kellene, hogy legyen.

Akarsz szavazni arról, hogy 0^0=1 vagy sem? Vagy simán közlöd a saját álláspontodat, és azt tekintsük végső igazságnak? https://en.wikipedia.org/wiki/Zero_to_the_power_of_zero

if (p==0) return 1.0;

Azt nem én írtam. Feltűnt a pluszjel a harmadik sor elején?

Ahogy ezt az amúgy minden esetet helyesn lekezelő sima pow() meg is old. Tök felesleges ez a körítés a pow() köré.

Ne nekem mondd.

Szerkesztve: 2021. 04. 14., sze – 20:54

mypow(-4.0, 1, 2) az miért NaN, miért nem -2.0?

Ok leesett, my bad. :)

Viszont a mypow(-4.0, 1, 1)-nek már biztosan -4-et kellene adnia, de arra is NaN-t ad. TCH "optimizációjával" (vagyis bugfixével) jó.

Írtam 290 tesztesetet ami rávilágított még egy potenciális hibára. Ha negatív szám van tört hatványra emelve, akkor az implementáció eredményt ad, miközben ekkor az biztosan nem valós szám. Illetve egy apróság, ha -0 van páratlan hatványra emelve az eredmény +0 lesz.

A tesztek segítségével implementáltam egy saját megoldást (batpow.c) és az összehasonlítás kedvéért a "naiv" megoldást is (stdpow.c).

Ööö... az a gond, hogy az q-adik gyökvonás (illetve (p/q)-adik hatványozás, ahol p és q relatív prím) igazából egy halmazt eredményez, aminek q eleme van [ha az alapszám nem nulla], ezek közül 0, 1 vagy 2 lehet valós. Ez az implementáció azt célozta meg, hogy vagy az egyetlen valós vagy a pozitív valós gyököt adja vissza.

Jahogy mi ennek az alkalmazási területe? Mondjuk az van a gyerek házifeladatában, hogy x^3 = -64

Valljuk be, nem lenne a számítástudomány fénypontja, ha erre én azt mondanám, hogy 'sajnos ezt programmal nem tudom kiszámolni, mert negatív számból nem tudok köbgyököt vonni'

Valóban, ezt nem teljesen gondoltam végig. Amikor ellenőriztem, csak megnéztem WolframAlpha-val, ami első eredménynek egy komplex számot adott. Ha picit lejjebb görgettem volna ott felsorolja az összes gyököt, köztük a valós eredményt is. Illetve gondolhattam volna, hogy ha (-2)^3 = -8, akkor az inverze is értelmezhető, vagyis lehet negatív számnak valós köbgyöke.

Ha negatív szám van tört hatványra emelve, akkor az implementáció eredményt ad, miközben ekkor az biztosan nem valós szám.

Ez definíció kérdése. Ha szigorúan veszed, hogy mi egy függvény definíciója [f:A→B függvény egy f=(A,B, Rel(f)) hármas, ahol Relf(f) regy függvényszerű reláció, azaz Rel(f)⊆A྾B (tehát (x,y) párok halmaza), ahol nem szerepelhet benne (x,y) és (x,z), y≠z], akkor nem lehet a függvény értéke halmaz. Pl. az 1/2-dik hatvány esetén: a szokásos mondat, hogy "nem egy gyöke van" arra utal, hogy nem egy olyan függvényt lehet megadni, amire f(x)^2=x. A komplex függvénytanban szoktak is különböző gyökfüggvényeket bevezetni, és pl. egy olyat, hogy a gyök főága, ami egy a exp(i ɸ) komplex számhoz (a≥, 0≤ɸ<π) hozzárendeli a exp(i ɸ/2)-t.

A fenti gyököknek vágása van: egy a 0-t a végtelennel összekötő görbe mentén szinguláris (ugrása van). Lehet továbbmenni: ha nem a komplex számokon értelmezed a gyököt, hanem egy Riemann-felületen, akkor a függvénynek nem lesz vágása (tehát most a függvény sqrt(z,n) lesz, ami megadja a z gyökei közül az n-ediket (nem az n-edig gyökét!). Az egy külön kérdés, hogy a kérdezőnek mire van szüksége: egy tetszőleges gyökfüggvényre, vagy az összes gyökre. (Tört hatványokra hasonló lesz, csak nem 2 érték lesz, hanem a kitevő nevezőjényi, egyszerűsítés után).

Illetve egy apróság, ha -0 van páratlan hatványra emelve az eredmény +0 lesz.

Azok úgyis egyenlőek, nem?

A "sírva, rossz helyen verem magamra a faszom" tipikus esete. :D

Vagy egy sunyibb lehetőség, ha a látszat csal: "Van egy feladatom. Nem tudom/akarom egyedül megoldani. Figyelemfelkeltő/önsajnáltató/önfényező/fröcsögő posztban felteszem a HUP-ra, hogy mások oldják meg nekem." :D

:)

Azt nem tudjuk egyszerűbben is? Nem tudom mit kell tudnia, de ilyen egyszerű dolgokat tud szerintem ez is:


#include<stdio.h>
#include<math.h>

double mypow(double x, int a, int b){
    if(x>=0) return pow(x, ((double)a)/((double)b));
    if( b % 2 == 1) return -pow(-x, ((double)a)/((double)b));
    return NAN;
}

int main(){
    double x = -8.0;
    int a = 1;
    int b = 3;
    double y = mypow(x, a, b);
    printf("%f ^ (%d / %d) = %f \n", x, a, b, y);

    b = 2;
    y = mypow(x, a, b);
    printf("%f ^ (%d / %d) = %f \n", x, a, b, y);

    return 0;
}

ha nem felejtettem ki valami fontosat.

Hát, szvsz. ez egy meglehetősen rossz megoldás. Egyrészt nem szerencsés inputfüggetlenül epsilonnal összehasonlítani (mi van, ha az 'a' nagyrágrendje mondjuk 10^20? Kb. soha se lesz igaz. Mi van, az a 'a' nagyságrendje 10^-20? Akkor meg mindig igaz lesz). De mégha ki is ütöd ezt belőle, nem szerencsés elvégezni az inverz műveletet csak azért, hogy összehasonlítsd újra az alappal. Elvégezni 4 db. pow-ot csak azért, hogy kitaláld, hogy melyik előjelű megoldás a jó, nem túl szép megoldás. Az eredeti kód ennél sokkal logikusabban működik (bár azon is lenne mit csiszolni). Remélem nem értettem semmit félre a kódban, Rustban nem igazán vagyok otthon.

Epsilon helyesbítést köszi. Mondhatni végre egy ember, aki figyelt. :)
Tempó: igen, az mindig kérdés részemről is, hogy milyen terhelhetőségre kell. Ahogy elnézem, powf-ből egy 8 évvel ezelőtti Core i3-as laptop proci 26 millió/mp tempóval tud magonként végezni.
Ha másodpercenként például 10.000 alatt kell ilyen számítás, akkor az átláthatóbb forráskód lehet hangsúlyosabb. Ha másodpercenként milliós nagyságrendben lesz használva, akkor viszont a gyorsabb futású megoldás lesz hangsúlyos.

Nem csak a futási sebesség miatt nem annyira jó ez a megoldás. De nem is annyira robosztus, valszeg egy csomó corner case-t nem kezel le helyesen (mi van, ha a bemeneten ez-az 0, vagy épp nagyon nagy/nagyon pici szám, overflow/underflow lehet, stb.). Az eredeti kód ebből a szempontból jobb, csak az meg túl van bonyolítva. Valszeg felesleges az lnko bele, elég lenne csak a 2-s faktorokat kiütni (ami ugye csak bitművelet a sok osztás helyett). Ill., a sok speciális eset kezelése se biztos, hogy kell bele. Valószínűleg meg lehet fogalmazni egyszerűbben. De amúgy tényleg jól rávilágít arra, hogy egy ilyen egyszerűbb matek problémához is kell megfelelő tudás jól lekódolni. De mondjuk ez szinte minden floating point-os rutinra igaz. Még egy egyszerű másodfokú megoldóképletnél is észnél kell lenni, ha robosztus rutint akar az ember írni.

Igen, érdemes átgondolni :) lehet kiderül, hogy nincs a bemenettől ilyen módon függő corner case. De nem igazán szeretnék ebbe energiát tenni, tekintve hogy nem én biztosan nem így csinálnám meg, meg ugye ez a kód így ebben a formában biztosan hibás, pl. a mypow(-1e-24, 1, 3)-re pozitív eredményt ad a negatív helyett (a korábban írt epsilon probléma miatt). Ha ezt kijavítod valahogy (pl. relatív hibát nézve az abszolut helyett, vagy valami ilyesmi), akkor is valszeg lehet olyan esetet találni, amit helytelenül kezel le (persze ez függ attól, hogyan javítod ki). A pow ugye nem ad garanciát arra nézve, hogy a számított érték a lehető legpontosabb érték. Szóval macera ez alapján belátni, hogy a számított pow érték pontosan milyen közel lesz az 'a'-hoz.

Ehhez képest az eredeti megoldás egyszerűbb, könnyebben be lehet látni, hogy helyesen működik (ha helyesen működik, igaziból azt se volt kedvem komolyabban átnézni), mivel nincs benne semmi floating point "mágia".

(Nem néztem át amúgy a kódodat tüzetesen, szóval nem tudom pontosan miért van benne 4-5 pow, meg pontosan miért is játszol az előjelekkel)

Szerintem nem helyes így gondolkodni erről. A mostani verzió a bemenetek egy igen jelentős részén hibázhat. Ha az 'a' kisebb, mint -1/epsilon, vagy 0 és -epsilon között van, akkor a hiba lehetősége fennáll (amennyire jól értem a dolgot, nem tettem bele különösen sok energiát). Ez pedig ugye a double egy igen jelentős tartománya, a negatív számok több, mint 90%-a. Csak azokra a számokra működik jól alapból, amik egyáltalán nincsenek a tartomány szélén semmilyen szempontból. Szóval először ezt ki kell javítani, és utána lehet arról beszélni, hogy vajon hol hibázhat. De csak hogy egy példát mondjak, hogy pl. hol mehet félre: ha az "a = -DBL_MAX", és mondjuk köbgyököt vonunk. Lehet, hogy a köbgyök eredményét újra köbre emelve már inf-et kapsz vissza a kerekítések miatt. Amit a mostani egyenlőségvizsgálat nem kezel le jól, de kérdéses, hogy ahogy kijavítod az epsilon hibát, az vajon ezt az esetet lekezeli-e. lll., hogy egyáltalán szép-e, hogy ilyet lekezelünk. Számomra eléggé egyértelmű, hogy nem.

Bocsánat, kimaradt a tesztelés, és az eredménye.


int main (void)
{
    int base, p, q;

    base= -4, p=  1, q= 2; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -4, p=  1, q= 1; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));

    base= -8, p=  1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p= -1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p= -2, q=-3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p= -1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p= -2, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p=  2, q=-3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base= -8, p= -1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=-16, p=  1, q= 2; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=-16, p=  2, q= 2; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=-16, p= -3, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));

    base=  8, p=  1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  8, p= -2, q=-3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  8, p= -1, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  8, p= -2, q= 3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  8, p=  2, q=-3; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));

    base=  0, p=  1, q= 0; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  0, p=  0, q= 1; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  0, p=  1, q= 1; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    base=  0, p= -1, q= 1; printf ("(%g)^(%d/%d)=%g\n", (double)base, p, q, mypow(base,p,q));
    return 0;
}

(-4)^(1/2)=nan
(-4)^(1/1)=-4
(-8)^(1/3)=-2
(-8)^(-1/3)=-0.5
(-8)^(-2/-3)=4
(-8)^(-1/3)=-0.5
(-8)^(-2/3)=0.25
(-8)^(2/-3)=0.25
(-8)^(-1/3)=-0.5
(-16)^(1/2)=nan
(-16)^(2/2)=-16
(-16)^(-3/3)=-0.0625
(8)^(1/3)=2
(8)^(-2/-3)=4
(8)^(-1/3)=0.5
(8)^(-2/3)=0.25
(8)^(2/-3)=0.25
(0)^(1/0)=nan
(0)^(0/1)=1
(0)^(1/1)=0
(0)^(-1/1)=nan