/*                                  */ /* IBM形式単精度浮動少数点シミュレーション */ /* */ /* 符号フラグ(1 bit) */ /*   | */ /* | 指数部(16^nとし負は2の補数:7bits) */ /* | */ /*   | 仮数部(24 bits) */ /*   ↓ */ /* +-+−−−−−−−+−−−−−−−−−−−−−−−−−−−+ */ /* | | 指数部    | 仮数部(24 bits/ 3 bytes) | */ /* +-+−−−−−−−+−−−−−−−−−−−−−−−−−−−+ */ /* */ #include "stdafx.h" #include "math.h" #include "stdlib.h" typedef unsigned long myFloat; /*-----------------------------------------------------------------------------*/ /* ■フィールド分割/結合 */ void sepField(myFloat DT, int *Flag, int *Exp, myFloat *Coeff){ *Coeff = DT & 0X00FFFFFF; int ExpDT = DT >> 24; if((ExpDT & 0X40) ==0) *Exp = ExpDT & 0x3F; else *Exp = ExpDT | 0xFFFFFFC0; *Flag= ExpDT>>7; } myFloat comField(int Exp, myFloat Coeff){/*符号以外をセット*/ return (((long)(Exp & 0x7F))<<24) | (Coeff & 0xFFFFFF); } /*-----------------------------------------------------------------------------*/ /* ■共通ルーチン */ void normalize(myFloat *C, int *exp){/*正規化(仮数部先頭4ビット≠0にする) */ if((*C)!=0) while((*C & 0xF00000)==0) { *C<<=4; (*exp)--;} } myFloat myAbs(myFloat A) { return A & 0x7FFFFFFF;}/*絶対値*/ myFloat myMinuss(myFloat A){ return A ^ 0x80000000;}/*マイナス処理*/ /*-----------------------------------------------------------------------------*/ /* ■絶対値の大きさの比較       */ /*       */ /*  @ 仮数部0のとき(値0), 仮数部0でない方を大きいとする。           */ /*  A 両者ともに 0 のとき, 等しいとする。             */ /*  B 両者ともに非ゼロのとき, 指数部が同じであれば仮数部で比較。 */ /*   仮数部の場合,減算結果がshort intに収まらない可能性があることに注意。 */ /*   指数部が異なれば指数部で比較する(減算する)。       */ /*                                       */ int absCmp(myFloat A, myFloat B){/*絶対値の大きさ比較 */ int flagA,flagB, expA, expB ;myFloat cA, cB; sepField(A,&flagA,& expA,&cA); sepField(B,&flagB,& expB,&cB); if(cA==0){if(cB==0) return 0;else return -1;}else if(cB==0) return 1; normalize(&cA, &expA); normalize(&cB,&expB);/*正規化*/ return ((expA==expB) ? ((cA>cB)? 1: ((cA < cB) ? -1: 0)) : expA-expB); } /*-----------------------------------------------------------------------------*/ /* ■絶対値加算 */ /* +−+−+−+−+−+−+ */ /* | | | | | | |   A仮数部 */ /* +−+−+−+−+−+−+ */ /* | +−+−+−+−+−+−+ */ /* | | | | | |…|…| B仮数部(…は無視) */ /* | +−+−+−+−+−+−+ */ /* ・---→|        pow(2,4) =16           */ /* | |(Aの指数部−Bの指数部)×4ビット右シフトして加算 */ /* */ /* 加算すると上位ビットがあふれる可能性がある */ /* +−+−+−+−+−+−+−+ (例) 0x0F+0x0F=0x1E */ /* | 1| 2| 3| 4| 5| 6| 7|   あふれた場合,右に4ビット */ /* +−+−+−+−+−+−+−+ シフトして指数を1増やす */ /* \ \ \ \ \ \ \ \            */ /* +−+−+−+−+−+−+−+          */ /* | 1| 2| 3| 4| 5| 6|…|   (…は無視) */ /* +−+−+−+−+−+−+−+ */ /* なお4ビットの桁落ちを少なくするため16進7桁にして計算する。 */ /*                                  */ myFloat absAdd(myFloat A, myFloat B){/*絶対値加算*/ int flagA,flagB,expA, expB, expC, dE; myFloat cA, cB, cC; sepField(A,&flagA,& expA,&cA); sepField(B,&flagB,& expB,&cB); if(cA==0) return B; else if(cB==0) return A; cA<<=4;cB<<=4;/*16進7桁にする*/ dE = expA - expB; expC = expA;/*桁数を合わせる*/ if(dE >0) cB>>=(dE<<2); else if(dE<0){cA>>=((-dE)<<2); expC=expB;} /*加算と桁あふれのとき4bits右シフト*/ cC = cA+cB; if((cC & 0xF0000000)!=0){expC++; cC>>=4;} return comField(expC,cC>>4);/*16進6桁に戻して設定*/ } /*-----------------------------------------------------------------------------*/ /* ■絶対値減算 */ /* +−+−+−+−+−+−+ */ /* | | | | | | |   A仮数部 */ /* +−+−+−+−+−+−+ */ /* | +−+−+−+−+−+−+ */ /* | | | | | |…|…| B仮数部(…は無視) */ /* | +−+−+−+−+−+−+ */ /* ・---→|        pow(2,4) =16           */ /* | |(Aの指数部−Bの指数部)×4ビット右シフトして減算 */ /* なお,桁落ち対策用に16進8桁にして調整している。 */ /* 減算した結果上位が0になる可能性がある */ /* +−+−+−+−+−+−+ */ /* | 0| 1| 2| 3| 4| 5| 0の個数×4ビット左シフトして */ /* +−+−+−+−+−+−+ シフトして指数を1減らす */ /* / / / / / /            */ /* +−+−+−+−+−+−+          */ /* | 1| 2| 3| 4| 5| 0|   (最後に0が詰まる) */ /* +−+−+−+−+−+−+ */ /*                                  */ myFloat absSub(myFloat A, myFloat B){/*絶対値減算*/ int flagA,flagB,expA, expB, expC, dE, cmp; myFloat cA, cB, cC; if((cmp=absCmp(A,B))==0) return 0;/*|A|>|B|であることを保証する*/ if(cmp<0){ cC = A; A = B; B = cC;} sepField(A,&flagA,& expA,&cA); sepField(B,&flagB,& expB,&cB); if(cB==0) return A;else if(cA==0) return B; normalize(&cA, &expA); normalize(&cB,&expB);/*正規化*/ cA<<=8; cB<<=8;/*桁落ちを少なくするために*/ dE=expA-expB;expC=expA; cB>>=(dE<<2); /*桁合わせ。|A|>|B|なので,必ずdE≧0*/ cC=(cA-cB); if ((cC & 0xFF000000)==0) expC-=2;/*桁落ち対策*/ else if((cC & 0xF0000000)==0) { expC--;cC>>=4;} else cC>>=8; normalize(&cC,&expC); return comField(expC,cC); } /*-------------------------------------------------------------------------------*/ /* ■絶対値乗算 */ /* +−+−+−+−+−+−+ */ /* | 1| 2| 3| 4| 5| 6|   仮数部を2バイトずつに分ける */ /* +−+−+−+−+−+−+ */ /* / / / \ \ \ \ \            */ /* +−+−+−+−+−+−+ +−+−+−+−+−+−+ */ /* | 0| 0| 0| 0| 1| 2| | 0| 0| 3| 4| 5| 6|                 */ /* +−+−+−+−+−+−+ +−+−+−+−+−+−+ */ /* |←-------------→| |←-------------→|            */ /* | BB[1]    | |   BB[0] | B */ /* | AA[1] | | AA[0] | A */ /* */ /* @ CC[0], CC[1], CC[2]をクリアする。 */ /* A i=0, 1に対してB〜Dを行う。       */ /*   B 下位あふれ値(以下ではAB)を0にする。 */ /* C j=0, 1に対してDを行う。       */ /*   D AB += AA[i]*BB[j]+CC[i+j]とし, */ /* 新しいCC[i+j]はABの下位16ビットとし, */ /*   次のあふれ値ABは上位16ビットとする(右に16ビットシフト)。 */ /* E 最後に残ったあふれ値ABは必ず0となる(少数部同士の乗算は1未満)。 */ /* F CC[2]の16ビットとCC[1]の8ビットとすることができるが, */ /* CC[2]の16ビット先頭4ビットが0となる場合があるので,        */ /* そのとき4ビット左にシフトし,指数部を1減らす。 */ /*                                      */ /*  初期値    CC[0],CC[1],CC[2],CC[3]=0                 */ /*                                      */ /*  i=0 のとき CC[0]=AA[0]*BB[0]+CC[0] の下位16ビット            */ /*   CC[1]=AA[0]*BB[1]+CC[1]+(CC[0]計算時のあふれ) の下位16ビット */ /*   CC[2]=(CC[1]計算時のあふれ) (注)0でないことがありうる. */ /*                                       */ /*  i=1 のとき CC[1]=AA[1]*BB[0]+CC[1] の下位16ビット        */ /*   CC[2]=AA[1]*BB[1]+CC[2]+(CC[1]計算時のあふれ) の下位16ビット */ /*         CC[3]=(CC[2]計算時のあふれ)(注)ただし,これは必ず0.     */ /*                                       */ myFloat absMult(myFloat A, myFloat B){/*絶対値乗算*/ int flagA,flagB,expA, expB, expC; myFloat cA, cB, cC; long AA[2],BB[2],CC[4],AB; sepField(A,&flagA,& expA,&cA); if(cA==0) return 0; sepField(B,&flagB,& expB,&cB); if(cB==0) return 0; normalize(&cA, &expA); normalize(&cB,&expB);/*正規化*/ AA[0]=cA & 0xFFFF; AA[1]=cA>>16; BB[0]=cB & 0xFFFF; BB[1]=cB>>16; for(int i=0;i<4;i++) CC[i]=0; for(int i=0;i<2;i++){ /*乗算実行 あふれがABに残る*/ AB=0; for(int j=0; j<2; AB+=AA[i]*BB[j]+CC[i+j], CC[i+j]=AB & 0xFFFF, AB>>=16, j++); CC[i+2]=AB;/*CC[3]=0となるが,i=0のときのC[2]≠0となることがある*/ } if((CC[2] &0xF000)!=0){cC=(CC[2]<< 8)|(CC[1]>> 8); expC=expA+expB ;} else {cC=(CC[2]<<12)|(CC[1]>> 4); expC=expA+expB-1;} return comField(expC,cC); } /*-------------------------------------------------------------------------------*/ /* ■絶対値除算 */ /* 正規化してから実行する。 */ /* 16進数6桁であるが,先頭が0となる可能性があるので7回計算を繰り返す。 */ /* (正規化して先頭桁は必ず非ゼロなので,先頭2桁共に0となることはありえない) */ /* @ 剰余=Aとする。 */ /* A i=0〜6として, CC[i]=剰余/B, 新しい剰余=(剰余 % B)*16 を繰り返す。 */ /* B CC[0]が非ゼロのとき割り算結果は1以上となるので */ /* 指数を1増やしてCC[0]〜CC[5]を結果とする(16進数として詰めなおす)。 */ /* C CC[0]がゼロのとき割り算結果は1未満となるので */ /* CC[1]〜CC[6]を結果とする(16進数として詰めなおす)。         */ /*                                       */ myFloat absDiv(myFloat A, myFloat B){/*絶対値除算*/ int flagA,flagB, expA, expB, expC, ist,ied,CC[7]; myFloat cA, cB, cC; sepField(A,&flagA,& expA,&cA); if(cA==0) return 0;/* 0/B=0 */ sepField(B,&flagB,& expB,&cB); if(cB==0) return 0x3FFFFFFF;/*Zero Divided */ normalize(&cA, &expA); normalize(&cB,&expB);/*正規化*/ cC=cA; for(int i=0;i<7; CC[i]=cC/cB, cC=(cC % cB)<<4, i++);/*除算実行*/ if(CC[0] !=0 ){expC=expA-expB+1;ist=0;ied=6;}/*除算結果が1以上のとき*/ else {expC=expA-expB ;ist=1;ied=7;}/*除算結果が1未満のとき*/ cC=0; for(int i=ist,SH=20;i=1; A/=16, exp++); return (flag<<31) | comField(exp,(long int)(A * 0x1000000+0.5));/* A*16^6*/ } /*----------------------------------------------------------------------------*/ double myToCFloat(myFloat A){/*シミュレータ用形式からCの浮動小数点に変換*/ int flag,exp;myFloat C; sepField(A,&flag,& exp,&C); double D=(double)C*pow(16.0,(double)exp-6.0); return (flag? -D : D); } /*----------------------------------------------------------------------------*/ /* ■符号を考慮した演算 */ /*符号同一のとき加算(符号同一), 異なるとき減算(符号は引かれる側の符号) */ myFloat myAdd(myFloat A, myFloat B){/*浮動小数点加算*/ if((A & 0x80000000)== (B & 0x80000000))return absAdd(A,B)|(A & 0x80000000); int cmp;if((cmp=absCmp(A,B))==0) return 0; if(cmp>0) return absSub(A,B)| (A & 0x80000000); else return absSub(B,A)| (B & 0x80000000); } /*----------------------------------------------------------------------------*/ /*符号異なるとき加算(Aの符号と同じ),同一のとき減算(符号はAと同一またはBの逆)*/ myFloat mySub(myFloat A, myFloat B){/*浮動小数点減算*/ if((A & 0x80000000)!=(B & 0x80000000))return absAdd(A,B)|(A & 0x80000000); int cmp;if((cmp=absCmp(A,B))==0) return 0; if(cmp>0) return absSub(A,B)| (A & 0x80000000); else return absSub(B,A)| ((B ^ 0x80000000) & 0x80000000); } /*----------------------------------------------------------------------------*/ /*符号同一のとき正, 異なるとき負となる*/ myFloat myMult(myFloat A, myFloat B){/*浮動小数点乗算*/ if((A & 0x80000000)!=(B & 0x80000000)) return 0x80000000 | absMult(A,B); else return 0x7FFFFFFF & absMult(A,B); } /*----------------------------------------------------------------------------*/ /*符号同一のとき正, 異なるとき負となる*/ myFloat myDiv(myFloat A, myFloat B){/*浮動小数点除算*/ if((A & 0x80000000)!=(B & 0x80000000)) return 0x80000000 | absDiv(A,B); else return 0x7FFFFFFF & absDiv(A,B); } /*----------------------------------------------------------------------------*/ /* ■longまたはintとの変換 */ long signed int myLint(myFloat A){/*myFloat to long signed integer*/ int flag, exp; myFloat C;long unsigned D; sepField(A,&flag,&exp,&C); if(C==0) return 0; normalize(&C, &exp);/*正規化(仮数部先頭4ビット≠0にする) */ if (exp <= 0) return 0; else if(exp == 6) D=C; else if(exp > 6) D = ((long unsigned) C<< ((exp-6)<<2)) & 0x7FFFFFFF; else D = ((long unsigned) C>> (24-(exp<<2))) & 0x7FFFFFFF; return (flag ? -(long signed)D:(long signed)D); } /*----------------------------------------------------------------------------*/ short signed int myInt(myFloat A){/*myFloat to short signed integer*/ int flag, exp; myFloat C;long unsigned D;sepField(A,&flag,&exp,&C); if(C==0) return 0; normalize(&C, &exp);/*正規化(仮数部先頭4ビット≠0にする) */ if (exp <= 0) return 0; else if(exp == 6) D=C; else if(exp > 6) D = ((long unsigned) C<< ((exp-6)<<2)) & 0x7FFF; else D = ((long unsigned) C>> (24-(exp<<2))) & 0x7FFF; return (flag ? -(short signed)D:(short signed)D); } /*----------------------------------------------------------------------------*/ myFloat intToMyFloat(int A){/*intをmyFloatに変換*/ if(A==0) return 0; int exp=6,flag=(A<0);myFloat C=abs(A); normalize(&C, &exp); return (flag<<31) | comField(exp,C); } /*----------------------------------------------------------------------------*/ /*seperate myFloat to long and fraction*/ void sepLongFrac(myFloat A, int *flag, long unsigned int *D, myFloat *FR){ int exp, SH; myFloat C; sepField(A,flag,&exp,&C); if(C==0) { *flag = 0; *D = 0; *FR = 0; return;} normalize(&C, &exp);/*正規化(仮数部先頭4ビット≠0にする) */ if (exp <= 0) { *D = 0; *FR=A & 0x7FFFFFFF;} else if(exp == 6) { *D = C; *FR=0;} else if(exp > 6) { *D = ((long unsigned ) C<< ((exp-6)<<2)) ;*FR=0;} else { SH = 24-(exp<<2); *D = (long unsigned) (C>> SH); C -= ((*D)<=10; LLD/=10,R*=10); if(flag) str[i++]='-'; while(R>=1){ D = LD / R; str[i++]=D+'0'; LD %= R; R /= 10;} return i; } /*----------------------------------------------------------------------------*/ int myFloatToString(myFloat A,char str[],int M){ long unsigned LDT;int flag; myFloat mF, F10=0x01A00000; sepLongFrac(A, &flag, &LDT, &mF);int i=longToString(flag, LDT, str); str[i++]='.'; for(int j=0;j=0; LLD=absDiv(LLD,F10),R=absMult(R,F10)); while(absCmp(R,F1)>0){ D = myInt(absDiv(AA,R)); str[i++]=D+'0'; AA= absSub(AA, absMult(R, intToMyFloat(D))); R=absDiv(R,F10); } D = myInt(AA); str[i++]=D+'0'; AA= absSub(AA, absMult(R, intToMyFloat(D))); str[i++]='.'; for(int j=0;j='0' && str[ip]<='9'; mF = mF*10+(str[ip++]-'0')); if(mF!=0){normalize(&mF, &exp);mF=comField(exp, mF);} if(str[ip]=='.'){ ip++; while(str[ip]>='0' && str[ip]<='9'){ F1=absDiv(F1,F10); mF=absAdd(mF,absMult(F1,intToMyFloat(str[ip++]-'0'))); } } return ((flag)? mF | 0x80000000 : mF); } int intToStringN(int flag, int LD, char str[],int M){ int i=0, D, R=1; for( int k=1;k=1){ D = LD / R; str[i++]=D+'0'; LD %= R; R /= 10;} return i; } int myFloatToExpString(myFloat DT, char ch[], int Np){ int Num=0, E, NNP=Np; myFloat A; A=myAbs(DT); myFloat C10=0x01A00000,C1=0x01100000; if((DT & 0x80000000)!=0)ch[Num++]='-'; E=0; while(absCmp(A,C10)>0) { E++; A=myDiv (A,C10);} while(absCmp(A,C1 )<0) {E--; A=absMult(A,C10);} int N =myFloatToString(A,&ch[Num],Np); Num+=(N-1); ch[Num++]='E';Num+=intToStringN(E<0, abs(E), &ch[Num],3); ch[Num++]=0; return Num; } /*----------------------------------------------------------------------------*/ /* ■テストメイン */ int _tmain(int argc, _TCHAR* argv[]){/*テストメイン*/ long unsigned LDT;int flag; myFloat mF; char str[2048]; int tmp; tmp= myFloatToExpString(cToMyFloat(123.567), str,7); printf("\n myFloatToExpString %s ", str); tmp= myFloatToExpString(0x7FFFFFFF, str,7); printf("\n myFloatToExpString %s ", str); tmp= myFloatToStringB(0x3FFFFFFF, str,7); printf("\n myFloatToStringB %s\n tmp=%d", str,tmp); tmp= myFloatToExpString (0x3FFFFFFF, str,7); printf("\n myFloatToExpString %s ", str); tmp= myFloatToString (0x3FFFFFFF, str,0); printf("\n myFloatToString %s ", str); printf("\n myMult %lf = 15.15* 15.15",myToCFloat(myMult(cToMyFloat(15.15),cToMyFloat(15.15)))); printf("\n myMult %lf = -12.5* -0.25",myToCFloat(myMult(cToMyFloat(-12.5),cToMyFloat(-0.25)))); printf("\n myMult %lf = 0.12345* 10",myToCFloat(myMult(cToMyFloat(0.12345),cToMyFloat(10)))); printf("\n myDiv %lf = 1.5/ 2",myToCFloat(myDiv(cToMyFloat(1.5),cToMyFloat(2)))); printf("\n myDiv %lf = 10.5/ -0.5",myToCFloat(myDiv(cToMyFloat(10.5),cToMyFloat(-0.5)))); printf("\n myAdd %lf = 15.5-10.24",myToCFloat(myAdd(cToMyFloat(15.5),cToMyFloat(-10.24)))); printf("\n mySub %lf = -10.5-(-15.24)",myToCFloat(mySub(cToMyFloat(-10.5),cToMyFloat(-15.24)))); printf("\n mySub %lf = -10.5-0",myToCFloat(mySub(cToMyFloat(-10.5),cToMyFloat(0)))); printf("\n mySub %lf = 123.5-123",myToCFloat(mySub(cToMyFloat(123.5),cToMyFloat(123)))); printf("\n Dec 10 = %08X",myInt(cToMyFloat(10))); sepLongFrac(cToMyFloat(-10.456), &flag, &LDT, &mF); printf("\n -10.456 flag = %d, LDT=%08X, mF=%lf",flag, LDT, myToCFloat(mF)); printf("\n intToMyFloat %lf %08X",myToCFloat(intToMyFloat(789)),intToMyFloat(789)); int i = myFloatToString(cToMyFloat(-12.12345),str,6); printf("\n -12.12345 String [%s]",str); myFloatToString(stringToMyFloat("-13.890"),str,6); printf("\n -13.890 string to myFloat %s",str); getchar(); return 0; }