Commit a7004f12 by ulrich_y

### Simplified bernoulli polynomials

`(cherry picked from commit 3b557d56)`
parent 0f92e2f1
 ... ... @@ -253,44 +253,55 @@ CONTAINS END FUNCTION trilog FUNCTION BERNOULLI_POLYNOMIAL(n, x) result(res) integer, parameter :: maxn = 15 integer n complex(kind=prec) :: x, res select case(n) case(1) res = -1/2. + x case(2) res = 1/6. - x + x**2 case(3) res = x/2. - 3*x**2/2. + x**3 case(4) res = -1/30. + x**2 - 2*x**3 + x**4 case(5) res = -x/6. + 5*x**3/3 - 5*x**4/2 + x**5 case(6) res = 1/42. - x**2/2 + 5*x**4/2 - 3*x**5 + x**6 case(7) res = x/6. - 7*x**3/6 + 7*x**5/2 - 7*x**6/2 + x**7 case(8) res = -1/30. + 2*x**2/3 - 7*x**4/3 + 14*x**6/3 - 4*x**7 + x**8 case(9) res = -3*x/10 + 2*x**3 - 21*x**5/5 + 6*x**7 - 9*x**8/2 + x**9 case(10) res = 5/66. - 3*x**2/2 + 5*x**4 - 7*x**6 + 15*x**8/2 - 5*x**9 + x**10 case(11) res = 5*x/6 - 11*x**3/2 + 11*x**5 - 11*x**7 + 55*x**9/6 - 11*x**10/2 + x**11 case(12) res = -691/2730. + 5*x**2 - 33*x**4/2 + 22*x**6 - 33*x**8/2 + 11*x**10 - 6*x**11 + x**12 case(13) res = -691*x/210 + 65*x**3/3 - 429*x**5/10 + 286*x**7/7 - 143*x**9/6 + 13*x**11 - 13*x**12/2 + x**13 case(14) res = 7/6. - 691*x**2/30 + 455*x**4/6 - 1001*x**6/10 + 143*x**8/2 - 1001*x**10/30 + 91*x**12/6 - 7*x**13 + x**14 case(15) res = 35*x/2 - 691*x**3/6 + 455*x**5/2 - 429*x**7/2 + 715*x**9/6 - 91*x**11/2 + 35*x**13/2 - 15*x**14/2 + x**15 case default print*,"Bernoulli beyond 15 is not implemented" stop end select restype :: x, res restype :: xpow(maxn+1) integer, parameter :: coeffN(maxn+1, maxn) = reshape((/ & - 1, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, - 1, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, + 1, - 3, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 0, + 1, - 2, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, - 1, 0, + 5, - 5, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, & + 1, 0, - 1, 0, + 5, - 3, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, & 0, + 1, 0, - 7, 0, + 7, - 7, + 1, 0, 0, 0, 0, 0, 0, 0, 0, & - 1, 0, + 2, 0, - 7, 0, + 14, - 4, + 1, 0, 0, 0, 0, 0, 0, 0, & 0, - 3, 0, + 2, 0, - 21, 0, + 6, - 9, + 1, 0, 0, 0, 0, 0, 0, & + 5, 0, - 3, 0, + 5, 0, - 7, 0, + 15, - 5, + 1, 0, 0, 0, 0, 0, & 0, + 5, 0, - 11, 0, + 11, 0, - 11, 0, + 55, - 11, + 1, 0, 0, 0, 0, & - 691, 0, + 5, 0, - 33, 0, + 22, 0, - 33, 0, + 11, - 6, + 1, 0, 0, 0, & 0, - 691, 0, + 65, 0, - 429, 0, + 286, 0, - 143, 0, + 13, - 13, + 1, 0, 0, & + 7, 0, - 691, 0, + 455, 0, -1001, 0, + 143, 0, -1001, 0, + 91, - 7, + 1, 0, & 0, + 35, 0, - 691, 0, + 455, 0, - 429, 0, + 715, 0, - 91, 0, + 35, - 15, + 1 /), & (/maxn+1, maxn/)) integer, parameter :: coeffD(maxn+1, maxn) = reshape((/ & + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 6, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 2, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 30, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 6, + 1, + 3, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 42, + 1, + 2, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 6, + 1, + 6, + 1, + 2, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 30, + 1, + 3, + 1, + 3, + 1, + 3, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 10, + 1, + 1, + 1, + 5, + 1, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 66, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 6, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 6, + 2, + 1, + 1, + 1, + 1, + 1, & +2730, + 1, + 1, + 1, + 2, + 1, + 1, + 1, + 2, + 1, + 1, + 1, + 1, + 1, + 1, + 1, & + 1, + 210, + 1, + 3, + 1, + 10, + 1, + 7, + 1, + 6, + 1, + 1, + 2, + 1, + 1, + 1, & + 6, + 1, + 30, + 1, + 6, + 1, + 10, + 1, + 2, + 1, + 30, + 1, + 6, + 1, + 1, + 1, & + 1, + 2, + 1, + 6, + 1, + 2, + 1, + 2, + 1, + 6, + 1, + 2, + 1, + 2, + 2, + 1 /), & (/maxn+1, maxn/)) real(kind=prec), parameter :: coeff(maxn+1,maxn) = coeffN/real(coeffD) integer i if (n>maxn) then print*,"Bernoulli beyond 15 is not implemented" stop endif xpow(1:n+1) = (/ ( x**i, i = 0, n ) /) res = sum( xpow(1:n+1) * coeff(1:n+1,n) ) END FUNCTION ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!