EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2026-06-02 08:17:09
0001 // 0002 /// @file lisk.hpp 0003 // LiSK 0004 // 0005 // Created by Sebastian on 08/03/16. 0006 // Copyright © 2016 Sebastian Kirchner. All rights reserved. 0007 // 0008 // This file is part of LiSK. 0009 // 0010 // LiSK is free software: you can redistribute it and/or modify 0011 // it under the terms of the GNU General Public License as published by 0012 // the Free Software Foundation, either version 3 of the License, or 0013 // (at your option) any later version. 0014 // 0015 // LiSK is distributed in the hope that it will be useful, 0016 // but WITHOUT ANY WARRANTY; without even the implied warranty of 0017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0018 // GNU General Public License for more details. 0019 // 0020 // You should have received a copy of the GNU General Public License 0021 // along with LiSK. If not, see <http://www.gnu.org/licenses/>. 0022 0023 #ifndef lisk_hpp 0024 #define lisk_hpp 0025 0026 #include <iostream> 0027 #include <complex> 0028 #include <vector> 0029 #include <string> 0030 #include <unordered_map> 0031 #include <cln/cln.h> 0032 0033 /** LiSK namespace */ 0034 namespace LiSK { 0035 0036 /** 0037 * Type defs 0038 * @note Define subsequently used abbreviations 0039 */ 0040 typedef std::complex<double> cp_d; 0041 0042 /** 0043 * LiSK - Library for Evaluating Classical Polylogarithms and Li_22 0044 * @note A lightweight library for evaluating classical polylogarithms 0045 * Li_n and the special function Li_{22} based on arXiv:1601.02649. 0046 * Currently a double precision as well as an arbitrary precision 0047 * implementation is available; the latter through the CLN library by 0048 * Bruno Haible ( http://www.ginac.de/CLN/ ). 0049 * CLN is an essential part of this library and is also used in the 0050 * initialisation phase of the double prescision implementation. 0051 * @note Although this library is header-only the actual implementation 0052 * is stored in lisk.cpp which is included at the end of lisk.hpp. 0053 * 0054 */ 0055 template <typename numtype> 0056 class LiSK { 0057 0058 private: 0059 0060 // The user can change the setup as follows: 0061 // 0062 // 1.) Set the step size _step for extending the vectors of constants. 0063 // A small value reduces the initialisation time but requires the 0064 // initialisation more often. Larger values vice versa. (default = 5) 0065 // 0066 // 2.) Vary _offset to determine the order of magnitude of the 'small' 0067 // imaginary part, i.e. iep = i*10^(offset-prec). With prec=17 0068 // for double approximation. Otherwise prec is set during 0069 // construction of the LiSK object. (default = 2) 0070 // 0071 // 3.) Enable/disable the conversion to floating point numbers of the 0072 // stored constants, i.e. the floating point approximation of rational 0073 // numbers are stored in the look-up tables instead of the exact value. 0074 // IMPORTANT: Although enabling (true) the conversion leads to 0075 // performance improvements, it does not work well for all arguments 0076 // in the complex plane. Especially real input parameters seem to pose 0077 // a sincere problem. Handle with care! (default = false) 0078 // 0079 // 4.) Influence LiSKs initialisation phase by raising the number of pre- 0080 // calculated constants in _user_max_values. Higher values lead to 0081 // a longer initialisation time but less re-computing during the main 0082 // run. Default minimum values are set during the construction nevertheless. 0083 0084 /** Step size 0085 * @note Set step size used during computation of constants 0086 * in case an extension is required. 0087 */ 0088 const unsigned _step = 5; 0089 0090 /** Small positive imaginary part 0091 * @note Small positive imaginary part which will be assigned as x->x-iep. 0092 * For the double precision implementation we set _iep_d= i*10^(offset-17) 0093 * and in the arbitrary precision case as _iep_clN = i*10^(offset-prec), 0094 * with prec set in constructor. The offset can be variied, but the user 0095 * has to know what he's doing (default = 2). 0096 * @{ 0097 */ 0098 const int _offset = 2; 0099 const cp_d _iep_d; 0100 const cln::cl_N _iep_clN; 0101 /** @} */ 0102 0103 /** Enable/disable floating point conversion of constants 0104 * @note If this flag is set true all constants in the look-up 0105 * tables are given by their floating point approximation, 0106 * respecting the required precision. However, problems can occur 0107 * when using real input variables. Be careful! (default = false) 0108 */ 0109 const bool _cln_tofloat = false; 0110 0111 /** User defined upper limits for various constants. 0112 * @note The limits, which determine how many constants are pre-calculated 0113 * can be set here. The keys should not be changed. Note: A minimum will be 0114 * assumed during construction and the user choice may be overwritten. But 0115 * higher values will be respected. 0116 */ 0117 const std::unordered_map<std::string, size_t> _user_max_values { 0118 {"nHarmNum",3}, // Harmonic numbers H(n) 0119 {"nPosZeta",4}, // Zeta(n) for n>=0 0120 {"nNegZeta",100}, // Zeta(n) for n<0 0121 {"nFactorials",100}, // Factorials n! 0122 {"nBernNum",100}, // Bernoulli numbers B(n) 0123 {"nLiCij",200}, // Cij/(j+1)! from eq.(5.10) 0124 {"nLiBn_eq59",100} // B_{2n}/(2n*(2n+m)!) from eq.(5.9) 0125 }; 0126 0127 0128 // ### Do NOT change anything from here ### 0129 0130 /** Precision value 0131 * @note Only relevant in arbitrary precision case. 0132 * Minimum is double precision with _prec=17 0133 */ 0134 const size_t _prec; 0135 0136 /** ErrorEstimate 0137 * @note Try to estimate how far we have to expand to get the required 0138 * precision. Use practical approach as also used in GiNaC, i.e. 0139 * check if new result differs from previous one. Due to included 0140 * zeros in the constants we need to check for them. 0141 * @param res Latest computed result 0142 * @param prev_res In previous step computed result 0143 * @param current Check current if zero 0144 * @return True, if res!=prev_res || tmp==0. False otherwise 0145 */ 0146 bool _ErrorEstimate(const numtype& res, const numtype& prev_res, const numtype& current) const; 0147 0148 /** Constants 0149 * @note Collection of constants which are precomputed during 0150 * initialisation and stored. If additional constants are 0151 * required during runtime the existing set will be extended. 0152 * Initial maximum values can be set in _init. 0153 * @{ 0154 */ 0155 template <typename T> 0156 struct constants{ 0157 /** Harmonic numbers H(n) for n=1,2,3,... */ 0158 std::vector<T> HarmNum; 0159 /** zeta(n) for n=0,(1),2,3,4. n=1 is kept but set to 0 */ 0160 std::vector<T> PosZeta; 0161 /** zeta(n) for n<0 and n odd. All negative even n give zero */ 0162 std::vector<T> NegZeta; 0163 /** Factorials n! (n>170 exceeds double precision). */ 0164 std::vector<T> Factorials; 0165 /** Bernoulli numbers of first kind, i.e. B_1 = -1/2. The vector is 0166 * {B_0,B_1,B_2n} with n=1,2,3,...*/ 0167 std::vector<T> BernNum; 0168 /** Li constants Cij/(j+1)! from eq.(5.10) */ 0169 std::vector<std::vector<T>> LiCij; 0170 /** Li constants B_{2n}/(2n*(2n+m)!) from eq.(5.9) */ 0171 std::vector<std::vector<T>> LiBn_eq59; 0172 /** Max values for available constants. The keys are named following the 0173 * convention "n"+"vector name", e.g. nHarmNum. In addition "nweight" is set 0174 * and determines the highest weigth for which the constants are computed. 0175 */ 0176 std::unordered_map<std::string, size_t> max_values; 0177 0178 /** C'stor 0179 * @note Pre-compute constants needed for the evaluation of Li_n and Li_22. 0180 * @param Li_n Set expected weight of classical polylogarithms 0181 */ 0182 constants(const size_t Li_n); 0183 }; 0184 constants<numtype> _constants; 0185 constants<cln::cl_R> _constants_clR; // Real numbers because this respects the 0186 // different nature of rationals and floats 0187 /** @} */ 0188 0189 /** Initialise constants 0190 * @param max_values Map setting the maximum for each constant indicated by its key value. See constants::max_value 0191 */ 0192 void _init(const std::unordered_map<std::string, size_t> max_values); 0193 0194 /** Convert CLN to desired type 0195 * @param x Value to be converted 0196 * @return x in desired type 0197 */ 0198 numtype convert(const cln::cl_R& x) const; 0199 0200 /** Maximum of vector entries 0201 * @param v Vector from which maximal value is determined, e.g. {x1,x2} 0202 * @return Maximal value of v, e.g. x1 if x1>x2 0203 */ 0204 template<typename T> 0205 T _vec_max(std::vector<T> v) const; 0206 0207 /** Extend Harmonic numbers 0208 * @note Extends HarmNum vector up to nHarmNum, starting with H_0=0. 0209 */ 0210 void _HarmNum_extend(); 0211 0212 /** Extend positve Zeta values 0213 * @note Extends PosZeta vector with zeta(n), n=0,1,2,... and zeta(1)=0. 0214 */ 0215 void _PosZeta_extend(); 0216 0217 /** Extend factorials 0218 * @note Extends Factorials vector n! with n=0,1,2,... 0219 */ 0220 void _Factorials_extend(); 0221 0222 /** Extend Bernoulli numbers 0223 * @note Extends BernNum vector B_n with n=0,1,2,3,... . B_n=0 for odd n>1. 0224 * Zeros are kept for now. If that poses a performance or memory problem 0225 * later on it could be changed. 0226 */ 0227 void _BernNum_extend(); 0228 0229 /** Negative zeta values 0230 * @note Extends NegZeta vector with zeta(-n), n=1,2,3.... zeta(-n)=0 for even n>0. 0231 * Zeros are kept for now. 0232 */ 0233 void _NegZeta_extend(); 0234 0235 /** Extend Li constants LiCij 0236 * @note Compute the coefficients tCij = LiCij/(j+1)! from eq.(5.10). Extends 0237 * the multi-dim vector LiCij containing tCij with i=1,2,3,4 and j=0,1,2,... 0238 */ 0239 void _LiCij_extend(); 0240 0241 /** Extend Li constants LiBn_eq59 0242 * @note Compute the coefficients tBnm = B_{2n}/(2n*(2n+m)!) for n>0 and m=1,2,3 0243 */ 0244 void _LiBn_eq59_extend(); 0245 0246 /** Li_n(1-exp(-alpha)) 0247 * @note Implementation of eq.(5.10) up to desired accuracy. 0248 * @param n Weight of Li_n 0249 * @param x Point of evaluation, i.e. x=1-exp(-alpha) 0250 * @return Li_n(x) up to desired accuracy. 0251 */ 0252 numtype _Lin_1mexpal(const size_t n, const numtype& x); 0253 0254 /** Li sums 0255 * @note Compute Li_n(~al) = precal + sign * sum_{m=1} B_{2m}/(2m*(2m+n-1)!)*al^{2m+n-1} as needed 0256 * for the classical polylogarithms (eq.(5.9)). 0257 * @param precal Terms to be added to sum 0258 * @param al Accordingly mapped argument of Li_n(x). See eq.(5.9) and eq.(5.10) 0259 * @param n Weight specifier Li_n 0260 * @param fac_2n If true, each summand will be multiplied by 2m 0261 * @param sign If true, the entire sum is multiplied by -1 0262 * @return Li_n(~al) 0263 */ 0264 numtype _Li_sumBn(const numtype &precal,const numtype &al, const size_t n, const bool fac_2n=false, const bool sign=false); 0265 0266 /** General Li_n(x) 0267 * @param x Point of evaluation 0268 * @param n Weight of Li_n 0269 * @return Li_n(x) 0270 */ 0271 numtype _Lin_basis(const size_t n, const numtype& x); 0272 0273 /** Inversion formula for Li(n,x) 0274 * @note Inversion formula from eq.(5.6). Used for Li_n with n>4. 0275 * @param x Point of evaluation 0276 * @param n Weight of Li_n 0277 * @return Li_n(x) where |x|>1 0278 */ 0279 numtype _Lin_inverse(const size_t n, const numtype& x); 0280 0281 /** All order Li_n(exp(-alpha)) 0282 * @note Implementation of eq.(5.8) for all n>0. 0283 * @param n Weight specifier Li_n 0284 * @param x Point of evaluation, i.e. x=exp(-alpha) 0285 * @return Li_n(x) up to desired accuracy 0286 */ 0287 numtype _Lin_expal(const size_t n, const numtype& x); 0288 0289 /** Power function 0290 * @note Power function for exponents of integer type. 0291 * Nothing more is needed here. 0292 * @param x Basis 0293 * @param n Exponent 0294 * @return x^n 0295 */ 0296 numtype Pow(const numtype& x, const int n) const; 0297 0298 /** Real part 0299 * @param x Complex number 0300 * @return Re(x) 0301 * @{ 0302 */ 0303 inline 0304 cln::cl_R Re(const cln::cl_N& x) const{ 0305 return cln::realpart(x); 0306 }; 0307 0308 inline 0309 double Re(const std::complex<double>& x) const{ 0310 return x.real(); 0311 }; 0312 /** @} */ 0313 0314 /** Is_zero 0315 * @note Using CLNs zerop() to determine if number 0316 * is zero. 0317 * @param x Complex number 0318 * @return True, if x=0. False, otherwise 0319 */ 0320 bool is_zero(const numtype& x) const; 0321 0322 /** Add imaginary part 0323 * @note Add small imaginary part depending on precision. 0324 * We follow the widely used convention for mathematical 0325 * software, e.g. "The C Standard", Mathematica, GiNaC, etc. 0326 * @param x Complex value 0327 * @return x-iep 0328 */ 0329 numtype _add_iep(const numtype& x) const; 0330 0331 /** Adjust digits 0332 * @note Adjust digits of input variables such that they fullfil 0333 * desired precision. Only relevant for arbitrary precision part. 0334 * @param x Complex value 0335 * @return x adjusted to prec digits 0336 */ 0337 inline 0338 cln::cl_N _adjustDigits(const cln::cl_N& x) const; 0339 0340 /** Li(n,x) 0341 * @note Classical polylogarithm of weight n at point x but no 0342 * imaginary part is assigned 0343 * @param n Weight of classical polylogarithm 0344 * @param x Point of evaluation 0345 * @return Li(n,x) 0346 * @{ 0347 */ 0348 numtype _Li1(const numtype& x); 0349 numtype _Li2(const numtype& x); 0350 numtype _Li3(const numtype& x); 0351 numtype _Li4(const numtype& x); 0352 /** @} */ 0353 0354 /** Li22(x,y) 0355 * @note No imaginary parts will be assigned. 0356 * @param x Point of evaluation 0357 * @param y Point of evaluation 0358 * @return Li22(x,y) 0359 */ 0360 numtype _Li22(const numtype& x, const numtype& y); 0361 0362 /** Original Li22(x,y) 0363 * @note Original definition of Li22(x,y) from (6.1) if 0364 * |x|<1 and |xy|<1. 0365 * @param x Point of evaluation 0366 * @param y Point of evaluation 0367 * @return Li22(x,y) 0368 */ 0369 numtype _Li22_orig(const numtype& x, const numtype& y) const; 0370 0371 /** Stuffle relation for Li22(x,y) 0372 * @param x Point of evaluation 0373 * @param y Point of evaluation 0374 * @return -Li22(y,x) - Li4(x*y) + Li2(x)*Li2(y) 0375 */ 0376 inline 0377 numtype _Li22_stuffle(const numtype& x, const numtype& y); 0378 0379 /** Inversion relation for Li22(x,y) 0380 * @param x Point of evaluation 0381 * @param y Point of evaluation 0382 * @return Inversion relation from eq.(6.4) 0383 */ 0384 numtype _Li22_inversion(const numtype& x, const numtype& y); 0385 0386 public: 0387 0388 /** C'stor 0389 * @note During initialisation (construction) constants for polylogarithms up to 0390 * weight n are computed. Occurences of polylogarithms with higher weights during 0391 * runtime are possible but might lead to longer evaluation times due to adaption 0392 * of the constants. If weights n>4 are expected to occur n should be specified 0393 * accordingly to reduce evaluation time during the actual computation. 0394 * If CLN is used for arbitrary precision then the desired precision is set with prec. 0395 * @param n Prepare for polylogs of weight n 0396 * @param prec Set precision when using CLN 0397 */ 0398 LiSK(const size_t n=4, const size_t prec=34); 0399 0400 /** Li(n,x) 0401 * @note Classical polylogarithm of weight n at point x. Public wrapper 0402 * for the actual private implementation. A small imaginary part is assigned, 0403 * i.e. x->x+iep. For the n<=4 the abbreviations Lin(x)=Li(n,x) are supported. 0404 * @param n Weight of classical polylogarithm 0405 * @param x Point of evaluation 0406 * @return Li(n,x) 0407 * @{ 0408 */ 0409 numtype Li(const size_t n, const numtype x); 0410 numtype Li1(const numtype x); 0411 numtype Li2(const numtype x); 0412 numtype Li3(const numtype x); 0413 numtype Li4(const numtype x); 0414 /** @} */ 0415 0416 /** Li22(x,y) 0417 * @note Public wrapper for Li22(x,y). A small imaginary part is assigned to x 0418 * and y, i.e. x->x+iep and y->y+iep. 0419 * @param x Point of evaluation 0420 * @param y Point of evaluation 0421 * @return Li22(x,y) 0422 */ 0423 numtype Li22(const numtype x, const numtype y); 0424 0425 }; 0426 0427 /************************************************************************* 0428 * Implementation of LiSK * 0429 *************************************************************************/ 0430 #include "lisk.ipp" 0431 0432 } // End LiSK namespace 0433 0434 #endif /* lisk_hpp */
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|