プチコン 非公式コミュニティ
プレイ日記
こういち
◆ou0jbJnEJ0Kb
2025/7/21 15:49
2025/7/21 15:49
このサイト
https://rohaki.main.jp/mokuji.html
を参考に、多項式の根を求めるプログラムを書いてみました。
実装難易度高いわりにそれなりの頻度で使うんですよね。
最近はサブルーチンを細かく分割するのにハマってます(IF/FOR/WHILEの入れ子が3段ぐらいを目安に分割してる)。
'最大1000次式
DIM B[1000] 'aをpx+qで割った商
DIM C[1000] 'bをpおよびqで偏微分した値
'二次方程式を解く
'a,b,c…係数(ax^2+bx+c)
'rp,jp…正の解(実部,虚部)
'rn,jn…負の解(実部,虚部)
DEF SOLVEQUAD A,B,C OUT RP,JP,RN,JN
VAR D#=B*B-4*A*C '判別式
IF D#>0 THEN '実数解を持つ
IF B>0 THEN 'Bが正なら、正の解で分子を有理化
RP=2*C/(-B-SQR(D#)) '分子を有理化した解の公式
JP=0
RN=(-B+SQR(D#))/2/A '解の公式
JN=0
ELSE 'Bが負なら、負の解で分子を有理化
RP=(-B+SQR(D#))/2/A '解の公式
JP=0
RN=2*C/(-B+SQR(D#)) '分子を有理化した解の公式
JN=0
ENDIF
ELSE '複素数解を持つ
RP=-B/2/A
JP=SQR(-D#)/2/A
RN=-B/2/A
JN=-SQR(D#)/2/A
ENDIF
END
'aをpx+qで割った商を求める
DEF BAIRSTOW_DIV A[],B[],P,Q,N
VAR I%
B[0]=A[0]:B[1]=A[1]-P*B[0]
FOR I%=2 TO N-1
B[I%]=A[I%]-P*B[I%-1]-Q*B[I%-2]
NEXT I%
END
//bをqで偏微分した値,bをqで偏微分した値を求める
DEF BAIRSTOW_DIF B[],C[],P,Q,N
VAR I%
C[0]=B[0]:C[1]=B[1]-P*C[0]
FOR I%=2 TO N-1
C[I%]=B[I%]-P*C[I%-1]-Q*C[I%-2]
NEXT I%
END
(続く)
https://rohaki.main.jp/mokuji.html
を参考に、多項式の根を求めるプログラムを書いてみました。
実装難易度高いわりにそれなりの頻度で使うんですよね。
最近はサブルーチンを細かく分割するのにハマってます(IF/FOR/WHILEの入れ子が3段ぐらいを目安に分割してる)。
'最大1000次式
DIM B[1000] 'aをpx+qで割った商
DIM C[1000] 'bをpおよびqで偏微分した値
'二次方程式を解く
'a,b,c…係数(ax^2+bx+c)
'rp,jp…正の解(実部,虚部)
'rn,jn…負の解(実部,虚部)
DEF SOLVEQUAD A,B,C OUT RP,JP,RN,JN
VAR D#=B*B-4*A*C '判別式
IF D#>0 THEN '実数解を持つ
IF B>0 THEN 'Bが正なら、正の解で分子を有理化
RP=2*C/(-B-SQR(D#)) '分子を有理化した解の公式
JP=0
RN=(-B+SQR(D#))/2/A '解の公式
JN=0
ELSE 'Bが負なら、負の解で分子を有理化
RP=(-B+SQR(D#))/2/A '解の公式
JP=0
RN=2*C/(-B+SQR(D#)) '分子を有理化した解の公式
JN=0
ENDIF
ELSE '複素数解を持つ
RP=-B/2/A
JP=SQR(-D#)/2/A
RN=-B/2/A
JN=-SQR(D#)/2/A
ENDIF
END
'aをpx+qで割った商を求める
DEF BAIRSTOW_DIV A[],B[],P,Q,N
VAR I%
B[0]=A[0]:B[1]=A[1]-P*B[0]
FOR I%=2 TO N-1
B[I%]=A[I%]-P*B[I%-1]-Q*B[I%-2]
NEXT I%
END
//bをqで偏微分した値,bをqで偏微分した値を求める
DEF BAIRSTOW_DIF B[],C[],P,Q,N
VAR I%
C[0]=B[0]:C[1]=B[1]-P*C[0]
FOR I%=2 TO N-1
C[I%]=B[I%]-P*C[I%-1]-Q*C[I%-2]
NEXT I%
END
(続く)
コメント
こういち
2025/7/21 15:50
◆ou0jbJnEJ0Kb
'aを割り切るような二次式(x^2+px+q)を求める
DEF BAIRSTOW_STEP A[],N OUT P,Q
VAR L%=1000000 '反復回数(実行遅いなら3〜4桁ぐらい減らして良い)
VAR I%
VAR DP#=1,DQ#=1
REPEAD
P=RND(500)/25-10:Q=RND(500)/25-10 '-10〜10の乱数
FOR I%=0 TO L-1
BAIRSTOW_DIV A,B,P,Q,N
BAIRSTOW_DIF B,C,P,Q,N
'連立方程式を解く
VAR R#=C[N-3]*C[N-3]-(C[N-2]-b[N-2])*C[N-4]
DP#=(B[N-2]*C[N-3]-B[N-1]*C[N-4])/R#
DQ#=(-B[N-2]*(C[N-2]-B[N-2])+B[N-1]*C[N-3])/R#
P=P+DP#:Q=Q+DQ#
NEXT I%
UNTIL ABS(DP#)<1e-5&&ABS(DQ#)<1e-5 '誤差が大きければ初期値を変えてやり直し
END
DEF BAIRSTOW(A[],RO[],JO[],N)
VAR CNT%=0
VAR I%
'前処理
WHILE 0==A[0]
FOR I%=0 TO N-2
A[I%]=A[I%+1]
NEXT I%
A[N-1]=0
N=N-1
WEND
IF ABS(A[0]-1)>qe-5 THEN
VAR TMP#=A[0]
FOR I%=0 TO N-1
A[I%]=A[I%]/TMP#
NEXT I%
ENDIF
VAR P#,Q#
WHILE N>3
'aを割り切るような二次式(x^2+px+q)を求める
BAIRSTOW_STEP A,N OUT P#,Q#
'二次方程式を解く
VAR RP#,JP# '正の解
VAR RN#,JN# '負の解
SOLVEQUAD 1,P#,Q# OUT RP#,JP#,RN#,JN#
RO[CNT%]=RP#:JO[CNT%]=JP#
RO[CNT%+1]=RN#:JO[CNT%+1]=JN#
CNT%=CNT%+2
'新しいaを求める
BAIRSTOW_DIV A,A,P#,Q#,N
N=N-2
WEND
IF 3==N THEN '残りが2次式
'二次方程式を解く
SOLVEQUAD 1,P#,Q# OUT RP#,JP#,RN#,JN#
RO[CNT%]=RP#:JO[CNT%]=JP#
RO[CNT%+1]=RN#:JO[CNT%+1]=JN#
CNT%=CNT%+2
ELSE '残りが1次式
RO[CNT%]=-A[1]
JO[CNT%]=0
CNT%=CNT%+1
ENDIF
RETURN CNT%
END
DEF BAIRSTOW_STEP A[],N OUT P,Q
VAR L%=1000000 '反復回数(実行遅いなら3〜4桁ぐらい減らして良い)
VAR I%
VAR DP#=1,DQ#=1
REPEAD
P=RND(500)/25-10:Q=RND(500)/25-10 '-10〜10の乱数
FOR I%=0 TO L-1
BAIRSTOW_DIV A,B,P,Q,N
BAIRSTOW_DIF B,C,P,Q,N
'連立方程式を解く
VAR R#=C[N-3]*C[N-3]-(C[N-2]-b[N-2])*C[N-4]
DP#=(B[N-2]*C[N-3]-B[N-1]*C[N-4])/R#
DQ#=(-B[N-2]*(C[N-2]-B[N-2])+B[N-1]*C[N-3])/R#
P=P+DP#:Q=Q+DQ#
NEXT I%
UNTIL ABS(DP#)<1e-5&&ABS(DQ#)<1e-5 '誤差が大きければ初期値を変えてやり直し
END
DEF BAIRSTOW(A[],RO[],JO[],N)
VAR CNT%=0
VAR I%
'前処理
WHILE 0==A[0]
FOR I%=0 TO N-2
A[I%]=A[I%+1]
NEXT I%
A[N-1]=0
N=N-1
WEND
IF ABS(A[0]-1)>qe-5 THEN
VAR TMP#=A[0]
FOR I%=0 TO N-1
A[I%]=A[I%]/TMP#
NEXT I%
ENDIF
VAR P#,Q#
WHILE N>3
'aを割り切るような二次式(x^2+px+q)を求める
BAIRSTOW_STEP A,N OUT P#,Q#
'二次方程式を解く
VAR RP#,JP# '正の解
VAR RN#,JN# '負の解
SOLVEQUAD 1,P#,Q# OUT RP#,JP#,RN#,JN#
RO[CNT%]=RP#:JO[CNT%]=JP#
RO[CNT%+1]=RN#:JO[CNT%+1]=JN#
CNT%=CNT%+2
'新しいaを求める
BAIRSTOW_DIV A,A,P#,Q#,N
N=N-2
WEND
IF 3==N THEN '残りが2次式
'二次方程式を解く
SOLVEQUAD 1,P#,Q# OUT RP#,JP#,RN#,JN#
RO[CNT%]=RP#:JO[CNT%]=JP#
RO[CNT%+1]=RN#:JO[CNT%+1]=JN#
CNT%=CNT%+2
ELSE '残りが1次式
RO[CNT%]=-A[1]
JO[CNT%]=0
CNT%=CNT%+1
ENDIF
RETURN CNT%
END
コメントを書く
- こちらは「プチコン3号」「プチコンBIG」など、プチコンシリーズに関する話題を扱ったコミュニティです
- プチコンシリーズにまったく関係ない書き込みはご遠慮下さい。削除の対象となります
- こちらにはその他のゲームや雑談のコミュニティはなく、作る予定もありません (ひとりで管理できないため)。ごめんなさい
- ユーザー登録なしで書き込みができます
- 秘密の合い言葉は成りすましの防止 (トリップ機能)、書き込みの編集時の本人認証に使用します
- 秘密の合い言葉に他人に推測されやすい言葉、他サービスと同じパスワードは入力しないでください。
- 書き込むと、投稿時に入力したお名前と秘密の暗号が記憶され、ログイン状態になります