整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月 2日(土)10時21分14秒
  問題
nを任意の正の整数とする。
a≧b≧0、a+2b=n
を満たす整数の組(a,b)のそれぞれに対して、積abを考える。
その積のすべての和をf(n)とすると、f(n)をnの式で表せ。

答え
a≧bなので、n=a+2b≧(1+2)bより、1≦b≦[n/3]である。
Σab
=Σ[B=1,[N/3]](N-2B)B
=NΣ[B=1,[N/3]]B - 2Σ[B=1,[N/3]]B^2
=N{x(x+1)/2} - 2{x(x+1)(2x+1)/6}　ただし、x=[N/3]
=x(x+1)(3N-2-4x)/6
または
=( -4x^3+3(N-2)x^2+(3N-2)x )/6
（終り）


FOR N=1 TO 50
   PRINT "N="; N

   LET S=0 !a≧b≧1として考える
   FOR B=1 TO INT( N/(1+2) ) !n=a+2b≧(1+2)bより
      LET A=N-2*B
      IF A< B THEN STOP !論理エラー

      PRINT A;B !題意を満たす
      LET S=S+A*B

   NEXT B
   PRINT "積="; S


   LET X=INT(N/3) !f(n)
   !!PRINT X*(X+1)*(3*N-2-4*X)/6
   PRINT ( -4*X^3+3*(N-2)*X^2+(3*N-2)*X )/6
NEXT N

END



問題
nを任意の正の整数とする。
a≧b≧c≧0、a+2b+3c=n
を満たす整数の組(a,b,c)のそれぞれに対して、積abcを考える。
その積のすべての和をf(n)とすると、f(n)をnの式で表せ。

考察
類題として問題文が展開ができるので、いくつか先を見ていくと、

漸化式を用いて、

FOR n=1 TO 60
   PRINT n;":"; F2(n,INT(n/3)); F3(n,INT(n/6)); F4(n,INT(n/10)); F5(n,INT(n/15))
NEXT n
END

!a+2b=n、a≧b≧0、Σab=Σ[b=1,[n/3]]b*(n-2b)
EXTERNAL FUNCTION F2(n,k) !Σ[i=1,k]i(n-2i)
LET S=0
FOR i=1 TO k
   LET S=S+i*(n-2*i)
NEXT i
LET F2=S
END FUNCTION

!a+2b+3c=n、a≧b≧c≧0、Σabc=Σ[c=1,[n/6]] c*{ Σ[b=c,[(n-3c)/3]] b*((n-3c)-2b) }
EXTERNAL FUNCTION F3(n,k)
LET S=0
FOR i=1 TO k
   LET x=n-3*i
   LET S=S+i*(F2(x,INT(x/3)) - F2(x,i-1))
NEXT i
LET F3=S
END FUNCTION

EXTERNAL FUNCTION F4(n,k)
LET S=0
FOR i=1 TO k
   LET x=n-4*i
   LET S=S+i*(F3(x,INT(x/6)) - F3(x,i-1))
NEXT i
LET F4=S
END FUNCTION

EXTERNAL FUNCTION F5(n,k)
LET S=0
FOR i=1 TO k
   LET x=n-5*i
   LET S=S+i*(F4(x,INT(x/10)) - F4(x,i-1))
NEXT i
LET F5=S
END FUNCTION

（終り）


Σabc=Σ[c=1,[n/6]] c*{ Σ[b=c,[(n-3c)/3]] b*((n-3c)-2b) } なので、
Σ[b=c,[(n-3c)/3]] b*((n-3c)-2b) の部分を、
　a+2b=n の式 f(N)=( -4x^3+3(N-2)x^2+(3N-2)x )/6
のxに、[n/3]-c と c-1 を代入して求めると、


FOR N=1 TO 60
   PRINT "N="; N

   LET S=0 !a≧b≧c≧1として考える
   FOR C=1 TO INT( N/(1+2+3) ) !n=a+2b+3c≧(1+2+3)cより
      FOR B=C TO INT( (N-3*C)/(1+2) ) !n-3c=a+2b≧(1+2)bより
         LET A=N-3*C-2*B
         IF A< B THEN STOP !論理エラー

         PRINT A;B;C !題意を満たす
         LET S=S+A*B*C

      NEXT B
   NEXT C
   PRINT "積="; S


   LET y=INT(N/6) !f(n)
   LET x=INT(N/3)
   PRINT y*(y+1)* ( 16*y^3  +3*(5*x-2)*y^2  +(10*x^2 +5*(5-4*N)*x -14)*y  -5*x*(4*x^2 +(5-3*N)*x -N+1) +4 )/60


NEXT N

END


a+2b+3c+4d=n、a≧b≧c≧d≧0 となると、手計算では難しい。。。

 		 

Re: 整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月 4日(月)13時15分20秒
  > No.2993[元記事へ]

> 問題
> nを任意の正の整数とする。
> a≧b≧0、a+2b=n
> を満たす整数の組(a,b)のそれぞれに対して、積abを考える。
> その積のすべての和をf(n)とすると、f(n)をnの式で表せ。

答え
A=a-bとする。a≧bより、A≧0
また、A=(n-2b)-b=n-3b　∴A+3b=n
Σ[a+2b=n,a≧b≧0]ab
=Σ[A+3b=n,A≧0,b≧0](A+b)b
3b=n-A≦nより、b≦n/3なので、
=Σ[1≦b≦[n/3]](A+b)b　ただし、A=n-3b
（終り）


FOR N=1 TO 50
   PRINT "N="; N

   LET S=0
   FOR b=1 TO INT(N/3) !b≧1として考える
      LET A=N-3*b
      IF A<0 THEN STOP !論理エラー

      PRINT A;b !題意を満たす
      LET S=S+(A+b)*b

   NEXT b
   PRINT "積=";S
NEXT N

END



別解　母関数
(A+b)b=Ab+b^2なので、項Abと項b^2について考える。
項Abのf(n)に対する寄与をf(n:Ab)と表す。
f(n:Ab)の母関数Σ[k=0,∞]f(k:Ab)*x^kについて
　f(n:Ab)は、
　　p[x]=0*1+1*x+2*x^2+3*x^3+4*x^4+5*x^5+6*x^6+7*x^7+ …  =x/(1-x)^2
　　q[x]=0*1+1*x^3+2*x^6+3*x^9+4*x^12+5*x^15+6*x^18+7*x^21+ …  =x^3/(1-x^3)^2
　として、積pqを展開したときのx^nの係数になる。
　よって、x/(1-x)^2 * x^3/(1-x^3)^2 をマクローリン展開したときのx^nの係数になる。

同様に
　f(n:b^2)は、
　　p[x]=1+x+x^2+x^3+x^4+x^5+x^6+x^7+ … =1/(1-x)
　　q[x]=0^2*1+1^2*x^3+2^2*x^6+3^2*x^9+4^2*x^12+5^2*x^15+6^2*x^18+7^2*x^21+ …  =(x^6+x^3)/(1-x^3)^3
　として、積pqを展開したときのx^nの係数になる。
　よって、1/(1-x) * (x^6+x^3)/(1-x^3)^3 をマクローリン展開したときのx^nの係数になる。

したがって、f(n)の母関数は、
Σ[k=0,∞]f(k)*x^k
=x^4/{(1-x)^2(1-x^3)^2} + (x^6+x^3)/{(1-x)(1-x^3)^3}
=x^3(2x^3+x^2+x+1)/{(1-x)(1-x^3)^3}
これをマクローリン展開したときのx^nの係数になる。、
（終り）


!マクローリン展開
!　P(x)=P(0)+{P'(0)/1!}x+{P''(0)/2!}x^2+{P'''(0)/3!}x^3+ …
!例
!P(x)=1/(1-x)=1+x+x^2+x^3+x^4+ …
!
!考察
!P(0)は、x=0を代入して、1/(1-0)=1
!P'(0)は、
!　P(x)=f(x)/g(x)より、P'=(f'g-fg')/g^2なので、
!　1'(1-x)-1(1-x)'/(1-x)^2=1/(1-x)^2
!　x=0を代入して、1
!P''(0)は、
!　P''(x)=(P'(x))'なので、
!　上記の結果をあらためて、f(x)=f'g-fg'、g(x)=g^2と考える。上記と同じ議論を繰り返す。
!（終り）

OPTION ARITHMETIC RATIONAL !有理数モード

PUBLIC NUMERIC N !次数
LET N=50

!変数xの多項式 Σ[k=0,n]a(k)*x^k=a(n)*x^n+a(n-1)*x^(n-1)+ … +a(1)*x+a(0)
!「dim a(0 to n) !係数」で定義する
DIM P(0 TO N)
DIM F(0 TO N),G(0 TO N) !p=f/g
MAT F=ZER
MAT G=ZER

DATA 6 !次数　f=x^3(2x^3+x^2+x+1)=2*x^6+x^5+x^4+x^3　※
DATA 0,0,0,1,1,1,2 !係数　※展開して次数が小さい方から
READ R
FOR i=0 TO R
   READ F(i)
NEXT i
!!!MAT PRINT F; !debug
CALL poly_disp(F) !多項式を表示する
PRINT

DATA 10 !次数　g=(1-x)(1-x^3)^3=x^10-x^9-3x^7+3x^6+3x^4-3x^3-x+1　※
DATA 1,-1,0,-3,3,0,3,-3,0,-1,1 !係数　※展開して次数が小さい方から
READ R
FOR i=0 TO R
   READ G(i)
NEXT i
!!!MAT PRINT G; !debug
CALL poly_disp(G) !多項式を表示する
PRINT


LET P(0)=F(0)/G(0)/FACT(0) !定数項
PRINT 0; P(0) !debug

DIM F1(0 TO N),G1(0 TO N), W1(0 TO N),W2(0 TO N),W3(0 TO N) !作業用
FOR i=1 TO N !i階微分

   CALL poly_diff(F,W3) !f'g
   CALL poly_mul(W3,G, W1)
   !!!MAT PRINT W1; !debug

   CALL poly_diff(G,W3) !fg'
   CALL poly_mul(F,W3, W2)
   !!!MAT PRINT W2; !debug

   CALL poly_sub(W1,W2, F1) !分子 f'g-fg'
   !!!MAT PRINT F1; !debug


   CALL poly_mul(G,G, G1) !分母 g^2
   !!!MAT PRINT G1; !debug


   LET P(i)=F1(0)/G1(0)/FACT(i) !x^iの係数
   PRINT i; P(i) !debug


   MAT F=F1 !次へ
   MAT G=G1
NEXT i

CALL poly_disp(P) !多項式を表示する
PRINT

END


!変数xの多項式 Σ[k=0,n]a(k)*x^k=a(n)*x^n+a(n-1)*x^(n-1)+ … +a(1)*x+a(0)
!「dim a(0 to n) !係数」で定義する

!演算関連

EXTERNAL SUB poly_add(v1(),v2(), v()) !加算 v=v1+v2
OPTION ARITHMETIC RATIONAL !有理数モード
MAT v=v1+v2
END SUB

EXTERNAL SUB poly_sub(v1(),v2(), v()) !減算 v=v1-v2
OPTION ARITHMETIC RATIONAL !有理数モード
MAT v=v1-v2
END SUB

EXTERNAL SUB poly_mul(v1(),v2(), v()) !乗算 v=v1*v2
OPTION ARITHMETIC RATIONAL !有理数モード
DIM w(0 TO 2*N) !桁数は2倍になる
MAT w=ZER
FOR i=0 TO N !係数
   FOR j=0 TO N
      LET w(i+j)=w(i+j)+v1(i)*v2(j) !畳み込み
   NEXT j
NEXT i
FOR i=0 TO N !※下n桁をコピーする　※オーバーフローは考慮していない
   LET v(i)=w(i)
NEXT i
END SUB

EXTERNAL SUB poly_diff(v1(), v()) !微分 v=v1'
OPTION ARITHMETIC RATIONAL !有理数モード
FOR i=1 TO N
   LET v(i-1)=v1(i)*i
NEXT i
LET v(N)=0
END SUB


!表示関連

EXTERNAL SUB poly_disp(A()) !多項式を表示する　a(X)=ΣAkX^k=AnX^n+An-1X^n-1+…+A1X+A0
OPTION ARITHMETIC RATIONAL !有理数モード
CALL mono_disp(A(N),N)
FOR i=N-1 TO 0 STEP -1 !次項
   LET w=A(i)
   IF w>0 THEN PRINT "+";
   IF w<>0 OR N=0 THEN CALL mono_disp(w,i)
NEXT i
END SUB

EXTERNAL SUB mono_disp(ak,k) !単項式を表示する　Ak*X^k
OPTION ARITHMETIC RATIONAL !有理数モード
IF k<>0 THEN !x^nで
   IF ak=0 OR ak=1 THEN !係数が0,1なら
   ELSEIF ak=-1 THEN !係数が－１なら
      PRINT "-"; !符号
   ELSE
      PRINT STR$(ak);"*";
   END IF
END IF
IF k=0 THEN !次数が０なら
   PRINT STR$(ak);
ELSEIF k=1 THEN !次数が１なら
   PRINT "X";
ELSE
   IF ak<>0 THEN PRINT "X^";STR$(k); !係数が0以外なら
END IF
END SUB

 		 

Re: 整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月 5日(火)13時15分2秒
  > No.3009[元記事へ]

> 問題
> nを任意の正の整数とする。
> a≧b≧c≧0、a+2b+3c=n
> を満たす整数の組(a,b,c)のそれぞれに対して、積abcを考える。
> その積のすべての和をf(n)とすると、f(n)をnの式で表せ。

答え
A=a-b、B=b-cとする。a≧b≧cより、A≧0、B≧0
また、A=(n-3c-2b)-b=n-3c-3b=n-3(b-c)-6c=n-3B-6c　∴A+3B+6c=n
Σ[a+2b+3c=n,a≧b≧c≧0]abc
=Σ[A+3B+6c=n,A≧0,B≧0,c≧0](A+B+c)(B+c)c
（終り）


FOR N=1 TO 50
   PRINT "N="; N

   LET S=0
   FOR c=1 TO INT(N/6) !c≧1として考える
      FOR B=0 TO INT((N-6*c)/3) !B≧0として考える
         LET A=N-6*c-3*B
         IF A<0 THEN STOP !論理エラー

         PRINT A;B;c !題意を満たす
         LET S=S+(A+B+c)*(B+c)*c

      NEXT B
   NEXT c
   PRINT "積=";S
NEXT N

END


別解　母関数
(A+B+c)(B+c)c=(AB+B^2)c+(A+2B)c^2+c^3である。
項ABcのf(n)に対する寄与をf(n:ABc)と表す。
変数1a,A,1b,B,B^2,c,c^2,c^3に相当するxの多項式は、
　　1a =1+x+x^2+x^3+x^4+x^5+x^6+x^7+ … =1/(1-x)
　　A =0*1+1x+2x^2+3x^3+4x^4+5x^5+6x^6+7x^7+ …  =x/(1-x)^2
　　1b =1+x^3+x^6+x^9+x^12+x^15+x^18+x^21+ … =1/(1-x^3)
　　B =0*1+1x^3+2x^6+3x^9+4x^12+5x^15+6x^18+7x^21+ …  =x^3/(1-x^3)^2
　　B^2 =0^2*1+1^2x^3+2^2x^6+3^2x^9+4^2x^12+5^2x^15+6^2x^18+7^2x^21+ …  =(x^6+x^3)/(1-x^3)^3
　　c =0*1+1x^6+2x^12+3x^18+4x^24+5x^30+6x^36+7x^42+ …  =x^6/(1-x^6)^2
　　c^2 =0^2*1+1^2x^6+2^2x^12+3^2x^18+4^2x^24+5^2x^30+6^2x^36+7^2x^42+ …  =(x^12+x^6)/(1-x^6)^3
　　c^3 =0^3*1+1^3x^6+2^3x^12+3^3x^18+4^3x^24+5^3x^30+6^3x^36+7^3x^42+ …  =(x^18+4x^12+x^6)/(1-x^6)^4
である。
f(n:ABc)の母関数Σ[k=0,∞]f(k:ABc)*x^kについて
　f(n:ABc)は、積A*B*cをマクローリン展開したときのx^nの係数になる。
同様に
　f(n:B^2c)は、積1a*B^2*c
　f(n:Ac^2)は、積A*1b*c^2
　f(n:Bc^2)は、積1a*B*c^2
　f(n:c^3)は、積1a*1b*c^3

したがって、f(n)の母関数は、
Σ[k=0,∞]f(k)*x^k
=f(n:ABc) + f(n:B^2c) + f(n:Ac^2) + 2*f(n:Bc^2) + f(n:c^3)
= x/(1-x)^2 * x^3/(1-x^3)^2 * x^6/(1-x^6)^2
　+ 1/(1-x) * (x^6+x^3)/(1-x^3)^3 * x^6/(1-x^6)^2
　+ x/(1-x)^2 * 1/(1-x^3) * (x^12+x^6)/(1-x^6)^3
　+ 2 * { 1/(1-x) * x^3/(1-x^3)^2 * (x^12+x^6)/(1-x^6)^3 }
　+ 1/(1-x) * 1/(1-x^3) * (x^18+4x^12+x^6)/(1-x^6)^4
= (6x^20+8x^19+10x^18+12x^17+13x^16+14x^15+17x^14+16x^13+15x^12+8x^11+7x^10+6x^9+3x^8+2x^7+x^6 )
　/ { (1-x^3)^2 * (1-x^6)^4 }
これをマクローリン展開したときのx^nの係数になる。、
（終り）

先のマクローリン展開のプログラムでは、DATA文の箇所を、

DATA 20 !次数　f=6x^20+8x^19+10x^18+12x^17+13x^16+14x^15+17x^14+16x^13+15x^12+8x^11+7x^10+6x^9+3x^8+2x^7+x^6
DATA 0,0,0,0,0,0,1,2,3,6,7,8,15,16,17,14,13,12,10,8,6 !係数　※展開して次数が小さい方から


DATA 30 !次数　g=x^30-2x^27-3x^24+8x^21+2x^18-12x^15+2x^12+8x^9-3x^6-2x^3+1
DATA 1,0,0,-2,0,0,-3,0,0,8,0,0,2,0,0,-12,0,0,2,0,0,8,0,0,-3,0,0,-2,0,0,1 !係数　※展開して次数が小さい方から

とする。
 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月10日(日)09時40分8秒
  > No.3009[元記事へ]

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

山中さんが作成された下記「マクローリン展開」プログラムを、任意関数P(x)のマクローリン展開に利用したのですが、任意関数を

DEF F(x)=....
DEF G(x)=....
DEF P(x)=F(x)/G(x)

と設定して簡単に使える様にできないでしょうか？

> !マクローリン展開
> !　P(x)=P(0)+{P'(0)/1!}x+{P''(0)/2!}x^2+{P'''(0)/3!}x^3+ …
> !例
> !P(x)=1/(1-x)=1+x+x^2+x^3+x^4+ …
> !
> !考察
> !P(0)は、x=0を代入して、1/(1-0)=1
> !P'(0)は、
> !　P(x)=f(x)/g(x)より、P'=(f'g-fg')/g^2なので、
> !　1'(1-x)-1(1-x)'/(1-x)^2=1/(1-x)^2
> !　x=0を代入して、1
> !P''(0)は、
> !　P''(x)=(P'(x))'なので、
> !　上記の結果をあらためて、f(x)=f'g-fg'、g(x)=g^2と考える。上記と同じ議論を繰り返す。
> !（終り）
>
> OPTION ARITHMETIC RATIONAL !有理数モード
>
> PUBLIC NUMERIC N !次数
> LET N=50
>
> !変数xの多項式 Σ[k=0,n]a(k)*x^k=a(n)*x^n+a(n-1)*x^(n-1)+ … +a(1)*x+a(0)
> !「dim a(0 to n) !係数」で定義する
> DIM P(0 TO N)
> DIM F(0 TO N),G(0 TO N) !p=f/g
> MAT F=ZER
> MAT G=ZER
>
> DATA 6 !次数　f=x^3(2x^3+x^2+x+1)=2*x^6+x^5+x^4+x^3　※
> DATA 0,0,0,1,1,1,2 !係数　※展開して次数が小さい方から
> READ R
> FOR i=0 TO R
>    READ F(i)
> NEXT i
> !!!MAT PRINT F; !debug
> CALL poly_disp(F) !多項式を表示する
> PRINT
>
> DATA 10 !次数　g=(1-x)(1-x^3)^3=x^10-x^9-3x^7+3x^6+3x^4-3x^3-x+1　※
> DATA 1,-1,0,-3,3,0,3,-3,0,-1,1 !係数　※展開して次数が小さい方から
> READ R
> FOR i=0 TO R
>    READ G(i)
> NEXT i
> !!!MAT PRINT G; !debug
> CALL poly_disp(G) !多項式を表示する
> PRINT
>
>
> LET P(0)=F(0)/G(0)/FACT(0) !定数項
> PRINT 0; P(0) !debug
>
> DIM F1(0 TO N),G1(0 TO N), W1(0 TO N),W2(0 TO N),W3(0 TO N) !作業用
> FOR i=1 TO N !i階微分
>
>    CALL poly_diff(F,W3) !f'g
>    CALL poly_mul(W3,G, W1)
>    !!!MAT PRINT W1; !debug
>
>    CALL poly_diff(G,W3) !fg'
>    CALL poly_mul(F,W3, W2)
>    !!!MAT PRINT W2; !debug
>
>    CALL poly_sub(W1,W2, F1) !分子 f'g-fg'
>    !!!MAT PRINT F1; !debug
>
>
>    CALL poly_mul(G,G, G1) !分母 g^2
>    !!!MAT PRINT G1; !debug
>
>
>    LET P(i)=F1(0)/G1(0)/FACT(i) !x^iの係数
>    PRINT i; P(i) !debug
>
>
>    MAT F=F1 !次へ
>    MAT G=G1
> NEXT i
>
> CALL poly_disp(P) !多項式を表示する
> PRINT
>
> END
>
>
> !変数xの多項式 Σ[k=0,n]a(k)*x^k=a(n)*x^n+a(n-1)*x^(n-1)+ … +a(1)*x+a(0)
> !「dim a(0 to n) !係数」で定義する
>
> !演算関連
>
> EXTERNAL SUB poly_add(v1(),v2(), v()) !加算 v=v1+v2
> OPTION ARITHMETIC RATIONAL !有理数モード
> MAT v=v1+v2
> END SUB
>
> EXTERNAL SUB poly_sub(v1(),v2(), v()) !減算 v=v1-v2
> OPTION ARITHMETIC RATIONAL !有理数モード
> MAT v=v1-v2
> END SUB
>
> EXTERNAL SUB poly_mul(v1(),v2(), v()) !乗算 v=v1*v2
> OPTION ARITHMETIC RATIONAL !有理数モード
> DIM w(0 TO 2*N) !桁数は2倍になる
> MAT w=ZER
> FOR i=0 TO N !係数
>    FOR j=0 TO N
>       LET w(i+j)=w(i+j)+v1(i)*v2(j) !畳み込み
>    NEXT j
> NEXT i
> FOR i=0 TO N !※下n桁をコピーする　※オーバーフローは考慮していない
>    LET v(i)=w(i)
> NEXT i
> END SUB
>
> EXTERNAL SUB poly_diff(v1(), v()) !微分 v=v1'
> OPTION ARITHMETIC RATIONAL !有理数モード
> FOR i=1 TO N
>    LET v(i-1)=v1(i)*i
> NEXT i
> LET v(N)=0
> END SUB
>
>
> !表示関連
>
> EXTERNAL SUB poly_disp(A()) !多項式を表示する　a(X)=ΣAkX^k=AnX^n+An-1X^n-1+…+A1X+A0
> OPTION ARITHMETIC RATIONAL !有理数モード
> CALL mono_disp(A(N),N)
> FOR i=N-1 TO 0 STEP -1 !次項
>    LET w=A(i)
>    IF w>0 THEN PRINT "+";
>    IF w<>0 OR N=0 THEN CALL mono_disp(w,i)
> NEXT i
> END SUB
>
> EXTERNAL SUB mono_disp(ak,k) !単項式を表示する　Ak*X^k
> OPTION ARITHMETIC RATIONAL !有理数モード
> IF k<>0 THEN !x^nで
>    IF ak=0 OR ak=1 THEN !係数が0,1なら
>    ELSEIF ak=-1 THEN !係数が－１なら
>       PRINT "-"; !符号
>    ELSE
>       PRINT STR$(ak);"*";
>    END IF
> END IF
> IF k=0 THEN !次数が０なら
>    PRINT STR$(ak);
> ELSEIF k=1 THEN !次数が１なら
>    PRINT "X";
> ELSE
>    IF ak<>0 THEN PRINT "X^";STR$(k); !係数が0以外なら
> END IF
> END SUB
>
> 		 

Re: 整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月10日(日)12時28分33秒
  > No.3012[元記事へ]

島村1243さんへのお返事です。

> DEF F(x)=....
> DEF G(x)=....
> DEF P(x)=F(x)/G(x)
> と設定して簡単に使える様にできないでしょうか？

提示したプログラムは、有理式に特化した手法です。
途中で算出される分数は、既約分数ではありませんが、厳密に微分した式になります。
よって、x=0を代入した値は誤差はありません。

一般の関数では、個々の関数に応じて微分しないといけませんので、プログラムは複雑になります。
汎用的な計算では、微分係数などで微分を近似します。したがって、誤差が累積していきます。
中点法、前進法、後退法の1,2,3次などを使えばある程度の精度は保証されますが、期待通りにはなりません。

F(X)=X/(1-X)^2の場合、係数 0,1,2,3,4,5,…

!マクローリン展開にみる数値微分の誤差と計算量 - 再帰呼び出しによるn階微分

LET A=0 !X=0　※マクローリン展開
LET N=9 !N次近似

DIM c(0 TO N) !係数
LET c(0)=DF(A,0) !f(a)
PRINT c(0)

FOR i=1 TO N
   LET c(i)=DF(A,i)/FACT(i)
   PRINT c(i);"X^";STR$(i)
NEXT i

END


EXTERNAL FUNCTION F(X) !関数
!LET F=SIN(X)
!LET F=COS(X)
!LET F=EXP(X)
LET F=X/(1-X)^2 !0,1,2,3,4,5,…
!LET F=1/(1-X) !1,1,1,1,…
!LET F=(1-X)^4 !1,-4,6,-4,1,0,0,…
END FUNCTION

EXTERNAL FUNCTION DF(x,n) !f(x)のn階微分
IF n>0 THEN !微分係数から求める
   LET h=1/1024
   LET DF=(DF(x+h,n-1)-DF(x,n-1))/h
ELSE
   LET DF=F(x)
END IF
END FUNCTION


実行結果

0
1.00195598975279 X^1
2.00881583539961 X^2
3.02355721412682 X^3
4.04920235221334 X^4
5.08883308900731 X^5
6.14561086632151 X^6
7.22195861264053 X^7
8.4220216127928 X^8
7.07402933721543 X^9



初等関数のように、n階微分がnの一般式で定義できる場合は、精度はかなり保証されます。

F(X)=SIN(X)の場合

!マクローリン展開にみる数値微分の誤差と計算量 - 再帰呼び出しによるn階微分

LET A=0 !X=0　※マクローリン展開
LET N=9 !N次近似

DIM c(0 TO N) !係数
LET c(0)=DF(A,0) !f(a)
PRINT c(0)

FOR i=1 TO N
   LET c(i)=DF(A,i)/FACT(i)
   PRINT c(i);"X^";STR$(i)
NEXT i

END


EXTERNAL FUNCTION F(X) !関数
LET F=SIN(X)
!LET F=COS(X)
!LET F=EXP(X)
END FUNCTION

EXTERNAL FUNCTION DF(x,n) !f(x)のn階微分
LET DF=SIN(x+n*PI/2) !f(x)=SIN(x)
!LET DF=COS(x+n*PI/2) !f(x)=COS(x)
!LET DF=EXP(x) !f(x)=EXP(x)
END FUNCTION


実行結果

0
1 X^1
2.31321691639751E-19 X^2
-.166666666666667 X^3
-3.85536152732919E-20 X^4
8.33333333333333E-3 X^5
1.9276807636646E-21 X^6
-1.98412698412698E-4 X^7
-4.58971610396332E-23 X^8
2.75573192239859E-6 X^9

 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月10日(日)13時14分8秒
  山中和義さんへのお返事です。

山中さん、早速のご教示有難う御座いました。

山中さんがご教示くださった
「汎用的な計算では、微分係数などで微分を近似する」方法のプログラムにF(x)=sin(x)と設定し出力したマクローリン係数と

「初等関数のように、n階微分がnの一般式で定義できる場合」の方法で出力されたF(x)=sin(x)のマクローリン係数

とを見比べると、x^5までの係数は概ね一致ですが、x^6以上の項では大きな差が出ており、任意関数を数値計算でマクローリン展開する場合は誤差の累積に注意が必要とのことが良く分かりました。 		 

Re: 整数の組による積の和

 投稿者：しばっち  投稿日：2013年 3月12日(火)21時50分7秒
  > No.3013[元記事へ]

再帰を使った方法より、下記の方が若干精度がいいようです

EXTERNAL  FUNCTION DF(X,K)
LET H=1/1024
FOR J=0 TO K
   LET  S=S+(-1)^J*COMB(K,J)*F(X+(K/2-J)*H)
NEXT J
LET  DF=S/(H^K)
END FUNCTION

また1000桁モードを使用して、下記のようにすればある程度は使えそうです

OPTION ARITHMETIC DECIMAL_HIGH
PUBLIC NUMERIC H
LET X=0
LET N=9
LET EPS=1E-10 !'計算精度
FOR I=0 TO N
   LET H=1/128
   DO
      LET H=H/2
      LET A=DF1(X,I)
      LET B=DF2(X,I)
   LOOP UNTIL ABS(A - B) < EPS
   LET DF=(A+B)/2
   PRINT DF/FACT(I)
NEXT I
END

EXTERNAL FUNCTION F(X)
OPTION ARITHMETIC DECIMAL_HIGH
LET F=X/(1-X)^2
END FUNCTION

EXTERNAL FUNCTION DF1(X,N)
OPTION ARITHMETIC DECIMAL_HIGH
IF N>0 THEN
   LET DF1=(DF1(X-2*H,N-1)-4*DF1(X-H,N-1)+3*DF1(X,N-1))/(2*H) !'3点前進法
   !'LET DF1=(3*DF1(X-4*H,N-1)-16*DF1(X-3*H,N-1)+36*DF1(X-2*H,N-1)-48*DF1(X-H,N-1)+25*DF1(X,N-1))/(12*H) !'5点前進法
ELSE
   LET DF1=F(X)
END IF
END FUNCTION

EXTERNAL FUNCTION DF2(X,N)
OPTION ARITHMETIC DECIMAL_HIGH
IF N>0 THEN
   LET DF2=(-3*DF2(X,N-1)+4*DF2(X+H,N-1)-DF2(X+2*H,N-1))/(2*H) !'3点後退法
   !'LET DF2=(-25*DF2(X,N-1)+48*DF2(X+H,N-1)-36*DF2(X+2*H,N-1)+16*DF2(X+3*H,N-1)-3*DF2(X+4*H,N-1))/(12*H) !'5点後退法
ELSE
   LET DF2=F(X)
END IF
END FUNCTION
 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月13日(水)09時25分4秒
  > No.3015[元記事へ]

しばっちさんへのお返事です。

しばっちさん、ご教示有難う御座いました。
関数　F(X)=X/(1-X)^2
についてご提示のあった二つのプロ（理屈は理解不能）を実行してみました。
---------------------------------------------------------
> 再帰を使った方法より、下記の方が若干精度がいいようです
---------------------------------------------------------
このプロの場合は結果が
0
1.00000071525602 X^1
2.00000381470259 X^2
3.00001192084544 X^3
4.00002806202375 X^4
5.00015669041249 X^5
6.01465745988539 X^6
24.1746473103899 X^7
-2263.61524295701 X^8
303480.118194826 X^9
となり、X^6の係数までは精度が高かったですが、X^7以降の係数は当てにならない結果でした。演算時間は速いので６次（必ずそう言えるか不明な点が気になるが）以下の係数には使えるようです。

----------------------------------------------------------------------
> また1000桁モードを使用して、下記のようにすればある程度は使えそうです
----------------------------------------------------------------------
上記２番目のプロは、９次までの係数が極めて高い精度で得られたので、幾らの次数まで正しい値かということを意識せず？に使えそうです。
ただし、私のコンピュータは遅い（AMD-AthlonXP2200+ 2.2GHz相当)ので８次以上の結果が出るまで５～６分必要としました。 		 

Re: 整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月13日(水)19時05分11秒
  > No.3014[元記事へ]

島村1243さんへのお返事です。

5次程度なら、次の連立方程式を解けば求まります。


!マクローリン展開
!f(x)=f(0)+{f'(0)/1!}x+{f''(0)/2!}x^2+{f'''(0)/3!}x^3+ …
!    =f(0)+a[1]x+a[2]x^2+a[3]x^3+a[4]x^4+ …

!DEF F(X)=SIN(X)
!DEF F(X)=COS(X)
!DEF F(X)=EXP(X)
DEF F(X)=X/(1-X)^2 !0,1,2,3,4,5,…
!DEF F(X)=1/(1-X) !1,1,1,1,1,…
!DEF F(x)=(1-X)^4 !1,-4,6,-4,1,0,0,…


LET N=9 !次数

!A=
!　 x   x^2   x^3    x^4　…
!　2x  4x^2  8x^3  16x^4　…
!　3x  9x^2 27x^3  81x^4　…
!　4x 16x^2 64x^3 256x^4　…
!　　　：

LET x=1E-3 !x≒0

DIM A(N,N)
FOR i=1 TO N
   FOR J=1 TO N
      LET A(i,J)=(i*x)^J
   NEXT J
NEXT i
!!!MAT PRINT A; !debug

!b=
!　f( x)-f(0)
!　f(2x)-f(0)
!　f(3x)-f(0)
!　f(4x)-f(0)
!　　　：
DIM b(N)
LET F0=F(0)
FOR i=1 TO N
   LET b(i)=F(i*x)-F0
NEXT i

DIM p(N),iA(N,N) !連立方程式Ap=bを解く
MAT iA=INV(A)
MAT p=iA*b

PRINT F(0) !各係数　定数項
MAT PRINT p; !x,x^2,x^3,…の項

END


実行結果

0
1.00000000000017  1.99999999953607  3.00000058623451  3.999024242  5.03640518 -18.75867688  1689.07931105014 -130559.084358461  2609878


 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月14日(木)13時16分45秒
  山中和義さんへのお返事です。

山中さん、度重なるご教示有難う御座いました。

> 5次程度なら、次の連立方程式を解けば求まります。

上記プログラムの考え方も記載して頂き、私の頭でも「なるほどな～」と良く分かる計算法で感心するばかりでした。

コード中に「ｘ=1E-3 !x≒0」の箇所が有りましたが、ｘ=1E-2、ｘ=1E-4、等と変えてみましたが、ｘを小さくすれば精度が上がると言う事でもない様ですね。
（この「ｘ=1E-3 !x≒0」は「行列Aの高次係数が大きくなり過ぎない様にする為」という理由でしょうか？）

さてこのプログラムで「F(X)=X/(1-X)^2」を、ｘ=1E-3のまま1000桁モードで計算させたら瞬時に９次までの高精度結果が得られました。難解な数値計算理論を使っていないのに、この方法は凄いですね。





 		 

Re: 整数の組による積の和

 投稿者：山中和義  投稿日：2013年 3月14日(木)14時20分0秒
  > No.3018[元記事へ]

島村1243さんへのお返事です。

> コード中に「ｘ=1E-3 !x≒0」の箇所が有りましたが、ｘ=1E-2、ｘ=1E-4、等と変えてみましたが、ｘを小さくすれば精度が上がると言う事でもない様ですね。

1E-3あたりが実用値でしょう。


> （この「ｘ=1E-3 !x≒0」は「行列Aの高次係数が大きくなり過ぎない様にする為」という理由でしょか？）

提示したプログラムは、X=x,2x,3x,…を代入して、これを満たす補間多項式を得ます。
x=0は、0を代入するのですが、x^kがすべて0になって、行列Aが成り立ちません。

Xは、テイラー展開におけるX=aと考えることもできます。
このことは、プログラムの最後に、

!F(x)の値
LET S=F(0)
FOR i=1 TO N
   LET S=S+p(i)*x^i
NEXT i
PRINT S; F(x)

を追加して、
先頭部分の関数定義とXの値を（別の関数と値でよいのですが）、

　DEF F(X)=SIN(X)

　LET x=PI/3

と変更して、確認できます。
 		 

Re: 整数の組による積の和

 投稿者：しばっち  投稿日：2013年 3月14日(木)21時36分13秒
  高階微分の計算に時間がかかるのは、再帰ルーチン内で多重呼び出しを行っているためです
これを直接(10階微分まで)求めるようにすれば速くなります
(精度については注意を払う必要がある)

OPTION ARITHMETIC DECIMAL_HIGH
PUBLIC NUMERIC H
LET X=0
LET EPS=1E-20 !'精度
FOR I=1 TO 10
   LET H=1/2^10
   DO
      LET H=H/2
      SELECT CASE I
      CASE 1
         LET A=(-25*F(X)+48*F(X+H)-36*F(X+2*H)+16*F(X+3*H)-3*F(X+4*H))/(12*H)
         LET B=(3*F(X-4*H)-16*F(X-3*H)+36*F(X-2*H)-48*F(X-H)+25*F(X))/(12*H)
      CASE 2
         LET A=(F(X)-2*F(X+H)+F(X+2*H))/(H^2)
         LET B=(F(X-2*H)-2*F(X-H)+F(X))/(H^2)
      CASE 3
         LET A=(-49*F(X)+232*F(X+H)-461*F(X+2*H)+496*F(X+3*H)-307*F(X+4*H)+104*F(X+5*H)-15*F(X+6*H))/(8*H^3)
         LET B=(15*F(X-6*H)-104*F(X-5*H)+307*F(X-4*H)-496*F(X-3*H)+461*F(X-2*H)-232*F(X-H)+49*F(X))/(8*H^3)
      CASE 4
         LET A=(F(X)-4*F(X+H)+6*F(X+2*H)-4*F(X+3*H)+F(X+4*H))/(H^4)
         LET B=(F(X-4*H)-4*F(X-3*H)+6*F(X-2*H)-4*F(X-H)+F(X))/(H^4)
      CASE 5
         LET A=(-81*F(X)+575*F(X+H)-1790*F(X+2*H)+3195*F(X+3*H)-3580*F(X+4*H)+2581*F(X+5*H)-1170*F(X+6*H)+305*F(X+7*H)-35*F(X+8*H))/(6*H^5)
         LET B=(35*F(X-8*H)-305*F(X-7*H)+1170*F(X-6*H)-2581*F(X-5*H)+3580*F(X-4*H)-3195*F(X-3*H)+1790*F(X-2*H)-575*F(X-H)+81*F(X))/(6*H^5)
      CASE 6
         LET A=(F(X)-6*F(X+H)+15*F(X+2*H)-20*F(X+3*H)+15*F(X+4*H)-6*F(X+5*H)+F(X+6*H))/(H^6)
         LET B=(F(X-6*H)-6*F(X-5*H)+15*F(X-4*H)-20*F(X-3*H)+15*F(X-2*H)-6*F(X-H)+F(X))/(H^6)
      CASE 7
         LET A=(-605*F(X)+5628*F(X+H)-23583*F(X+2*H)+58632*F(X+3*H)-95802*F(X+4*H)+107520*F(X+5*H)-83958*F(X+6*H)+45048*F(X+7*H)-15897*F(X+8*H)+3332*F(X+9*H)-315*F(X+10*H))/(24*H^7)
         LET B=(315*F(X-10*H)-3332*F(X-9*H)+15897*F(X-8*H)-45048*F(X-7*H)+83958*F(X-6*H)-107520*F(X-5*H)+95802*F(X-4*H)-58632*F(X-3*H)+23583*F(X-2*H)-5628*F(X-H)+605*F(X))/(24*H^7)
      CASE 8
         LET A=(F(X)-8*F(X+H)+28*F(X+2*H)-56*F(X+3*H)+70*F(X+4*H)-56*F(X+5*H)+28*F(X+6*H)-8*F(X+7*H)+F(X+8*H))/(H^8)
         LET B=(F(X-8*H)-8*F(X-7*H)+28*F(X-6*H)-56*F(X-5*H)+70*F(X-4*H)-56*F(X-3*H)+28*F(X-2*H)-8*F(X-H)+F(X))/(H^8)
      CASE 9
         LET A=(-169*F(X)+1932*F(X+H)-10128*F(X+2*H)+32196*F(X+3*H)-69129*F(X+4*H)+105624*F(X+5*H)-117768*F(X+6*H)+96552*F(X+7*H)-57771*F(X+8*H)+24604*F(X+9*H)-7080*F(X+10*H)+1236*F(X+11*H)-99*F(X+12*H))/(4*H^9)
         LET B=(99*F(X-12*H)-1236*F(X-11*H)+7080*F(X-10*H)-24604*F(X-9*H)+57771*F(X-8*H)-96552*F(X-7*H)+117768*F(X-6*H)-105624*F(X-5*H)+69129*F(X-4*H)-32196*F(X-3*H)+10128*F(X-2*H)-1932*F(X-H)+169*F(X))/(4*H^9)
      CASE 10
         LET A=(F(X)-10*F(X+H)+45*F(X+2*H)-120*F(X+3*H)+210*F(X+4*H)-252*F(X+5*H)+210*F(X+6*H)-120*F(X+7*H)+45*F(X+8*H)-10*F(X+9*H)+F(X+10*H))/(H^10)
         LET B=(F(X-10*H)-10*F(X-9*H)+45*F(X-8*H)-120*F(X-7*H)+210*F(X-6*H)-252*F(X-5*H)+210*F(X-4*H)-120*F(X-3*H)+45*F(X-2*H)-10*F(X-H)+F(X))/(H^10)
      END SELECT
   LOOP UNTIL ABS(A - B) < EPS
   LET DF=(A+B)/2
   PRINT DF/FACT(I)
NEXT I
END

EXTERNAL FUNCTION F(X)
OPTION ARITHMETIC DECIMAL_HIGH
LET F=X/(1-X)^2
END FUNCTION
 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月15日(金)16時25分44秒
  > No.3019[元記事へ]

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

山中さん、ご教示有難うございます。

> Xは、テイラー展開におけるX=aと考えることもできます。

ご指示頂いた様に、x=1E-3ではなく、x=pi/3と設定しF(x)=sin(x)について確認しましたら、ご教示のとおりになりました。
上記「Xはテーラー展開における・・・」の意味は下記で解釈すると納得出来たのですが誤りでしょうか？

X=aの近傍X=a+hに於けるF(a+h)はテイラー展開によって
  F(a+h)=F(a)+F'(a)*h+F''(a)*h^2/2!+・・・
だから、これを
  F(a+h)-F(a)=F'(a)*h+F''(a)*h^2/2!+・・・
と変形して「a=0,h=pi/3」と置けば、「多元連立式[A][p]=[b]を[p]について解きマクローリン係数を求めるプログラム」を利用してF(pi/3)が得られる。 		 

Re: 整数の組による積の和

 投稿者：島村1243  投稿日：2013年 3月16日(土)14時48分27秒
  > No.3020[元記事へ]

しばっちさんへのお返事です。

> 高階微分の計算に時間がかかるのは、再帰ルーチン内で多重呼び出しを行っているためです
> これを直接(10階微分まで)求めるようにすれば速くなります
> (精度については注意を払う必要がある)

ご教示頂いたコードを走らせて、瞬時に結果が得られることを確認しました。
どうも有難うございました。 		 

戻る