倍精度2進浮動小数点数の10進変換 (2012/08/02)

Raspberry Pi は、普通にLinuxとしてできることはなんでもプログラミングできます。 Jun's homepage では Raspberry Pi のハードウェアに特化した部分から試してみたいと思います。候補は GPUを使った OpenGL ES2のプログラミングか、ARM11のコプロセッサであるVFPv2を使った浮動小数点演算をアセンブラで試すのが面白そうです。 VFPv2 はベクトル演算も可能なので色々実験するネタには困りません。

が、何をするにも浮動小数点演算の結果を表示して確認する必要があるので、最初のネタとして倍精度2進浮動小数点数を10進変換して標準出力に表示するサブルーチン(DoubleString)を作成してみました。 DoubleString ルーチンは d0 に設定された倍精度2進浮動小数点数を10進浮動小数点数の文字列に変換して、その文字列(AsciiZ)の先頭アドレスを r0 に返します。標準出力への表示のために、このページで作成したeabi 版の stdio.sを使用しています。

とにかく、Raspberry Pi での私の最初のプログラミングの成果として置いておきます。

@-------------------------------------------------------------------------
@  file : d2str.s 16digit version
@  2012/08/02 Jun Mizutani (https://www.mztn.org/)
@-------------------------------------------------------------------------

.include        "stdio.s"
.text
.global _start
_start:
        adr     r1, D001
    1:
        fldd    d0, [r1]                @ mem --> d0
        add     r1, r1, #8
        bl      DoubleString
        bl      OutAsciiZ
        bl      NewLine
        fcmpezd d0                      @ exit, if d0==0.0
        fmstat
        bne     1b                      @ next
        bl      Exit

@------------------------------------
D001:   .double 1.2345678987654321
        .double -1.2345678987654321
        .double 9.87654E33
        .double -9.87654E-33
        .double 1.2345678987654321E299
        .double -1.2345678987654321E-299
        .double 12345678987654
        .double 12345678987653
        .double 12345678987652
        .double -123456789876543
        .double -123456789876542
        .double -123456789876541
        .double 923456789876543
        .double 923456789876542
        .double 923456789876541
        .double 1.000
        .double 0.1
        .double 0

@-------------------------------------------------------------------------
@ float(d0) to asciiz string
@   r0 returns string address
@-------------------------------------------------------------------------
DoubleString:
        stmfd   sp!, {r1-r6, lr}        @ push
        vpush   {d0-d5}                 @ fstmfdd sp!, {d1}

        mov     r4, #0                  @ plus
        fcmpezd d0                      @ check 0
        fmstat
        beq     6f                      @ d0 = 0.0

        fabsdlt d0, d0                  @ if negative, d0 = abs(d0); r4=1
        movlt   r4, #1

        ldr     ip, numBuf
        cmp     r4, #0                  @ if negative, print '-'s
        mov     r0, #' '
        movne   r0, #'-'
        strb    r0, [ip], #2            @ additional space for decimal point

        @------------------------------
        mov     r6, #0                  @ exponent
        mov     r3, #TblSize
        mov     r2, #512
        fldd    d3, DOne                @ 1.00 --> d3
        fldd    d4, DTen                @ 10.0 --> d4
        fcmpd   d0, d3
        fmstat
        bge     2f                      @ goto 2f, if d0 >= 1.0

    1:  @ negative exponent (f<1)
        sub     r3, r3, #8
        mov     r2, r2, ASR #1          @ start from w=256
        cmp     r2, #0
        beq     4f                      @ range01_1

        adr     r5, NegTable
        add     r5, r5, r3
        fldd    d5, [r5]
        fcmpd   d0, d5                  @ f > 10^-w
        fmstat
        sublt   r6, r6, r2              @ exponent = exponent - w
        fdivdlt d0, d0, d5              @ f = f / 10^-w = f * 10^w
        b       1b

    2:  @ positive (f>=1)
        sub     r3, r3, #8
        mov     r2, r2, ASR #1          @ w
        cmp     r2, #0
        beq     3f                      @ div 10

        adr     r5, PosTable
        add     r5, r5, r3
        fldd    d5, [r5]
        fcmpd   d0, d5                  @ f <= 10^n
        fmstat
        addge   r6, r6, r2              @ exponent = exponent + w
        fdivdge d0, d0, d5              @ f = f / 10^w
        b       2b

    3:  @ div10
        fdivd   d0, d0, d4              @ f = f / 10
        add     r6, r6, #1              @ w = w + 1

        @------------------------------
    4:  @ mantissa should be 0.1 <= f < 1

        bl      VFPRoundTowardZero
        fldd    d2, Exp08
        fmuld   d0, d0, d2
        ftosid  s30, d0                 @ d0 --> integer into s2
        fmrs    r4, s30                 @ set upper 8 digit

        fsitod  d1, s30
        fsubd   d0, d0, d1              @ d0 = d0 - (upper 8)
        fmuld   d0, d0, d2              @ d0 = d0 * 10^8

        bl      VFPRoundToNearest
        ftosid  s30, d0                 @ d0 --> integer into s30
        fmrs    r5, s30                 @ set lower 8 digit
        ldr     r2, digit8
        cmp     r5, r2
        addgt   r4, r4, #1              @ add carry to upper 8
        mov     r0, r4
        bl      _put_uint               @ upper 8 digits
        mov     r0, r5
        bl      _put_uint               @ lower 8 digits

        sub     r0, r6, #1              @ r0 = r6-1
        cmp     r0, #0
        movne   r1, #'E'
        strneb  r1, [ip], #1
        blne    _put_int                @ exponent

    5:
        ldr     r0, numBuf
        ldrb    r1, [r0, +#2]           @ shift left 1 digit for decimal point
        strb    r1, [r0, +#1]
        mov     r1, #'.'
        strb    r1, [r0, +#2]
        vpop    {d0-d5}                 @ fldmfdd sp!, {d1}
        ldmfd   sp!, {r1-r6, pc}        @ pop & return

    6:  @ zero
        ldr     ip, numBuf
        mov     r0, #' '
        strb    r0, [ip], #2            @ additional space for decimal point
        mov     r0, #0
        bl      _put_uint               @ upper 8 digits
        bl      _put_uint               @ lower 8 digits
        b       5b

@==============================================================
numBuf: .long   num_buffer
digit8: .long   99999999

        .align  3
DOne:   .double 1.0
DTen:   .double 10.0
Exp08:  .double 1.0e8

NegTable:
        .double 1.0e-1
        .double 1.0e-2
        .double 1.0e-4
        .double 1.0e-8
        .double 1.0e-16
        .double 1.0e-32
        .double 1.0e-64
        .double 1.0e-128
        .double 1.0e-256
PosTable:
        .double 1.0e+1
        .double 1.0e+2
        .double 1.0e+4
        .double 1.0e+8
        .double 1.0e+16
        .double 1.0e+32
        .double 1.0e+64
        .double 1.0e+128
        .double 1.0e+256
TblSize = . - PosTable

@------------------------------------
@ Output 8 digit unsigned number to buffer
@ r0 : number
@ ip : buffer
_put_uint:
        stmfd   sp!, {r0-r3, lr}        @ push
        mov     r3, #8                  @
        mov     r2, #0                  @ counter
    1:  mov     r1, #10                 @ r1 = 10
        bl      udiv                    @ division by 10
        add     r2, r2, #1              @ counter++
        stmfd   sp!, {r1}               @ least digit (reminder)
        cmp     r2, r3
        bne     1b                      @ done ?
    2:  ldmfd   sp!, {r0}               @ most digit
        add     r0, r0, #'0'            @ ASCII
        strb    r0, [ip], #1            @ store a char into buffer
        subs    r2, r2, #1              @ counter--
        bne     2b
        strb    r2, [ip]                @ store 0 into buffer
        ldmfd   sp!, {r0-r3, pc}        @ pop & return

@------------------------------------
@ Output signed number to buffer
@ r0 : number
@ ip : buffer
_put_int:
        stmfd   sp!, {r0-r2, lr}        @ push
        mov     r2, #0                  @ counter
        cmp     r0, #0
        sublt   r0, r2, r0
        movlt   r1, #'-'
        strltb  r1, [ip], #1
    1:  mov     r1, #10                 @ r1 = 10
        bl      udiv                    @ division by 10
        add     r2, r2, #1              @ counter++
        stmfd   sp!, {r1}               @ least digit (reminder)
        cmp     r0, #0
        bne     1b                      @ done ?
    2:  ldmfd   sp!, {r0}               @ most digit
        add     r0, r0, #'0'            @ ASCII
        strb    r0, [ip], #1            @ store a char into buffer
        subs    r2, r2, #1              @ counter--
        bne     2b
        strb    r2, [ip]                @ store 0 into buffer
        ldmfd   sp!, {r0-r2, pc}        @ pop & return

@-------------------------------------------------------------------------
@ Set VFP rounding mode
@-------------------------------------------------------------------------
VFPRoundTowardZero:
        stmfd   sp!, {r4, r5, lr}
        mov     r4, #0x00c00000         @ rounding mode(11)
        fmrx    r5, FPSCR               @ fpscr --> r5
        orr     r5, r4                  @ Round toward zero(RZ) mode
        fmxr    FPSCR, r5               @ r5    --> FPSCR
        ldmfd   sp!, {r4, r5, pc}

VFPRoundToNearest:
        stmfd   sp!, {r4, r5, lr}
        mov     r4, #0xff3fffff         @ rounding mode(00)
        fmrx    r5, FPSCR               @ fpscr --> r5
        and     r5, r4                  @ Round to nearest(RN) mode
        fmxr    FPSCR, r5               @ r5    --> FPSCR
        ldmfd   sp!, {r4, r5, pc}
@==============================================================
.bss
                .align  2
num_buffer:     .skip   32

実行例

ソース中に埋め込んだ浮動小数点数が正しく変換されるか確認できます。十分な精度で変換できていると(私は)思います。

jun@raspberrypi ~/pi_asm/d2str $ make
as -mcpu=arm1176jzf-s -mfpu=vfpv2 -o d2str.o d2str.s
ld -o d2str d2str.o
jun@raspberrypi ~/pi_asm/d2str $ ./d2str
 1.234567898765432
-1.234567898765432
 9.876540000000000E33
-9.876539999999999E-33
 1.234567898765432E299
-1.234567898765432E-299
 1.234567898765400E13
 1.234567898765300E13
 1.234567898765200E13
-1.234567898765430E14
-1.234567898765420E14
-1.234567898765410E14
 9.234567898765430E14
 9.234567898765421E14
 9.234567898765408E14
 1.000000000000000
 1.000000000000000E-1
 0.000000000000000

Makefile

アセンブルに使ったMakefileです。TABコードに注意。

# Makefile for vfp
#     2012/08/01  Jun Mizutani


OBJS   = d2str.o
TARGET = d2str

ARCH=-mcpu=arm1176jzf-s -mfpu=vfpv2

AS=as
LD=ld
STRIP=strip

.s.o :
        ${AS} ${ARCH} -o $@ $<

all : ${TARGET}

${TARGET} : ${OBJS}
        ${LD} -o $@ ${OBJS}

strip : ${TARGET}
        ${STRIP} ${TARGET}

dist :
        rm *.o

clean :
        rm *.o

「Raspberry Pi で軽く三次元」に続く...