二次方程式を解きたいのですが。

 投稿者:nuメール  投稿日:2012年 1月25日(水)22時40分32秒
  !以下のプログラムで二次方程式を解きたいのですが、どうにも手詰まりです。
!特にa=12,b=23,c=34,d=45,e=56,f=67の場合、如何ともし難いので、
!お分かりになられる方、ご教授を下さるようお願いいたします。
!----------------------------------------------------------------------
!二次方程式を解くためのプログラム
!根はわざと開かないような仕様にいたしました。
!----------------------------------------------------------------------
00100 LET t0=TIME
00500 DECLARE EXTERNAL FUNCTION GCD
01000 !【1】配列
01100 DIM x(2)
01200 DIM y(3)
01300 GOTO 02000 !【2】代入へ

02000 !【2】代入
02100 PRINT "「二次方程式の解を計算します。"
02110 PRINT "  a/d*x^2 + b/e*x + c/f = 0 の形式で、"
02120 PRINT "  a,b,c,d,e,fの順序で入力してください。」"
02130 PRINT "  <例> a*x^2 + b*x + c = 0 のとき、"
02140 PRINT "       a,b,c,1,1,1と入力します。"
02300 INPUT a,b,c,d,e,f
02400 LET y(1)=a*e*f !置き換え
02500 LET y(2)=b*d*f
02600 LET y(3)=c*d*e
02700 LET G=GCD(GCD(y(1),y(2)),y(3))
02750 LET a=y(1)/G
02800 LET b=y(2)/G
02850 LET c=y(3)/G
02900 GOTO 03000 !【3】判別式分岐へ

03000 !【3】判別式分岐
03100 LET D=b^2-4*a*c !D=平方数
03200 IF D>=0 THEN GOTO 04000 !【4】実数解のときの平方数分岐
03300 IF D<0 THEN GOTO 03400
03400 LET D=-D
03500 GOTO 08000!【8】虚数解のときの分岐

04000 !【4】実数解のときの平方数分岐
04100 IF INT(SQR(D))=SQR(D) THEN GOTO 05000!【5】!実数解のときのD=平方数のときの解
04200 IF INT(SQR(D))><SQR(D) THEN GOTO 06000!【6】!実数解のときのD≠平方数のときの分岐

05000 !【5】実数解のときのD=平方数のときの解
05100 LET x(1)=(-b-SQR(D))/(2*a)
05200 LET x(2)=(-b+SQR(D))/(2*a)
05300 PRINT "x(1)=";x(1)
05400 PRINT "x(2)=";x(1)
05500 GOTO 20000 !ENDへ

06000 !【6】実数解のときのD≠平方数のときの分岐
06100 FOR i=2 TO INT(SQR(D))
06300    IF MOD(D/(i^2)-INT(D/(i^2)),1)><0 THEN GOTO 06900 !iの探索を行うループ
06400    IF MOD(D/(i^2)-INT(D/(i^2)),1)=0 THEN GOTO 06500 !除数分岐へ
06500    IF MOD((i^2)-INT((i^2)),1)><0 THEN GOTO 06900
06600    IF MOD((i^2)-INT((i^2)),1)=0 THEN GOTO 06620 !実数解のときの-D≠平方数のときの解
06620    IF i/GCD(2*a,i)=1 THEN GOTO 06820
06640    IF i/GCD(2*a,i)><1 THEN GOTO 06700
06700    PRINT "x(1)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"+";i/GCD(2*a,i);"√";D/(i^2);"/";2*a/GCD(2*a,i)
06800    PRINT "x(2)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"-";i/GCD(2*a,i);"√";D/(i^2);"/";2*a/GCD(2*a,i)
06805    PRINT "⇔"
06810    GOTO 06900
06820    PRINT "x(1)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"+";"√";D;"/";2*a
06840    PRINT "x(2)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"-";"√";D;"/";2*a
06850    PRINT "⇔"
06900 NEXT i
07000 PRINT"これが解です。"
07100 GOTO 20000 !ENDへ

08000 !【8】虚数解のときの平方数分岐
08100 IF INT(SQR(D))=SQR(D) THEN GOTO 09000!【9】!虚数解のときの-D=平方数のときの解
08200 IF INT(SQR(D))><SQR(D) THEN GOTO 10000!【10】!虚数解のときの-D≠平方数のときの分岐

09000 !【9】虚数解のときの-D=平方数のときの解
09100 LET x(1)=(-b-SQR(D))/(2*a)
09200 LET x(2)=(-b+SQR(D))/(2*a)
09300 PRINT "x(1)=";x(1)
09400 PRINT "x(2)=";x(2)
09500 GOTO 20000 !ENDへ

10000 !【10】虚数解のときの-D≠平方数のときの分岐
10100 FOR i=2 TO INT(SQR(D))
10300    IF MOD(D/(i^2)-INT(D/(i^2)),1)><0 THEN GOTO 10900 !iの探索を行うループ
10400    IF MOD(D/(i^2)-INT(D/(i^2)),1)=0 THEN GOTO 10500 !除数分岐へ
10500    IF MOD((i^2)-INT((i^2)),1)><0 THEN GOTO 10900
10600    IF MOD((i^2)-INT((i^2)),1)=0 THEN GOTO 10620 !実数解のときの-D≠平方数のときの解
10620    IF i/GCD(2*a,i)=1 THEN GOTO 10820
10640    IF i/GCD(2*a,i)><1 THEN GOTO 10700
10700    PRINT "x(1)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"-";i/GCD(2*a,i);"√";D/(i^2);"i/";2*a/GCD(2*a,i)
10800    PRINT "x(2)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"+";i/GCD(2*a,i);"√";D/(i^2);"i/";2*a/GCD(2*a,i)
10805    PRINT "⇔"
10810    GOTO 10900
10820    PRINT "x(1)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"+";"√";D;"i/";2*a
10840    PRINT "x(2)=";-b/GCD(2*a,b);"/";2*a/GCD(2*a,b);"-";"√";D;"i/";2*a
10850    PRINT "⇔"
10900 NEXT i
11100 PRINT"これが解です。"
11150 PRINT TIME-t0;"秒かかりました。"
11200 GOTO 20000 !ENDへ
20000 END

30000 EXTERNAL FUNCTION GCD(a,b)
30100 DO
30200    LET r=MOD(a,b)
30300    IF r=0 THEN EXIT DO
30400    LET a=b
30500    LET b=r
30600 LOOP
30700 LET GCD=b
30800 END FUNCTION
 

二次方程式を解きたいのですが。(訂正あり)

 投稿者:nuメール  投稿日:2012年 1月25日(水)23時25分49秒
  > No.1731[元記事へ]

訂正です。
  (誤)05400 PRINT "x(2)=";x(1)
→(正)05400 PRINT "x(2)=";x(2)
 

Re: 二次方程式を解きたいのですが。

 投稿者:SECOND  投稿日:2012年 1月26日(木)09時46分34秒
  > No.1731[元記事へ]

06100 FOR i=2 TO INT(SQR(D))
10100 FOR i=2 TO INT(SQR(D))  →  FOR i=1 TO INT(SQR(D)) 両方とも、変えてみる。

√ の中の約数 i を探す所を、
2から始めると、見つからない場合の処理が、必要になりますが、文中にそれは有りません。
1から始めれば、最初に1つは、見つかるので、見つからない場合のブランクを、防げます。

<追記>
FOR i=INT(SQR(D)) TO 1 STEP -1   !1へ向って、逆回しにする。
   IF FP( D/(i^2))=0 THEN        ! 小数部が0なら、
      (
       )         !最初に見つかる最大の約数iで PRINT し、
      EXIT FOR   !for~next から出る。最悪でも i=1 で 漏れない。
   END IF
NEXT i
 

Re: 二次方程式を解きたいのですが。

 投稿者:山中和義  投稿日:2012年 1月26日(木)13時28分21秒
  > No.1731[元記事へ]

nuさんへのお返事です。

有理数(分数)の計算ができるので、こちらでまず考えてみました。


!2次方程式Ax^2+Bx+C=0すると、
! 判別式D=B^2-4AC
! 解の公式より、x=(-B/(2A))±√{(B^2-4AC)}/(2A)
!となる。

OPTION ARITHMETIC RATIONAL !有理数モード

LET A=12/45
LET B=23/56
LET C=34/67


LET D=B^2-4*A*C
IF D=0 THEN !重根
   PRINT "x="; -B/(2*A); "(重根)"

ELSE
   LET S=NUMER(D) !分子 ※符号なし
   LET T=DENOM(D) !分母 ※符号なし
   !!!PRINT D;S;T !debug
   CALL SqNormalize(S,P,Q) !√S=P√Qと変形する
   CALL SqNormalize(T,X,Y)
   !!!PRINT P;Q; X;Y !debug
   IF Q=1 AND Y=1 THEN !平方根の中が平方数なら、m形
      PRINT "x="; (-B+P/X)/(2*A); ","; (-B-P/X)/(2*A)

   ELSE !m±√n形
      LET W=-B/(2*A) !有理数部(実部)
      IF W<>0 THEN PRINT "x="; W;

      PRINT "±";

      IF D<0 THEN PRINT "i * "; !虚数単位

      LET W=(P/X)/(2*A*Y) !平方根部(虚部) ※有理化
      LET Z=NUMER(W)
      IF Z>1 THEN PRINT Z;
      PRINT "√"; Q*Y;
      LET Z=DENOM(W)
      IF Z>1 THEN PRINT "/"; Z

   END IF

END IF

END

EXTERNAL SUB SqNormalize(n, p,q) !平方根の中をできるだけ小さな正の整数に直す
OPTION ARITHMETIC RATIONAL !有理数モード
!※n=p^2*q、n,p,q≧0とすると、SQR(n)=p*SQR(q)と変形できる。
LET q=1 !※SQR(0)=0*SQR(1)とする ※n=0なら、1行下のFOR文でp=0は設定される
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



実行結果
x=-345/448 ±i * √ 1180210695 / 30016
 

Re: 二次方程式を解きたいのですが。

 投稿者:nuメール  投稿日:2012年 1月26日(木)16時50分17秒
  > No.1733[元記事へ]

SECONDさんへのお返事です。

SECOND様、ご回答いただきありがとうございました。
大変勉強になりました。
早くなることを期待して、
FOR i=2 TO INT(SQR(D))にしたことが主な原因だったようです。
FOR i=1 TO INT(SQR(D))に元どおり手直しします。
あとで配列を変えて、最後まで根を簡単にするようにしたいと思います。
参考にさせていただきます。

nu
 

Re: 二次方程式を解きたいのですが。

 投稿者:nuメール  投稿日:2012年 1月26日(木)16時52分27秒
  > No.1734[元記事へ]

山中和義さんへのお返事です。

山中様、ご回答いただきありがとうございました。
大変勉強になりました。

有理数モードの有効性に驚かされました。
おそらく、DATAからREADの流れで、一気に複数の問題を解くのに向いていると思います。
参考にさせていただきます。

un
 

戻る