Skip to content

GitLab

Commit 3b557d56 authored by ulrich_y's avatar ulrich_y
Browse files

Simplified bernoulli polynomials

parent bedbe8d8
Hide whitespace changes
Inline Side-by-side
Showing with 47 additions and 36 deletions
src/maths_functions.f90
View file @ 3b557d56
......@@ -266,44 +266,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
......
  • ulrich_y @ulrich_y ·
    Owner

    This was generated in Mathematica using

    coeff = CoefficientList[BernoulliB[#, x], x, 16] & /@ Range[15];


ClearAll[MyToString]
MyToString[n_Integer] /; n < 0 := "-" <> MyToString[-n, 1]
MyToString[n_Integer] /; n > 0 := "+" <> MyToString[n, 1]
MyToString[0] := " " <> MyToString[0, 1]
MyToString[n_Integer, 1] := StringPadLeft[ToString[n], 4]

MakeMatrix[mat_, name_: ""] := StringJoin[
  "    integer, parameter :: coeff" <> name,
    "(" ,
      StringRiffle[ToString /@ {Length[mat[[1]]], Length[mat]},", "],
    ") = reshape((/ &\n" ,
      StringRiffle[Map[
        "        " <> StringRiffle[#, ", "] &,
        Map[MyToString, mat, {2}]
      ], ", &\n"],
    " /), &\n            (/",
       StringRiffle[ToString /@ {Length[mat[[1]]], Length[mat]},", "],
     "/))"]


MakeMatrix[Numerator[coeff], "N"] <> "\n" <> MakeMatrix[Denominator[coeff], "D"]
  • ulrich_y @ulrich_y

    mentioned in commit a7004f12

     ·

    mentioned in commit a7004f12

    Toggle commit list
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment