ab.com コミュニティ

ActiveBasicを通したコミュニケーション
現在時刻 - 2017年12月16日(土) 00:41

All times are UTC+09:00




新しいトピックを投稿する  トピックへ返信する  [ 1 件の記事 ] 
作成者 メッセージ
投稿記事Posted: 2008年10月12日(日) 19:37 
オフライン

登録日時: 2005年7月25日(月) 13:27
記事: 893
住所: 埼玉県東松山市
暇だったので作ってみました。[hide=定義とか]
コード:
'1024bit浮動小数点演算ライブラリ
Const B=1024/32'4byte単位で何ワード使うか。変更可能。ただしあまり大きな数にするとalloc-free周りで落ちる。

'計算をする必要がないことを示すフラグ
Const FLAGS_ZERO = &H00000001'0
Const FLAGS_NAN  = &H00000002'非数
Const FLAGS_INF  = &H00000004'無限大

'内部表現
Type T1
	sign As Long'0=正、非0=負
	exp As Long'バイアスはかかってない。基数は2
	flags As Long
	float As *DWord
End Type

Function CreateT1() As *T1'=new T1
	Dim ret As *T1
	ret=malloc(SizeOf(T1))
	ret->float=calloc(B*SizeOf(DWord))
	ret->flags=FLAGS_ZERO
	ret->exp=0
	ret->sign=0
	CreateT1=ret
End Function

Function CreateCopyT1(src As *T1) As *T1'=new T1(src)
	If src=0 then CreateCopyT1=0:Exit Function
	Dim ret As *T1
	ret=CreateT1()
	CopyT1(ret,src)
	CreateCopyT1=ret
End Function

Sub CopyT1(dst As *T1,src As *T1)'dst=src
	If src=0 then Exit Sub
	If dst=0 then Exit Sub
	dst->exp=src->exp
	dst->flags=src->flags
	dst->sign=src->sign
	memcpy(dst->float,src->float,B*SizeOf(DWord))
End Sub

Sub DeleteT1(ByRef t As *T1)'Delete t
	If t then
		free(t->float)
		free(t)
		t=0
	End If
End Sub

'桁あわせをする。桁落ちの危険性がある。主に加算・減算で使う。
Sub f1(r1 As *T1,r2 As *T1)
	If r1=0 then Exit Sub
	If r2=0 then Exit Sub
	If r1->flags then Exit Sub
	If r1->exp=r2->exp then Exit Sub
	If r1->exp>r2->exp then
		Shr(r2->float,r1->exp-r2->exp,B)
		r2->exp=r1->exp
	Else
		Shr(r1->float,r2->exp-r1->exp,B)
		r1->exp=r2->exp
	End If
End Sub

'正規化する。
Sub f2(r1 As *T1)
	If r1=0 then Exit Sub
	If r1->flags then Exit Sub
	If r1->float[0]=1 then
		Exit Sub
	ElseIf r1->float[0]>1 then
		While r1->float[0]<>1
			Shr(r1->float,1,B)
			r1->exp=r1->exp+1
			If r1->exp=&H80000000 then r1->flags=r1->flags Or FLAGS_INF
		Wend
	Else
		Dim i As Long
		For i=0 To B-1
			If r1->float[i] then Exit For
		Next
		If i=B then
			r1->flags=r1->flags Or FLAGS_ZERO
		Else
			While r1->float[0]<>1
				Shl(r1->float,1,B)
				r1->exp=r1->exp-1
				If r1->exp=&H7fffffff then r1->flags=r1->flags Or FLAGS_ZERO
			Wend
		End If
	End If
End Sub

'lengthは4バイト単位での指定。sbitで指定された数だけvp配列全体を右シフト
Sub Shr(vp As *DWord,sbit As Long,length As Long)
	If vp=0 then Exit Sub
	Dim td As Long,tm As Long
	Dim i As Long
	td=sbit\32
	tm=sbit mod 32
	If td then
		For i=length-td-1 To 0 Step -1
			vp[i+td]=vp[i]
		Next
		For i=0 To td-1
			vp[i]=0
		Next
	End If
	If tm then
		For i=length-1 To 1 Step -1
			vp[i]=(vp[i]>>tm)+(vp[i-1]<<(32-tm))
		Next
		vp[0]=vp[0]>>tm
	End If
End Sub

'lengthは4バイト単位での指定。sbitで指定された数だけvp配列全体を左シフト
Sub Shl(vp As *DWord,sbit As Long,length As Long)
	If vp=0 then Exit Sub
	Dim td As Long,tm As Long
	Dim i As Long
	td=sbit\32
	tm=sbit mod 32
	If td then
		For i=0 To length-td-1
			vp[i]=vp[i+td]
		Next
		For i=length-td-1 To length-1
			vp[i]=0
		Next
	End If
	If tm then
		For i=0 To length-2
			vp[i]=(vp[i]<<tm)+(vp[i+1]>>(32-tm))
		Next
		vp[length-1]=vp[length-1]<<tm
	End If
End Sub

Sub plus(r1 As *T1,r2 As *T1)'加算。r1=r1+r2
	If r1=0 then Exit Sub
	If r2=0 then Exit Sub
	Dim t As *T1
	If (r1->flags Or r2->flags) And (FLAGS_NAN Or FLAGS_INF) then
		r1->flags=r1->flags Or r2->flags
	ElseIf r1->flags And FLAGS_ZERO then
		CopyT1(r1,r2)
		r1->flags=(r1->flags And not FLAGS_ZERO) Or (r2->flags And FLAGS_ZERO)
	ElseIf r2->flags And FLAGS_ZERO then
		'nop
	Else
		If (r1->sign=0 And r2->sign) Or (r1->sign And r2->sign=0) then
			t=CreateCopyT1(r2)
			t->sign=t->sign=0
			minus(r1,t)
			DeleteT1(t)
		Else
			t=CreateCopyT1(r2)
			f1(r1,t)
			Dim i As Long,c=0 As Long
			For i=B-1 To 0 Step -1
				If c then
					c=(not r1->float[i])<=t->float[i]
					r1->float[i]=r1->float[i]+t->float[i]+1
				Else
					c=(not r1->float[i])<t->float[i]
					r1->float[i]=r1->float[i]+t->float[i]
				End If
			Next
			f2(r1)
			DeleteT1(t)
		End If
	End If
End Sub

Sub minus(r1 As *T1,r2 As *T1)'減算。r1=r1-r2
	If r1=0 then Exit Sub
	If r2=0 then Exit Sub
	Dim t As *T1
	If (r1->flags Or r2->flags) And (FLAGS_NAN Or FLAGS_INF) then
		r1->flags=r1->flags Or r2->flags
	ElseIf r1->flags And FLAGS_ZERO then
		CopyT1(r1,r2)
		r1->sign=r2->sign=0
		r1->flags=(r1->flags And not FLAGS_ZERO) Or (r2->flags And FLAGS_ZERO)
	ElseIf r2->flags And FLAGS_ZERO then
		'nop
	Else
		If (r1->sign=0 And r2->sign) Or (r1->sign And r2->sign=0) then
			t=CreateCopyT1(r2)
			t->sign=t->sign=0
			plus(r1,t)
			DeleteT1(t)
		Else
			t=CreateCopyT1(r2)
			f1(r1,t)
			Dim i As Long,c=0 As Long
			For i=B-1 To 0 Step -1
				If c then
					c=r1->float[i]<=t->float[i]
					r1->float[i]=r1->float[i]-t->float[i]-1
				Else
					c=r1->float[i]<t->float[i]
					r1->float[i]=r1->float[i]-t->float[i]
				End If
			Next
			If c then
				r1->sign=r1->sign=0
				For i=0 To B-1
					r1->float[i]=not r1->float[i]
				Next
				c=r1->float[B-1]=-1
				For i=B-1 To 0 Step -1
					If c then
						c=r1->float[i]=-1
						r1->float[i]=r1->float[i]+1
					Else
						Exit For
					End If
				Next
			End If
			f2(r1)
			DeleteT1(t)
		End If
	End If
End Sub

Sub multiply(r1 As *T1,r2 As *T1)'乗算。r1=r1*r2
	If r1=0 then Exit Sub
	If r2=0 then Exit Sub
	Dim exp As Long,flags As Long,sign As Long
	Dim ml As DWord,mh As DWord
	Dim tmp[B] As DWord,res[2*B-1] As DWord
	Dim i As Long,j As Long,c=0 As Long
	f2(r1)
	If (r1->sign And r2->sign=0) Or (r1->sign=0 And r2->sign) then
		sign=-1
	Else
		sign=0
	End If
	exp=r1->exp+r2->exp
	flags=r1->flags Or r2->flags
	If flags then
		r1->flags=flags
	Else
		FillMemory(res,2*B*SizeOf(DWord),0)
		For i=B-1 To 0 Step -1
			FillMemory(tmp,(B+1)*SizeOf(DWord),0)
			For j=B-1 To 0 Step -1
				ml=r1->float[i]*r2->float[j]
				mh=(r1->float[i] As QWord*r2->float[j])>>32
				c=(not tmp[j+1])<ml
				tmp[j+1]=tmp[j+1]+ml
				tmp[j]=mh-c
			Next
			c=0
			For j=B To 0 Step -1
				If c then
					c=(not res[i+j])<=tmp[j]
					res[i+j]=res[i+j]+tmp[j]+1
				Else
					c=(not res[i+j])<tmp[j]
					res[i+j]=res[i+j]+tmp[j]
				End If
			Next
			Do
				If i+j<0 then Exit Do
				If c then
					c=res[i+j]=-1
					res[i+j]=res[i+j]+1
				Else
					Exit Do
				End If
				j=j-1
			Loop
		Next
		memcpy(r1->float,res+4,B*SizeOf(DWord))
		r1->exp=exp
		r1->sign=sign
		f2(r1)
	End If
End Sub

Sub divid(r1 As *T1,r2 As *T1)'除算。r1=r1/r2
	If r1=0 then Exit Sub
	If r2=0 then Exit Sub
	If r2->flags then Exit Sub
	Dim t As *T1
	t=CreateT1()
	ReciprocalT1(t,r2)
	multiply(r1,t)
	DeleteT1(t)
End Sub

Sub ReciprocalT1(dst As *T1,src As *T1)'逆数を取る。dst=1/src
	If src=0 then Exit Sub
	If src->flags then Exit Sub
	Dim x As *T1,xd As *T1,i As Long
	Dim _2 As *T1
	xd=CreateT1()
	x=CreateT1()
	If src->float[1]<&H6A09E668 then
		x->exp=-src->exp
	Else
		x->exp=-src->exp-1
	End If
	x->sign=src->sign
	x->float[0]=1
	x->flags=0
	CopyT1(xd,x)
	_2=CreateT1FromDouble(2)
	'X[n+1]=X[n]*(2-a*X[n])→lim[n→∞]X[n]=1/a
	For i=1 To 10'1024bitの場合は10。4096bitの場合は12程度。
		multiply(xd,src)
		minus(xd,_2)
		xd->sign=xd->sign=0
		multiply(xd,x)
		CopyT1(x,xd)
	Next
	CopyT1(dst,x)
	DeleteT1(x)
	DeleteT1(xd)
	DeleteT1(_2)
End Sub

Function CreateT1FromSingle(src As Single) As *T1'IEEE 754
	's-e-f=1-8-23 b=127
	Dim _src As *DWord,exp As Long
	Dim dst As *T1
	_src=VarPtr(src)
	dst=CreateT1()
	dst->flags=0
	dst->sign=_src[0]>>31
	exp=(_src[0]>>23) And &Hff
	If exp=0 then
		If _src[0] and &H007fffff then
			dst->exp=exp-126
			dst->float[0]=_src[0] and &H007fffff
			Shr(dst->float,23,3)
			f2(dst)
		Else
			dst->flags=FLAGS_ZERO
			dst->float[0]=0
		End If
	ElseIf exp=255 then
		If _src[0] and &H007fffff then
			dst->flags=FLAGS_NAN
		Else
			dst->flags=FLAGS_INF
		End If
	Else
		dst->exp=exp-127
		dst->float[0]=(_src[0] and &H007fffff)+&H00800000
		Shr(dst->float,23,3)
	End If
	CreateT1FromSingle=dst
End Function

Function GetSingleFromT1(src As *T1) As Single'IEEE 754
	If src=0 then GetSingleFromT1=0:Exit Function
	Dim dst As *DWord,exp As Long
	dst=VarPtr(GetSingleFromT1)
	If src->flags And FLAGS_INF then
		dst[0]=&H7f800000
	ElseIf src->flags And FLAGS_NAN then
		dst[0]=&H7fffffff
	ElseIf src->flags And FLAGS_ZERO then
		dst[0]=&H00000000
	Else
		exp=src->exp+127
		If exp>=255 then'無限大
			dst[0]=&H7f800000
		ElseIf exp=<0 then
			exp=src->exp+126
			If exp>=-23 then
				dst[0]=&H00400000+(src->float[1]>>9)
				Do
					exp=exp+1
					dst[0]=dst[0]>>1
					If exp=0 then Exit Do
				Loop
			Else
				dst[0]=&H00000000
			End If
		Else
			dst[0]=(exp<<23)+(src->float[1]>>9)
		End If
	End If
	dst[0]=dst[0] Or ((src->sign<>0)<<31)
End Function

Function CreateT1FromDouble(src As Double) As *T1'IEEE 754
	's-e-f=1-11-52 b=1023
	Dim _src As *DWord,exp As Long
	Dim dst As *T1
	_src=VarPtr(src)
	dst=CreateT1()
	dst->flags=0
	dst->sign=_src[1]>>31
	exp=(_src[1]>>20) And &H7ff
	If exp=0 then
		If (_src[1] and &H000fffff) or _src[0] then
			dst->exp=-1022
			dst->float[0]=_src[1] and &H000fffff
			dst->float[1]=_src[0]
			Shr(dst->float,20,3)
			f2(dst)
		Else
			dst->flags=FLAGS_ZERO
			dst->float[0]=0
		End If
	ElseIf exp=2047 then
		If (_src[1] and &H000fffff) or _src[0] then
			dst->flags=FLAGS_NAN
		Else
			dst->flags=FLAGS_INF
		End If
	Else
		dst->exp=exp-1023
		dst->float[0]=(_src[1] and &H000fffff)+&H00100000
		dst->float[1]=_src[0]
		Shr(dst->float,20,3)
	End If
	CreateT1FromDouble=dst
End Function

Function GetDoubleFromT1(src As *T1) As Double'IEEE 754
	If src=0 then GetDoubleFromT1=0:Exit Function
	Dim dst As *DWord,exp As Long
	dst=VarPtr(GetDoubleFromT1)
	If src->flags And FLAGS_INF then
		dst[1]=&H7ff00000
		dst[0]=&H00000000
	ElseIf src->flags And FLAGS_NAN then
		dst[1]=&H7fffffff
		dst[0]=&Hffffffff
	ElseIf src->flags And FLAGS_ZERO then
		dst[1]=&H00000000
		dst[0]=&H00000000
	Else
		exp=src->exp+1023
		If exp>=2047 then'無限大
			dst[1]=&H7ff00000
			dst[0]=&H00000000
		ElseIf exp=<0 then
			exp=src->exp+1022
			If exp>=-52 then
				dst[1]=&H00100000+(src->float[1]>>12)
				dst[0]=(src->float[1]<<20)+(src->float[2]>>12)
				Do
					exp=exp+1
					dst[0]=(dst[1]<<31)+(dst[0]>>1)
					dst[1]=dst[1]>>1
					If exp=0 then Exit Do
				Loop
			Else
				dst[1]=&H00000000
				dst[0]=&H00000000
			End If
		Else
			dst[1]=(exp<<20)+(src->float[1]>>12)
			dst[0]=(src->float[1]<<20)+(src->float[2]>>12)
		End If
	End If
	dst[1]=dst[1] Or ((src->sign<>0)<<31)
End Function
[/hide][hide=使用例]
コード:
Dim a As *T1,b As *T1
Dim v As Double

a=CreateT1FromDouble(1.7)
b=CreateT1FromDouble(0.4)
plus(a,b)
v=GetDoubleFromT1(a)
DeleteT1(a)
DeleteT1(b)

a=CreateT1FromDouble(1.7)
b=CreateT1FromDouble(0.4)
minus(a,b)
v=GetDoubleFromT1(a)
DeleteT1(a)
DeleteT1(b)

a=CreateT1FromDouble(1.7)
b=CreateT1FromDouble(0.4)
multiply(a,b)
v=GetDoubleFromT1(a)
DeleteT1(a)
DeleteT1(b)

a=CreateT1FromSingle(1.7)
v=GetSingleFromT1(a)
DeleteT1(a)

a=CreateT1FromDouble(1.7)
b=CreateT1FromDouble(0.4)
divid(a,b)
v=GetDoubleFromT1(a)
DeleteT1(a)
DeleteT1(b)
[/hide]型名他で適当な名前がついてるのは気にしない事。

#それにしても久しぶりに書いたなぁ・・・

_________________
Website→http://web1.nazca.co.jp/himajinn13sei/top.html
ここ以外の場所では「暇人13世」というHNを主として使用。

に署名を書き換えて欲しいと言われたので暇だしやってみるテスト。


通報する
ページトップ
期間内表示:  ソート  
新しいトピックを投稿する  トピックへ返信する  [ 1 件の記事 ] 

All times are UTC+09:00


オンラインデータ

このフォーラムを閲覧中のユーザー: なし & ゲスト[1人]


トピック投稿:  可
返信投稿:  可
記事編集: 不可
記事削除: 不可
ファイル添付: 不可

検索:
ページ移動:  
cron
Powered by phpBB® Forum Software © phpBB Limited
Japanese translation principally by KONISHI Yohsuke