円周率の計算 (8)

円周率

8 月 20 日の記事で検討した、spigot アルゴリズムを base = 100、つまり 100 進数として扱う場合に 8 ビット演算におさまる条件は、100 進に限らず、他の場合でも有効であることが分かりました。
それは、

  1. num[] の初期値を base にかかわらず「2」とする
  2. 変数 out の初期値も「2」にする
  3. i に関するループは i = 2 で抜ける

ことにより、base が 1000 でも 10000 でも、変数 temp の値が無駄に大きくなることを防げます。
1. の条件から、何進数で扱おうとも整数部は、ひと桁の数「3」として得られることになります。
そこで、今度は「16 ビット縛り」として、PC 上の x86 用 16 ビット DOS アプリケーションとして、spigot アルゴリズムによるプログラムをアセンブリ言語で作成してみました。
COM 形式の実行ファイルの大きさは 228 バイト、配列 num[] のために 64 K バイトのセグメントを使い、小数点以下 9860 桁まで正しい値を出力します。
COM 形式のファイルとするために、16 ビット・セグメント、tiny モデルに設定したアセンブラのソースを下に示します。
TASM と TLINK を使って作成しましたが、「MASM32」パッケージの「ML」と、16 ビット・リンカ「LINK16」を使ってもアセンブルできます。
(8 月 25 日追記: オブジェクトサイズを 247 バイトから 228 バイトに削減したバージョンのプログラムに差し替えました)
spigot86.asm

;******************************************************
;* spigot86.asm -- pi calculation by spigot algorithm *
;*                 using 32 Kword (64 Kbyte) array    *
;*                 with 16-bit arithmetics            *
;*                 9860 decimal digits accuracy       *
;*                                                    *
;* 2011-08-25: object size 247 byte --> 228 byte      *
;* 2011-08-24: Created by pcm1723                     *
;******************************************************
	page	80
	.286		; for PUSHA/POPA
	.model	tiny	; single segment for .COM file
	dosseg		; DOS-style segment order
;
; equates
;
nterms	=	32768	; number of terms
n_step	=	13	; calc. interval
n_lim	=	(54*n_step)	; lower limit for n
base	=	10000	; base 10000 number
;
;  DGROUP Data segment here
;  
	.data
;
dig_cnt	dw	9999h	; digit counter (packed BCD)
;	
;	.stack	; no stack seg for tiny model
;
	.code
;
Start:			; Code start from here
	.startup
; calc segment value for num[] array 
	mov	ax,ds	; AX <-- DS
	add	ax,1000h	; room for 64K .COM
	mov	es,ax	; set ES to point free memory
	call	spigot
DosExit:
	.exit
;	
banner	db	0dh,0ah,'pi = ','$'
;
; Program code
;
spigot	proc
;
	mov	dx,offset banner	; 'pi = '
	mov	ah,9	; DOS string output
	int	21h
;
; register assignments
;
; ES = segment for num[] array
; DS = segment for DGROUP (CODE/DATA/STACK)
;
; AX = low (temp), mul/div work etc.
; BX = (n), low (temp), work
; CX = (i), (denom)
; DX = mul/div work
;
; BP = const. 10000 (base)
; SI = (out), high (temp)
; DI = (num[]) pointer  
;
	std		; set direction backward
	mov	bp,base		; BP = base (10000)
	mov	bx,nterms	; BX = n <-- nterms
	call	IniNump	; initalize num[i] pointer
	mov	ax,2	; initial value
	rep	stosw	; repeat store 
	xchg	ax,si	; SI <-- 2
pi_002:			; for (n) loop
	call	IniNump	; i <-- (n-1)
	push	bx	; save (n)
	push	si	; save (out)
	xor	ax,ax	; clear low (temp)
pi_003:			; for (i) loop
	mul	cx	; DX:AX = low (temp) * i
	xchg	ax,bx	; BX <-- AX, (low temp)
	mov	si,dx	; SI <-- DX, (high temp)
	mov	ax,es:[di]	; AX <- num[i]
	mul	bp	; DX:AX = num[i] * base
	add	ax,bx
	adc	dx,si	; DX:AX += (si:bx); (temp)
;	
	push	cx	; save i (CX)
	shl	cx,1	; i * 2
	dec	cx	; CX = denom <-- (i * 2) - 1
	div	cx	; DX = DX:AX % denom
			; AX = DX:AX / denom
	pop	cx	; restore i (CX)
	xchg	ax,dx	; AX <-- (temp % denom)
	stosw		; num[i] <-- AX
	xchg	ax,dx	; AX <-- low (temp / denom)
	cmp	cx,2	; test for (i == 2)?
	loopne	pi_003	; jump if (i != 2)
;	
	pop	si	; restore (out)
	xor	dx,dx	; DX <-- 0
	div	bp	; DX = DX:AX % base
			; AX = DX:AX / base
	add	ax,si	; AX = out + temp / base
	mov	si,dx	; SI = out = temp % base 
	call	Out4Dec	; output number in AX
;	
	pop	bx	; restore (n)
	sub	bx,n_step	; n -= n_step
	cmp	bx,n_lim	; lower limit ?
	jae	pi_002	; jump if (n >= n_lim)
	ret
spigot	endp	
;
; (i) and num[] pointer initialize from (n)
;
; BX = n (argument)
; CX = i <-- (n - 1)
; DI <-- &(num[i])
; 
IniNump	proc	near
	mov	cx,bx
	dec	cx	; CX = i <-- (n-1)    
	mov	di,cx	; DI <-- i (word index)
	shl	di,1	; convert to byte index
	ret
IniNump	endp
;
; output 4 decimal digits in AX (0000..9999)
;
; volatile: BX
;
Out4Dec	proc	near
	cmp	byte ptr[dig_cnt+1],99h	; the very first digit?
	je	Out1Dec	; output 1 digit (integer part) only
	push	ax	; save AX
	mov	bl,100	; BL <-- 100
	div	bl	; AH <-- AX % 100
			; AL <-- AX / 100
	call	Out2Dec	; output upper 2 digits
	xchg	al,ah	; AL <-- AH
	jmp	short Out2D	; lower 2 digits
Out4Dec	endp
;
; output 2 decimal digits in AL (00..99)
;
Out2Dec	proc	near
	push	ax	; save AX
Out2D	label	near
	aam		; AH <-- AL / 10
			; AL <-- AL % 10
	xchg	al,ah
	call	Out1Dec	; 10's digit
	xchg	al,ah
	jmp	short Out1D	; 1's digit
Out2Dec	endp
;
; output 1 decimal digit in AL (0..9)
;
; volatile: BX
;
Out1Dec	proc	near
	push	ax	; save AX
Out1D	label	near
	add	al,'0'	; conv to ASCII
	call	outc_AL	; output a char
; increment digit count
	mov	ax,[dig_cnt]	; AX <-- digit count
	add	al,1	; increment low byte
	daa		; BCD adjust
	xchg	al,ah	; AL <--> AH
	adc	al,0	; increment high byte
	daa		; BCD adjust	
	xchg	al,ah	; AL <--> AH
	mov	[dig_cnt],ax	; save it
;			; test for delimiter output
	test	al,0fh	; dig_cnt == 'xxx0' ?
	jne	Out1D_4	; jump if not
	mov	bl,' '	; assume delim = ' '
	cmp	ax,0	; the very first digit?
	jne	Out1D_1	; jump if not
;			; integer part	
	inc	ah	; AX = 0100h for 1 newline
	mov	bl,'.'	; decimal point
Out1D_1:	
	xchg	ax,bx	; save AX, AL <-- delimiter char
	call	outc_AL	; output delimiter
	xchg	ax,bx	; restore AX
;			; test for newline output	
	cmp	al,50h	; dig_cnt == 'xx50' ?
	je	Out1D_3	; jump if so
	cmp	al,00h	; dig_cnt == 'xx00' ?
	jne	Out1D_4	; jump if not
	test	ah,0fh	; dig_cnt == 'x000' ?
	jne	Out1D_3	; jump if not
	call	out_nl	; output newline
Out1D_3:	
	call	out_nl	; output newline
Out1D_4:	
	pop	ax	; restore AX
	ret
Out1Dec	endp
;
; output newline (CR / LF)
;
out_nl	proc	near
	pusha
	mov	al,0dh	; CR
	call	outc_AL
	mov	al,0ah	; LF
	jmp	short outc_n
out_nl	endp
;
; output a char in AL
;
outc_AL	proc	near
	pusha
outc_n	label	near	
	mov	dl,al	; DL <-- AL
	mov	ah,2	; DOS out char
	int	21h
	popa
	ret
outc_AL	endp
;
	end

出来上がった COM ファイルを intel HEX 形式に変換したものを下に示します。
spigot86.hex

:100100008CD80500108EC0E80C00B44CCD210D0A2F
:100110007069203D2024BA0E01B409CD21FDBD1027
:1001200027BB0080E84100B80200F3AB96E8380036
:10013000535633C0F7E1938BF2268B05F7E503C3E3
:1001400013D651D1E149F7F15992AB9283F902E00C
:10015000E35E33D2F7F503C68BF2E813005B83EB63
:100160000D81FBBE0273C6C38BCB498BF9D1E7C3AC
:10017000803EE30199741850B364F6F3E8040086F6
:10018000C4EB0150D40A86C4E8040086C4EB0150D5
:100190000430E84300A1E20104012786C4140027CB
:1001A00086C4A3E201A80F7523B32083F800750469
:1001B000FEC4B32E93E82000933C50740C3C0075B1
:1001C0000BF6C40F7503E80500E8020058C360B0E1
:1001D0000DE80400B00AEB01608AD0B402CD2161C1
:0401E000C300999926
:0400000300000100F8
:00000001FF

適当な HEX ファイル・ユーティリティー、たとえば GNU objcopy を使って、

objcopy -I ihex -O binary spigot86.hex spigot86.com

とすれば、バイナリ・ファイルを得ることができます。
出来上がったバイナリは、Windows の「コマンド・プロンプト」上で実行できます。
Windows XP と Windows Vista では確認しましたが、Windows 7 は持っていないので未確認です。
実行時間は、出力を「NUL」にリダイレクトした状態で数秒といったところで、通常の出力を表示させる状態では、出力に要する時間が支配的となります。
実際には問題となる可能性は低いのですが、プログラムをロードして実行可能なコンベンショナル・メモリが 128 K バイト以上ないと正常に実行されません。
配列 num[] のために、16 ビットセグメントの限度いっぱいの 64 K バイトのメモリを使用しているのですが、簡単のため、DOS に対する正式な手続きを取らずに、手抜きの方法で実現しています。
出力結果を (途中省略して) 下に示します。
最後の行の「2162484391」までが正しく、「2199」からは正しく求まっていません。

pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 
5820974944 5923078164 0628620899 8628034825 3421170679 
8214808651 3282306647 0938446095 5058223172 5359408128 
4811174502 8410270193 8521105559 6446229489 5493038196 
4428810975 6659334461 2847564823 3786783165 2712019091 
4564856692 3460348610 4543266482 1339360726 0249141273 
7245870066 0631558817 4881520920 9628292540 9171536436 
7892590360 0113305305 4882046652 1384146951 9415116094 
3305727036 5759591953 0921861173 8193261179 3105118548 
0744623799 6274956735 1885752724 8912279381 8301194912 
9833673362 4406566430 8602139494 6395224737 1907021798 
6094370277 0539217176 2931767523 8467481846 7669405132 
0005681271 4526356082 7785771342 7577896091 7363717872 
1468440901 2249534301 4654958537 1050792279 6892589235 
4201995611 2129021960 8640344181 5981362977 4771309960 
5187072113 4999999837 2978049951 0597317328 1609631859 
5024459455 3469083026 4252230825 3344685035 2619311881 
7101000313 7838752886 5875332083 8142061717 7669147303 
5982534904 2875546873 1159562863 8823537875 9375195778 
1857780532 1712268066 1300192787 6611195909 2164201989 


    <<<<<< 1001 桁目から 9000 桁目まで省略 >>>>>>
    
    <<<<<< 以下は 9001 桁目から 9864 桁まで >>>>>>

9588970695 3653494060 3402166544 3755890045 6328822505 
4525564056 4482465151 8754711962 1844396582 5337543885 
6909411303 1509526179 3780029741 2076651479 3942590298 
9695946995 5657612186 5619673378 6236256125 2163208628 
6922210327 4889218654 3648022967 8070576561 5144632046 
9279068212 0738837781 4233562823 6089632080 6822246801 
2248261177 1858963814 0918390367 3672220888 3215137556 
0037279839 4004152970 0287830766 7094447456 0134556417 
2543709069 7939612257 1429894671 5435784687 8861444581 
2314593571 9849225284 7160504922 1242470141 2147805734 
5510500801 9086996033 0276347870 8108175450 1193071412 
2339086639 3833952942 5786905076 4310063835 1983438934 
1596131854 3475464955 6978103829 3097164651 4384070070 
7360411237 3599843452 2516105070 2705623526 6012764848 
3084076118 3013052793 2054274628 6540360367 4532865105 
7065874882 2569815793 6789766974 2205750596 8344086973 
5020141020 6723585020 0724522563 2651341055 9240190274 
2162484391 2199