ab.com コミュニティ https://www.activebasic.com/forum/ |
|
[AB4]浮動小数点演算ライブラリ https://www.activebasic.com/forum/viewtopic.php?t=2484 |
ページ 1 / 1 |
作成者: | konisi [ 2008年10月12日(日) 19:37 ] |
記事の件名: | [AB4]浮動小数点演算ライブラリ |
暇だったので作ってみました。 定義とか [ここをクリックすると内容が表示されます] コード: '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) #それにしても久しぶりに書いたなぁ・・・ |
ページ 1 / 1 | 全ての表示時間は UTC+09:00 です |
Powered by phpBB® Forum Software © phpBB Limited https://www.phpbb.com/ |