ab.com コミュニティ

コード: 全て選択

'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 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
		Next
		For i=0 To td-1
			vp=0
		Next
	End If
	If tm then
		For i=length-1 To 1 Step -1
			vp=(vp>>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=vp[i+td]
		Next
		For i=length-td-1 To length-1
			vp=0
		Next
	End If
	If tm then
		For i=0 To length-2
			vp=(vp<<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)<=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 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

コード: 全て選択

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)