[IA32 assembler] Machtsverheffen

[Aaargh!]

Hier lijkt geen instructie voor ze zijn (ik maak gebruik van de FPU) iemand enig idee of dit efficienter te doen is dan een lusje waarin je X keer een getal met zichzelf vermenigvuldigd ?

[Rukapul]

Wat voor getallen wil je vermenigvuldigen? x^7, waarbij y geen integer is kan bijvoorbeeld niet op de manier die jij beschrijft. Indien y wel een integer is dan kun je optimaliseren met volgens mij log x vermenigvuldigingen. Kan het op Google niet direct vinden, maar het wordt oa behandeld in Applied Cryptography en maakt RSA gevoelig voor poweranalysis.

Hier staat volgens mij het algoritme (geleend van Openssl) met het belangrijkste stukje in de FOR loop:

code:

/* this one works - simple but works */
int BN_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
    {
    int i,bits,ret=0;
    BIGNUM *v,*rr;

    BN_CTX_start(ctx);
    if ((r == a) || (r == p))
        rr = BN_CTX_get(ctx);
    else
        rr = r;
    if ((v = BN_CTX_get(ctx)) == NULL) goto err;

    if (BN_copy(v,a) == NULL) goto err;
    bits=BN_num_bits(p);

    if (BN_is_odd(p))
        { if (BN_copy(rr,a) == NULL) goto err; }
    else    { if (!BN_one(rr)) goto err; }

    for (i=1; i<bits; i++)
        {
        if (!BN_sqr(v,v,ctx)) goto err;
        if (BN_is_bit_set(p,i))
            {
            if (!BN_mul(rr,rr,v,ctx)) goto err;
            }
        }
    ret=1;
err:
    if (r != rr) BN_copy(r,rr);
    BN_CTX_end(ctx);
    return(ret);
    }

[ Voor 57% gewijzigd door Rukapul op 10-02-2003 19:58 ]

[hobbit_be]

het lusje unrollen? en misschien opsplitsen in meerdere paths zodat je FPU lekker druk draait. x^4 = (x*x)*(x*x) die tussen haakjes kun je dan elke pipe zetten en op het laatst vermenigvuldigen (SSE ?). wel gek dat er geen asm opcode voor is?

[madwizard]

In het boek 'the art of assembly language' (online te lezen webster.cs.ucr.edu), chapter 14 (van de oude 16-bits versie maar maakt niet uit, fpu is hetzelfde gebleven) wordt alles over de FPU instructies uitgelegd en bovendien staan aan het eind een aantal functies die niet direct in de FPU zitten maar wel gemaakt kunnen worden (zoals machtsverheffen).
Het komt er op neer:

code:

(xlog is log met grondgetal x, lg is log met grondgetal 2)
x^y = z
xlog(z) = y
lg(z)/lg(x) = y (want xlog(y) is lg(y)/lg(x))
lg(z) = lg(x) * y
z = 2 ^ (lg(x)*y))

De laatste kun je wel uitrekenen met FPU instructies.

Wat ook wel handig is is Agner Fog's optimisation manual, zit bijvoorbeeld in het MASM32 pakket, daar staan tips in over code optimalisatie, onder andere voor bepaalde algoritmes zoals e^x.

edit:
Hier is een geoptimaliseerde snippet voor x ^ y.

[ Voor 8% gewijzigd door madwizard op 11-02-2003 11:56 ]

[Verwijderd]

hobbit_be schreef op 10 February 2003 @ 19:57:
het lusje unrollen? en misschien opsplitsen in meerdere paths zodat je FPU lekker druk draait. x^4 = (x*x)*(x*x) die tussen haakjes kun je dan elke pipe zetten en op het laatst vermenigvuldigen (SSE ?). wel gek dat er geen asm opcode voor is?

Een loopnz gebruiken en de code binnen de loop zo kort mogelijk houden zal veel sneller zijn dan het unrollen, de processor kan dan immers de loop in zijn pipeline houden...

meerder paden is leuk, maar (x*x) == (x*x), dus het heeft geen zin dat meerder keren te gaan uitrekenen...

De oplossing die Rukapul geeft is volgens mij de enige efficiente manier, waarbij de gedacht er achter is, dat

x^2*k == x^k * x^k


en dat je machtsverheffingen kunt opslitsen zoals in:

x¹¹ == x⁸ * x² * x¹

Wat overeenkomt met de binaire representatie:

x^0b1011 == x^0b1000 * x^0b0010 * x^0b0001


Als je dezelfde x blijft, maar je tot variabel machten verheft (bijvoorbeeld in het uitwerken van een taylor-reeks), kun je zelfs het beste de (recursieve) kwadraten van x ergens bijhouden, dan wordt het allemaal nog sneller...

[ Voor 8% gewijzigd door Verwijderd op 11-02-2003 13:17 ]

[EXX]

Mijn oude 8 bitter rekent x^y uit als e^( y * ln(x)) (oftwel EXP(Y*LOG(X)). Misschien is dat wat.

Overigens zit daar natuurlijk geen FPU in, de EXP en LOG worden uitgerekend m.b.v taylor reeksen, die weer als tabellen zijn opgeslagen. Anders duurt het uitrekenen bij zo een functie een eeuwigheid bij een 2 MHz CPU.

Ik kan je de Assembly code listing wel geven

[Verwijderd]

voor een machtverheffing met een positief geheel 32-bits getal, zou je het volgend kunnen gebruiken:

C:

double pwr_dbl_ul(double x ,unsigned long power)
{
    double  retval  = 1;
    double  current_value  = x;

    while (power)
    {
        if (power & 1)
            retval *= current_value;

        power >>= 1;

        current_value *= current_value;
    }
    return retval;
}

Wat in intel-syntax-assembly-volgens-gcc weer zou kunnen worden uitgeschreven als:

GAS:

pwr_dbl_ul:
    mov %eax, DWORD PTR [%esp+12]
    test    %eax, %eax
    fld1
    fldl    QWORD PTR [%esp+4]
    je  _exit

_while:
    test    %eax, 1
    je  _skip
    fmul    %st(1), %st

_skip:
    shr %eax, 1
    fmul    %st, %st(0)
    jne _while

_exit:
    fstp    %st(0)
    ret

[Verwijderd]

EXX schreef op 11 februari 2003 @ 12:07:
Mijn oude 8 bitter rekent x^y uit als e^( y * ln(x)) (oftwel EXP(Y*LOG(X)). Misschien is dat wat.

Overigens zit daar natuurlijk geen FPU in, de EXP en LOG worden uitgerekend m.b.v taylor reeksen, die weer als tabellen zijn opgeslagen. Anders duurt het uitrekenen bij zo een functie een eeuwigheid bij een 2 MHz CPU.

Ik kan je de Assembly code listing wel geven

Dat is heeeeel traag om het op die manier te doen... "Met behulp van taylor reeksen", was dat niet:
a₀ + a₁ * x¹ + a₂ * x² + a₃ * x³ ... + a_k * x^k

Het algorithme waar Rukapul mee kwam is echt veel sneller, geloof mij maar...

[roelio]

Die Taylorreeksen zitten in tabellen, dus wordne niet volledig uitgerekend. Immers, dan zit je weer met machten

[Verwijderd]

limoentje schreef op 11 February 2003 @ 13:15:
Die Taylorreeksen zitten in tabellen, dus wordne niet volledig uitgerekend. Immers, dan zit je weer met machten

Ja en DUS kun je in jouw 8 bit geval zelfs net zo goed een tabel maken voor de machtsverheffing van 7 keer 256 items - in een keer...
(da's een tabel van 256 mogelijke waarden voor x, met de machten x²; x⁴; x⁸; x¹⁶; x³²; x⁶⁴; x¹²⁸. Uiteraard overbodig om een x¹ lookup te doen )

Lookuptable zijn leuk voor <= 16 bit integers, al helemaal niet praktisch realiseerbaar voor floats/doubles...

[EXX]

Verwijderd schreef op 11 February 2003 @ 13:36:
[...]

Ja en DUS kun je in jouw 8 bit geval zelfs net zo goed een tabel maken voor de machtsverheffing van 7 keer 256 items - in een keer...
(da's een tabel van 256 mogelijke waarden voor x, met de machten x²; x⁴; x⁸; x¹⁶; x³²; x⁶⁴; x¹²⁸. Uiteraard overbodig om een x¹ lookup te doen )

Lookuptable zijn leuk voor <= 16 bit integers, al helemaal niet praktisch realiseerbaar voor floats/doubles...

Die Taylorreeksen werden algemeen gebruikt voor machtverheffen, logaritmes, etc. De reeksen zelf zijn tabellen, maar in voor single precision floating point. Zo een tabel ziet dan zo uit:

code:

; Coefficients for EXP
; '!' means faculty of the number (3! = 1 * 2 * 3)

  01479 08              DEFB    08H             ;8 coefficients

  0147A 40              DEFB    40H             ;-1.41361E-4 = -1/7076
  0147B 2E              DEFB    2EH             ;approx. -1/5040 = -1/7!
  0147C 94              DEFB    94H
  0147D 74              DEFB    74H

  0147E 70              DEFB    70H             ;1.32988E-3 = 1/752
  0147F 4F              DEFB    4FH             ;approx. 1/720 = 1/6!
  01480 2E              DEFB    2EH
  01481 77              DEFB    77H

  01482 6E              DEFB    6EH             ;=8.30136E-3 = -1/120
  01483 02              DEFB    02H                          = -1/5!
  01484 88              DEFB    88H
  01485 7A              DEFB    7AH

  01486 E6              DEFB    E6H             ;0.0416574 = 1/24
  01486 A0              DEFB    A0H                        = 1/4!
  01488 2A              DEFB    2AH
  01489 7C              DEFB    7CH

  0148A 50              DEFB    50H             ;-0.166665 = -1/6
  0148B AA              DEFB    AAH                        = -1/3!
  0148C AA              DEFB    AAH
  0148D 7E              DEFB    7EH

  0148E FF              DEFB    FFH             ;0.5 = 1/2
  0148F FF              DEFB    FFH                  = 1/2!
  01490 7F              DEFB    7FH
  01491 7F              DEFB    7FH

  01492 00              DEFB    00H             ;-1 = -1/1!
  01493 00              DEFB    00H
  01494 80              DEFB    80H
  01495 81              DEFB    81H

  01496 00              DEFB    00H             ;1
  01497 00              DEFB    00H
  01498 00              DEFB    00H
  01499 81              DEFB    81H

Berekenen vd reeks gaat dan bv. als volgt:

code:

; Calculate Taylor-series of the following form:
;
; y = k1 + k2*x + k3*x*x*x... (k1, k2, k3 are coefficients)
;
; I: HL -> coefficients table
;          The first byte of the table indicates the number of coefficients
;          in the table. Then the coefficients follow, in an inverted
;          order (k1 last)
;    X = factor in series (x in the example)
; O: X = result of series calculation (y in the example)

  014A9 CDA409          CALL    09A4H           ;(SP) = X
  014AC 7E              LD      A,(HL)          ;A = number of coefficients
  014AD 23              INC     HL              ;HL -> 1st number (last
                                                ;coefficient in the series)
  014AE CDB109          CALL    09B1H           ;X = BCDE = (HL) = 1st coeff.
  014B1 06F1            LD      B,0F1H          ;--
* 014B1   F1            POP     AF              ;Restore counter
  014B3 C1              POP     BC              ;BCDE = X
  014B4 D1              POP     DE
  014B5 3D              DEC     A               ;Any more coefficients ?
  014B6 C8              RET     Z               ;No: done, return

  014B7 D5              PUSH    DE              ;Save BCDE
  014B8 C5              PUSH    BC
  014B9 F5              PUSH    AF              ;Save counter
  014BA E5              PUSH    HL              ;Save pointer
  014BB CD4708          CALL    0847H           ;Result = result + BCDE
  014BE E1              POP     HL              ;Restore pointer
  014BF CDC209          CALL    09C2H           ;BCDE = (HL) = coefficient
  014C2 E5              PUSH    HL              ;Save pointer
  014C3 CD1607          CALL    0716H           ;Result = result + BCDE
  014C6 E1              POP     HL              ;Restore pointer
  014C7 18E9            JR      14B2H           ;Calculate next term

Naar de huidige maatstaven met 32 bit FPU's is dit natuurlijk hopeloos onnauwkeurig, maar het was voor 8-bits computers de enige zinvolle manier.

[Verwijderd]

EXX schreef op 11 februari 2003 @ 13:49:
[...]

Die Taylorreeksen werden algemeen gebruikt voor machtverheffen, logaritmes, etc. De reeksen zelf zijn tabellen, maar in voor single precision floating point. Zo een tabel ziet dan zo uit:

<knip>

Naar de huidige maatstaven met 32 bit FPU's is dit natuurlijk hopeloos onnauwkeurig, maar het was voor 8-bits computers de enige zinvolle manier.

32 bit FPU's

Dit zijn de coefficienten die jij geeft. Dat is hetgeen wat waarmee het resultaat van x^k mee wordt gewogen. De machtverheffing zelf gebeurt alsnog ergens anders in een van die calls - waarschijnlijk volgens het algoritme waar Rukapul mee kwam...

Maar die snippet van madwizard is idd het handigs voor een <double/float> tot-de-macht <double/float>.

[EXX]

Verwijderd schreef op 11 February 2003 @ 14:21:
[...]

32 bit FPU's

Eeeh, zijn de FPU's van tegenwoordig dan niet 32 bit?

Dit zijn de coefficienten die jij geeft. Dat is hetgeen wat waarmee het resultaat van x^k mee wordt gewogen. De machtverheffing zelf gebeurt alsnog ergens anders in een van die calls - waarschijnlijk volgens het algoritme waar Rukapul mee kwam...

De machtverheffing wordt als volgt uitgevoerd:

Je hebt X^Y. Eerste berekent ie LOG(X), vervolgens Y*LOG(X) en dan EXP(Y*LOG(X)). Bij zowel de EXP als de LOG functie maakt ie dan gebruik van de Taylor reeks tabellen.

Maar die snippet van madwizard is idd het handigs voor een <double/float> tot-de-macht <double/float>.

Ongetwijfeld

[Verwijderd]

EXX schreef op 11 February 2003 @ 14:59:
Je hebt X^Y. Eerste berekent ie LOG(X), vervolgens Y*LOG(X) en dan EXP(Y*LOG(X)). Bij zowel de EXP als de LOG functie maakt ie dan gebruik van de Taylor reeks tabellen.

Misschien komt dit als een schok voor jou, maar de taylor-reeks gebruikt daarbij een reeks machten...

Namelijk:
LOG(x) => c_l0 + c_l1 * x + c_l2 * x² + c_l3 * x³ ...

EXP(x) => c_e0 + c_e1 * x + c_e2 * x² + c_e3 * x³ ...

Dus jij voorkomt ook niet dat er machten geheven moeten worden, het gaat nu om de manier hoe dit zo efficient mogelijk te doen...Zeer goed mogelijk dat de TS dit nodig heeft juist t.b.v. een taylor-reeks berekening...

[EXX]

Oei, ik begin het een beetje te begrijpen, we praten langs elkaar heen

Nou, dan maar helemaal tot op de bodem:

De LOG functie:

code:

; X = LOG ( X )
; -------------
; Calculates the natural logarithm of X
;
; I: X = numerical value (<> 0)
; O: X = LOG (numerical value)

  00809 CD5509          CALL    0955H           ;TEST2
  0080C B7              OR      A               ;Argument = 0 ?
  0080D EA4A1E          JP      PE,1E4AH        ;Yes: ?FC Error

  00810 212441          LD      HL,4124H        ;HL -> Exp (argument)
  00813 7E              LD      A,(HL)          ;A = Exp (argument)
  00814 013580          LD      BC,8035H        ;BCDE = 0.707107 = SQR(2)/2
  00817 11F304          LD      DE,04F3H
                                                ;Conversion Arg = x * 2 ^ n
  0081A 90              SUB     B               ;Exp (Arg) - 128 = n (B is 128)
  0081B F5              PUSH    AF              ;save n
  0081C 70              LD      (HL),B          ;X = x (set Exp to 80H)
  0081D D5              PUSH    DE              ;Save BCDE
  0081E C5              PUSH    BC
  0081F CD1607          CALL    0716H           ;X = X + BCDE = X+1/2*SQR(2)
  00822 C1              POP     BC              ;Restore BCDE
  00823 D1              POP     DE
  00824 04              INC     B               ;Exp(BCDE) + 1
                                                ;BCDE = BCDE * 2 + SQR(2)
  00825 CDA208          CALL    08A2H           ;X = BCDE / X = SQR(2) / X
  00828 21F807          LD      HL,07F8H        ;HL -> 1.0
  0082B CD1007          CALL    0710H           ;X = (HL) - X = 1.0 - X
  0082E 21FC07          LD      HL,07FCH        ;HL -> numeric table
  00831 CD9A14          CALL    149AH           ;Compute row1 (Taylor!)
  00834 018080          LD      BC,8080H        ;BCDE = -0.5
  00837 110000          LD      DE,0000H
  0083A CD1607          CALL    0716H           ;X = BCDE + X = X - 0.5
  0083D F1              POP     AF              ;Restore n
  0083E CD890F          CALL    0F89H           ;X = X + n
                                                ;and multiply with LOG(2)

; X = X * LOG(2)

  00841 013180          LD      BC,8031H        ;BCDE = 0.693147 = LOG(2)
  00844 111872          LD      DE,7218H

; SMUL: X = BCDE * X
; Multiply two single precision numbers
;
; I: BCDE = 1st factor
;    X    = 2nd factor
; O: X    = product

  00847 CD5509          CALL    0955H           ;TEST2
  0084A C8              RET     Z               ;X = 0: result = 0

  0084B 2E00            LD      L,00H           ;Flag = 0 (MUL indication)
  0084D CD1409          CALL    0914H           ;Process exponent
  00850 79              LD      A,C             ;Save CDE (1st factor) in
  00851 324F41          LD      (414FH),A       ;system RAM from 414FH onwards
  00854 EB              EX      DE,HL
  00855 225041          LD      (4150H),HL
  00858 010000          LD      BC,0000H        ;BCDE = 00000000H
  0085B 50              LD      D,B
  0085C 58              LD      E,B
  0085D 216507          LD      HL,0765H        ;Put new return address
  00860 E5              PUSH    HL              ;to 0765H
  00861 216908          LD      HL,0869H        ;Put new return address twice
  00864 E5              PUSH    HL              ;to 0869H
  00865 E5              PUSH    HL
  00866 212141          LD      HL,4121H        ;HL -> 2nd factor
  00869 7E              LD      A,(HL)          ;A = next byte of mantissa of
                                                ;the 2nd factor
  0086A 23              INC     HL              ;Pointer + 1
  0086B B7              OR      A               ;Byte = 00H ?
  0086C 2824            JR      Z,0892H         ;Yes: continue at 0892H

  0086E E5              PUSH    HL              ;Save pointer
  0086F 2E08            LD      L,08H           ;L = counter for 8 bits
  00871 1F              RRA                     ;Shift next bit into C-flag
  00872 67              LD      H,A             ;Save byte in H
  00873 79              LD      A,C             ;A = MSB
                                                ;Bit set by last shift ?
  00874 300B            JR      NC,0881H        ;No: continue at 0881H

  00876 E5              PUSH    HL              ;Save HL
  00877 2A5041          LD      HL,(4150H)      ;CDE = CDE + 1st factor
  0087A 19              ADD     HL,DE
  0087B EB              EX      DE,HL
  0087C E1              POP     HL              ;Restore HL
  0087D 3A4F41          LD      A,(414FH)
  00880 89              ADC     A,C             ;A = MSB
  00881 1F              RRA                     ;CDEB one bit to the right
  00882 4F              LD      C,A
  00883 7A              LD      A,D             ;Shift D
  00884 1F              RRA
  00885 57              LD      D,A
  00886 7B              LD      A,E             ;Shift E
  00887 1F              RRA
  00888 5F              LD      E,A
  00889 78              LD      A,B             ;Shift B
  0088A 1F              RRA
  0088B 47              LD      B,A
  0088C 2D              DEC     L               ;Counter - 1
  0088D 7C              LD      A,H             ;Byte back into A
  0088E 20E1            JR      NZ,0871H        ;Check next bit

  00890 E1              POP     HL              ;HL -> X
  00891 C9              RET                     ;RET twice to 0896H and
                                                ;once to 0765H

De Taylor berekening die wordt aangeroepen.

code:

; Series calculation 1
; Calculate Taylor-series of the following form:
;
; y = k1*x + k2*x*x*x + k3*x*x*x*x*x... (k1, k2, k3 are coefficients)
;
; I: HL -> coefficients table
;          The first byte of the table indicates the number of coefficients
;          in the table. Then the coefficients follow, in an inverted
;          order (k1 last)
;    X = factor in series (x in the example)
; O: X = result of series calculation (y in the example)

  0149A CDA409          CALL    09A4H           ;(SP) = X
  0149D 11320C          LD      DE,0C32H        ;DE = address of X = X * (SP)
  014A0 D5              PUSH    DE              ;Save as new RET address
  014A1 E5              PUSH    HL              ;Save table pointer
  014A2 CDBF09          CALL    09BFH           ;BCDE = X
  014A5 CD4708          CALL    0847H           ;X = BCDE * X = X * X
  014A8 E1              POP     HL              ;Restore table pointer and
                                                ;Series calulation 2 using
                                                ;X * X and multiply result
                                                ;again with X
; Series calculation 2
; Calculate Taylor-series of the following form:
;
; y = k1 + k2*x + k3*x*x*x... (k1, k2, k3 are coefficients)
;
; I: HL -> coefficients table
;          The first byte of the table indicates the number of coefficients
;          in the table. Then the coefficients follow, in an inverted
;          order (k1 last)
;    X = factor in series (x in the example)
; O: X = result of series calculation (y in the example)

  014A9 CDA409          CALL    09A4H           ;(SP) = X
  014AC 7E              LD      A,(HL)          ;A = number of coefficients
  014AD 23              INC     HL              ;HL -> 1st number (last
                                                ;coefficient in the series)
  014AE CDB109          CALL    09B1H           ;X = BCDE = (HL) = 1st coeff.
  014B1 06F1            LD      B,0F1H          ;--
* 014B1   F1            POP     AF              ;Restore counter
  014B3 C1              POP     BC              ;BCDE = X
  014B4 D1              POP     DE
  014B5 3D              DEC     A               ;Any more coefficients ?
  014B6 C8              RET     Z               ;No: done, return

  014B7 D5              PUSH    DE              ;Save BCDE
  014B8 C5              PUSH    BC
  014B9 F5              PUSH    AF              ;Save counter
  014BA E5              PUSH    HL              ;Save pointer
  014BB CD4708          CALL    0847H           ;Result = result + BCDE
  014BE E1              POP     HL              ;Restore pointer
  014BF CDC209          CALL    09C2H           ;BCDE = (HL) = coefficient
  014C2 E5              PUSH    HL              ;Save pointer
  014C3 CD1607          CALL    0716H           ;Result = result + BCDE
  014C6 E1              POP     HL              ;Restore pointer
  014C7 18E9            JR      14B2H           ;Calculate next term

De gebruikte tabellen met coefficienten:

code:

; Table of SNG coefficients for LOG function

  007FC 03              DEFB    03H             ;3 Coefficients

  007FD AA              DEFB    AAH             ;0.598979 approx.
  007FE 56              DEFB    56H             ;2 / ( 5*LOG(2) )
  007FF 19              DEFB    19H
  00800 80              DEFB    80H

  00801 F1              DEFB    F1H             ;0.961471 approx.
  00802 22              DEFB    22H             ;2 / ( 3*LOG(2) )
  00803 76              DEFB    76H
  00804 80              DEFB    80H

  00805 45              DEFB    45H             ;2.88539  approx.
  00806 AA              DEFB    AAH             ;2 / ( 1*LOG(2) )
  00807 38              DEFB    38H
  00808 82              DEFB    82H

Note: De programmafragmenten zijn afkomstig uit de (TRS80) Level II BASIC, een Microsoft BASIC variant.

In deze hele handel zie ik geen "echt" machtverheffen. Er wordt een hoop vermenigvuldigd, gedeeld, opgeteld, afgetrokken en geschoven, maar dat zijn "normale" basisbewerkingen.

[.oisyn]

EXX schreef op 11 February 2003 @ 14:59:
[...]
Eeeh, zijn de FPU's van tegenwoordig dan niet 32 bit?

nee, die zijn al sinds het begin van de x87 80 bits

[EXX]

Aha, weer wat geleerd

Idd, als je met floating point moet gaan rekenen is 32 bits wel erg weinig. Dat had ik zelf ook wel kunnen bedenken

[ Voor 72% gewijzigd door EXX op 11-02-2003 15:56 ]

[.oisyn]

Valt wel mee, 32 bits floats worden vaak gebruikt als opslag (24 bits mantissa, 8 bits exponent en een sign-bit. Heej, da's 33 bits Klopt, de eerste bit van de mantissa is altijd 1 en wordt derhalve niet opgeslagen ).

[Aaargh!]

Verwijderd schreef op 11 februari 2003 @ 15:15:
Dus jij voorkomt ook niet dat er machten geheven moeten worden, het gaat nu om de manier hoe dit zo efficient mogelijk te doen...Zeer goed mogelijk dat de TS dit nodig heeft juist t.b.v. een taylor-reeks berekening...

Het is hier voor nodig.

[Verwijderd]

EXX schreef op 11 February 2003 @ 15:40:
Oei, ik begin het een beetje te begrijpen, we praten langs elkaar heen

Nou, dan maar helemaal tot op de bodem:

Jouw multiply maakt overigens gebruik van eenzelfde principe voor het mantissa gedeelte:

a * b = a* (b & 1) + X*(b & 2) + X * (b & 4) ...

Wat weer herschreven kan worden naar optelsommen van shift-lefts zoals:

C:

unsigned long mul_ab(unsigned short b, unsigned short a)
{
    unsigned long   retval  = 0;
    unsigned long   current_value  = b;

    while (a)
    {
        if (a & 1  == 1)
            retval += current_value;

        a >>= 1;
        current_value <<= 1;
    }
    return retval;
}

[ Voor 27% gewijzigd door Verwijderd op 11-02-2003 16:18 ]

[Verwijderd]

Aaargh! schreef op 11 februari 2003 @ 16:10:
[...]

Het is hier voor nodig.

Dan zou ik zeker voor de snippet gaan, aangezien er heel waarschijnlijk veel verschillende waarden tot een variabele macht moeten worden geheven...

[EXX]

Verwijderd schreef op 11 februari 2003 @ 16:11:
[...]

Jouw multiply maakt overigens gebruik van eenzelfde principe voor het mantissa gedeelte:

a * b = a* (b & 1) + X*(b & 2) + X * (b & 4) ...

Wat weer herschreven kan worden naar optelsommen van shift-lefts zoals:

[stukje c-code]

Dat klopt, maar in het hele verhaal komt dus geen (complex) machtsverheffen as such voor. De hele berekening wordt met basis-rekenkunde en veel bitgeschuifel (vooral dit laatste is erg efficient/snel voor CPU's , zeker als die geen MULT of DIV instructies hebben) gedaan, samen met een tabel-gerelateerde Taylor reeks en een paar constantes. Dus als je rap moet zijn is zulk een algoritme misschien nog niet zo slecht.

edit:
Idd, mee eens, het algoritme uit de snippet is zeker het overwegen waard

Overigens kan mijn stukje voorbeeldcode op een moderne x86 CPU waarschijnlijk stukken efficienter, omdat die MUL en DIV instructies hebben , en veel bredere registers.

[ Voor 50% gewijzigd door EXX op 11-02-2003 16:26 ]

[Verwijderd]

EXX schreef op 11 februari 2003 @ 16:21:
[...]
Dat klopt, maar in het hele verhaal komt dus geen (complex) machtsverheffen as such voor.

Omdat de z80 dat ook niet kan nee.

Dus als je rap moet zijn is zulk een algoritme misschien nog niet zo slecht.

FOUT, het algoritme wat een MUL-circuit gebruikt in een cpu gebruikt hetzelfde algoritme als ik beschreef met mul_ab alleen dan hardwarematig waardoor ge-optimaliseerd sneller, vaak volledig met een hardwarematige lookuptabel, of in elk geval gewoon in 1 (effectieve) cycle, sneller kan niet.

De floating point exp en log instructies in een FPU gebruikt vaak weer een coeff. lookuptabel met hardwarematige taylor-polynomen. Omdat dit in hardware is gebouwd zal een hardwarevariant altijd minstens zo snel zijn als hetzelfde algorithme in software.

edit:
Anders worden dat soort HW instructies ook niet gebruikt en hebben de chipontwerpers wel wat beters te doen met hun ruimte

[ Voor 9% gewijzigd door Verwijderd op 11-02-2003 16:51 ]

[EXX]

Verwijderd schreef op 11 February 2003 @ 16:42:
FOUT, het algoritme wat een MUL-circuit gebruikt in een cpu gebruikt hetzelfde algoritme als ik beschreef met mul_ab alleen dan hardwarematig waardoor ge-optimaliseerd sneller, vaak volledig met een hardwarematige lookuptabel, of in elk geval gewoon in 1 (effectieve) cycle, sneller kan niet.

De floating point exp in een FPU gebruikt vaak weer een coeff. lookuptabel met hardwarematige taylor-polynomen. Omdat dit in hardware is gebouwd zal een hardwarevariant altijd minstens zo snel zijn als hetzelfde algorithme in software.

Yikes, ik zat even niet op te letten. Natuurlijk is een FPU waar dit soort zaken zijn ingebakken rapper. Ik was weer teveel "CPU only" aan het denken.

edit:
Is er geen machtverhef-functie in de huidige FPU? Ik heb hier wat docu van de antieke 8087 en die had al een FYL2X = Y * log2(X) (X en Y op de stack).
Zou je die niet kunnen gebruiken want x^y = 2^ (y * log2(x))

[ Voor 16% gewijzigd door EXX op 11-02-2003 16:58 ]

[Verwijderd]

EXX schreef op 11 februari 2003 @ 16:47:
[...]
Yikes, ik zat even niet op te letten. Natuurlijk is een FPU waar dit soort zaken zijn ingebakken rapper. Ik was weer teveel "CPU only" aan het denken.

D'r zijn ook genoeg CPU's die naast een ALU/AGU ook een FPU hebben, die dus "deel" uitmaakt van de "CPU". (Volgens mij was de 386 de eerste en laaste x86 processor met een FPU co-processor, maar misschien had de 486 die ook nog wel, maar voor de rest is de FPU gewoon onderdeelgeworden van de CPU zelf.)

De Z180 was overigens meteen veel populairder (en duurder ) omdat deze wel over een multiply instructie beschikte. Daar ben ik toen snel naar gevlucht...

[ Voor 21% gewijzigd door Verwijderd op 11-02-2003 17:03 ]

[.oisyn]

late reactie, maar:

Verwijderd schreef op 11 februari 2003 @ 13:02:
Wat in intel-syntax-assembly-volgens-gcc

dat noemen ze AT&T syntax

[Verwijderd]

.oisyn schreef op 11 February 2003 @ 17:01:
late reactie, maar:
[...]

dat noemen ze AT&T syntax

Nee hoor - AT&T is de syntax die gcc geeft zonder de optie "-mintel-syntax", die van links naar rechts leest. De assembly die ik gaf was met de optie... Maar grappig dat je die regel opviel (uiteindelijk)...

[Olaf van der Spek]

EXX schreef op 11 februari 2003 @ 14:59:
Eeeh, zijn de FPU's van tegenwoordig dan niet 32 bit?

Nee, de FPU is sinds de XT 80-bit (veel te laat).
Veel apps werken echter met 32-bit floats of 64-bit doubles.

[ Voor 21% gewijzigd door Olaf van der Spek op 11-02-2003 19:19 ]

[EXX]

Verwijderd schreef op 11 februari 2003 @ 17:00:
[...]

D'r zijn ook genoeg CPU's die naast een ALU/AGU ook een FPU hebben, die dus "deel" uitmaakt van de "CPU". (Volgens mij was de 386 de eerste en laaste x86 processor met een FPU co-processor, maar misschien had de 486 die ook nog wel, maar voor de rest is de FPU gewoon onderdeelgeworden van de CPU zelf.)

Ook de 8086/88 had een FPU, de 8087. Voor bij de 80186/88 had je de 80187 en voor de 80286 ( je raad het al) de 80287).

madwizard schreef op 11 februari 2003 @ 11:45:
In het boek 'the art of assembly language' (online te lezen webster.cs.ucr.edu), chapter 14 (van de oude 16-bits versie maar maakt niet uit, fpu is hetzelfde gebleven) wordt alles over de FPU instructies uitgelegd en bovendien staan aan het eind een aantal functies die niet direct in de FPU zitten maar wel gemaakt kunnen worden (zoals machtsverheffen).
Het komt er op neer:

code:

(xlog is log met grondgetal x, lg is log met grondgetal 2)
x^y = z
xlog(z) = y
lg(z)/lg(x) = y (want xlog(y) is lg(y)/lg(x))
lg(z) = lg(x) * y
z = 2 ^ (lg(x)*y))

De laatste kun je wel uitrekenen met FPU instructies.
Hier is een geoptimaliseerde snippet voor x ^ y.

Dit lijkt me idd de meest optimale oplossing.

Ik weet niet of het zin heft, maar de SSE / SSE2 instructieset bevatten ook rekenkundige bewerkingen. IIRC is de SSE / SSE2 unit behoorlijk rapper dan de FPU.

edit:
Dan mot je natuurlijk wel een CPU gebruiken met SSE / SSE2 aan boord

[ Voor 5% gewijzigd door EXX op 12-02-2003 10:00 ]

[.oisyn]

Verwijderd schreef op 11 February 2003 @ 17:28:
[...]

Nee hoor - AT&T is de syntax die gcc geeft zonder de optie "-mintel-syntax", die van links naar rechts leest. De assembly die ik gaf was met de optie... Maar grappig dat je die regel opviel (uiteindelijk)...

hmm idd, je hebt gelijk... ik wist niet eens dat die optie bestond
Ik had er in mijn DJGPP tijdperk altijd zo'n hekel aan, vandaar dat ik mijn asm altijd mbv nasm meecompilede

[RobIII]

code:

; #########################################################################
;
;                             FpuXexpY
;
;##########################################################################

  ; -----------------------------------------------------------------------
  ; This procedure was written by Raymond Filiatreault, December 2002
  ;
  ;            Src1^Src2 = antilog2[ log2(Src1) * Src2 ] -> Dest
  ;
  ; This FpuXexpY function raises a number (Src1) to a power (Src2)
  ; with the FPU and returns the result as an 80-bit REAL number at the
  ; specified destination (the FPU itself or a memory location), unless an
  ; invalid operation is reported by the FPU or the definition of the
  ; parameters (with uID) is invalid.
  ;
  ; Either the number or the power can be an 80-bit REAL number from the
  ; FPU itself or from memory, an immediate DWORD value or one in memory,
  ; or one of the FPU constants. The base number (Src1) must be positive.
  ;
  ; The sources are not checked for validity. This is the programmer's
  ; responsibility.
  ;
  ; Only EAX is used to return error or success. All other registers are
  ; preserved.
  ;
  ; -----------------------------------------------------------------------

    .486
    .model flat, stdcall  ; 32 bit memory model
    option casemap :none  ; case sensitive

    include Fpu.inc

    .data
    
    tempdw  dd    0    
    stword  dw    0

    .code

; #########################################################################

FpuXexpY proc public lpSrc1:DWORD, lpSrc2:DWORD, lpDest:DWORD, uID:DWORD
        
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
;
; Because a library is assembled before its functions are called, all
; references to external memory data must be qualified for the expected
; size of that data so that the proper code is generated.
;
;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

;the FPU will be initialized only if neither of the two source parameters
;is taken from the FPU itself.

      fclex                   ;clear exception flags on FPU
      test  uID,SRC1_FPU or SRC2_FPU     ;is data taken from FPU?
      jnz   @F
      finit

;----------------------------------------
;check source for Src1 and load it to FPU
;----------------------------------------

@@:
      mov   eax,lpSrc1
      test  uID,SRC1_FPU      ;is Src1 taken from FPU?
      jnz   src2              ;check next parameter Src2 for exponent
      
      test  uID,SRC1_REAL     ;is Src1 an 80-bit REAL in memory?
      jz    @F
      fld   tbyte ptr [eax]
      jmp   src2              ;check next parameter Src2 for exponent
@@:
      test  uID,SRC1_DMEM     ;is Src1 a 32-bit integer in memory?
      jz    @F
      fild  dword ptr [eax]
      jmp   src2              ;check next parameter Src2 for exponent
@@:
      test  uID,SRC1_DIMM     ;is Src1 an immediate 32-bit integer?
      jz    @F
      mov   tempdw,eax
      fild  tempdw
      jmp   src2              ;check next parameter Src2 for exponent
@@:
      test  uID,SRC1_CONST    ;is Src1 one of the FPU constants?
      jnz   @F                ;otherwise no correct flag for Src1

srcerr:
      finit
      xor   eax,eax
      ret

@@:
      test  eax,FPU_PI
      jz    @F
      fldpi                   ;load pi (3.14159...) on FPU
      jmp   src2              ;check next parameter Src2 for exponent
@@:
      test  eax,FPU_NAPIER
      jz    srcerr            ;no correct CONST flag for Src1
      fld1
      fldl2e
      fsub  st,st(1)
      f2xm1
      fadd  st,st(1)
      fscale
      fxch
      fstp  st(0)

;----------------------------------------
;check source for Src2 and load it to FPU
;----------------------------------------

src2:
      mov   eax,lpSrc2
      test  uID,SRC2_FPU      ;is Src2 taken from FPU?
      jz    @F                ;check next source for Src2
      test  uID,SRC1_FPU      ;is Src1 the same source
      jz    src21
      fld   st(0)             ;copy itself on FPU
src21:
      fxch
      jmp   dest0             ;go complete process

@@:
      test  uID,SRC2_REAL     ;is Src2 an 80-bit REAL in memory?
      jz    @F
      fld   tbyte ptr [eax]
      jmp   dest0             ;go complete process
@@:
      test  uID,SRC2_DMEM     ;is Src2 a 32-bit integer in memory?
      jz    @F
      fild  dword ptr [eax]
      jmp   dest0             ;go complete process
@@:
      test  uID,SRC2_DIMM     ;is Src2 an immediate 32-bit integer?
      jz    @F
      mov   tempdw,eax
      fild  tempdw
      jmp   dest0             ;go complete process
@@:
      test  uID,SRC2_CONST    ;is Src2 one of the FPU constants?
      jz    srcerr            ;no correct flag for Src2
      test  eax,FPU_PI
      jz    @F
      fldpi                   ;load pi (3.14159...) on FPU
      jmp   dest0             ;go complete process
@@:
      test  eax,FPU_NAPIER
      jz    srcerr            ;no correct CONST flag for Src1
      fld1
      fldl2e
      fsub  st,st(1)
      f2xm1
      fadd  st,st(1)
      fscale
      fxch
      fstp  st(0)

dest0:
      fxch                    ;set up FPU registers for next operation
      fyl2x                   ;->log2(Src1)*exponent
      
;the FPU can compute the antilog only with the mantissa
;the characteristic of the logarithm must thus be removed
      
      fld   st(0)             ;copy the logarithm
      frndint                 ;keep only the characteristic
      fsub  st(1),st          ;keeps only the mantissa
      fxch                    ;get the mantissa on top

      f2xm1                   ;->2^(mantissa)-1
      fld1
      fadd                    ;add 1 back

;the number must now be readjusted for the characteristic of the logarithm

      fscale                  ;scale it with the characteristic
      
      fstsw stword            ;retrieve exception flags from FPU
      fwait
      test  stword,1          ;test for invalid operation
      jnz   srcerr            ;clean-up and return if error
      
;the characteristic is still on the FPU and must be removed

      fxch                    ;get the characteristic on top
      fstp  st(0)             ;"pop" it

      mov   eax,lpDest
      test  uID,DEST_FPU      ;check where result should be stored
      jnz   @F                ;leave result on FPU if so indicated
      fstp  tbyte ptr [eax]   ;store result at specified address

@@:
      or    al,1              ;to insure EAX!=0
      ret
    
FpuXexpY endp

; #########################################################################

end

Uit Masm32...