diff --git a/GX-AnalyticContinuation/CMakeLists.txt b/GX-AnalyticContinuation/CMakeLists.txt index 9faa8837..4aadeead 100644 --- a/GX-AnalyticContinuation/CMakeLists.txt +++ b/GX-AnalyticContinuation/CMakeLists.txt @@ -1,3 +1,9 @@ +# *************************************************************************************************** +# Copyright (C) 2020-2024 GreenX library +# This file is distributed under the terms of the APACHE2 License. +# +# *************************************************************************************************** + # ----------------------------------------------- # CMakeLists for Analytic Continuation Library # ----------------------------------------------- diff --git a/GX-AnalyticContinuation/LICENSE-2.0.txt b/GX-AnalyticContinuation/LICENSE-2.0.txt new file mode 100644 index 00000000..d6456956 --- /dev/null +++ b/GX-AnalyticContinuation/LICENSE-2.0.txt @@ -0,0 +1,202 @@ + + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/GX-AnalyticContinuation/api/gx_ac.F90 b/GX-AnalyticContinuation/api/gx_ac.F90 index 07113113..1b806007 100644 --- a/GX-AnalyticContinuation/api/gx_ac.F90 +++ b/GX-AnalyticContinuation/api/gx_ac.F90 @@ -18,8 +18,8 @@ module gx_ac free_params, & arbitrary_precision_available - !> brief store the parameters of the tiehle pade model - !> (potentially in abitrary precision floats using GMP) + !> @brief store the parameters of the tiehle pade model + !!! (potentially in abitrary precision floats using GMP) type :: params logical :: initialized = .false. integer :: n_par @@ -43,7 +43,9 @@ module gx_ac #ifdef GMPXX_FOUND interface - !> auxiliary function to compute Thiele-Pade parameters using arbitrary precision numbers + !> @brief auxiliary function to compute Thiele-Pade parameters using + !! arbitrary precision numbers + !! !! @param[in] n_par - order of the interpolant !! @param[in] x_ref - array of the reference points !! @param[in] y_ref - array of the reference function values @@ -59,7 +61,9 @@ function thiele_pade_mp_aux(n_par, x_ref, y_ref, do_greedy, precision) bind(C, n type(c_ptr) :: thiele_pade_mp_aux end function thiele_pade_mp_aux - !> auxiliary function to evaluate the Thiele-Pade parameters using arbitrary precision numbers + !> @brief auxiliary function to evaluate the Thiele-Pade parameters using + !! arbitrary precision numbers + !! !! @param[in] x - point where the function needs to be evaluated !! @param[in] params_ptr - pointer to abstract type to store all parameters !! @return - interpolated function value @@ -70,7 +74,8 @@ function evaluate_thiele_pade_mp_aux(x, params_ptr) bind(C, name="evaluate_thiel complex(c_double_complex) :: evaluate_thiele_pade_mp_aux end function evaluate_thiele_pade_mp_aux - !> Frees the C++ pointers used in the Fortran type + !> @brief Frees the C++ pointers used in the Fortran type + !! !! @param[in] params_ptr - the C++ pointer subroutine free_pade_model(params_ptr) bind(C, name="free_pade_model") import :: c_ptr @@ -82,7 +87,9 @@ end subroutine free_pade_model contains - !> API function to compute Thiele-Pade approximations of a meromorphic function + !> @brief API function to compute Thiele-Pade approximations of a meromorphic + !! function + !! !! @param[in] n_par - order of the interpolant !! @param[in] x_ref - array of the reference points !! @param[in] y_ref - array of the reference function values @@ -116,8 +123,8 @@ end subroutine thiele_pade_api - !> API function to compute Thiele-Pade parameters - !! (potentially using arbitrary precision arithmetic) + !> @brief API function to compute Thiele-Pade parameters + !! (potentially using arbitrary precision arithmetic) !! !! @param[in] n_par - order of the interpolant !! @param[in] x - array of the reference points @@ -196,8 +203,8 @@ end function create_thiele_pade - !> API function to evaluate the Thiele-Pade approximation - !! (potentially using arbitrary precision numbers) + !> @brief API function to evaluate the Thiele-Pade approximation + !! (potentially using arbitrary precision numbers) !! !! @param[in] x - point where the function is evaluated !! @param[in] params - abstract type to store all parameters @@ -237,7 +244,7 @@ end function evaluate_thiele_pade_at - !> deallocate the pade model + !> @brief deallocate the pade model !! !! @param[inout] par - pade model struct subroutine free_params(par) diff --git a/GX-AnalyticContinuation/src/pade_approximant.f90 b/GX-AnalyticContinuation/src/pade_approximant.f90 index 7c458f1c..dacca819 100644 --- a/GX-AnalyticContinuation/src/pade_approximant.f90 +++ b/GX-AnalyticContinuation/src/pade_approximant.f90 @@ -1,5 +1,5 @@ ! *************************************************************************************************** -! Copyright (C) 2020-2023 GreenX library +! Copyright (C) 2020-2024 GreenX library ! This file is distributed under the terms of the APACHE2 License. ! ! *************************************************************************************************** diff --git a/GX-AnalyticContinuation/src/pade_mp.cpp b/GX-AnalyticContinuation/src/pade_mp.cpp index b89eb5db..978bca02 100644 --- a/GX-AnalyticContinuation/src/pade_mp.cpp +++ b/GX-AnalyticContinuation/src/pade_mp.cpp @@ -1,4 +1,11 @@ +// *************************************************************************************************** +// Copyright (C) 2020-2024 GreenX library +// This file is distributed under the terms of the APACHE2 License. +// +// *************************************************************************************************** + #include "pade_mp.h" + #include #include #include @@ -15,14 +22,14 @@ void thiele_pade_gcoeff_mp(const ComplexGMP *x, const std::vector &y, std::vector> &g_func, int n) { - g_func[n][0] = y[n]; - if (n == 0) - return; - - for (int idx = 1; idx <= n; idx++) { - g_func[n][idx] = (g_func[idx - 1][idx - 1] - g_func[n][idx - 1]) / - ((x[n] - x[idx - 1]) * g_func[n][idx - 1]); - } + g_func[n][0] = y[n]; + if (n == 0) + return; + + for (int idx = 1; idx <= n; idx++) { + g_func[n][idx] = (g_func[idx - 1][idx - 1] - g_func[n][idx - 1]) / + ((x[n] - x[idx - 1]) * g_func[n][idx - 1]); + } } /// @brief Gets the value of the Wmn matrices using the previously computed Pade @@ -38,20 +45,19 @@ void evaluate_thiele_pade_tab_mp(int n_par, const any_complex *x_ref, const any_complex &x, const any_complex *a_par, std::vector &acoef, std::vector &bcoef) { - - // Define constants and variables - const ComplexGMP c_one(std::complex(1.0, 0.0)); - any_complex delta; - - // Wallis' method iteration - delta = a_par[n_par + 1] * (x - x_ref[n_par]); - if (n_par == 0) { - acoef[1] = acoef[0]; - bcoef[1] = c_one + delta; - return; - } - acoef[n_par + 1] = acoef[n_par] + delta * acoef[n_par - 1]; - bcoef[n_par + 1] = bcoef[n_par] + delta * bcoef[n_par - 1]; + // Define constants and variables + const ComplexGMP c_one(std::complex(1.0, 0.0)); + any_complex delta; + + // Wallis' method iteration + delta = a_par[n_par + 1] * (x - x_ref[n_par]); + if (n_par == 0) { + acoef[1] = acoef[0]; + bcoef[1] = c_one + delta; + return; + } + acoef[n_par + 1] = acoef[n_par] + delta * acoef[n_par - 1]; + bcoef[n_par + 1] = bcoef[n_par] + delta * bcoef[n_par - 1]; } /// @brief Evaluates a Pade approximant constructed with Thiele's reciprocal @@ -61,34 +67,34 @@ void evaluate_thiele_pade_tab_mp(int n_par, const any_complex *x_ref, /// @return the value of the interpolant at x std::complex evaluate_thiele_pade_mp(const std::complex x, pade_model *params) { - mpf_set_default_prec(params->precision); - - // conversion from double to gmp - ComplexGMP x_mp(x); - - // Define constants - const ComplexGMP c_one(std::complex(1.0, 0.0)); - const mpf_class tol("1", params->precision); - - // Wallis' method coefficients - std::vector acoef(params->n_par + 1), bcoef(params->n_par + 1); - acoef[0] = params->a_par[0]; - bcoef[0] = c_one; - - // Compute continued fraction - for (int i_par = 0; i_par < params->n_par - 1; ++i_par) { - evaluate_thiele_pade_tab_mp(i_par, params->xref, x_mp, - params->a_par, acoef, bcoef); - if (bcoef[i_par + 1].abs() > tol) { - acoef[i_par + 1] = acoef[i_par + 1] / bcoef[i_par + 1]; - acoef[i_par] = acoef[i_par] / bcoef[i_par + 1]; - bcoef[i_par] = bcoef[i_par] / bcoef[i_par + 1]; - bcoef[i_par + 1] = c_one; + mpf_set_default_prec(params->precision); + + // conversion from double to gmp + ComplexGMP x_mp(x); + + // Define constants + const ComplexGMP c_one(std::complex(1.0, 0.0)); + const mpf_class tol("1", params->precision); + + // Wallis' method coefficients + std::vector acoef(params->n_par + 1), bcoef(params->n_par + 1); + acoef[0] = params->a_par[0]; + bcoef[0] = c_one; + + // Compute continued fraction + for (int i_par = 0; i_par < params->n_par - 1; ++i_par) { + evaluate_thiele_pade_tab_mp(i_par, params->xref, x_mp, + params->a_par, acoef, bcoef); + if (bcoef[i_par + 1].abs() > tol) { + acoef[i_par + 1] = acoef[i_par + 1] / bcoef[i_par + 1]; + acoef[i_par] = acoef[i_par] / bcoef[i_par + 1]; + bcoef[i_par] = bcoef[i_par] / bcoef[i_par + 1]; + bcoef[i_par + 1] = c_one; + } } - } - return std::complex(acoef[params->n_par - 1] / - bcoef[params->n_par - 1]); + return std::complex(acoef[params->n_par - 1] / + bcoef[params->n_par - 1]); } /// @brief Gets the Pade approximant of a meromorphic function F @@ -100,169 +106,169 @@ std::complex evaluate_thiele_pade_mp(const std::complex x, /// @return a pointer to the struct holding all model parameters pade_model *thiele_pade_mp(int n_par, const std::complex *x_ref, const std::complex *y_ref, int do_greedy, int precision) { - // set floating point precision of GMP lib - mpf_set_default_prec(precision); - - // Define constants - const ComplexGMP c_one(std::complex(1.0, 0.0)); - const ComplexGMP c_zero(std::complex(0.0, 0.0)); - const mpf_class tol("1e-6", precision); - - // auxiliary variables - int kdx, n_rem; - std::vector n_rem_idx; - mpf_class deltap, pval; - ComplexGMP pval_in, x_in, y_in, acoef_in, bcoef_in; - std::vector> g_func; - ComplexGMP *a_par_mp = new ComplexGMP[n_par]; - ComplexGMP *x_ref_mp = new ComplexGMP[n_par]; - ComplexGMP *xtmp = new ComplexGMP[n_par]; - std::vector acoef(n_par + 1), bcoef(n_par + 1); - std::vector y_ref_mp; - std::vector x(n_par), ytmp(n_par); - - // initialize variables - for (int i_par = 0; i_par < n_par; i_par++) { - x_ref_mp[i_par] = x_ref[i_par]; - x[i_par] = x_ref[i_par]; - y_ref_mp.push_back(y_ref[i_par]); - g_func.push_back(std::vector()); - a_par_mp[i_par] = ComplexGMP(); - n_rem_idx.push_back(i_par); - for (int j_par = 0; j_par < n_par; j_par++) { - g_func[i_par].push_back(ComplexGMP()); - } - } - - if (do_greedy == 1) { + // set floating point precision of GMP lib + mpf_set_default_prec(precision); + + // Define constants + const ComplexGMP c_one(std::complex(1.0, 0.0)); + const ComplexGMP c_zero(std::complex(0.0, 0.0)); + const mpf_class tol("1e-6", precision); + + // auxiliary variables + int kdx, n_rem; + std::vector n_rem_idx; + mpf_class deltap, pval; + ComplexGMP pval_in, x_in, y_in, acoef_in, bcoef_in; + std::vector> g_func; + ComplexGMP *a_par_mp = new ComplexGMP[n_par]; + ComplexGMP *x_ref_mp = new ComplexGMP[n_par]; + ComplexGMP *xtmp = new ComplexGMP[n_par]; + std::vector acoef(n_par + 1), bcoef(n_par + 1); + std::vector y_ref_mp; + std::vector x(n_par), ytmp(n_par); + + // initialize variables for (int i_par = 0; i_par < n_par; i_par++) { - x_ref_mp[i_par] = ComplexGMP(); - } - - /********************************** - a_0 coefficient - ***********************************/ - - // Finding the index that maximizes |F| - auto it_zero = - std::max_element(y_ref_mp.begin(), y_ref_mp.end(), - [](const ComplexGMP &a, const ComplexGMP &b) { - return a.abs() < b.abs(); - }); - int kdx_zero = std::distance(y_ref_mp.begin(), it_zero); - - // Add slected point - xtmp[0] = x[kdx_zero]; - ytmp[0] = y_ref_mp[kdx_zero]; - x_ref_mp[0] = x[kdx_zero]; - - // Compute the generating function - thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, 0); - a_par_mp[0] = g_func[0][0]; - - /********************************** - a_1 coefficient - ***********************************/ - - // Finding the index that maximizes abs(x-x_0) while excluding the first_kdx - int kdx_one = kdx_zero + 1; - ComplexGMP max_one; - for (int i = 0; i < x.size(); ++i) { - if (i == kdx_zero) - continue; - - ComplexGMP current_diff = x[i] - x_ref_mp[0]; - if (i == 0 || current_diff.abs() > max_one.abs()) { - kdx_one = i; - max_one = current_diff; - } + x_ref_mp[i_par] = x_ref[i_par]; + x[i_par] = x_ref[i_par]; + y_ref_mp.push_back(y_ref[i_par]); + g_func.push_back(std::vector()); + a_par_mp[i_par] = ComplexGMP(); + n_rem_idx.push_back(i_par); + for (int j_par = 0; j_par < n_par; j_par++) { + g_func[i_par].push_back(ComplexGMP()); + } } - // Add selected point - xtmp[1] = x[kdx_one]; - ytmp[1] = y_ref_mp[kdx_one]; - x_ref_mp[1] = x[kdx_one]; - - // Compute the generating function - thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, 1); - a_par_mp[1] = g_func[1][1]; - - /********************************** - Wallis' method coefficients - ***********************************/ + if (do_greedy == 1) { + for (int i_par = 0; i_par < n_par; i_par++) { + x_ref_mp[i_par] = ComplexGMP(); + } - acoef[0] = a_par_mp[0]; - bcoef[0] = c_one; + /********************************** + a_0 coefficient + ***********************************/ + + // Finding the index that maximizes |F| + auto it_zero = + std::max_element(y_ref_mp.begin(), y_ref_mp.end(), + [](const ComplexGMP &a, const ComplexGMP &b) { + return a.abs() < b.abs(); + }); + int kdx_zero = std::distance(y_ref_mp.begin(), it_zero); + + // Add slected point + xtmp[0] = x[kdx_zero]; + ytmp[0] = y_ref_mp[kdx_zero]; + x_ref_mp[0] = x[kdx_zero]; + + // Compute the generating function + thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, 0); + a_par_mp[0] = g_func[0][0]; + + /********************************** + a_1 coefficient + ***********************************/ + + // Finding the index that maximizes abs(x-x_0) while excluding the first_kdx + int kdx_one = kdx_zero + 1; + ComplexGMP max_one; + for (int i = 0; i < x.size(); ++i) { + if (i == kdx_zero) + continue; + + ComplexGMP current_diff = x[i] - x_ref_mp[0]; + if (i == 0 || current_diff.abs() > max_one.abs()) { + kdx_one = i; + max_one = current_diff; + } + } - // Update indexes of non-visited points - if (kdx_zero > kdx_one) - std::swap(kdx_zero, kdx_one); - - n_rem_idx.erase(n_rem_idx.begin() + kdx_one); - n_rem_idx.erase(n_rem_idx.begin() + kdx_zero); - n_rem = n_par - 2; - - // Add remaining points ensuring min |P_i(x_{1+1}) - F(x_{i+1})| - for (int idx = 2; idx < n_par; ++idx) { - pval = mpf_class("1e9999", precision); // Huge value - - for (int jdx = 0; jdx < n_rem; ++jdx) { - // Compute next convergent P_i(x_{i+1}) - evaluate_thiele_pade_tab_mp( - idx - 2, xtmp, x[n_rem_idx[jdx]], a_par_mp, acoef, bcoef); - pval_in = acoef[idx - 1] / bcoef[idx - 1]; - deltap = (pval_in - y_ref_mp[n_rem_idx[jdx]]).abs(); - - if (deltap < pval) { - pval = deltap; - x_in = x[n_rem_idx[jdx]]; - y_in = y_ref_mp[n_rem_idx[jdx]]; - acoef_in = acoef[idx - 1]; - bcoef_in = bcoef[idx - 1]; - kdx = jdx; + // Add selected point + xtmp[1] = x[kdx_one]; + ytmp[1] = y_ref_mp[kdx_one]; + x_ref_mp[1] = x[kdx_one]; + + // Compute the generating function + thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, 1); + a_par_mp[1] = g_func[1][1]; + + /********************************** + Wallis' method coefficients + ***********************************/ + + acoef[0] = a_par_mp[0]; + bcoef[0] = c_one; + + // Update indexes of non-visited points + if (kdx_zero > kdx_one) + std::swap(kdx_zero, kdx_one); + + n_rem_idx.erase(n_rem_idx.begin() + kdx_one); + n_rem_idx.erase(n_rem_idx.begin() + kdx_zero); + n_rem = n_par - 2; + + // Add remaining points ensuring min |P_i(x_{1+1}) - F(x_{i+1})| + for (int idx = 2; idx < n_par; ++idx) { + pval = mpf_class("1e9999", precision); // Huge value + + for (int jdx = 0; jdx < n_rem; ++jdx) { + // Compute next convergent P_i(x_{i+1}) + evaluate_thiele_pade_tab_mp( + idx - 2, xtmp, x[n_rem_idx[jdx]], a_par_mp, acoef, bcoef); + pval_in = acoef[idx - 1] / bcoef[idx - 1]; + deltap = (pval_in - y_ref_mp[n_rem_idx[jdx]]).abs(); + + if (deltap < pval) { + pval = deltap; + x_in = x[n_rem_idx[jdx]]; + y_in = y_ref_mp[n_rem_idx[jdx]]; + acoef_in = acoef[idx - 1]; + bcoef_in = bcoef[idx - 1]; + kdx = jdx; + } + } + + // Update indexes of non-visited points + n_rem_idx.erase(n_rem_idx.begin() + kdx); + n_rem -= 1; + + // Add the winning point and recompute generating function + x_ref_mp[idx] = x_in; + xtmp[idx] = x_in; + ytmp[idx] = y_in; + + // Rescale Wallis coefficients to avoid overflow + acoef[idx - 1] = acoef_in; + bcoef[idx - 1] = bcoef_in; + if (bcoef_in.abs() > tol) { + acoef[idx - 1] = acoef[idx - 1] / bcoef_in; + acoef[idx - 2] = acoef[idx - 2] / bcoef_in; + bcoef[idx - 2] = bcoef[idx - 2] / bcoef_in; + bcoef[idx - 1] = c_one; + } + + thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, idx); + + // Unpack parameters a_i = g_i(w_i) + a_par_mp[idx] = g_func[idx][idx]; + } + } else { + for (int i_par = 0; i_par < n_par; i_par++) { + thiele_pade_gcoeff_mp(x_ref_mp, y_ref_mp, g_func, i_par); + a_par_mp[i_par] = g_func[i_par][i_par]; } - } - - // Update indexes of non-visited points - n_rem_idx.erase(n_rem_idx.begin() + kdx); - n_rem -= 1; - - // Add the winning point and recompute generating function - x_ref_mp[idx] = x_in; - xtmp[idx] = x_in; - ytmp[idx] = y_in; - - // Rescale Wallis coefficients to avoid overflow - acoef[idx - 1] = acoef_in; - bcoef[idx - 1] = bcoef_in; - if (bcoef_in.abs() > tol) { - acoef[idx - 1] = acoef[idx - 1] / bcoef_in; - acoef[idx - 2] = acoef[idx - 2] / bcoef_in; - bcoef[idx - 2] = bcoef[idx - 2] / bcoef_in; - bcoef[idx - 1] = c_one; - } - - thiele_pade_gcoeff_mp(xtmp, ytmp, g_func, idx); - - // Unpack parameters a_i = g_i(w_i) - a_par_mp[idx] = g_func[idx][idx]; - } - } else { - for (int i_par = 0; i_par < n_par; i_par++) { - thiele_pade_gcoeff_mp(x_ref_mp, y_ref_mp, g_func, i_par); - a_par_mp[i_par] = g_func[i_par][i_par]; } - } - // Create pointer to the parameter struct - pade_model *params = new pade_model; - params->a_par = a_par_mp; - params->precision = precision; - params->n_par = n_par; - params->xref = x_ref_mp; + // Create pointer to the parameter struct + pade_model *params = new pade_model; + params->a_par = a_par_mp; + params->precision = precision; + params->n_par = n_par; + params->xref = x_ref_mp; - // Clean-up - delete[] xtmp; + // Clean-up + delete[] xtmp; - return params; + return params; } diff --git a/GX-AnalyticContinuation/src/pade_mp.h b/GX-AnalyticContinuation/src/pade_mp.h index 6f2394e8..15c97699 100644 --- a/GX-AnalyticContinuation/src/pade_mp.h +++ b/GX-AnalyticContinuation/src/pade_mp.h @@ -1,92 +1,98 @@ -#include +// *************************************************************************************************** +// Copyright (C) 2020-2024 GreenX library +// This file is distributed under the terms of the APACHE2 License. +// +// *************************************************************************************************** + #include +#include /// @brief datatype of a complex number using arbitrary precision numbers for /// real and imaginary part class ComplexGMP { -public: - /// real part - mpf_class real_mp; - /// imag part - mpf_class imag_mp; - - /// @brief initialize a complex number 0 + 0i - ComplexGMP() : real_mp(0), imag_mp(0) {} - - /// @brief initialize a complex number from real and imag part using arbitrary - /// precision numbers - /// @param real_mp real part arbitrary precision - /// @param imag_mp imaginary part arbitrary precision - ComplexGMP(const mpf_class &real_mp, const mpf_class &imag_mp) - : real_mp(real_mp), imag_mp(imag_mp) {} - - /// @brief initialize a complex number from a complex double - /// @param z complex number - ComplexGMP(std::complex z) - : real_mp(real(z)), imag_mp(imag(z)) {} - - /// @brief addition of two arbitrary precision complex numbers - /// @param rhs right summand - /// @return arbitrary precision comoplex number - ComplexGMP operator+(const ComplexGMP &rhs) const { - return ComplexGMP(real_mp + rhs.real_mp, imag_mp + rhs.imag_mp); - } - - /// @brief substraction of two arbitrary precision complex numbers - /// @param rhs subtrahend - /// @return arbitrary precision comoplex number - ComplexGMP operator-(const ComplexGMP &rhs) const { - return ComplexGMP(real_mp - rhs.real_mp, imag_mp - rhs.imag_mp); - } - - /// @brief multiplication of two arbitrary precision complex numbers - /// @param rhs right factor - /// @return arbitrary precision comoplex number - ComplexGMP operator*(const ComplexGMP &rhs) const { - mpf_class newreal_mp = real_mp * rhs.real_mp - imag_mp * rhs.imag_mp; - mpf_class newimag_mp = real_mp * rhs.imag_mp + imag_mp * rhs.real_mp; - return ComplexGMP(newreal_mp, newimag_mp); - } - - /// @brief division of two arbitrary precision complex numbers - /// @param rhs denominator - /// @return arbitrary precision comoplex number - ComplexGMP operator/(const ComplexGMP &rhs) const { - mpf_class denominator = - rhs.real_mp * rhs.real_mp + rhs.imag_mp * rhs.imag_mp; - mpf_class newreal_mp = - (real_mp * rhs.real_mp + rhs.imag_mp * imag_mp) / denominator; - mpf_class newimag_mp = - (rhs.real_mp * imag_mp - real_mp * rhs.imag_mp) / denominator; - return ComplexGMP(newreal_mp, newimag_mp); - } - - /// @brief type casting to complex double - /// @return complex double - operator std::complex() const { - return std::complex(real_mp.get_d(), imag_mp.get_d()); - } - - /// @brief overload the abs function - /// @return ComplexGMP - mpf_class abs() const { - mpf_class re_square = real_mp * real_mp; - mpf_class im_square = imag_mp * imag_mp; - mpf_class rho_square = re_square + im_square; - - mpf_class result; - mpf_sqrt(result.get_mpf_t(), rho_square.get_mpf_t()); - - return result; - } + public: + /// real part + mpf_class real_mp; + /// imag part + mpf_class imag_mp; + + /// @brief initialize a complex number 0 + 0i + ComplexGMP() : real_mp(0), imag_mp(0) {} + + /// @brief initialize a complex number from real and imag part using arbitrary + /// precision numbers + /// @param real_mp real part arbitrary precision + /// @param imag_mp imaginary part arbitrary precision + ComplexGMP(const mpf_class &real_mp, const mpf_class &imag_mp) + : real_mp(real_mp), imag_mp(imag_mp) {} + + /// @brief initialize a complex number from a complex double + /// @param z complex number + ComplexGMP(std::complex z) + : real_mp(real(z)), imag_mp(imag(z)) {} + + /// @brief addition of two arbitrary precision complex numbers + /// @param rhs right summand + /// @return arbitrary precision comoplex number + ComplexGMP operator+(const ComplexGMP &rhs) const { + return ComplexGMP(real_mp + rhs.real_mp, imag_mp + rhs.imag_mp); + } + + /// @brief substraction of two arbitrary precision complex numbers + /// @param rhs subtrahend + /// @return arbitrary precision comoplex number + ComplexGMP operator-(const ComplexGMP &rhs) const { + return ComplexGMP(real_mp - rhs.real_mp, imag_mp - rhs.imag_mp); + } + + /// @brief multiplication of two arbitrary precision complex numbers + /// @param rhs right factor + /// @return arbitrary precision comoplex number + ComplexGMP operator*(const ComplexGMP &rhs) const { + mpf_class newreal_mp = real_mp * rhs.real_mp - imag_mp * rhs.imag_mp; + mpf_class newimag_mp = real_mp * rhs.imag_mp + imag_mp * rhs.real_mp; + return ComplexGMP(newreal_mp, newimag_mp); + } + + /// @brief division of two arbitrary precision complex numbers + /// @param rhs denominator + /// @return arbitrary precision comoplex number + ComplexGMP operator/(const ComplexGMP &rhs) const { + mpf_class denominator = + rhs.real_mp * rhs.real_mp + rhs.imag_mp * rhs.imag_mp; + mpf_class newreal_mp = + (real_mp * rhs.real_mp + rhs.imag_mp * imag_mp) / denominator; + mpf_class newimag_mp = + (rhs.real_mp * imag_mp - real_mp * rhs.imag_mp) / denominator; + return ComplexGMP(newreal_mp, newimag_mp); + } + + /// @brief type casting to complex double + /// @return complex double + operator std::complex() const { + return std::complex(real_mp.get_d(), imag_mp.get_d()); + } + + /// @brief overload the abs function + /// @return ComplexGMP + mpf_class abs() const { + mpf_class re_square = real_mp * real_mp; + mpf_class im_square = imag_mp * imag_mp; + mpf_class rho_square = re_square + im_square; + + mpf_class result; + mpf_sqrt(result.get_mpf_t(), rho_square.get_mpf_t()); + + return result; + } }; struct pade_model { - int n_par; - int precision; - ComplexGMP *a_par; - ComplexGMP *xref; + int n_par; + int precision; + ComplexGMP *a_par; + ComplexGMP *xref; }; // function prototype for fortran interface @@ -96,10 +102,10 @@ std::complex evaluate_thiele_pade_mp(const std::complex x, pade_model *thiele_pade_mp(int n_par, const std::complex *x_ref, const std::complex *y_ref, int do_greedy, int precision); void free_pade_model(pade_model *model) { - if (model != nullptr) { - delete[] model->a_par; - delete[] model->xref; - delete model; - } + if (model != nullptr) { + delete[] model->a_par; + delete[] model->xref; + delete model; + } } } \ No newline at end of file diff --git a/GX-AnalyticContinuation/test/test_pade_approximant.f90 b/GX-AnalyticContinuation/test/test_pade_approximant.f90 index 49a8ae2b..a8516c07 100644 --- a/GX-AnalyticContinuation/test/test_pade_approximant.f90 +++ b/GX-AnalyticContinuation/test/test_pade_approximant.f90 @@ -26,38 +26,48 @@ module test_pade_approximant contains - ! Unfortunately Zofu asserts require real or double precision declarations - ! (etc), which are not equivalent to real(dp), complex(dp)... - ! so one needs a helper function such that a logical can be evaluated. - ! Sensible thing would be to fork the framework and modify, or open a PR. + + !> @brief helper function to assert if the difference between two complex + !! numbers is under a threshold + !! + !! Unfortunately Zofu asserts require real or double precision declarations + !! (etc), which are not equivalent to real(kind=dp), complex(kind=dp)... + !! so one needs a helper function such that a logical can be evaluated. + !! Sensible thing would be to fork the framework and modify, or open a PR. + !! + !! @param[in] a - complex number + !! @param[in] b - complex number + !! @param[in] tol - threshold + !! @return (abs(a - b) <= tol) logical function is_close_complex_dp(a, b, tol) - complex(dp), intent(in) :: a - complex(dp), intent(in) :: b - real(dp), optional, intent(in) :: tol + complex(kind=dp), intent(in) :: a + complex(kind=dp), intent(in) :: b + real(kind=dp), optional, intent(in) :: tol ! abs() evaluates to real, hence tolerance is real - real(dp) :: tolerance = 1.-8_dp + real(kind=dp) :: tolerance = 1.-8_dp if (present(tol)) tolerance = tol is_close_complex_dp = abs(a - b) <= tolerance end function is_close_complex_dp - !> Test the Pade interpolant against the function -1 / (x - x0) + !> @brief Test the Pade interpolant against the function -1 / (x - x0) + !! + !! @param test - test object subroutine test_pade(test) - !> Test object class(unit_test_type), intent(inout) :: test !> N sampling points integer, parameter :: n = 100 !> Variable and function, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) !> Pade approximant of f, and the its reference value - complex(dp) :: f_approx, ref + complex(kind=dp) :: f_approx, ref !> Some function center - complex(dp), parameter :: x0 = cmplx(2.0_dp, 2.0_dp, kind=dp) - complex(dp), parameter :: xx = cmplx(1.0_dp, 1.0_dp, kind=dp) + complex(kind=dp), parameter :: x0 = cmplx(2.0_dp, 2.0_dp, kind=dp) + complex(kind=dp), parameter :: xx = cmplx(1.0_dp, 1.0_dp, kind=dp) !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i !> Test setup @@ -81,7 +91,9 @@ end subroutine test_pade - !> Test the GMP Pade interpolant against the function -1 / (x - x0) + !> @brief Test the GMP Pade interpolant against the function -1 / (x - x0) + !! + !! @param test - test object subroutine test_pade_mp(test) !> Test object class(unit_test_type), intent(inout) :: test @@ -90,15 +102,15 @@ subroutine test_pade_mp(test) integer, parameter :: n = 100 type(params) :: params_thiele !> Variable and function, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) !> Pade approximant of f, and its reference value - complex(dp), dimension(:), allocatable :: xx - complex(dp), dimension(:), allocatable :: f_approx - complex(dp) :: ref + complex(kind=dp), dimension(:), allocatable :: xx + complex(kind=dp), dimension(:), allocatable :: f_approx + complex(kind=dp) :: ref !> Some function center - complex(dp), parameter :: x0 = cmplx(2.0_dp, 2.0_dp, kind=dp) + complex(kind=dp), parameter :: x0 = cmplx(2.0_dp, 2.0_dp, kind=dp) !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i !> Test setup @@ -125,21 +137,24 @@ end subroutine test_pade_mp - !> Test the Thiele-Pade interpolant against the function 1 / (-x^2 + 1) which has poles + !> @brief Test the Thiele-Pade interpolant against the function 1 / (-x^2 + 1) + !! which has poles + !! + !! @param test - test object subroutine test_thiele_pade_poles(test) class(unit_test_type), intent(inout) :: test !> N sampling points integer, parameter :: n = 100 !> Variable, function, and parameters, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) !> Pade approximant of f, and its reference value - complex(dp) :: ref - complex(dp), dimension(1) :: f_approx + complex(kind=dp) :: ref + complex(kind=dp), dimension(1) :: f_approx !> Test point - complex(dp), dimension(1), parameter :: xx = cmplx(1.0_dp, 3.0_dp, kind=dp) + complex(kind=dp), dimension(1), parameter :: xx = cmplx(1.0_dp, 3.0_dp, kind=dp) !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i !> Test setup @@ -163,21 +178,24 @@ end subroutine test_thiele_pade_poles - !> Test the GMP Thiele-Pade interpolant against the function 1 / (-x^2 + 1) which has poles + !> @brief Test the GMP Thiele-Pade interpolant against the function + !! 1 / (-x^2 + 1) which has poles + !! + !! @param test - test object subroutine test_thiele_pade_poles_mp(test) class(unit_test_type), intent(inout) :: test !> N sampling points integer, parameter :: n = 100 !> Variable, function, and parameters, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) type(params) :: params_thiele !> Pade approximant of f, and its reference value - complex(dp), dimension(:), allocatable :: xx - complex(dp), dimension(:), allocatable :: f_approx - complex(dp) :: ref + complex(kind=dp), dimension(:), allocatable :: xx + complex(kind=dp), dimension(:), allocatable :: f_approx + complex(kind=dp) :: ref !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i !> Test setup @@ -203,24 +221,27 @@ end subroutine test_thiele_pade_poles_mp - !> Test the Thiele-Pade interpolant against the function |x| which has a branch point + !> @brief Test the Thiele-Pade interpolant against the function |x| which has + !! a branch point + !! + !! @param test - test object subroutine test_thiele_pade_abs(test) class(unit_test_type), intent(inout) :: test !> N sampling points integer, parameter :: n = 100 !> Newman grid constant - real(dp), parameter :: eta = exp(-1.0_dp / sqrt(dble(n))) - real(dp), parameter :: delta_eta = 0.0005_dp + real(kind=dp), parameter :: eta = exp(-1.0_dp / sqrt(dble(n))) + real(kind=dp), parameter :: delta_eta = 0.0005_dp !> Variable, function, and parameters, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) !> Pade approximant of f, and its reference value - complex(dp) :: ref - complex(dp), dimension(1) :: f_approx + complex(kind=dp) :: ref + complex(kind=dp), dimension(1) :: f_approx !> Test point - complex(dp), dimension(1), parameter :: xx = cmplx(0.7_dp, 0.0_dp, kind=dp) + complex(kind=dp), dimension(1), parameter :: xx = cmplx(0.7_dp, 0.0_dp, kind=dp) !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i, npar !> Test setup @@ -249,25 +270,28 @@ end subroutine test_thiele_pade_abs - !> Test the GMP Thiele-Pade interpolant against the function |x| which has a branch point + !> @brief Test the GMP Thiele-Pade interpolant against the function |x| which + !! has a branch point + !! + !! @param test - test object subroutine test_thiele_pade_abs_mp(test) class(unit_test_type), intent(inout) :: test !> N sampling points integer, parameter :: n = 100 !> Newman grid constant - real(dp), parameter :: eta = exp(-1.0_dp / sqrt(dble(n))) - real(dp), parameter :: delta_eta = 0.0005_dp + real(kind=dp), parameter :: eta = exp(-1.0_dp / sqrt(dble(n))) + real(kind=dp), parameter :: delta_eta = 0.0005_dp !> Variable, function, and parameters, respectively - complex(dp), allocatable :: x(:), f(:) + complex(kind=dp), allocatable :: x(:), f(:) type(params) :: params_thiele !> Pade approximant of f, and its reference value - complex(dp), dimension(:), allocatable :: xx - complex(dp), dimension(:), allocatable :: f_approx - complex(dp) :: ref + complex(kind=dp), dimension(:), allocatable :: xx + complex(kind=dp), dimension(:), allocatable :: f_approx + complex(kind=dp) :: ref !> Test point !> Tolerance - real(dp) :: tol = 1.e-7_dp + real(kind=dp) :: tol = 1.e-7_dp integer :: i, npar !> Test setup