]>
Commit | Line | Data |
---|---|---|
20effc67 TL |
1 | // Jacobi theta functions |
2 | // Copyright Evan Miller 2020 | |
3 | // | |
4 | // Use, modification and distribution are subject to the | |
5 | // Boost Software License, Version 1.0. (See accompanying file | |
6 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
7 | // | |
8 | // Four main theta functions with various flavors of parameterization, | |
9 | // floating-point policies, and bonus "minus 1" versions of functions 3 and 4 | |
10 | // designed to preserve accuracy for small q. Twenty-four C++ functions are | |
11 | // provided in all. | |
12 | // | |
13 | // The functions take a real argument z and a parameter known as q, or its close | |
14 | // relative tau. | |
15 | // | |
16 | // The mathematical functions are best understood in terms of their Fourier | |
1e59de90 | 17 | // series. Using the q parameterization, and summing from n = 0 to INF: |
20effc67 | 18 | // |
1e59de90 TL |
19 | // theta_1(z,q) = 2 SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z) |
20 | // theta_2(z,q) = 2 SUM q^(n+1/2)^2 * cos((2n+1)z) | |
21 | // theta_3(z,q) = 1 + 2 SUM q^n^2 * cos(2nz) | |
22 | // theta_4(z,q) = 1 + 2 SUM (-1)^n * q^n^2 * cos(2nz) | |
20effc67 TL |
23 | // |
24 | // Appropriately multiplied and divided, these four theta functions can be used | |
25 | // to implement the famous Jacabi elliptic functions - but this is not really | |
26 | // recommended, as the existing Boost implementations are likely faster and | |
27 | // more accurate. More saliently, setting z = 0 on the fourth theta function | |
28 | // will produce the limiting CDF of the Kolmogorov-Smirnov distribution, which | |
1e59de90 | 29 | // is this particular implementation's raison d'etre. |
20effc67 TL |
30 | // |
31 | // Separate C++ functions are provided for q and for tau. The main q functions are: | |
32 | // | |
33 | // template <class T> inline T jacobi_theta1(T z, T q); | |
34 | // template <class T> inline T jacobi_theta2(T z, T q); | |
35 | // template <class T> inline T jacobi_theta3(T z, T q); | |
36 | // template <class T> inline T jacobi_theta4(T z, T q); | |
37 | // | |
38 | // The parameter q, also known as the nome, is restricted to the domain (0, 1), | |
39 | // and will throw a domain error otherwise. | |
40 | // | |
41 | // The equivalent functions that use tau instead of q are: | |
42 | // | |
43 | // template <class T> inline T jacobi_theta1tau(T z, T tau); | |
44 | // template <class T> inline T jacobi_theta2tau(T z, T tau); | |
45 | // template <class T> inline T jacobi_theta3tau(T z, T tau); | |
46 | // template <class T> inline T jacobi_theta4tau(T z, T tau); | |
47 | // | |
1e59de90 | 48 | // Mathematically, q and tau are related by: |
20effc67 | 49 | // |
1e59de90 | 50 | // q = exp(i PI*Tau) |
20effc67 | 51 | // |
1e59de90 TL |
52 | // However, the tau in the equation above is *not* identical to the tau in the function |
53 | // signature. Instead, `tau` is the imaginary component of tau. Mathematically, tau can | |
54 | // be complex - but practically, most applications call for a purely imaginary tau. | |
20effc67 TL |
55 | // Rather than provide a full complex-number API, the author decided to treat the |
56 | // parameter `tau` as an imaginary number. So in computational terms, the | |
57 | // relationship between `q` and `tau` is given by: | |
58 | // | |
59 | // q = exp(-constants::pi<T>() * tau) | |
60 | // | |
61 | // The tau versions are provided for the sake of accuracy, as well as conformance | |
62 | // with common notation. If your q is an exponential, you are better off using | |
63 | // the tau versions, e.g. | |
64 | // | |
65 | // jacobi_theta1(z, exp(-a)); // rather poor accuracy | |
66 | // jacobi_theta1tau(z, a / constants::pi<T>()); // better accuracy | |
67 | // | |
68 | // Similarly, if you have a precise (small positive) value for the complement | |
69 | // of q, you can obtain a more precise answer overall by passing the result of | |
70 | // `log1p` to the tau parameter: | |
71 | // | |
72 | // jacobi_theta1(z, 1-q_complement); // precision lost in subtraction | |
73 | // jacobi_theta1tau(z, -log1p(-q_complement) / constants::pi<T>()); // better! | |
74 | // | |
75 | // A third quartet of functions are provided for improving accuracy in cases | |
1e59de90 | 76 | // where q is small, specifically |q| < exp(-PI) = 0.04 (or, equivalently, tau |
20effc67 TL |
77 | // greater than unity). In this domain of q values, the third and fourth theta |
78 | // functions always return values close to 1. So the following "m1" functions | |
79 | // are provided, similar in spirit to `expm1`, which return one less than their | |
80 | // regular counterparts: | |
81 | // | |
82 | // template <class T> inline T jacobi_theta3m1(T z, T q); | |
83 | // template <class T> inline T jacobi_theta4m1(T z, T q); | |
84 | // template <class T> inline T jacobi_theta3m1tau(T z, T tau); | |
85 | // template <class T> inline T jacobi_theta4m1tau(T z, T tau); | |
86 | // | |
87 | // Note that "m1" versions of the first and second theta would not be useful, | |
88 | // as their ranges are not confined to a neighborhood around 1 (see the Fourier | |
89 | // transform representations above). | |
90 | // | |
91 | // Finally, the twelve functions above are each available with a third Policy | |
92 | // argument, which can be used to define a custom epsilon value. These Policy | |
93 | // versions bring the total number of functions provided by jacobi_theta.hpp | |
94 | // to twenty-four. | |
95 | // | |
96 | // See: | |
97 | // https://mathworld.wolfram.com/JacobiThetaFunctions.html | |
98 | // https://dlmf.nist.gov/20 | |
99 | ||
100 | #ifndef BOOST_MATH_JACOBI_THETA_HPP | |
101 | #define BOOST_MATH_JACOBI_THETA_HPP | |
102 | ||
103 | #include <boost/math/tools/complex.hpp> | |
104 | #include <boost/math/tools/precision.hpp> | |
105 | #include <boost/math/tools/promotion.hpp> | |
106 | #include <boost/math/policies/error_handling.hpp> | |
107 | #include <boost/math/constants/constants.hpp> | |
108 | ||
109 | namespace boost{ namespace math{ | |
110 | ||
111 | // Simple functions - parameterized by q | |
112 | template <class T, class U> | |
113 | inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q); | |
114 | template <class T, class U> | |
115 | inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q); | |
116 | template <class T, class U> | |
117 | inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q); | |
118 | template <class T, class U> | |
119 | inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q); | |
120 | ||
121 | // Simple functions - parameterized by tau (assumed imaginary) | |
1e59de90 TL |
122 | // q = exp(i*PI*TAU) |
123 | // tau = -log(q)/PI | |
20effc67 TL |
124 | template <class T, class U> |
125 | inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau); | |
126 | template <class T, class U> | |
127 | inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau); | |
128 | template <class T, class U> | |
129 | inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau); | |
130 | template <class T, class U> | |
131 | inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau); | |
132 | ||
133 | // Minus one versions for small q / large tau | |
134 | template <class T, class U> | |
135 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q); | |
136 | template <class T, class U> | |
137 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q); | |
138 | template <class T, class U> | |
139 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau); | |
140 | template <class T, class U> | |
141 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau); | |
142 | ||
143 | // Policied versions - parameterized by q | |
144 | template <class T, class U, class Policy> | |
145 | inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy& pol); | |
146 | template <class T, class U, class Policy> | |
147 | inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy& pol); | |
148 | template <class T, class U, class Policy> | |
149 | inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy& pol); | |
150 | template <class T, class U, class Policy> | |
151 | inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy& pol); | |
152 | ||
153 | // Policied versions - parameterized by tau | |
154 | template <class T, class U, class Policy> | |
155 | inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy& pol); | |
156 | template <class T, class U, class Policy> | |
157 | inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy& pol); | |
158 | template <class T, class U, class Policy> | |
159 | inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy& pol); | |
160 | template <class T, class U, class Policy> | |
161 | inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy& pol); | |
162 | ||
163 | // Policied m1 functions | |
164 | template <class T, class U, class Policy> | |
165 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy& pol); | |
166 | template <class T, class U, class Policy> | |
167 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy& pol); | |
168 | template <class T, class U, class Policy> | |
169 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy& pol); | |
170 | template <class T, class U, class Policy> | |
171 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy& pol); | |
172 | ||
173 | // Compare the non-oscillating component of the delta to the previous delta. | |
174 | // Both are assumed to be non-negative. | |
175 | template <class RealType> | |
176 | inline bool | |
177 | _jacobi_theta_converged(RealType last_delta, RealType delta, RealType eps) { | |
178 | return delta == 0.0 || delta < eps*last_delta; | |
179 | } | |
180 | ||
181 | template <class RealType> | |
182 | inline RealType | |
183 | _jacobi_theta_sum(RealType tau, RealType z_n, RealType z_increment, RealType eps) { | |
184 | BOOST_MATH_STD_USING | |
185 | RealType delta = 0, partial_result = 0; | |
186 | RealType last_delta = 0; | |
187 | ||
188 | do { | |
189 | last_delta = delta; | |
190 | delta = exp(-tau*z_n*z_n/constants::pi<RealType>()); | |
191 | partial_result += delta; | |
192 | z_n += z_increment; | |
193 | } while (!_jacobi_theta_converged(last_delta, delta, eps)); | |
194 | ||
195 | return partial_result; | |
196 | } | |
197 | ||
198 | // The following _IMAGINARY theta functions assume imaginary z and are for | |
199 | // internal use only. They are designed to increase accuracy and reduce the | |
200 | // number of iterations required for convergence for large |q|. The z argument | |
201 | // is scaled by tau, and the summations are rewritten to be double-sided | |
202 | // following DLMF 20.13.4 and 20.13.5. The return values are scaled by | |
1e59de90 | 203 | // exp(-tau*z^2/Pi)/sqrt(tau). |
20effc67 | 204 | // |
1e59de90 | 205 | // These functions are triggered when tau < 1, i.e. |q| > exp(-Pi) = 0.043 |
20effc67 TL |
206 | // |
207 | // Note that jacobi_theta4 uses the imaginary version of jacobi_theta2 (and | |
208 | // vice-versa). jacobi_theta1 and jacobi_theta3 use the imaginary versions of | |
209 | // themselves, following DLMF 20.7.30 - 20.7.33. | |
210 | template <class RealType, class Policy> | |
211 | inline RealType | |
212 | _IMAGINARY_jacobi_theta1tau(RealType z, RealType tau, const Policy&) { | |
213 | BOOST_MATH_STD_USING | |
214 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
215 | RealType result = RealType(0); | |
216 | ||
217 | // n>=0 even | |
218 | result -= _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::two_pi<RealType>(), eps); | |
219 | // n>0 odd | |
220 | result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>() + constants::pi<RealType>()), constants::two_pi<RealType>(), eps); | |
221 | // n<0 odd | |
222 | result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps); | |
223 | // n<0 even | |
224 | result -= _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>() - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps); | |
225 | ||
226 | return result * sqrt(tau); | |
227 | } | |
228 | ||
229 | template <class RealType, class Policy> | |
230 | inline RealType | |
231 | _IMAGINARY_jacobi_theta2tau(RealType z, RealType tau, const Policy&) { | |
232 | BOOST_MATH_STD_USING | |
233 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
234 | RealType result = RealType(0); | |
235 | ||
236 | // n>=0 | |
237 | result += _jacobi_theta_sum(tau, RealType(z + constants::half_pi<RealType>()), constants::pi<RealType>(), eps); | |
238 | // n<0 | |
239 | result += _jacobi_theta_sum(tau, RealType(z - constants::half_pi<RealType>()), RealType (-constants::pi<RealType>()), eps); | |
240 | ||
241 | return result * sqrt(tau); | |
242 | } | |
243 | ||
244 | template <class RealType, class Policy> | |
245 | inline RealType | |
246 | _IMAGINARY_jacobi_theta3tau(RealType z, RealType tau, const Policy&) { | |
247 | BOOST_MATH_STD_USING | |
248 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
249 | RealType result = 0; | |
250 | ||
251 | // n=0 | |
252 | result += exp(-z*z*tau/constants::pi<RealType>()); | |
253 | // n>0 | |
254 | result += _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::pi<RealType>(), eps); | |
255 | // n<0 | |
256 | result += _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType(-constants::pi<RealType>()), eps); | |
257 | ||
258 | return result * sqrt(tau); | |
259 | } | |
260 | ||
261 | template <class RealType, class Policy> | |
262 | inline RealType | |
263 | _IMAGINARY_jacobi_theta4tau(RealType z, RealType tau, const Policy&) { | |
264 | BOOST_MATH_STD_USING | |
265 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
266 | RealType result = 0; | |
267 | ||
268 | // n = 0 | |
269 | result += exp(-z*z*tau/constants::pi<RealType>()); | |
270 | ||
271 | // n > 0 odd | |
272 | result -= _jacobi_theta_sum(tau, RealType(z + constants::pi<RealType>()), constants::two_pi<RealType>(), eps); | |
273 | // n < 0 odd | |
274 | result -= _jacobi_theta_sum(tau, RealType(z - constants::pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps); | |
275 | // n > 0 even | |
276 | result += _jacobi_theta_sum(tau, RealType(z + constants::two_pi<RealType>()), constants::two_pi<RealType>(), eps); | |
277 | // n < 0 even | |
278 | result += _jacobi_theta_sum(tau, RealType(z - constants::two_pi<RealType>()), RealType (-constants::two_pi<RealType>()), eps); | |
279 | ||
280 | return result * sqrt(tau); | |
281 | } | |
282 | ||
283 | // First Jacobi theta function (Parameterized by tau - assumed imaginary) | |
1e59de90 | 284 | // = 2 * SUM (-1)^n * exp(i*Pi*Tau*(n+1/2)^2) * sin((2n+1)z) |
20effc67 TL |
285 | template <class RealType, class Policy> |
286 | inline RealType | |
287 | jacobi_theta1tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) | |
288 | { | |
289 | BOOST_MATH_STD_USING | |
290 | unsigned n = 0; | |
291 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
292 | RealType q_n = 0, last_q_n, delta, result = 0; | |
293 | ||
294 | if (tau <= 0.0) | |
295 | return policies::raise_domain_error<RealType>(function, | |
296 | "tau must be greater than 0 but got %1%.", tau, pol); | |
297 | ||
298 | if (abs(z) == 0.0) | |
299 | return result; | |
300 | ||
301 | if (tau < 1.0) { | |
302 | z = fmod(z, constants::two_pi<RealType>()); | |
303 | while (z > constants::pi<RealType>()) { | |
304 | z -= constants::two_pi<RealType>(); | |
305 | } | |
306 | while (z < -constants::pi<RealType>()) { | |
307 | z += constants::two_pi<RealType>(); | |
308 | } | |
309 | ||
310 | return _IMAGINARY_jacobi_theta1tau(z, RealType(1/tau), pol); | |
311 | } | |
312 | ||
313 | do { | |
314 | last_q_n = q_n; | |
315 | q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5) ); | |
316 | delta = q_n * sin(RealType(2*n+1)*z); | |
317 | if (n%2) | |
318 | delta = -delta; | |
319 | ||
320 | result += delta + delta; | |
321 | n++; | |
322 | } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); | |
323 | ||
324 | return result; | |
325 | } | |
326 | ||
327 | // First Jacobi theta function (Parameterized by q) | |
1e59de90 | 328 | // = 2 * SUM (-1)^n * q^(n+1/2)^2 * sin((2n+1)z) |
20effc67 TL |
329 | template <class RealType, class Policy> |
330 | inline RealType | |
331 | jacobi_theta1_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
332 | BOOST_MATH_STD_USING | |
333 | if (q <= 0.0 || q >= 1.0) { | |
334 | return policies::raise_domain_error<RealType>(function, | |
335 | "q must be greater than 0 and less than 1 but got %1%.", q, pol); | |
336 | } | |
337 | return jacobi_theta1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function); | |
338 | } | |
339 | ||
340 | // Second Jacobi theta function (Parameterized by tau - assumed imaginary) | |
1e59de90 | 341 | // = 2 * SUM exp(i*Pi*Tau*(n+1/2)^2) * cos((2n+1)z) |
20effc67 TL |
342 | template <class RealType, class Policy> |
343 | inline RealType | |
344 | jacobi_theta2tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) | |
345 | { | |
346 | BOOST_MATH_STD_USING | |
347 | unsigned n = 0; | |
348 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
349 | RealType q_n = 0, last_q_n, delta, result = 0; | |
350 | ||
351 | if (tau <= 0.0) { | |
352 | return policies::raise_domain_error<RealType>(function, | |
353 | "tau must be greater than 0 but got %1%.", tau, pol); | |
354 | } else if (tau < 1.0 && abs(z) == 0.0) { | |
355 | return jacobi_theta4tau(z, 1/tau, pol) / sqrt(tau); | |
356 | } else if (tau < 1.0) { // DLMF 20.7.31 | |
357 | z = fmod(z, constants::two_pi<RealType>()); | |
358 | while (z > constants::pi<RealType>()) { | |
359 | z -= constants::two_pi<RealType>(); | |
360 | } | |
361 | while (z < -constants::pi<RealType>()) { | |
362 | z += constants::two_pi<RealType>(); | |
363 | } | |
364 | ||
365 | return _IMAGINARY_jacobi_theta4tau(z, RealType(1/tau), pol); | |
366 | } | |
367 | ||
368 | do { | |
369 | last_q_n = q_n; | |
370 | q_n = exp(-tau * constants::pi<RealType>() * RealType(n + 0.5)*RealType(n + 0.5)); | |
371 | delta = q_n * cos(RealType(2*n+1)*z); | |
372 | result += delta + delta; | |
373 | n++; | |
374 | } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); | |
375 | ||
376 | return result; | |
377 | } | |
378 | ||
379 | // Second Jacobi theta function, parameterized by q | |
1e59de90 | 380 | // = 2 * SUM q^(n+1/2)^2 * cos((2n+1)z) |
20effc67 TL |
381 | template <class RealType, class Policy> |
382 | inline RealType | |
383 | jacobi_theta2_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
384 | BOOST_MATH_STD_USING | |
385 | if (q <= 0.0 || q >= 1.0) { | |
386 | return policies::raise_domain_error<RealType>(function, | |
387 | "q must be greater than 0 and less than 1 but got %1%.", q, pol); | |
388 | } | |
389 | return jacobi_theta2tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function); | |
390 | } | |
391 | ||
392 | // Third Jacobi theta function, minus one (Parameterized by tau - assumed imaginary) | |
1e59de90 | 393 | // This function preserves accuracy for small values of q (i.e. |q| < exp(-Pi) = 0.043) |
20effc67 | 394 | // For larger values of q, the minus one version usually won't help. |
1e59de90 | 395 | // = 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz) |
20effc67 TL |
396 | template <class RealType, class Policy> |
397 | inline RealType | |
398 | jacobi_theta3m1tau_imp(RealType z, RealType tau, const Policy& pol) | |
399 | { | |
400 | BOOST_MATH_STD_USING | |
401 | ||
402 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
403 | RealType q_n = 0, last_q_n, delta, result = 0; | |
404 | unsigned n = 1; | |
405 | ||
406 | if (tau < 1.0) | |
407 | return jacobi_theta3tau(z, tau, pol) - RealType(1); | |
408 | ||
409 | do { | |
410 | last_q_n = q_n; | |
411 | q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n)); | |
412 | delta = q_n * cos(RealType(2*n)*z); | |
413 | result += delta + delta; | |
414 | n++; | |
415 | } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); | |
416 | ||
417 | return result; | |
418 | } | |
419 | ||
420 | // Third Jacobi theta function, parameterized by tau | |
1e59de90 | 421 | // = 1 + 2 * SUM exp(i*Pi*Tau*(n)^2) * cos(2nz) |
20effc67 TL |
422 | template <class RealType, class Policy> |
423 | inline RealType | |
424 | jacobi_theta3tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) | |
425 | { | |
426 | BOOST_MATH_STD_USING | |
427 | if (tau <= 0.0) { | |
428 | return policies::raise_domain_error<RealType>(function, | |
429 | "tau must be greater than 0 but got %1%.", tau, pol); | |
430 | } else if (tau < 1.0 && abs(z) == 0.0) { | |
431 | return jacobi_theta3tau(z, RealType(1/tau), pol) / sqrt(tau); | |
432 | } else if (tau < 1.0) { // DLMF 20.7.32 | |
433 | z = fmod(z, constants::pi<RealType>()); | |
434 | while (z > constants::half_pi<RealType>()) { | |
435 | z -= constants::pi<RealType>(); | |
436 | } | |
437 | while (z < -constants::half_pi<RealType>()) { | |
438 | z += constants::pi<RealType>(); | |
439 | } | |
440 | return _IMAGINARY_jacobi_theta3tau(z, RealType(1/tau), pol); | |
441 | } | |
442 | return RealType(1) + jacobi_theta3m1tau_imp(z, tau, pol); | |
443 | } | |
444 | ||
445 | // Third Jacobi theta function, minus one (parameterized by q) | |
1e59de90 | 446 | // = 2 * SUM q^n^2 * cos(2nz) |
20effc67 TL |
447 | template <class RealType, class Policy> |
448 | inline RealType | |
449 | jacobi_theta3m1_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
450 | BOOST_MATH_STD_USING | |
451 | if (q <= 0.0 || q >= 1.0) { | |
452 | return policies::raise_domain_error<RealType>(function, | |
453 | "q must be greater than 0 and less than 1 but got %1%.", q, pol); | |
454 | } | |
455 | return jacobi_theta3m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol); | |
456 | } | |
457 | ||
458 | // Third Jacobi theta function (parameterized by q) | |
1e59de90 | 459 | // = 1 + 2 * SUM q^n^2 * cos(2nz) |
20effc67 TL |
460 | template <class RealType, class Policy> |
461 | inline RealType | |
462 | jacobi_theta3_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
463 | BOOST_MATH_STD_USING | |
464 | if (q <= 0.0 || q >= 1.0) { | |
465 | return policies::raise_domain_error<RealType>(function, | |
466 | "q must be greater than 0 and less than 1 but got %1%.", q, pol); | |
467 | } | |
468 | return jacobi_theta3tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol, function); | |
469 | } | |
470 | ||
471 | // Fourth Jacobi theta function, minus one (Parameterized by tau) | |
472 | // This function preserves accuracy for small values of q (i.e. tau > 1) | |
1e59de90 | 473 | // = 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz) |
20effc67 TL |
474 | template <class RealType, class Policy> |
475 | inline RealType | |
476 | jacobi_theta4m1tau_imp(RealType z, RealType tau, const Policy& pol) | |
477 | { | |
478 | BOOST_MATH_STD_USING | |
479 | ||
480 | RealType eps = policies::get_epsilon<RealType, Policy>(); | |
481 | RealType q_n = 0, last_q_n, delta, result = 0; | |
482 | unsigned n = 1; | |
483 | ||
484 | if (tau < 1.0) | |
485 | return jacobi_theta4tau(z, tau, pol) - RealType(1); | |
486 | ||
487 | do { | |
488 | last_q_n = q_n; | |
489 | q_n = exp(-tau * constants::pi<RealType>() * RealType(n)*RealType(n)); | |
490 | delta = q_n * cos(RealType(2*n)*z); | |
491 | if (n%2) | |
492 | delta = -delta; | |
493 | ||
494 | result += delta + delta; | |
495 | n++; | |
496 | } while (!_jacobi_theta_converged(last_q_n, q_n, eps)); | |
497 | ||
498 | return result; | |
499 | } | |
500 | ||
501 | // Fourth Jacobi theta function (Parameterized by tau) | |
1e59de90 | 502 | // = 1 + 2 * SUM (-1)^n exp(i*Pi*Tau*(n)^2) * cos(2nz) |
20effc67 TL |
503 | template <class RealType, class Policy> |
504 | inline RealType | |
505 | jacobi_theta4tau_imp(RealType z, RealType tau, const Policy& pol, const char *function) | |
506 | { | |
507 | BOOST_MATH_STD_USING | |
508 | if (tau <= 0.0) { | |
509 | return policies::raise_domain_error<RealType>(function, | |
510 | "tau must be greater than 0 but got %1%.", tau, pol); | |
511 | } else if (tau < 1.0 && abs(z) == 0.0) { | |
512 | return jacobi_theta2tau(z, 1/tau, pol) / sqrt(tau); | |
513 | } else if (tau < 1.0) { // DLMF 20.7.33 | |
514 | z = fmod(z, constants::pi<RealType>()); | |
515 | while (z > constants::half_pi<RealType>()) { | |
516 | z -= constants::pi<RealType>(); | |
517 | } | |
518 | while (z < -constants::half_pi<RealType>()) { | |
519 | z += constants::pi<RealType>(); | |
520 | } | |
521 | return _IMAGINARY_jacobi_theta2tau(z, RealType(1/tau), pol); | |
522 | } | |
523 | ||
524 | return RealType(1) + jacobi_theta4m1tau_imp(z, tau, pol); | |
525 | } | |
526 | ||
527 | // Fourth Jacobi theta function, minus one (Parameterized by q) | |
528 | // This function preserves accuracy for small values of q | |
1e59de90 | 529 | // = 2 * SUM q^n^2 * cos(2nz) |
20effc67 TL |
530 | template <class RealType, class Policy> |
531 | inline RealType | |
532 | jacobi_theta4m1_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
533 | BOOST_MATH_STD_USING | |
534 | if (q <= 0.0 || q >= 1.0) { | |
535 | return policies::raise_domain_error<RealType>(function, | |
536 | "q must be greater than 0 and less than 1 but got %1%.", q, pol); | |
537 | } | |
538 | return jacobi_theta4m1tau_imp(z, RealType (-log(q)/constants::pi<RealType>()), pol); | |
539 | } | |
540 | ||
541 | // Fourth Jacobi theta function, parameterized by q | |
1e59de90 | 542 | // = 1 + 2 * SUM q^n^2 * cos(2nz) |
20effc67 TL |
543 | template <class RealType, class Policy> |
544 | inline RealType | |
545 | jacobi_theta4_imp(RealType z, RealType q, const Policy& pol, const char *function) { | |
546 | BOOST_MATH_STD_USING | |
547 | if (q <= 0.0 || q >= 1.0) { | |
548 | return policies::raise_domain_error<RealType>(function, | |
549 | "|q| must be greater than zero and less than 1, but got %1%.", q, pol); | |
550 | } | |
551 | return jacobi_theta4tau_imp(z, RealType(-log(q)/constants::pi<RealType>()), pol, function); | |
552 | } | |
553 | ||
554 | // Begin public API | |
555 | ||
556 | template <class T, class U, class Policy> | |
557 | inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau, const Policy&) { | |
558 | BOOST_FPU_EXCEPTION_GUARD | |
559 | typedef typename tools::promote_args<T, U>::type result_type; | |
560 | typedef typename policies::normalise< | |
561 | Policy, | |
562 | policies::promote_float<false>, | |
563 | policies::promote_double<false>, | |
564 | policies::discrete_quantile<>, | |
565 | policies::assert_undefined<> >::type forwarding_policy; | |
566 | ||
567 | static const char* function = "boost::math::jacobi_theta1tau<%1%>(%1%)"; | |
568 | ||
569 | return policies::checked_narrowing_cast<result_type, Policy>( | |
570 | jacobi_theta1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
571 | forwarding_policy(), function), function); | |
572 | } | |
573 | ||
574 | template <class T, class U> | |
575 | inline typename tools::promote_args<T, U>::type jacobi_theta1tau(T z, U tau) { | |
576 | return jacobi_theta1tau(z, tau, policies::policy<>()); | |
577 | } | |
578 | ||
579 | template <class T, class U, class Policy> | |
580 | inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q, const Policy&) { | |
581 | BOOST_FPU_EXCEPTION_GUARD | |
582 | typedef typename tools::promote_args<T, U>::type result_type; | |
583 | typedef typename policies::normalise< | |
584 | Policy, | |
585 | policies::promote_float<false>, | |
586 | policies::promote_double<false>, | |
587 | policies::discrete_quantile<>, | |
588 | policies::assert_undefined<> >::type forwarding_policy; | |
589 | ||
590 | static const char* function = "boost::math::jacobi_theta1<%1%>(%1%)"; | |
591 | ||
592 | return policies::checked_narrowing_cast<result_type, Policy>( | |
593 | jacobi_theta1_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
594 | forwarding_policy(), function), function); | |
595 | } | |
596 | ||
597 | template <class T, class U> | |
598 | inline typename tools::promote_args<T, U>::type jacobi_theta1(T z, U q) { | |
599 | return jacobi_theta1(z, q, policies::policy<>()); | |
600 | } | |
601 | ||
602 | template <class T, class U, class Policy> | |
603 | inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau, const Policy&) { | |
604 | BOOST_FPU_EXCEPTION_GUARD | |
605 | typedef typename tools::promote_args<T, U>::type result_type; | |
606 | typedef typename policies::normalise< | |
607 | Policy, | |
608 | policies::promote_float<false>, | |
609 | policies::promote_double<false>, | |
610 | policies::discrete_quantile<>, | |
611 | policies::assert_undefined<> >::type forwarding_policy; | |
612 | ||
613 | static const char* function = "boost::math::jacobi_theta2tau<%1%>(%1%)"; | |
614 | ||
615 | return policies::checked_narrowing_cast<result_type, Policy>( | |
616 | jacobi_theta2tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
617 | forwarding_policy(), function), function); | |
618 | } | |
619 | ||
620 | template <class T, class U> | |
621 | inline typename tools::promote_args<T, U>::type jacobi_theta2tau(T z, U tau) { | |
622 | return jacobi_theta2tau(z, tau, policies::policy<>()); | |
623 | } | |
624 | ||
625 | template <class T, class U, class Policy> | |
626 | inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q, const Policy&) { | |
627 | BOOST_FPU_EXCEPTION_GUARD | |
628 | typedef typename tools::promote_args<T, U>::type result_type; | |
629 | typedef typename policies::normalise< | |
630 | Policy, | |
631 | policies::promote_float<false>, | |
632 | policies::promote_double<false>, | |
633 | policies::discrete_quantile<>, | |
634 | policies::assert_undefined<> >::type forwarding_policy; | |
635 | ||
636 | static const char* function = "boost::math::jacobi_theta2<%1%>(%1%)"; | |
637 | ||
638 | return policies::checked_narrowing_cast<result_type, Policy>( | |
639 | jacobi_theta2_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
640 | forwarding_policy(), function), function); | |
641 | } | |
642 | ||
643 | template <class T, class U> | |
644 | inline typename tools::promote_args<T, U>::type jacobi_theta2(T z, U q) { | |
645 | return jacobi_theta2(z, q, policies::policy<>()); | |
646 | } | |
647 | ||
648 | template <class T, class U, class Policy> | |
649 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau, const Policy&) { | |
650 | BOOST_FPU_EXCEPTION_GUARD | |
651 | typedef typename tools::promote_args<T, U>::type result_type; | |
652 | typedef typename policies::normalise< | |
653 | Policy, | |
654 | policies::promote_float<false>, | |
655 | policies::promote_double<false>, | |
656 | policies::discrete_quantile<>, | |
657 | policies::assert_undefined<> >::type forwarding_policy; | |
658 | ||
659 | static const char* function = "boost::math::jacobi_theta3m1tau<%1%>(%1%)"; | |
660 | ||
661 | return policies::checked_narrowing_cast<result_type, Policy>( | |
662 | jacobi_theta3m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
663 | forwarding_policy()), function); | |
664 | } | |
665 | ||
666 | template <class T, class U> | |
667 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1tau(T z, U tau) { | |
668 | return jacobi_theta3m1tau(z, tau, policies::policy<>()); | |
669 | } | |
670 | ||
671 | template <class T, class U, class Policy> | |
672 | inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau, const Policy&) { | |
673 | BOOST_FPU_EXCEPTION_GUARD | |
674 | typedef typename tools::promote_args<T, U>::type result_type; | |
675 | typedef typename policies::normalise< | |
676 | Policy, | |
677 | policies::promote_float<false>, | |
678 | policies::promote_double<false>, | |
679 | policies::discrete_quantile<>, | |
680 | policies::assert_undefined<> >::type forwarding_policy; | |
681 | ||
682 | static const char* function = "boost::math::jacobi_theta3tau<%1%>(%1%)"; | |
683 | ||
684 | return policies::checked_narrowing_cast<result_type, Policy>( | |
685 | jacobi_theta3tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
686 | forwarding_policy(), function), function); | |
687 | } | |
688 | ||
689 | template <class T, class U> | |
690 | inline typename tools::promote_args<T, U>::type jacobi_theta3tau(T z, U tau) { | |
691 | return jacobi_theta3tau(z, tau, policies::policy<>()); | |
692 | } | |
693 | ||
694 | ||
695 | template <class T, class U, class Policy> | |
696 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q, const Policy&) { | |
697 | BOOST_FPU_EXCEPTION_GUARD | |
698 | typedef typename tools::promote_args<T, U>::type result_type; | |
699 | typedef typename policies::normalise< | |
700 | Policy, | |
701 | policies::promote_float<false>, | |
702 | policies::promote_double<false>, | |
703 | policies::discrete_quantile<>, | |
704 | policies::assert_undefined<> >::type forwarding_policy; | |
705 | ||
706 | static const char* function = "boost::math::jacobi_theta3m1<%1%>(%1%)"; | |
707 | ||
708 | return policies::checked_narrowing_cast<result_type, Policy>( | |
709 | jacobi_theta3m1_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
710 | forwarding_policy(), function), function); | |
711 | } | |
712 | ||
713 | template <class T, class U> | |
714 | inline typename tools::promote_args<T, U>::type jacobi_theta3m1(T z, U q) { | |
715 | return jacobi_theta3m1(z, q, policies::policy<>()); | |
716 | } | |
717 | ||
718 | template <class T, class U, class Policy> | |
719 | inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q, const Policy&) { | |
720 | BOOST_FPU_EXCEPTION_GUARD | |
721 | typedef typename tools::promote_args<T, U>::type result_type; | |
722 | typedef typename policies::normalise< | |
723 | Policy, | |
724 | policies::promote_float<false>, | |
725 | policies::promote_double<false>, | |
726 | policies::discrete_quantile<>, | |
727 | policies::assert_undefined<> >::type forwarding_policy; | |
728 | ||
729 | static const char* function = "boost::math::jacobi_theta3<%1%>(%1%)"; | |
730 | ||
731 | return policies::checked_narrowing_cast<result_type, Policy>( | |
732 | jacobi_theta3_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
733 | forwarding_policy(), function), function); | |
734 | } | |
735 | ||
736 | template <class T, class U> | |
737 | inline typename tools::promote_args<T, U>::type jacobi_theta3(T z, U q) { | |
738 | return jacobi_theta3(z, q, policies::policy<>()); | |
739 | } | |
740 | ||
741 | template <class T, class U, class Policy> | |
742 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau, const Policy&) { | |
743 | BOOST_FPU_EXCEPTION_GUARD | |
744 | typedef typename tools::promote_args<T, U>::type result_type; | |
745 | typedef typename policies::normalise< | |
746 | Policy, | |
747 | policies::promote_float<false>, | |
748 | policies::promote_double<false>, | |
749 | policies::discrete_quantile<>, | |
750 | policies::assert_undefined<> >::type forwarding_policy; | |
751 | ||
752 | static const char* function = "boost::math::jacobi_theta4m1tau<%1%>(%1%)"; | |
753 | ||
754 | return policies::checked_narrowing_cast<result_type, Policy>( | |
755 | jacobi_theta4m1tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
756 | forwarding_policy()), function); | |
757 | } | |
758 | ||
759 | template <class T, class U> | |
760 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1tau(T z, U tau) { | |
761 | return jacobi_theta4m1tau(z, tau, policies::policy<>()); | |
762 | } | |
763 | ||
764 | template <class T, class U, class Policy> | |
765 | inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau, const Policy&) { | |
766 | BOOST_FPU_EXCEPTION_GUARD | |
767 | typedef typename tools::promote_args<T, U>::type result_type; | |
768 | typedef typename policies::normalise< | |
769 | Policy, | |
770 | policies::promote_float<false>, | |
771 | policies::promote_double<false>, | |
772 | policies::discrete_quantile<>, | |
773 | policies::assert_undefined<> >::type forwarding_policy; | |
774 | ||
775 | static const char* function = "boost::math::jacobi_theta4tau<%1%>(%1%)"; | |
776 | ||
777 | return policies::checked_narrowing_cast<result_type, Policy>( | |
778 | jacobi_theta4tau_imp(static_cast<result_type>(z), static_cast<result_type>(tau), | |
779 | forwarding_policy(), function), function); | |
780 | } | |
781 | ||
782 | template <class T, class U> | |
783 | inline typename tools::promote_args<T, U>::type jacobi_theta4tau(T z, U tau) { | |
784 | return jacobi_theta4tau(z, tau, policies::policy<>()); | |
785 | } | |
786 | ||
787 | template <class T, class U, class Policy> | |
788 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q, const Policy&) { | |
789 | BOOST_FPU_EXCEPTION_GUARD | |
790 | typedef typename tools::promote_args<T, U>::type result_type; | |
791 | typedef typename policies::normalise< | |
792 | Policy, | |
793 | policies::promote_float<false>, | |
794 | policies::promote_double<false>, | |
795 | policies::discrete_quantile<>, | |
796 | policies::assert_undefined<> >::type forwarding_policy; | |
797 | ||
798 | static const char* function = "boost::math::jacobi_theta4m1<%1%>(%1%)"; | |
799 | ||
800 | return policies::checked_narrowing_cast<result_type, Policy>( | |
801 | jacobi_theta4m1_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
802 | forwarding_policy(), function), function); | |
803 | } | |
804 | ||
805 | template <class T, class U> | |
806 | inline typename tools::promote_args<T, U>::type jacobi_theta4m1(T z, U q) { | |
807 | return jacobi_theta4m1(z, q, policies::policy<>()); | |
808 | } | |
809 | ||
810 | template <class T, class U, class Policy> | |
811 | inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q, const Policy&) { | |
812 | BOOST_FPU_EXCEPTION_GUARD | |
813 | typedef typename tools::promote_args<T, U>::type result_type; | |
814 | typedef typename policies::normalise< | |
815 | Policy, | |
816 | policies::promote_float<false>, | |
817 | policies::promote_double<false>, | |
818 | policies::discrete_quantile<>, | |
819 | policies::assert_undefined<> >::type forwarding_policy; | |
820 | ||
821 | static const char* function = "boost::math::jacobi_theta4<%1%>(%1%)"; | |
822 | ||
823 | return policies::checked_narrowing_cast<result_type, Policy>( | |
824 | jacobi_theta4_imp(static_cast<result_type>(z), static_cast<result_type>(q), | |
825 | forwarding_policy(), function), function); | |
826 | } | |
827 | ||
828 | template <class T, class U> | |
829 | inline typename tools::promote_args<T, U>::type jacobi_theta4(T z, U q) { | |
830 | return jacobi_theta4(z, q, policies::policy<>()); | |
831 | } | |
832 | ||
833 | }} | |
834 | ||
835 | #endif |