小さなねた

HSP を使う上でちょっと役に立つかもしれない技みたいなものを紹介します。
便利そうなものも、不必要そうなものもとりあえず紹介してみます。

誤りを見つけられた方は、BBS かメールで報告していただけるとありがたいです。

(注) このページの内容は、HSP 2.x を対象とした古めの内容になってきています。
HSP 3.0 向けの情報は、小さなねた 3 で扱っています。

2004/01/01 HSP でアドレス変数を作って使ってみよう
2003/07/05 アークタンジェントを求める
2002/08/07 10bit 精度の sin テーブル
2002/01/09 初版アップ

今回は、HSP でアドレス変数を作って使ってみようというねたです。

書いてるうちに長くなったので別ページにしました。

HSP でアドレス変数を作って使ってみよう

HSP の標準命令には、三角関数を求める命令がありません。
が、sin (cos) の値はマクローリン展開を使えばわりと簡単に求めることができます。

8bit 精度のものはどこかで見かけたので、ここでは小数部 10bit 精度の sin テーブルを作成するスクリプトを紹介します。

	dim sin, 1025 ; テーブル用配列

	sin.256 =  1024
	sin.768 = -1024

	repeat 256
		r = cnt << 23 / 166886 ; ラジアンに変換
		sin.cnt = r*r>>10*r>>13/3*r>>15*r>>17/5+r-(r*r>>10*r>>17/3)-(r*r>>10*r>>13/3*r>>15*r>>13/5*r>>13/6*r>>17/7)+4>>3

		; 位相をずらしてコピー
		a = 512 - cnt
		sin.a = sin.cnt
		a = 512 + cnt
		sin.a = -sin.cnt
		a = 1024 - cnt
		sin.a = -sin.cnt

		loop

	screen 0, 1024, 512, 1
	repeat 1024 ; 波形を表示
		line cnt, (sin.cnt + 1024) / 4
		loop
	stop

r = cnt << 23 / 166886

この式が、1024 で一周する角度単位を 13bit 固定小数のラジアンに変換するものです。

sin.cnt = r*r>>10*r>>13/3*r>>15*r>>17/5+r-(r*r>>10*r>>17/3)-(r*r>>10*r>>13/3*r>>15*r>>13/5*r>>13/6*r>>17/7)+4>>3

この式が、マクローリン展開の式です。32bit int で最高の精度で演算させようと最適化したため、わけの分からない式になってしまいました。
これで小数部 13bit 精度の sin を求め、それを四捨五入して 10bit にすることで滑らかな sin カーブを得ています。

マクローリン展開自体は簡単でも、なめらかな曲線を得るには演算時の精度を高める工夫が必要でした。

さて、先ほど求めた sin テーブルは、アークタンジェントの計算にも使えます。

以下が、実際のアークタンジェント計算ルーチンです。マウスでぐるぐる動きます。
まあ、入力が (y/x) ではなく (x, y) なので厳密には「アークタンジェント応用関数」ですが、ここではアークタンジェントと呼ぶことにします。

スクリプトは以下のとおり。
(長いので、コピーする場合は右クリックで すべて選択 を選んでください)

	screen 0, 512, 512

	; さっきの sin テーブル計算ルーチン
	dim sin, 1025
	sin.256 =  1024
	sin.768 = -1024
	repeat 256
		r = cnt << 23 / 166886 ; ラジアンに変換
		sin.cnt = r*r>>10*r>>13/3*r>>15*r>>17/5+r-(r*r>>10*r>>17/3)-(r*r>>10*r>>13/3*r>>15*r>>13/5*r>>13/6*r>>17/7)+4>>3
		; 位相をずらしてコピー
		a = 512 - cnt
		sin.a = sin.cnt
		a = 512 + cnt
		sin.a = -sin.cnt
		a = 1024 - cnt
		sin.a = -sin.cnt
		loop


	; line の中心座標
	cx = winx/2
	cy = winy/2

	repeat
		; マウス座標と中心座標の差
		v1 = mousex - cx
		v2 = mousey - cy

		; ここから アークタンジェントを求める
		; 入力 x, y = v1, v2
		; 出力 atan

		; 場合わけ
		if v1 + v2 > 0 {
			if v1 > v2 {
				x = v1
				y = v2
				a = 1
				c = 0
			} else {
				x = v2
				y = v1
				a = -1
				c = 256
			}
		} else {
			if v1 < v2 {
				x = -v1
				y = -v2
				a = 1
				c = 512
			} else {
				x = -v2
				y = -v1
				a = -1
				c = 768
			}
		}

		; 二分探索でアークタンジェントを求める
		if x {
			alp = 128
			bet = -128

			repeat 8
				ctr = alp + bet / 2

				tx = ctr + 256 & 1023
				ty = ctr & 1023

				if sin.ty<<10/sin.tx >= (y<<10/x) {
					alp = ctr
				} else {
					bet = ctr
				}
				loop

			atan = alp * a + c & 1023 ; 出力値
		}

		; 解が求められない場合
		else {
			atan = 0
		}

		; 描画
		redraw 0
		palcolor 255:boxf
		color

		x = atan + 256 & 1023
		y = atan

		line sin.x/4+cx, sin.y/4+cy, cx, cy

		redraw

		wait 1
		loop

これは、積分などを使わないで、二分探索でテーブルから数値を探す計算法です。

まず、入力 x, y による 4 パターンの場合分けを行い、探索ループへの入力を適切な値に変換し、出力とを加工するためのパラメータを設定します。
探索ループでは、sin テーブルを使い、試行中の角度 (ctr) の tan と入力 x, y の tan を比較します。8 回の比較で -128 〜 128 の解を求めます。それを先の場合分けによって加工し、出力値とします。

角度値の精度は、sin テーブルと同様一周を 1024 分割したものになります。
ang = (ANS * 2π / 1024) [rad]

HSP の標準命令には、平方根や指数の計算をする命令がありません。
正零整数の指数の計算はループを作って掛け算をさせれば簡単にできますが、ルートの計算はどうやればいいのでしょうか。

関数電卓で遊んでるときに思いつきました。

ある正の数 a の平方根 r を求めるのであれば、はじめに r の値を 0 とした後で r の最上位の桁から値を 9〜0 と変化させ、r*r <= a となったところでその桁の値を確定し、さらに小さい桁について同様の計算を行っていけば、ルートの値を絞り込んで求めることができるはずです。

それをビットシフト等を使って二進数で実行するスクリプトは、以下にようになります。

	a = 49 ; 平方根を求める数

	r = 0
	repeat 15 ; 平方根を求める ;出力の最大値:32767 (15bit)
		b = 16384 >> cnt | r
		if b * b <= a : r = b
		loop

	mes "√"+a+" = "+r ; 計算結果
	stop

上のスクリプトでは 15 回のループで 15 ビットの値を求めるため、出力される r の最大値は 32767 (111111111111111B) となります。
これは、符号付 32bit 整数値の最大値 2147483647 (7fffffffH) の平方根の値 46340 (.95…) よりも小さくなってしまいます。

ところが、上のスクリプトでループ回数を単に 16 回に増やしたのでは "if b*b<=a : r=b" の行で b*b の値が 32bit int の最大値を超える場合が出るため、計算結果がおかしくなります。

そこで、以下のスクリプトではその行に少し改良を加えて 16bit 対応にし、モジュール化して使い勝手を良くしました。出力される最大の値は 46340 です。

#module
; sqrt (結果が代入される変数), (ルートを求める正の数値)
#deffunc sqrt val, int
	mref v1, 16
	mref v2, 1

	v1 = 0
	repeat 16
		a = $8000 >> cnt | v1
		if a * a <= v2 & (a <= 46340) : v1 = a
		loop
	return
#global

	randomize

	repeat 24 * 4
		rnd a, 2000 : rnd a, a + 1

		sqrt b, a * 1000 * 1000

		c = b \ 1000 : str c, 3
		b = b / 1000
		pos cnt / 24 * 160, cnt \ 24 * 20
		mes "√"+a+" = "+b+"."+c
		loop
	stop

このスクリプトを応用すれば、三乗根の計算なども行うことができます。

; 三乗根
#deffunc cbrt val, int
	mref v01, 16
	mref v02, 1

	v01 = 0
	repeat 11
		a = 1024 >> cnt | v01
		if a * a * a <= v02 & (a <= 1290) {
			v01 = a
		}
		loop

	return

本を一冊ずつ積み上げてそれをまた上から取り出していくように、次々と入力されて蓄えられたデータが、入力された順番と逆の順番に取り出されるような記憶方式を、FILO (First In Last Out) もしくは LIFO (Last In First Out) 記憶というらしいです。
たとえばアセンブラで出てくる PUSH, POP 命令がそれなわけですが、その記憶方式を HSP で真似してみようと。

以下のスクリプトで、モジュール内の変数 stack がスタックインデックス兼スタック領域となっています。
stack.0 が現在のスタックのインデックスを表し、stack.1〜15 が入力されたデータの記憶に使われます。

#module

; push (入力する値)
#deffunc push int
	mref v1
	stack++
	stack.stack = v1
	return

; pop (値を取り出す変数)
#deffunc pop val
	mref v1, 16
	v1 = stack.stack
	stack--
	return
#global

	a = 100	; 元の値
	b = 200

	push a	; a の値をスタックに待避
	push b	; b の値をスタックに待避

	a = 0	; 値を破壊してみる
	b = 0

	pop b	; スタックに待避していた値を取り出す
	pop a

	mes a	; ほうら、元通り
	mes b

	stop

HSP スクリプトで使う場合、通常はこの記憶方式に有効な使い道はなさそうです。

しかし、再帰的なルーチンを組む場合、ローカル変数の使えない HSP ではこの FILO 記憶方式が役に立ちそうです。