対数の計算

 投稿者:山中和義  投稿日:2009年 6月29日(月)19時47分16秒
  ・同じ底での加減算と整数倍

未サポート
・異なる底での加減算と整数倍、乗除算
!対数の計算

!自然数nを、素因数分解 n=2^a*3^b*5^c* … とすると、
!LOG(n)=LOG(2^a*3^b*5^c* … )=a*LOG(2)+b*LOG(3)+c*LOG(5)+ … となる。
!真数が同じ素数のものを同類項としてまとめることができる。

OPTION ARITHMETIC RATIONAL

LET szLg=100 !扱う対数の範囲 ※必要に応じて変更のこと
LET cBASE=2 !仮の底 ※正の実数 a≠1 をとる

!演算関連
SUB LogSet(n, a()) !a=LOG(n)とする
   IF n<=0 OR INT(n)<>n THEN
      PRINT "非負の整数ではありません。"; n
      STOP
   END IF

   MAT a=ZER

   LET x=n !素因数分解 n=2^a*3^b*5^c* …
   LET m=2
   DO WHILE x>1 !1まで繰り返す
      LET K=0
      DO WHILE MOD(x,m)=0 !mで割り切れるなら(mは因数)
         LET x=x/m !割り切れるまで(m^Kの形)
         LET K=K+1
      LOOP

      LET a(m)=K !m^K ※ … + K*LOG(m) + …

      IF m>2 THEN !2,3,5,7,9,・・・で調べる
         LET m=m+2
      ELSE
         LET m=m+1
      END if
   LOOP
END SUB
DIM w0(szLg),w1(szLg) !作業用
SUB LogSetQ(x,y, a()) !a=LOG(x/y)、x:整数、y:正の整数 とする
   IF y<=0 OR INT(y)<>y THEN
      PRINT "正の整数を設定してください。"; y
      STOP
   END IF
   CALL LogSet(x, w0) !LOG(x/y)=LOG(x)-LOG(y)
   CALL LogSet(y, w1)
   CALL LogSub(w0,w1, a)
END SUB
SUB LogSetN(n, a()) !a=nとする
   MAT a=ZER
   LET a(cBASE)=n
END SUB

SUB LogAdd(a(),b(), c()) !加算 c=a+b
   MAT c=a+b
END SUB
SUB LogSub(a(),b(), c()) !減算 c=a-b
   MAT c=a-b
END SUB
SUB LogMulN(a(),N, c()) !乗算 c=a*N ※{Log(x)+Log(y)+ … }*N
   MAT c=(N)*a
END SUB
SUB LogDivN(a(),N, c()) !除算 c=a/N ※{Log(x)+Log(y)+ … }*N
   IF N=0 THEN
      PRINT "0では割れません。"
      STOP
   END IF
   MAT c=(1/N)*a
END SUB

!出力関連
SUB LogPrint(a()) !LOG(x)+LOG(y)+ … 形式で表示する
   LET flg=0
   FOR i=1 TO UBOUND(a) !小さい順に
      LET Ai=a(i) !係数
      IF Ai<>0 THEN
         IF flg=1 THEN PRINT " + "; !継続なら

         IF Ai<0 THEN
            PRINT "( ";Ai;") ";
         ELSE
            IF i=cBASE OR Ai<>1 THEN PRINT Ai; !係数が1以外なら
         END IF

         IF i<>cBASE THEN !真数が底以外なら
            IF Ai<>1 THEN PRINT "* "; !係数が1以外なら
            PRINT "Log";STR$(cBASE);"(";i;")";
         END IF
         LET flg=1
      END IF
   NEXT i
   IF flg=0 THEN PRINT " 0";
   PRINT
END SUB
SUB LogPrintQ(a()) !m/n+LOG(x/y) 形式で表示する
   CALL LogPack(a, K,x1,x2)

   IF K<>0 THEN PRINT K; !有理数の部分

   !無理数の部分
   IF x1/x2=1 THEN !真数=1の場合
      IF K=0 THEN PRINT "+ 0";
   ELSEIF x1=1 THEN !真数=1/x2の場合
      PRINT "- Log";STR$(cBASE);"(";x2;")";
   ELSE
      PRINT "+ Log";STR$(cBASE);"(";x1/x2;")";
   END IF
   PRINT
END SUB

!補助ルーチン
DIM b(szLg),c(szLg) !作業用
SUB LogYYY(a(),b(), K) !c=a-b^Kを調べる
   LET K=0 !引いた回数

   FOR i=1 TO UBOUND(b)
      IF b(i)<>0 THEN EXIT FOR !最小の素因数を探す
   NEXT i
   IF i>UBOUND(a) OR a(i)=0 THEN EXIT SUB !共通な素因数がない場合

   LET cSGN=1 !符号
   IF a(i)*b(i)<0 THEN LET cSGN=-1

   MAT c=a
   DO
      IF cSGN>0 THEN MAT c=c-b ELSE MAT c=c+b

      FOR i=1 TO UBOUND(c)
         IF c(i)*a(i)<0 THEN EXIT DO !引きすぎかどうか確認する
      NEXT i

      LET K=K+1
   LOOP
   LET K=cSGN*K

   IF cSGN>0 THEN MAT c=c+b ELSE MAT c=c-b !引き戻し
END SUB
SUB LogPack(a(), K,nn,mm) !1つの項にまとめる a*Log(x)+b*Log(y)+ … =Log(x^a*y^b* … )
   IF INT(cBASE)=cBASE THEN !整数なら
      CALL LogSet(cBASE, b) !底を素因数分解する
   ELSE
      CALL LogSetQ(NUMER(cBASE),DENOM(cBASE), b)
   END IF
   CALL LogYYY(a,b, K) !真数=底^K ?

   IF K=0 THEN
      CALL LogYYY(b,a, K) !底と真数をひっくり返す。底=真数^K ?
      IF k<>0 THEN LET K=1/K
   END IF

   IF K<>0 THEN MAT a=c


   LET nn=1 !分子
   LET mm=1 !分母
   FOR i=1 TO UBOUND(a) !真数を1つにまとめる
      IF a(i)<0 THEN !係数が負なら、分母へ
         LET mm=mm*i^ABS(a(i))
      ELSE
         LET nn=nn*i^a(i)
      END IF
   NEXT i


   IF nn<>1 AND mm=1 THEN !真数が1より大きな整数なら
      FOR i=1 TO UBOUND(b) !底
         IF b(i)<>0 THEN
            IF a(i)=0 THEN !共通素因子がなければ
               LET K1=0
               EXIT FOR
            END IF
            LET K1=a(i)/b(i) !底=真数^K1 ?
         END IF
      NEXT i
      IF K1<>0 THEN
         LET K=K+K1
         LET nn=1
      END IF
   END IF

   !!!PRINT "K=";K; "nn=";nn; "mm=";mm !debug
END SUB
!------------------------------ ここまでがサブルーチン

DIM T1(szLg),T2(szLg),T3(szLg),T4(szLg) !作業用


!●例1 LOG(2/3)+LOG(12/25)-LOG(8/15) = LOG(3)-LOG(5) = LOG(3/5) の計算

LET cBASE=2 !底

CALL LogSetQ(2,3, T1) !T1=LOG(2/3)
CALL LogSetQ(12,25, T2) !T1=LOG(12/25)
CALL LogAdd(T1,T2, T1)
CALL LogSetQ(8,15, T2) !T1=LOG(8/15)
CALL LogSub(T1,T2, T1)

!!!MAT PRINT T1; !debug
CALL LogPrint(T1)
CALL LogPrintQ(T1)



!●例2 LOG10(7/4)-LOG10(9)-2*LOG10(5/3)-LOG10(49)/2 = -2 の計算

LET cBASE=10 !底

CALL LogSetQ(7,4, T1) !T1=LOG(7/4)
CALL LogSet(9, T2) !T2=LOG(9)
CALL LogSub(T1,T2, T1)

CALL LogSetQ(5,3, T2) !T2=LOG(5/3)
CALL LogMulN(T2,2, T2)
CALL LogSub(T1,T2, T1)

CALL LogSet(49, T2) !T2=LOG(49)
CALL LogDivN(T2,2, T2)
CALL LogSub(T1,T2, T1)

!MAT PRINT T1; !debug
CALL LogPrint(T1)
CALL LogPrintQ(T1)



!●例3 LOG27(1/9) = -2/3 の計算

!LET cBASE=9
LET cBASE=27

!CALL LogSetQ(1,3, T1)
CALL LogSetQ(1,9, T1)
!MAT PRINT T1; !debug
CALL LogPrint(T1)
CALL LogPrintQ(T1)


END
 

戻る