定義とか [ここをクリックすると内容が表示されます]
コード: 全て選択
'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)＃それにしても久しぶりに書いたなぁ・・・