十進BASIC 第2掲示板 過去ログ

!平方根の計算

!p*SQR(q)を、「SQR(q)とその係数p」として、配列a(q)=pで表せる。
!∵n=p^2*q、n,p,q≧0とすると、SQR(n)=p*SQR(q)と変形できる。
!　これより、SQR(32)=4*SQR(2)となり、同類項をまとめる場合、共通項SQR(2)として扱える。
!　また、SQR()を配列として解釈すると、SQR(n)を一意に表現できる。

LET szRt=10 !扱う平方根の範囲　※必要に応じて変更のこと

!演算関連
SUB SqSet(n, a()) !a=SQR(n)とする
   IF INT(n)<>n THEN
      PRINT "整数を設定してください。"; n
      STOP
   END IF
   MAT a=ZER
   CALL SqNormalize(ABS(n), p,q)
   IF n<0 THEN LET q=-q
   LET a(q)=p
END SUB
SUB SqSetQ(x,y, a()) !a=SQR(x/y)、x:整数、y:正の整数 とする
   IF y<=0 OR INT(y)<>y THEN
      PRINT "正の整数を設定してください。"; y
      STOP
   END IF
   CALL SqSet(x*y, a)
   CALL SqDivN(a,y, a)
END SUB
SUB SqSetN(n, a()) !a=nとする
   MAT a=ZER
   LET a(1)=n
END SUB

SUB SqAdd(a(),b(), c()) !加算 c=a+b
   MAT c=a+b
END SUB
SUB SqSub(a(),b(), c()) !減算 c=a-b
   MAT c=a-b
END SUB
SUB SqMulN(a(),N, c()) !乗算 c=a*N　※{√(x)+√(y)+ … }*N
   MAT c=(N)*a
END SUB
SUB SqDivN(a(),N, c()) !除算 c=a/N　※{√(x)+√(y)+ … }/N
   IF N=0 THEN
      PRINT "０では割れません。"
      STOP
   END IF
   MAT c=(1/N)*a
END SUB
DIM w(-szRt TO szRt) !作業用
SUB SqMulS(a(),B, c()) !乗算 c=a*B　※{√(x)+√(y)+ … }*√(B)
   MAT w=ZER
   FOR i=LBOUND(a) TO UBOUND(a)
      IF a(i)<>0 THEN !係数が０以外なら
         IF i*B<0 THEN
            CALL SqNormalize(ABS(i*B), p,q) !※√(-a)*√(b)=√(-a*b)、√(a)*√(-b)=√(-a*b)
            LET q=-q
         ELSE
            CALL SqNormalize(i*B, p,q) !※√(a)*√(b)=√(a*b)
            IF B<0 THEN LET p=-p !※√(-a)*√(-b)=-√(a*b)
         END IF
         LET w(q)=w(q)+a(i)*p !√(a[])*√(B)
      END IF
   NEXT i
   MAT c=w
END SUB
SUB SqDivS(a(),B, c()) !除算 c=a/B　※{√(x)+√(y)+ … }/√(B)
   IF B=0 THEN
      PRINT "０では割れません。"
      STOP
   END IF
   CALL SqMulS(a,B, c) !√(i)*√(B)/B　※分母を有理化する
   CALL SqDivN(c,B, c)
END SUB

DIM w0(-szRt TO szRt),w1(-szRt TO szRt) !作業用
SUB SqMul(a(),b(), c()) !乗算 c=a*b　※{√(x)+√(y)+ … }*{√(u)+√(v)+ … }
   MAT w0=ZER
   FOR k=LBOUND(b) TO UBOUND(b)
      LET Bk=b(k) !√(a[])*√(b[k])
      IF Bk<>0 THEN !係数が０以外なら
         CALL SqMulS(a,k, w1) !平方根の中の部分
         CALL SqMulN(w1,Bk, w1) !係数の部分
         MAT w0=w0+w1
      END IF
   NEXT k
   MAT c=w0
END SUB
DIM w2(-szRt TO szRt),w3(-szRt TO szRt),w4(-szRt TO szRt),w5(-szRt TO szRt) !作業用
SUB SqDiv(a(),b(), c()) !除算 c=a/b　※{√(x)+√(y)+ … }/{√(u)+√(v)+ … }
   MAT w2=a
   MAT w3=b
   DO
      FOR k=LBOUND(b) TO UBOUND(b) !平方根を探す
         IF k<>1 AND w3(k)<>0 THEN EXIT FOR
      NEXT k
      IF k>UBOUND(b) THEN EXIT DO !平方根がなくなれば、有理化を終了する

      MAT w4=w3 !{√(u)+√(v)+ … }-√(k)をつくって、(s+t)*(s-t)=s^2-t^2の形へ
      LET w4(k)=-w3(k)

      CALL SqMul(w2,w4, w5) !分子側　{√(x)+√(y)+ … }*{√(u)+√(v)+ … -√(k)}
      MAT w2=w5
      CALL SqMul(w3,w4, w5) !分母側　{{√(u)+√(v)+ … }+√(k)}*{{√(u)+√(v)+ … }-√(k)}
      MAT w3=w5
   LOOP
   CALL SqDivN(w2,w3(1), c) !分子/分母
END SUB
SUB SqPow(a(),n, x()) !べき乗 x=a^n　※{√(x)+√(y)+ … }^n
   IF INT(n)<>n THEN
      PRINT "べき乗数が整数ではありません。"; n
      STOP
   END IF
   CALL SqSetN(1, w3) !x=1
   MAT w2=a !b=a
   LET n2=ABS(n)
   DO UNTIL n2=0
      IF MOD(n2,2)=1 THEN CALL SqMul(w3,w2, w3) !x=x*b
      CALL SqMul(w2,w2, w2) !b=b*b
      LET n2=INT(n2/2)
   LOOP
   IF n<0 THEN !負なら逆数にする
      CALL SqSetN(1, w2)
      CALL SqDiv(w2,w3, w3) !x=1/x
   END IF
   MAT x=w3
END SUB

SUB SqNormalize(n, p,q) !平方根の中をできるだけ小さな正の整数に直す
!※n=p^2*q、n,p,q≧0とすると、SQR(n)=p*SQR(q)と変形できる。
   LET q=1 !※SQR(0)=0*SQR(1)とする　※n=0なら、１行下のFOR文でp=0は設定される
   FOR p=INT(SQR(n)) TO 1 STEP -1 !約数p^2の候補を大きい方から
   !FOR p=INTSQR(n) TO 1 STEP -1 !約数p^2の候補を大きい方から　※有理数モードのとき
      LET q=n/p^2
      IF q=INT(q) THEN EXIT FOR !qは自然数より
   NEXT p
END SUB

!出力関連
SUB SqPrint(a()) !√(x)+√(y)+ … 形式で表示する
   LET flg=0
   FOR i=LBOUND(a) 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=1 OR Ai<>1 THEN PRINT Ai; !係数が１以外なら
         END IF
         IF i<>1 THEN !SQR(1)以外なら
            IF Ai<>1 THEN PRINT "* "; !係数が１以外なら
            PRINT "SQR(";i;")";
         END IF
         LET flg=1
      END IF
   NEXT i
   IF flg=0 THEN PRINT " 0";
   PRINT
END SUB
!------------------------------ ここまでがサブルーチン

DIM T1(-szRt TO szRt),T2(-szRt TO szRt),T3(-szRt TO szRt) !作業用



!●例１　SQR(32)-SQR(3)*{2*SQR(2)+SQR(6)}+12/SQR(6) の計算

DIM c2(-szRt TO szRt),c3(-szRt TO szRt),c6(-szRt TO szRt),c32(-szRt TO szRt) !定数
CALL SqSet(32, c32) !c32=SQR(32)
CALL SqSet(3, c3) !c3=SQR(3)
CALL SqSet(2, c2) !c2=SQR(2)
CALL SqSet(6, c6) !c6=SQR(6)

CALL SqMulN(c2,2, T1) !T1=2*SQR(2)
CALL SqAdd(T1,c6, T1) !T1=2*SQR(2)+SQR(6)
!CALL SqPrint(T1)
CALL SqMul(c3,T1, T3) !T3=SQR(3)*{2*SQR(2)+SQR(6)}
!!!または、CALL SqMulS(T1,3, T3) !T3=SQR(3)*{2*SQR(2)+SQR(6)}
!CALL SqPrint(T3)

CALL SqSub(c32,T3, T1) !T1=SQR(32)-SQR(3)*{2*SQR(2)+SQR(6)}

DIM n12(-szRt TO szRt) !定数
CALL SqSetN(12, n12) !n12=12

CALL SqDivS(n12,6, T2) !T2=12/SQR(6)
!CALL SqPrint(T2)
CALL SqAdd(T1,T2, T1) !T1=SQR(32)-SQR(3)*{2*SQR(2)+SQR(6)}+12/SQR(6)

CALL SqPrint(T1) !結果




!●例２ (-5+SQR(-9))/(3-SQR(-4)) の計算

DIM n3(-szRt TO szRt),nm5(-szRt TO szRt) !定数
CALL SqSetN(3, n3) !n3=3
CALL SqSetN(-5, nm5) !nm5=-5

DIM cm4(-szRt TO szRt),cm9(-szRt TO szRt) !定数
CALL SqSet(-4, cm4) !cm4=SQR(-4)
CALL SqSet(-9, cm9) !cm9=SQR(-9)

CALL SqAdd(nm5,cm9, T1) !T1=-5+SQR(-9)
CALL SqSub(n3,cm4, T2) !T2=3-SQR(-4)
CALL SqDiv(T1,T2, T3) !T3=(-5+SQR(-9))/(3-SQR(-4))

CALL SqPrint(T3) !結果


END