変数xの有理式
B(n-1)*x^(n-1)+B(n-2)*x^(n-2)+ … +B(1)*x+B(0)
---------------------------------------------------------
A(n)*x^n+A(n-1)*x^(n-1)+A(n-2)*x^(n-2)+ … +A(1)*x+A(0)
とする。分母の多項式が、
(x-α1)*(x-α2)* … *(x-αn) ※重根は含まない
と因数分解可能のとき
k1 k2 kn
------- + ------- + … + -------
(x-α1) (x-α2) (x-αn)
と部分分数に分解できる。
LET N=3 !次数
DIM A(0 TO N),B(0 TO N-1) !係数
DATA 6,-1,-3 !分子 6*x^2-x-3
DATA 1,0,-1,0 !分母 x^3-x=x*(x+1)*(x-1)
DATA 0,1 !根,重複度
DATA -1,1
DATA 1,1
!DATA 0,1,1 !分子 s+1
!DATA 1,-5,8,-4 !分母 s^3-5*s^2+8*s-4=(s-1)*(s-2)^2
!DATA 1,1 !根,重複度
!DATA 2,2
FOR i=N-1 TO 0 STEP -1
READ B(i)
NEXT i
FOR i=N TO 0 STEP -1
READ A(i)
NEXT i
!ヘヴィサイド（Heaviside）の展開定理
! F(x)=B(x)/A(x)とする。
! x-αが分母A(x)の重複度１の因数とすると
! 部分分数 C/(x-α) の係数は、C=(x-α)*F(x)│x=αとなる。
!
! x-αが分母A(x)の重複度mの因数とすると
! 部分分数 C1/(x-α)+C2/(x-α)^2+ … +Cm/(x-α)^m の係数は、
! Cm-k=1/k!*d^k/dx^k{(x-α)^m*F(x)}│x=α、k=0,1,2,…,m-1となる。
DIM Q(0 TO N),T1(0 TO 4*N),T2(0 TO 4*N) !作業用
DIM AA(0 TO 2*N),BB(0 TO 2*N),dA(0 TO 2*N),dB(0 TO 2*N)
LET Lp=N
DO WHILE Lp>0 !各根に対して
READ v,m !α、重複度
SELECT CASE m
CASE 1 !単根なら
CALL poly_divByLin(A,v, Q,R) !(x-α)*F(x)=B(x)/(A(x)/(x-α))
CALL PrintOut(B,Q,0,v)
CASE ELSE !重複根なら
MAT AA=A
MAT BB=B
FOR j=1 TO m !(x-α)^m*F(x)の０階微分
CALL poly_divByLin(AA,v, Q,R) !(x-α)^m*F(x)=B(x)/(A(x)/(x-α)^m)
MAT AA=Q
NEXT j
CALL PrintOut(BB,AA,m,v) !Cm
FOR k=1 TO m-1 !ｋ階微分
!微分の公式 (B/A)'=(B'*A-B*A')/A^2
CALL poly_diff(BB, dB) !B'
CALL poly_diff(AA, dA) !A'
CALL poly_mul(dB,AA, T1) !分子
CALL poly_mul(BB,dA, T2)
MAT T1=T1-T2
CALL poly_copy(T1, BB)
CALL poly_mul(AA,AA, T1) !分母
CALL poly_copy(T1, AA)
CALL PrintOut(BB,AA,m-k,v/fact(k)) !Cm-k
NEXT k
END SELECT
LET Lp=Lp-m !次へ
LOOP
SUB PrintOut(B(),A(),k,v) !分数を表示する
PRINT poly_val(B,v)/poly_val(A,v); !分子側=C
IF v>0 THEN !分母側=x-α
PRINT "/ ( x -";v;")";
ELSEIF v<0 THEN
PRINT "/ ( x +";ABS(v);")";
ELSE
PRINT "/ x ";
END IF
IF k>1 THEN PRINT " ^";k ELSE PRINT !べき乗
END SUB
END
!補助ルーチン
!変数xの多項式 Σ[k=0,n]a(k)*x^k=a(n)*x^n+a(n-1)*x^(n-1)+ … +a(1)*x+a(0)
!演算関連
EXTERNAL FUNCTION poly_degree(A()) !次数を得る
FOR i=UBOUND(A) TO 1 STEP -1
IF A(i)<>0 THEN EXIT FOR !係数が０でない最初の位置
NEXT i
LET poly_degree=i
END FUNCTION
EXTERNAL SUB poly_divByLin(A(),v, Q(),R) !多項式a(x)をx-αで割ったときの商q(x)と余りRを求める
MAT Q=ZER
LET aa=poly_degree(A) !次数
IF aa>0 THEN !１次式以上なら
LET qq=aa-1 !次数
LET Q(qq)=A(aa) !商 ※組立除法
FOR i=qq TO 1 STEP -1
LET Q(i-1)=A(i)+Q(i)*v
NEXT i
LET R=A(0)+Q(0)*v !余り
ELSE
LET Q(0)=0 !商
LET R=A(0) !余り
END IF
END SUB
EXTERNAL FUNCTION poly_val(A(),v) !多項式を計算する（数値代入）
LET k=0
FOR i=UBOUND(A) TO 0 STEP -1 !Hornerの方法
LET k=k*v+A(i) !(…(((An*X+An-1)*X+An-2)*X+An-3)*X+…+A1)*X+A0
NEXT i
LET poly_val=k
END FUNCTION
EXTERNAL SUB poly_copy(A(), S()) !代入 S=A
MAT S=ZER
FOR i=0 TO UBOUND(S)
LET S(i)=A(i)
NEXT i
END SUB
EXTERNAL SUB poly_mul(A(),B(), S()) !乗算 S=A*B ※S≠A、S≠B
MAT S=ZER
FOR i=UBOUND(A) TO 0 STEP -1
FOR j=UBOUND(B) TO 0 STEP -1
LET S(i+j)=S(i+j)+A(i)*B(j) !すべての係数をかける
NEXT j
NEXT i
END SUB
EXTERNAL SUB poly_diff(A(), S()) !微分 S=A'
FOR i=1 TO UBOUND(A)
LET S(i-1)=i*A(i)
NEXT i
END SUB