Cza­sa­mi przy­cho­dzi taka chwi­la w życiu ko­de­ra, gdy musi po­li­czyć pier­wia­stek. Koder sce­no­wy powie teraz „żaden pro­blem, wy­star­czy sta­bli­co­wać i od­czy­tać”, koder sys­te­mo­wy na­to­miast „otwie­ram math#?trans.library i po kło­po­cie”. Po­sia­dacz ko­pro­ce­so­ra ma­te­ma­tycz­ne­go stwier­dzi z wyż­szo­ścią „to prze­cież kwe­stia jed­nej in­struk­cji ko­pro­ca...”. Nie­ste­ty wszyst­kie po­dej­ścia mają swoje wady. Ta­bli­ca, ow­szem jest bar­dzo szyb­ka, ale zaj­mu­je tro­chę pa­mię­ci, zwłasz­cza gdy wy­ma­ga­na jest nieco więk­sza pre­cy­zja. Z dru­giej stro­ny wy­ko­rzy­sta­nie bi­blio­tek ma­te­ma­tycz­nych jest wy­god­ne i wy­ma­ga w miarę małej ilo­ści pa­mię­ci. Pier­wiast­ki są li­czo­ne bar­dzo do­kład­nie, ale też przy oka­zji bar­dzo wolno. Z kolei za­sto­so­wa­nie ko­pro­ce­so­ra ma­te­ma­tycz­ne­go ma pod­sta­wo­wą wadę, którą jest ko­niecz­ność po­sia­da­nia owego ko­pro­ce­so­ra :-). No i wresz­cie nie za­wsze mamy ocho­tę na prze­li­cza­nie licz­by cał­ko­wi­tej na for­mat zmien­no­prze­cin­ko­wy. W ta­kich przy­pad­kach nie po­zo­sta­je nam nic in­ne­go, jak za­siąść do kla­wia­tu­ry i na­pi­sać sa­me­mu taką pro­ce­dur­kę.

Star­si wy­ja­da­cze pa­mię­ta­ją za­pew­ne ar­ty­kuł na ten temat, który swego czasu uka­zał się w „Ami­gow­cu” Autor za­pro­po­no­wał tam wy­ko­rzy­sta­nie al­go­ryt­mu pier­wiast­ko­wa­nia pi­sem­ne­go. Kod źró­dło­wy zaj­mo­wał całą stro­nę ga­ze­ty, a i naj­szyb­szy też chyba nie był. Błąd po­le­gał na złym po­dej­ściu do pro­ble­mu. Pier­wiast­ko­wa­nie pi­sem­ne to al­go­rytm dla czło­wie­ka, który liczy w ukła­dzie dzie­sięt­nym. Kom­pu­ter liczy w sys­te­mie dwój­ko­wym i po­trzeb­ny jest al­go­rytm, który to weź­mie pod uwagę.

Za­sta­nów­my się teraz nad de­fi­ni­cją pier­wiast­ka kwa­dra­to­we­go. Jest to taka licz­ba, która pod­nie­sio­na do kwa­dra­tu daje nam licz­bę wej­ścio­wą. Pro­ste. Ale nie za­po­mi­naj­my, że dzia­ła­my na licz­bach cał­ko­wi­tych. Ob­ci­na­my wszyst­ko po prze­cin­ku i dla­te­go wynik pod­nie­sio­ny do kwa­dra­tu może być tro­chę mniej­szy niż licz­ba pier­wiast­ko­wa­na. Można więc po­wie­dzieć, że pier­wia­stek cał­ko­wi­ty to taka moż­li­wie naj­więk­sza licz­ba, któ­rej kwa­drat jest nie więk­szy niż licz­ba pier­wiast­ko­wa­na. Druga ważna spra­wa to tzw. prze­li­czal­ność zbio­ru liczb n-bi­to­wych. Co to ozna­cza? Je­że­li mamy re­jestr np. 16-bi­to­wy, to mo­że­my tam za­pi­sać skoń­czo­ną ilość róż­nych liczb (do­kład­nie 65536 i ani grama wię­cej). Czyli pro­blem pier­wiast­ko­wa­nia spro­wa­dza się do wy­bra­nia jed­nej od­po­wied­niej licz­by spo­śród skoń­czo­ne­go ich zbio­ru. Tylko jak wy­brać tę licz­bę? Jaki jest naj­szyb­szy al­go­rytm zna­le­zie­nia ta­kiej licz­by? Oczy­wi­ście po­ło­wie­nie prze­dzia­łu, jest on al­go­ryt­mem jak naj­bar­dziej „kom­pu­te­ro­wym” opar­tym na dwój­ko­wym sys­te­mie li­cze­nia. A więc ro­bi­my tak:

1. Wy­bie­ra­my licz­bę w po­ło­wie prze­dzia­łu.
2. Pod­no­si­my ją do kwa­dra­tu i po­rów­nu­je­my z licz­bą pier­wiast­ko­wa­ną.
3. Je­że­li jest mniej­sza to prze­dzia­łem staje się jego dolna po­ło­wa, w prze­ciw­nym wy­pad­ku górna.
4. Je­że­li prze­dział zwę­ził się do jed­nej licz­by, to ob­li­czy­li­śmy pier­wia­stek.
5. A jeśli nie to skok do punk­tu 1.

Jaki prze­dział wy­brać? Je­że­li pier­wiast­ku­je­my licz­bę n-bi­to­wą, pier­wia­stek na pe­wno zmie­ści się w n/2 bi­tach, oczy­wi­ście je­że­li pier­wiast­ku­je­my licz­bę o nie­pa­rzy­stej licz­bie bitów, to po­wsta­ją­ce pół bita na­le­ży za­okrą­glić w górę. Je­że­li ktoś nie wie­rzy, to pro­szę sobie to roz­pi­sać ko­rzy­sta­jąc z twier­dze­nia o po­tę­dze po­tę­gi: sqrt(2ⁿ) = (2ⁿ)^1/2 = 2^n/2. Teraz upro­ści­my sobie spra­wę jesz­cze bar­dziej i po­pra­cu­je­my nad samym po­ło­wie­niem prze­dzia­łu. Ob­la­ta­nym w elek­tro­ni­ce pod­po­wiem, że za­sto­su­je­my me­to­dę sto­so­wa­ną w aprok­sy­mu­ją­cych prze­twor­ni­kach C/A. Za­miast bawić się w dwie stop­nio­wo za­wę­ża­ją­ce się gra­ni­ce prze­dzia­łu bę­dzie­my je­chać jedną licz­bą jego środ­kiem. Dzię­ki temu wynik bę­dzie­my two­rzyć ka­su­jąc i usta­wia­jąc ko­lej­ne jego bity, za­czy­na­jąc od naj­star­sze­go.

Wy­ja­śnię to na przy­kła­dzie wy­cią­ga­nia pier­wiast­ka czte­ro­bi­to­we­go z licz­by ośmio­bi­to­wej. Prze­dział wy­ni­ku to 0–15. Usta­wia­jąc naj­star­szy bit wy­ni­ku otrzy­ma­my licz­bę 8, czyli śro­dek prze­dzia­łu. Pod­no­si­my ją do kwa­dra­tu i spraw­dza­my, czy nie jest więk­sza od licz­by pier­wiast­ko­wa­nej, je­że­li tak, to ka­su­je­my bit. Potem prze­cho­dzi­my do na­stęp­ne­go bitu i znowu spraw­dza­my. Je­że­li licz­ba jest za duża, to ka­su­je­my ten bit i prze­cho­dzi­my do na­stęp­ne­go itd. Oto przy­kład, li­czy­my pier­wia­stek z 174:

Al­go­rytm dla pier­wiast­ka z licz­by n-bi­to­wej wy­ma­ga wy­ko­na­nia n/2 mno­żeń liczb n/2-bi­to­wych, n/2 po­rów­nań liczb n-bi­to­wych, n/2 usta­wień bitu i od 0 do n/2 ska­so­wań bitu. Szyb­kość ob­li­czeń w małym stop­niu za­le­ży od licz­by pier­wiast­ko­wa­nej.

Od teo­rii pora więc przejść do prak­ty­ki. Oto pro­ce­du­ra w asem­ble­rze pier­wiast­ko­wa­nia licz­by 32-bi­to­wej:

; liczba w d0, wynik też w d0.
; pierwiastkowanie liczby 32-bitowej.

Root32:
    MOVEM.L d2/d3,-(sp)
    CLR.L;  d1
    MOVEQ   #15,d3
Root32Loop:
    BSET    d3,d1
    MOVE.W  d1,d2
    MULU    d2,d2
    CMP.L   d0,d2
    BLS.S   NoClear32
    BCLR    d3,d1
NoClear32:
    DBF     d3,Root32Loop
    MOVE.L  d1,d0
    MOVEM.L (sp)+,d2/d3
    RTS

Jak widać, z całej stro­ny wy­dru­ku zro­bi­ło się rap­tem 14 roz­ka­zów. Jed­nak na tym nie po­prze­sta­nie­my. Jak wia­do­mo, gdy po­trzeb­na jest więk­sza do­kład­ność ob­li­czeń, a licz­by zmien­no­prze­cin­ko­we na­le­ży omi­jać z da­le­ka, sto­su­je się ska­lo­wa­nie liczb. Mno­ży­my licz­by przed ob­li­cze­nia­mi, z re­gu­ły prze­su­wa­jąc je bi­to­wo w lewo o ileś tam bitów, a po skoń­czo­nych ob­li­cze­niach dzie­li­my je prze­su­wa­jąc w pra­wo. Dzię­ki temu błąd za­okrą­gle­nia po­peł­nia­my tylko raz. Pisał o tym swego czasu Mi­klesz w Ma­ga­zy­nie Amiga. Przy pier­wiast­ko­wa­niu, po­dob­nie jak przy mno­że­niu i dzie­le­niu zmie­nia się czyn­nik ska­lu­ją­cy. Je­że­li chce­my wynik pier­wiast­ko­wa­nia mieć prze­ska­lo­wa­ny o n po­zy­cji, po­win­ni­śmy licz­bę pier­wiast­ko­wa­ną prze­su­nąć o 2n miejsc dwój­ko­wych. Łatwo może się wtedy oka­zać, że za­brak­nie nam 32 bitów. Co robić? Nic prost­sze­go, za­sto­su­je­my pier­wiast­ko­wa­nie licz­by 64-bi­to­wej z wy­ni­kiem 32 bi­to­wym:

INCLUDE "utility/utility_lib.i"
XREF _UtilityBase

; pierwiastkowanie liczby 64-bitowej. 
; Liczba: starsze długie słowo w d1,
; młodsze w d0. Wynik w d0.

Root64:
    MOVEM.L d2-d5,-(sp)
    MOVE.L  d0,d2
    MOVE.L  d1,d3
    MOVEA.L _UtilityBase,a6
    CLR.L   d4
    MOVEQ   #31,d5
Root64Loop:
    BSET    d5,d4
    MOVE.L  d4,d0
    MOVE.L  d4,d1
    JSR     _LVOUMult64(a6)
    CMP.L   d1,d3
    BCS.S   Clear64
    BNE.S   NoClear64
    CMP.L   d0,d2
    BHI.S   NoClear64
Clear64:
    BCLR    d5,d4
NoClear64:
    DBF     d5,Root64Loop
    MOVE.L  d4,d0
    MOVEM.L (sp)+,d2-d5
    RTS

Jak widać po­słu­ży­łem się tu mno­że­niem 64-bi­to­wym z uti­li­ty.li­bra­ry. To umoż­li­wia dzia­ła­nie pro­ce­du­ry na wszyst­kich pro­ce­so­rach od 68000. Je­że­li ktoś nie ma ocho­ty na otwie­ra­nie bi­blio­te­ki, a pisze pod 020 lub 030, to może wsta­wić za­miast JSR-a „MULU.L d1,d1:d0”. Nie­ste­ty tego roz­ka­zu nie ma w pro­ce­so­rze 060 i chyba też w 040, tak więc UMul­t64() jest roz­wią­za­niem naj­bez­piecz­niej­szym. Ewen­tu­al­nie można roz­pi­sać ręcz­nie to mno­że­nie na roz­ka­zach MULU.W (choć­by prze­pi­su­jąc kod z bi­blio­te­ki), ale nie wiem, czy warto tak wo­jo­wać o jedną parę JSR – RTS. I to tyle tym razem.