Skip to content

Commit 3d23a9f

Browse files
ExpanderExpander
committed
addin Fortran implementation of Clausen function Cl2
1 parent aa9fd64 commit 3d23a9f

File tree

6 files changed

+145
-0
lines changed

6 files changed

+145
-0
lines changed

‎src/fortran/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@ if(CMAKE_Fortran_COMPILER)
22
add_library(polylog_fortran
33
fortran_wrappers.f90
44
fast_clog.f90
5+
Cl2.f90
56
Li2.f90
67
Li3.f90
78
Li4.f90

‎src/fortran/Cl2.f90

+89
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
!*********************************************************************
2+
! This file is part of Polylogarithm.
3+
!
4+
! Polylogarithm is licenced under the MIT License.
5+
!*********************************************************************
6+
7+
!*********************************************************************
8+
!> @brief Clausen function \f$\mathrm{Cl}_2(\theta) = \mathrm{Im}(\mathrm{Li}_2(e^{i\theta}))\f$
9+
!> @param x real angle
10+
!> @return \f$\mathrm{Cl}_2(\theta)\f$
11+
!> @author Alexander Voigt
12+
!> Implemented as economized Padé approximation.
13+
!*********************************************************************
14+
15+
double precision function dcl2(x)
16+
implicit none
17+
double precision :: x, y, z, z2, z4, p, q, p0, p1, h, sgn
18+
double precision, parameter :: PI = 3.14159265358979324D0
19+
double precision, parameter :: PI2 = 2*PI, PIH = PI/2, PI28 = PI*PI/8
20+
double precision, parameter :: cp(4) = (/ &
21+
2.7951565822419270D-2, &
22+
-8.8865360514541522D-4, &
23+
6.8282348222485902D-6, &
24+
-7.5276232403566808D-9 /)
25+
double precision, parameter :: cq(4) = (/ &
26+
1.0000000000000000D+0, &
27+
-3.6904397961160525D-2, &
28+
3.7342870576106476D-4, &
29+
-8.7460760866531179D-7 /)
30+
double precision, parameter :: cr(6) = (/ &
31+
6.4005702446195512D-1, &
32+
-2.0641655351338783D-1, &
33+
2.4175305223497718D-2, &
34+
-1.2355955287855728D-3, &
35+
2.5649833551291124D-5, &
36+
-1.4783829128773320D-7 /)
37+
double precision, parameter :: cs(6) = (/ &
38+
1.0000000000000000D+0, &
39+
-2.5299102015666356D-1, &
40+
2.2148751048467057D-2, &
41+
-7.8183920462457496D-4, &
42+
9.5432542196310670D-6, &
43+
-1.8184302880448247D-8 /)
44+
45+
sgn = 1
46+
47+
if (x .lt. 0) then
48+
x = -x
49+
sgn = -1
50+
endif
51+
52+
if (x .ge. PI2) then
53+
x = mod(x, PI2)
54+
endif
55+
56+
if (x .gt. PI) then
57+
p0 = 6.28125D0
58+
p1 = 0.0019353071795864769253D0
59+
x = (p0 - x) + p1
60+
sgn = -sgn
61+
endif
62+
63+
if (x .eq. 0 .or. x .eq. PI) then
64+
dcl2 = 0
65+
return
66+
endif
67+
68+
if (x .lt. PIH) then
69+
y = x*x
70+
z = y - PI28
71+
z2 = z*z
72+
p = cp(1) + z * cp(2) + z2 * (cp(3) + z * cp(4))
73+
q = cq(1) + z * cq(2) + z2 * (cq(3) + z * cq(4))
74+
h = x*(1 - log(x) + y*p/q/2)
75+
else
76+
y = PI - x
77+
z = y*y - PI28
78+
z2 = z*z
79+
z4 = z2*z2
80+
p = cr(1) + z * cr(2) + z2 * (cr(3) + z * cr(4)) + &
81+
z4 * (cr(5) + z * cr(6))
82+
q = cs(1) + z * cs(2) + z2 * (cs(3) + z * cs(4)) + &
83+
z4 * (cs(5) + z * cs(6))
84+
h = y*p/q
85+
endif
86+
87+
dcl2 = sgn*h
88+
89+
end function dcl2

‎src/fortran/fortran_wrappers.f90

+10
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,16 @@
55
!*********************************************************************
66

77

8+
subroutine cl2_fortran(x, res) bind(C)
9+
use, intrinsic :: iso_c_binding
10+
implicit none
11+
real(c_double), intent(in) :: x
12+
real(c_double), intent(out) :: res
13+
double precision dcl2
14+
res = dcl2(x)
15+
end subroutine cl2_fortran
16+
17+
818
subroutine li2_fortran(x, res) bind(C)
919
use, intrinsic :: iso_c_binding
1020
implicit none

‎src/fortran/fortran_wrappers.h

+3
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,9 @@ extern "C" {
66

77
#ifdef ENABLE_FORTRAN
88

9+
/** Clausen function with n=2, Fortran implementation */
10+
void cl2_fortran(const double* x, double* res);
11+
912
/** real polylogarithm with n=2 (dilogarithm), Fortran implementation */
1013
void li2_fortran(const double* x, double* res);
1114

‎test/bench_Cl.cpp

+16
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#include "alt.h"
22
#include "bench.hpp"
3+
#include "fortran_wrappers.h"
34
#include "Cl2.hpp"
45
#include "Cl3.hpp"
56
#include "Cl4.hpp"
@@ -15,6 +16,16 @@
1516

1617
namespace {
1718

19+
#ifdef ENABLE_FORTRAN
20+
21+
double poly_Cl2_fortran(double x) {
22+
double res{};
23+
cl2_fortran(&x, &res);
24+
return res;
25+
}
26+
27+
#endif
28+
1829
double Cl2_via_Li2(double x) noexcept
1930
{
2031
return std::imag(polylogarithm::Li2(std::polar(1.0, x)));
@@ -75,6 +86,11 @@ int main() {
7586
bench_fn([&](double x) { return polylogarithm::Cl2(x); }, values_d,
7687
"polylogarithm C++", "double");
7788

89+
#ifdef ENABLE_FORTRAN
90+
bench_fn([&](double x) { return poly_Cl2_fortran(x); }, values_d,
91+
"polylogarithm Fortran", "double");
92+
#endif
93+
7894
bench_fn([&](double x) { return Cl2_via_Li2(x); }, values_d,
7995
"via Li2 C++", "double");
8096

‎test/test_Cl2.cpp

+26
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22

33
#include "doctest.h"
44
#include "alt.h"
5+
#include "fortran_wrappers.h"
56
#include "Cl2.hpp"
67
#include "Li2.hpp"
78
#include "read_data.hpp"
@@ -19,6 +20,16 @@
1920

2021
namespace {
2122

23+
#ifdef ENABLE_FORTRAN
24+
25+
double poly_Cl2_fortran(double x) {
26+
double res{};
27+
cl2_fortran(&x, &res);
28+
return res;
29+
}
30+
31+
#endif
32+
2233
double Cl2_via_Li2(double x) noexcept
2334
{
2435
return std::imag(polylogarithm::Li2(std::polar(1.0, x)));
@@ -102,6 +113,9 @@ TEST_CASE("test_real_fixed_values")
102113
const auto cl64_koelbig = clausen_2_koelbig(x64);
103114
const auto cl64_pade = clausen_2_pade(x64);
104115
const auto cl64_poly = polylogarithm::Cl2(x64);
116+
#ifdef ENABLE_FORTRAN
117+
const auto cl64_poly_f = poly_Cl2_fortran(x64);
118+
#endif
105119
const auto cl64_li2 = Cl2_via_Li2(x64);
106120
const auto cl64_wu = clausen_2_wu(x64);
107121
const auto cl128_poly = polylogarithm::Cl2(x128);
@@ -115,6 +129,9 @@ TEST_CASE("test_real_fixed_values")
115129
INFO("x(64) = " << x64);
116130
INFO("Cl2(64) real = " << cl64_expected << " (expected)");
117131
INFO("Cl2(64) real = " << cl64_poly << " (polylogarithm C++)");
132+
#ifdef ENABLE_FORTRAN
133+
INFO("Cl2(64) real = " << cl64_poly_f << " (polylogarithm Fortran)");
134+
#endif
118135
INFO("Cl2(64) real = " << cl64_li2 << " (via Li2 C++)");
119136
#ifdef ENABLE_GSL
120137
INFO("Cl2(64) real = " << cl64_gsl << " (GSL)");
@@ -137,6 +154,15 @@ TEST_CASE("test_real_fixed_values")
137154
} else {
138155
CHECK_CLOSE(cl64_poly , cl64_expected , 100*eps64);
139156
}
157+
#ifdef ENABLE_FORTRAN
158+
if (std::abs(x64 - 2*pi64) > 1e-2) {
159+
CHECK_CLOSE(cl64_poly_f , cl64_expected , 2*eps64);
160+
} else if (std::abs(x64 - 2*pi64) > 1e-12) {
161+
CHECK_CLOSE(cl64_poly_f , cl64_expected , 10*eps64);
162+
} else {
163+
CHECK_CLOSE(cl64_poly_f , cl64_expected , 100*eps64);
164+
}
165+
#endif
140166
CHECK_CLOSE(cl64_li2 , cl64_expected , 10*eps64);
141167
if (std::abs(x64 - 2*pi64) > 1e-3) {
142168
CHECK_CLOSE(cl64_bernoulli, cl64_expected , 2*eps64);

0 commit comments

Comments
 (0)