diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/definitions.h | 23 | ||||
| -rw-r--r-- | src/double_util.c | 30 | ||||
| -rw-r--r-- | src/double_util.h | 43 | ||||
| -rw-r--r-- | src/hh_integral.c | 276 | ||||
| -rw-r--r-- | src/hh_integral.h | 63 | ||||
| -rw-r--r-- | src/hh_integral_double.c | 107 | ||||
| -rw-r--r-- | src/hh_integral_double.h | 54 | ||||
| -rw-r--r-- | src/hh_root.c | 144 | ||||
| -rw-r--r-- | src/hh_root.h | 39 | ||||
| -rw-r--r-- | src/hh_root_double.c | 75 | ||||
| -rw-r--r-- | src/hh_root_double.h | 35 | ||||
| -rw-r--r-- | src/hhtop.c | 597 | ||||
| -rw-r--r-- | src/parser.c | 274 | ||||
| -rw-r--r-- | src/parser.h | 37 | ||||
| -rw-r--r-- | src/ss_integral.c | 958 | ||||
| -rw-r--r-- | src/ss_integral.h | 200 | ||||
| -rw-r--r-- | src/ss_integral_double.c | 437 | ||||
| -rw-r--r-- | src/ss_integral_double.h | 170 | ||||
| -rw-r--r-- | src/types.h | 35 | ||||
| -rw-r--r-- | src/zz_integral.c | 237 | ||||
| -rw-r--r-- | src/zz_integral.h | 53 | ||||
| -rw-r--r-- | src/zz_integral_double.c | 89 | ||||
| -rw-r--r-- | src/zz_integral_double.h | 35 |
23 files changed, 4011 insertions, 0 deletions
diff --git a/src/definitions.h b/src/definitions.h new file mode 100644 index 0000000..16c640f --- /dev/null +++ b/src/definitions.h @@ -0,0 +1,23 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#ifndef DEFINITIONS_S +#define DEFINITIONS_S + +#define VERSION "1.0" + +#endif + diff --git a/src/double_util.c b/src/double_util.c new file mode 100644 index 0000000..a4b5949 --- /dev/null +++ b/src/double_util.c @@ -0,0 +1,30 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "double_util.h" +#include <mpfr.h> + +// convert hh_params +int hh_params_todouble(hh_params_double* params_double, hh_params params){ + params_double->omega=params.omega; + params_double->t1=mpfr_get_ld(params.t1, MPFR_RNDN); + params_double->t2=mpfr_get_ld(params.t2, MPFR_RNDN); + params_double->lambda=mpfr_get_ld(params.lambda, MPFR_RNDN); + params_double->sinphi=mpfr_get_ld(params.sinphi, MPFR_RNDN); + params_double->phi=mpfr_get_ld(params.phi, MPFR_RNDN); + params_double->W=mpfr_get_ld(params.W, MPFR_RNDN); + return(0); +} diff --git a/src/double_util.h b/src/double_util.h new file mode 100644 index 0000000..f903809 --- /dev/null +++ b/src/double_util.h @@ -0,0 +1,43 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + convert hh_params to double +*/ + +#ifndef DOUBLE_UTIL_H +#define DOUBLE_UTIL_H + +#include "types.h" + +// params +typedef struct hh_params_double { + int omega; + long double t1; + long double t2; + long double lambda; + long double sinphi; + long double phi; + long double W; +} hh_params_double; + +// format for I function (used to compute sunrise diagrams) +#define TYPE_I_DOUBLE long double (*I)(long double, long double, long double, long double, long double, long double, long double, long double, long double, long double) + +// convert hh_params +int hh_params_todouble(hh_params_double* params_double, hh_params params); + +#endif diff --git a/src/hh_integral.c b/src/hh_integral.c new file mode 100644 index 0000000..f8eda10 --- /dev/null +++ b/src/hh_integral.c @@ -0,0 +1,276 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "hh_integral.h" + +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> +#include <libinum.h> + +int hh_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights){ + mpfr_t lower, upper; + // arguments for first integral + hh_argsint1 args1; + int ret; + + args1.params=params; + args1.abcissa=abcissa; + args1.weights=weights; + + mpfr_inits(lower, upper, NULL); + + // compute the boundaries of the integral over theta + // pi/6 + mpfr_const_pi(upper, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + // -pi/6 + mpfr_neg(lower, upper, MPFR_RNDN); + + // integrate + ret=integrate_gauss_mpfr(out, &hh_integrand1, lower, upper, abcissa, weights, &args1); + + mpfr_clears(lower, upper, NULL); + + return(ret); +} + +// integrand of the integral over theta +int hh_integrand1(mpfr_t* out, mpfr_t theta, void* args){ + hh_argsint2 args2; + mpfr_t lower, upper; + int ret; + + mpfr_inits(upper, lower, NULL); + + // recover parameters + args2.params=((hh_argsint1*)args)->params; + args2.theta=theta; + + // boundaries of the integral over rho + // 1/cos(theta-pi/6) + mpfr_const_pi(upper, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + mpfr_sub(upper, theta, upper, MPFR_RNDN); + mpfr_cos(upper, upper, MPFR_RNDN); + mpfr_ui_div(upper, 1, upper, MPFR_RNDN); + + // 0 + mpfr_set_ui(lower, 0, MPFR_RNDN); + + ret=integrate_gauss_mpfr(out, &hh_integrand2, lower, upper, ((hh_argsint1*)args)->abcissa, ((hh_argsint1*)args)->weights, &args2); + + mpfr_clears(upper, lower, NULL); + + return(ret); +} + +// integrand of the integral over rho +int hh_integrand2(mpfr_t* out, mpfr_t rho, void* args){ + mpfr_t tmp1, tmp2, tmp3; + + mpfr_inits(tmp1, tmp2, tmp3, NULL); + + // out = Omega^2 + // tmp1 = alpha2 + hh_Omega2_alpha2(*out, tmp1, rho, ((hh_argsint2*)args)->theta, tmp2, tmp3); + // tmp1 = m + hh_m(tmp1, ((hh_argsint2*)args)->params.W, ((hh_argsint2*)args)->params.t2, ((hh_argsint2*)args)->params.sinphi, tmp1); + // out = xi^2 + hh_xi2(*out, *out, tmp1, ((hh_argsint2*)args)->params.t1, tmp2); + + // out = rho/sqrt(xi^2)*m + mpfr_sqrt(*out, *out, MPFR_RNDN); + mpfr_div(*out, rho, *out, MPFR_RNDN); + mpfr_mul(*out, *out, tmp1, MPFR_RNDN); + + mpfr_clears(tmp1, tmp2, tmp3, NULL); + + return(0); +} + +// derivative +int hh_d_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights){ + mpfr_t lower, upper; + // arguments for first integral + hh_argsint1 args1; + int ret; + + args1.params=params; + args1.abcissa=abcissa; + args1.weights=weights; + + mpfr_inits(lower, upper, NULL); + + // compute the boundaries of the integral over theta + // pi/6 + mpfr_const_pi(upper, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + // -pi/6 + mpfr_neg(lower, upper, MPFR_RNDN); + + // integrate + ret=integrate_gauss_mpfr(out, &hh_d_integrand1, lower, upper, abcissa, weights, &args1); + + mpfr_clears(lower, upper, NULL); + + return(ret); +} + +// derivative of the integrand of the integral over theta +int hh_d_integrand1(mpfr_t* out, mpfr_t theta, void* args){ + hh_argsint2 args2; + mpfr_t lower, upper; + int ret; + + mpfr_inits(lower, upper, NULL); + + // recover parameters + args2.params=((hh_argsint1*)args)->params; + args2.theta=theta; + + // boundaries of the integral over rho + // 1/cos(theta-pi/6) + mpfr_const_pi(upper, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + mpfr_sub(upper, theta, upper, MPFR_RNDN); + mpfr_cos(upper, upper, MPFR_RNDN); + mpfr_ui_div(upper, 1, upper, MPFR_RNDN); + + // 0 + mpfr_set_ui(lower, 0, MPFR_RNDN); + + ret=integrate_gauss_mpfr(out, &hh_d_integrand2, lower, upper, ((hh_argsint1*)args)->abcissa, ((hh_argsint1*)args)->weights, &args2); + + mpfr_clears(lower, upper, NULL); + + return(ret); +} + +// derivative of the integrand of the integral over rho +int hh_d_integrand2(mpfr_t* out, mpfr_t rho, void* args){ + mpfr_t tmp1, tmp2, tmp3; + + mpfr_inits(tmp1, tmp2, tmp3, NULL); + + // out = Omega^2 + // tmp1 = alpha2 + hh_Omega2_alpha2(*out, tmp1, rho, ((hh_argsint2*)args)->theta, tmp2, tmp3); + // tmp1 = m + hh_m(tmp1, ((hh_argsint2*)args)->params.W, ((hh_argsint2*)args)->params.t2, ((hh_argsint2*)args)->params.sinphi, tmp1); + // out = xi^2 + hh_xi2(*out, *out, tmp1, ((hh_argsint2*)args)->params.t1, tmp2); + + // tmp2 = 1-m^2/xi^2 + mpfr_pow_ui(tmp2, tmp1, 2,MPFR_RNDN); + mpfr_div(tmp2, tmp2, *out, MPFR_RNDN); + mpfr_ui_sub(tmp2, 1, tmp2, MPFR_RNDN); + + // out = rho/sqrt(xi^2)*(1-m^2/xi^2) + mpfr_sqrt(*out, *out, MPFR_RNDN); + mpfr_div(*out, rho, *out, MPFR_RNDN); + mpfr_mul(*out, *out, tmp2, MPFR_RNDN); + + mpfr_clears(tmp1, tmp2, tmp3, NULL); + + return(0); +} + +// Omega^2 and alpha_2 +// provide two initialized tmp mpfr_t's +// Omega2 and alpha2 must be initialized +int hh_Omega2_alpha2(mpfr_t Omega2, mpfr_t alpha2, mpfr_t rho, mpfr_t theta, mpfr_t tmp1, mpfr_t tmp2){ + // Omega2 and alpha2 will be used as tmp variables whenever possible + + // Omega2 = pi + mpfr_const_pi(Omega2, MPFR_RNDN); + + // tmp1 = pi/sqrt(3)*rho + mpfr_sqrt_ui(tmp1, 3, MPFR_RNDN); + mpfr_div(tmp1, Omega2, tmp1, MPFR_RNDN); + mpfr_mul(tmp1, tmp1, rho, MPFR_RNDN); + + // alpha2 = cos(theta) + mpfr_cos(alpha2, theta, MPFR_RNDN); + + // tmp1 = cos(pi/sqrt(3)*rho*cos(theta)) + mpfr_mul(tmp1, tmp1, alpha2, MPFR_RNDN); + //// alpha2 free + mpfr_cos(tmp1, tmp1, MPFR_RNDN); + + // alpha2 = pi/3*(1+rho*sin(theta)) + mpfr_sin(alpha2, theta, MPFR_RNDN); + mpfr_mul(alpha2, alpha2, rho, MPFR_RNDN); + mpfr_add_ui(alpha2, alpha2, 1, MPFR_RNDN); + mpfr_div_ui(alpha2, alpha2, 3, MPFR_RNDN); + mpfr_mul(alpha2, alpha2, Omega2, MPFR_RNDN); + //// Omega2 free + + // Omega2 = cos(pi/3*(1+rho*sin(theta))) + mpfr_cos(Omega2, alpha2, MPFR_RNDN); + + // tmp2 = sin(pi/3*(1+rho*sin(theta))) + mpfr_sin(tmp2, alpha2, MPFR_RNDN); + //// alpha2 free + + // alpha2 = -2*sin(pi/3*(1+rho*sin(theta)))*(cos(pi/3*(1+rho*sin(theta)))+cos(pi/sqrt(3)*rho*cos(theta))) + mpfr_add(alpha2, Omega2, tmp1, MPFR_RNDN); + mpfr_mul(alpha2, alpha2, tmp2, MPFR_RNDN); + //// tmp2 free + mpfr_mul_si(alpha2, alpha2, -2, MPFR_RNDN); + + // tmp1 = cos(pi/3*(1+rho*sin(theta)))-cos(pi/sqrt(3)*rho*cos(theta)) + mpfr_sub(tmp1, Omega2, tmp1, MPFR_RNDN); + // Omega2 = 1+4*cos(pi/3*(1+rho*sin(theta)))*(cos(pi/3*(1+rho*sin(theta)))-cos(pi/sqrt(3)*rho*cos(theta))) + mpfr_mul(Omega2, Omega2, tmp1, MPFR_RNDN); + //// tmp1 free + mpfr_mul_ui(Omega2, Omega2, 4, MPFR_RNDN); + mpfr_add_ui(Omega2, Omega2, 1, MPFR_RNDN); + + return(0); +} + +// m +// out must be initialized +// out and alpha2 can point to the same number +int hh_m(mpfr_t out, mpfr_t W, mpfr_t t2, mpfr_t sinphi, mpfr_t alpha2){ + // out = W-2*t2*sinphi*alpha2 + mpfr_mul(out, alpha2, sinphi, MPFR_RNDN); + mpfr_mul(out, out, t2, MPFR_RNDN); + mpfr_mul_ui(out, out, 2, MPFR_RNDN); + mpfr_sub(out, W, out, MPFR_RNDN); + return(0); +} + +// xi^2 +// provide one initialized tmp mpfr_t +// out must be initialized +// out and Omega2 can point to the same number +// tmp and m can point to the same number +int hh_xi2(mpfr_t out, mpfr_t Omega2, mpfr_t m, mpfr_t t1, mpfr_t tmp){ + // out = t1^2*Omega^2 + mpfr_mul(out, Omega2, t1, MPFR_RNDN); + mpfr_mul(out, out, t1, MPFR_RNDN); + + // tmp = m^2 + mpfr_pow_ui(tmp, m, 2, MPFR_RNDN); + + // out = m^2+t1^2*Omega^2 + mpfr_add(out, out, tmp, MPFR_RNDN); + + return(0); +} diff --git a/src/hh_integral.h b/src/hh_integral.h new file mode 100644 index 0000000..4db9bf3 --- /dev/null +++ b/src/hh_integral.h @@ -0,0 +1,63 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + integrals required for the computation +*/ + +#ifndef HH_INTEGRAL_H +#define HH_INTEGRAL_H + +#include "types.h" +#include <libinum.h> +#include <mpfr.h> + +// extra arguments for the integration over theta (includes the integration options) +typedef struct hh_argsint1 { + array_mpfr abcissa; + array_mpfr weights; + hh_params params; +} hh_argsint1; + +// extra arguments for integration over rho (includes theta) +typedef struct hh_argsint2 { + mpfr_ptr theta; + hh_params params; +} hh_argsint2; + +// the integral +int hh_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights); +// integrand of the first integration +int hh_integrand1(mpfr_t* out, mpfr_t theta, void* args); +// integrand of the second integration +int hh_integrand2(mpfr_t* out, mpfr_t rho, void* args); + +// derivative of the integral +int hh_d_integrate(mpfr_t* out, hh_params params, array_mpfr abcissa, array_mpfr weights); +// integrand of the first integration +int hh_d_integrand1(mpfr_t* out, mpfr_t k1, void* args); +// integrand of the second integration +int hh_d_integrand2(mpfr_t* out, mpfr_t k2, void* args); + +// functions +// Omega^2 and alpha_2 +int hh_Omega2_alpha2(mpfr_t Omega2, mpfr_t alpha2, mpfr_t rho, mpfr_t theta, mpfr_t tmp1, mpfr_t tmp2); +// m +int hh_m(mpfr_t out, mpfr_t W, mpfr_t t2, mpfr_t sinphi, mpfr_t alpha2); +// xi^2 +int hh_xi2(mpfr_t out, mpfr_t Omega2, mpfr_t m, mpfr_t t1, mpfr_t tmp); + +#endif diff --git a/src/hh_integral_double.c b/src/hh_integral_double.c new file mode 100644 index 0000000..2776849 --- /dev/null +++ b/src/hh_integral_double.c @@ -0,0 +1,107 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "hh_integral_double.h" + +#include <math.h> + +#define PI 3.1415926535897932385L + +// compute the integral +int hh_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights){ + hh_argsint1_double args; + int ret; + + args.params=params; + args.abcissa=abcissa; + args.weights=weights; + + ret=integrate_gauss_ldouble(out, &hh_integrand1_double, -PI/6, PI/6, abcissa, weights, &args); + + return(ret); +} + +// integrand of the integral over theta +int hh_integrand1_double(long double* out, long double theta, void* args){ + hh_argsint2_double nargs; + int ret; + hh_argsint1_double* argument=(hh_argsint1_double*)args; + + nargs.params=argument->params; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &hh_integrand2_double, 0, 1./cos(theta-PI/6), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over rho +int hh_integrand2_double(long double* out, long double rho, void* args){ + hh_argsint2_double* argument=(hh_argsint2_double*)args; + long double m, xi, alpha2, O2; + + alpha2=-2*sinl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))+cosl(PI/sqrtl(3.)*rho*cosl(argument->theta))); + O2=1+4*cosl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))-cosl(PI/sqrtl(3.)*rho*cosl(argument->theta))); + m=argument->params.W-2*argument->params.t2*argument->params.sinphi*alpha2; + xi=sqrtl(m*m+argument->params.t1*argument->params.t1*O2); + *out=rho*m/xi; + + return(0); +} + + +// derivative +int hh_d_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights){ + hh_argsint1_double args; + int ret; + + args.params=params; + args.abcissa=abcissa; + args.weights=weights; + + ret=integrate_gauss_ldouble(out, &hh_d_integrand1_double, -PI/6, PI/6, abcissa, weights, &args); + + return(ret); +} + +// derivative of the integrand of the integral over theta +int hh_d_integrand1_double(long double* out, long double theta, void* args){ + hh_argsint2_double nargs; + int ret; + hh_argsint1_double* argument=(hh_argsint1_double*)args; + + nargs.params=argument->params; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &hh_d_integrand2_double, 0, 1./cos(theta-PI/6), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// derivative of the integrand of the integral over rho +int hh_d_integrand2_double(long double* out, long double rho, void* args){ + hh_argsint2_double* argument=(hh_argsint2_double*)args; + long double m, xi, alpha2, O2; + + alpha2=-2*sinl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))+cosl(PI/sqrtl(3.)*rho*cosl(argument->theta))); + O2=1+4*cosl(PI/3*(1+rho*sinl(argument->theta)))*(cosl(PI/3*(1+rho*sinl(argument->theta)))-cosl(PI/sqrtl(3.)*rho*cosl(argument->theta))); + m=argument->params.W-2*argument->params.t2*argument->params.sinphi*alpha2; + xi=sqrtl(m*m+argument->params.t1*argument->params.t1*O2); + *out=rho/xi*(1-m*m/xi/xi); + + return(0); +} + diff --git a/src/hh_integral_double.h b/src/hh_integral_double.h new file mode 100644 index 0000000..9f631cb --- /dev/null +++ b/src/hh_integral_double.h @@ -0,0 +1,54 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#ifndef HH_INTEGRAL_DOUBLE_H +#define HH_INTEGRAL_DOUBLE_H + +#include "double_util.h" +#include "libinum.h" + +// arguments for the integral over theta +typedef struct hh_argsint1_double{ + hh_params_double params; + array_ldouble weights; + array_ldouble abcissa; +} hh_argsint1_double; +// arguments for the integral over rho +typedef struct hh_argsint2_double{ + hh_params_double params; + long double theta; +} hh_argsint2_double; + + +// compute the integral +int hh_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights); + +// integrand of the integral over theta +int hh_integrand1_double(long double* out, long double theta, void* args); + +// integrand of the integral over rho +int hh_integrand2_double(long double* out, long double rho, void* args); + +// derivative +int hh_d_integrate_double(long double* out, hh_params_double params, array_ldouble abcissa, array_ldouble weights); + +// derivative of the integrand of the integral over theta +int hh_d_integrand1_double(long double* out, long double theta, void* args); + +// derivative of the integrand of the integral over rho +int hh_d_integrand2_double(long double* out, long double rho, void* args); + +#endif diff --git a/src/hh_root.c b/src/hh_root.c new file mode 100644 index 0000000..6df4c82 --- /dev/null +++ b/src/hh_root.c @@ -0,0 +1,144 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "hh_root.h" + +#include <mpfr.h> + +#include "hh_integral.h" + +// wrapper for the integration function, used for the Newton scheme +int integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args){ + mpfr_t tmp; + hh_params params; + int ret; + + mpfr_init(tmp); + + mpfr_set(((args_integration*)extra_args)->params.W, in, MPFR_RNDN); + + // out = I+ + ret=hh_integrate(out, ((args_integration*)extra_args)->params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights); + if(ret<0){ + mpfr_clear(tmp); + return(ret); + } + + // clone params (to set sinphi=-sinphi in order to compute I_-) + mpfr_inits(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL); + mpfr_set(params.t1, ((args_integration*)extra_args)->params.t1, MPFR_RNDN); + mpfr_set(params.t2, ((args_integration*)extra_args)->params.t2, MPFR_RNDN); + mpfr_set(params.lambda, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN); + mpfr_set(params.W, ((args_integration*)extra_args)->params.W, MPFR_RNDN); + mpfr_neg(params.sinphi, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN); + + // tmp = I- + ret=hh_integrate(&tmp, params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights); + if(ret<0){ + mpfr_clear(tmp); + mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL); + return(ret); + } + + mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, NULL); + + // out=I+ + I- + mpfr_add(*out, *out, tmp, MPFR_RNDN); + //// tmp free + + // tmp = sqrt(3) + mpfr_sqrt_ui(tmp, 3, MPFR_RNDN); + + // out = W-sqrt(3)*lambda/6*(I+ + I-) + mpfr_mul(*out, *out, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN); + mpfr_div_ui(*out, *out, 6, MPFR_RNDN); + mpfr_mul(*out, *out, tmp, MPFR_RNDN); + mpfr_sub(*out, ((args_integration*)extra_args)->params.W, *out, MPFR_RNDN); + + // tmp = 3*sqrt(3)*t2*sin(phi) + mpfr_mul_ui(tmp, tmp, 3, MPFR_RNDN); + mpfr_mul(tmp, tmp, ((args_integration*)extra_args)->params.t2, MPFR_RNDN); + mpfr_mul(tmp, tmp, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN); + + // W+w*3*sqrt(3)*t2*sin(phi)-sqrt(3)*lambda/6*(I+ + I-) + if(((args_integration*)extra_args)->params.omega==1){ + mpfr_add(*out, *out, tmp, MPFR_RNDN); + } + else{ + mpfr_sub(*out, *out, tmp, MPFR_RNDN); + } + //// tmp free + + mpfr_clear(tmp); + return(0); +} + +// wrapper for the derivative of the integration function, used for the Newton scheme +int d_integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args){ + mpfr_t tmp; + hh_params params; + int ret; + + mpfr_init(tmp); + + mpfr_set(((args_integration*)extra_args)->params.W, in, MPFR_RNDN); + + // out = dI+ + ret=hh_d_integrate(out, ((args_integration*)extra_args)->params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights); + if(ret<0){ + mpfr_clear(tmp); + return(ret); + } + + // clone params (to set sinphi=-sinphi in order to compute dI_-) + mpfr_inits(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL); + mpfr_set(params.t1, ((args_integration*)extra_args)->params.t1, MPFR_RNDN); + mpfr_set(params.t2, ((args_integration*)extra_args)->params.t2, MPFR_RNDN); + mpfr_set(params.lambda, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN); + mpfr_set(params.W, ((args_integration*)extra_args)->params.W, MPFR_RNDN); + params.omega=((args_integration*)extra_args)->params.omega; + mpfr_neg(params.sinphi, ((args_integration*)extra_args)->params.sinphi, MPFR_RNDN); + mpfr_neg(params.phi, ((args_integration*)extra_args)->params.phi, MPFR_RNDN); + + // tmp = dI- + ret=hh_d_integrate(&tmp, params, ((args_integration*)extra_args)->abcissa, ((args_integration*)extra_args)->weights); + if(ret<0){ + mpfr_clear(tmp); + mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL); + return(ret); + } + + mpfr_clears(params.t1, params.t2, params.lambda, params.W, params.sinphi, params.phi, NULL); + + // out=dI+ + dI- + mpfr_add(*out, *out, tmp, MPFR_RNDN); + //// tmp free + + // tmp = sqrt(3) + mpfr_sqrt_ui(tmp, 3, MPFR_RNDN); + + // out = 1-sqrt(3)*lambda/6*(dI+ + dI-) + mpfr_mul(*out, *out, ((args_integration*)extra_args)->params.lambda, MPFR_RNDN); + mpfr_div_ui(*out, *out, 6, MPFR_RNDN); + mpfr_mul(*out, *out, tmp, MPFR_RNDN); + //// tmp free + mpfr_ui_sub(*out, 1, *out, MPFR_RNDN); + + mpfr_clear(tmp); + + return(0); +} + diff --git a/src/hh_root.h b/src/hh_root.h new file mode 100644 index 0000000..2546fa3 --- /dev/null +++ b/src/hh_root.h @@ -0,0 +1,39 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + wrappers for the functions used to find the desired roots +*/ + +#ifndef HH_ROOT_H +#define HH_ROOT_H + +#include "types.h" +#include <libinum.h> + +// the arguments for the integration wrapper +typedef struct args_integration{ + array_mpfr abcissa; + array_mpfr weights; + hh_params params; +} args_integration; + +// wrapper for the integration function, used for the Newton scheme +int integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args); +// wrapper for the derivative of the integration function, used for the Newton scheme +int d_integration_wrapper(mpfr_t* out, mpfr_t in, void* extra_args); + +#endif diff --git a/src/hh_root_double.c b/src/hh_root_double.c new file mode 100644 index 0000000..5474bcb --- /dev/null +++ b/src/hh_root_double.c @@ -0,0 +1,75 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "hh_root_double.h" +#include <math.h> +#include "hh_integral_double.h" + +// wrapper for the integration function, used for the Newton scheme +int hh_integration_wrapper_double(long double* out, long double in, void* args){ + hh_args_integration_double* argument=(hh_args_integration_double*)args; + long double val; + int ret; + + argument->params.W=in; + + ret=hh_integrate_double(&val, argument->params, argument->abcissa, argument->weights); + + if(ret<0){ + return(ret); + } + + // repeat with -sinphi + argument->params.sinphi=-argument->params.sinphi; + + ret=hh_integrate_double(out, argument->params, argument->abcissa, argument->weights); + + // reset sinphi + argument->params.sinphi=-argument->params.sinphi; + + *out=argument->params.W+3*sqrtl(3)*argument->params.omega*argument->params.t2*argument->params.sinphi-argument->params.lambda*sqrtl(3)/6*(*out+val); + + return(ret); +} + +// wrapper for the derivative of the integration function, used for the Newton scheme +int hh_d_integration_wrapper_double(long double* out, long double in, void* args){ + hh_args_integration_double* argument=(hh_args_integration_double*)args; + long double val; + int ret; + + argument->params.W=in; + + ret=hh_d_integrate_double(&val, argument->params, argument->abcissa, argument->weights); + + if(ret<0){ + return(ret); + } + + // repeat with -sinphi + argument->params.sinphi=-argument->params.sinphi; + argument->params.phi=-argument->params.phi; + + ret=hh_d_integrate_double(out, argument->params, argument->abcissa, argument->weights); + + // reset sinphi + argument->params.sinphi=-argument->params.sinphi; + argument->params.phi=-argument->params.phi; + + *out=1-argument->params.lambda*sqrtl(3)/6*(*out+val); + + return(ret); +} diff --git a/src/hh_root_double.h b/src/hh_root_double.h new file mode 100644 index 0000000..2201322 --- /dev/null +++ b/src/hh_root_double.h @@ -0,0 +1,35 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#ifndef HH_ROOT_DOUBLE_H +#define HH_ROOT_DOUBLE_H + +#include <libinum.h> +#include "double_util.h" + +typedef struct hh_args_integration_double{ + hh_params_double params; + array_ldouble abcissa; + array_ldouble weights; +} hh_args_integration_double; + +// wrapper for the integration function, used for the Newton scheme +int hh_integration_wrapper_double(long double* out, long double in, void* args); + +// wrapper for the derivative of the integration function, used for the Newton scheme +int hh_d_integration_wrapper_double(long double* out, long double in, void* args); + +#endif diff --git a/src/hhtop.c b/src/hhtop.c new file mode 100644 index 0000000..0889b76 --- /dev/null +++ b/src/hhtop.c @@ -0,0 +1,597 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include <stdio.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> +#include <math.h> +#include <libinum.h> + +#include "types.h" +#include "hh_integral.h" +#include "hh_root.h" +#include "hh_integral_double.h" +#include "hh_root_double.h" +#include "ss_integral_double.h" +#include "zz_integral.h" +#include "zz_integral_double.h" +#include "parser.h" +#include "definitions.h" + +// usage message +int print_usage(); +// read arguments +int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double); + +// compute the first loop correction to the topological phase diagram +int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order); +// using doubles instead of mpfr +int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order); + +// compute the sunrise diagram +int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads); +// using doubles instead of mpfr +int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads); + +// codes for possible computations +#define COMPUTATION_PHASE 1 +#define COMPUTATION_ZZ 2 +#define COMPUTATION_ZZZZ 3 + +int main(int argc, char* argv[]){ + hh_params params; + mpfr_t tolerance; + unsigned int maxiter; + unsigned int order; + int ret; + unsigned int computation_nr; + unsigned int threads; + unsigned int use_double; + + // default computation: phase diagram + computation_nr=COMPUTATION_PHASE; + + mpfr_inits(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL); + // read command line arguments + ret=read_args(argc, argv, ¶ms, tolerance, &maxiter, &order, &computation_nr, &threads, &use_double); + if(ret<0){ + mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL); + return(-1); + } + if(ret>0){ + mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL); + return(0); + } + + // phase diagram + if(computation_nr==COMPUTATION_PHASE){ + // compute the first-loop correction to the phase diagram + if(use_double==0){ + compute_hh(params, tolerance, maxiter, order); + } + else{ + compute_hh_double(params, tolerance, maxiter, order); + } + } + else if(computation_nr==COMPUTATION_ZZ){ + // compute second-order correction to z1-z2 + if(use_double==0){ + compute_ss(&zz_I, params, tolerance, maxiter, order, threads); + } + else{ + compute_ss_double(&zz_I_double, params, tolerance, maxiter, order, threads); + } + } + else if(computation_nr==COMPUTATION_ZZZZ){ + // compute second-order correction to z1+z2 + if(use_double==0){ + compute_ss(&ZZ_I, params, tolerance, maxiter, order, threads); + } + else{ + compute_ss_double(&ZZ_I_double, params, tolerance, maxiter, order, threads); + } + } + + mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL); + + return(0); +} + +// usage message +int print_usage(){ + fprintf(stderr, "usage:\n hhtop phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax]\n hhtop z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n hhtop z1+z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n\n hhtop -D phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter]\n hhtop -D z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-C threads]\n\n hhtop -V [-v]\n\n"); + return(0); +} + +// read command line arguments +#define CP_FLAG_PARAMS 1 +#define CP_FLAG_ORDER 2 +#define CP_FLAG_TOLERANCE 3 +#define CP_FLAG_MAXITER 4 +#define CP_FLAG_MPFR_PREC 5 +#define CP_FLAG_MPFR_EXP 6 +#define CP_FLAG_THREADS 7 +int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double){ + int i; + int ret; + // temporary long int + long int tmp_lint; + // temporary unsigned int + unsigned int tmp_uint; + // pointers + char* ptr; + // flag that indicates what argument is being read + int flag=0; + // pointer to various arguments + char* tolerance_str=NULL; + char* params_str=NULL; + // whether to print the variables after they are read + int print_vars=0; + // keep track of which flags were used (to check for incompatibilities) + unsigned char pflag=0; + unsigned char Oflag=0; + unsigned char tflag=0; + unsigned char Nflag=0; + unsigned char Pflag=0; + unsigned char Eflag=0; + unsigned char vflag=0; + unsigned char Cflag=0; + unsigned char Dflag=0; + unsigned char Vflag=0; + + // defaults + mpfr_set_d(params->t1, 1., MPFR_RNDN); + mpfr_set_d(params->t2, .1, MPFR_RNDN); + mpfr_set_d(params->lambda, .01, MPFR_RNDN); + mpfr_set_d(params->sinphi, 1., MPFR_RNDN); + mpfr_set_d(tolerance, 1e-11, MPFR_RNDN); + *maxiter=1000000; + *order=10; + params->omega=+1; + *threads=1; + *use_double=0; + + mpfr_const_pi(params->phi, MPFR_RNDN); + mpfr_div_ui(params->phi, params->phi, 2, MPFR_RNDN); + + // default W=-3*sqrt(3)*t2*sin(phi) + mpfr_sqrt_ui(params->W, 3, MPFR_RNDN); + mpfr_mul_si(params->W, params->W, -3, MPFR_RNDN); + mpfr_mul(params->W, params->W, params->sinphi, MPFR_RNDN); + mpfr_mul(params->W, params->W, params->t2, MPFR_RNDN); + + // loop over arguments + for(i=1;i<argc;i++){ + // flag + if(argv[i][0]=='-'){ + for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){ + switch(*ptr){ + // parameters + case 'p': + flag=CP_FLAG_PARAMS; + pflag=1; + break; + // order of the integration + case 'O': + flag=CP_FLAG_ORDER; + Oflag=1; + break; + // tolerance + case 't': + flag=CP_FLAG_TOLERANCE; + tflag=1; + break; + // maximal number of Newton steps + case 'N': + flag=CP_FLAG_MAXITER; + Nflag=1; + break; + // mpfr precision + case 'P': + flag=CP_FLAG_MPFR_PREC; + Pflag=1; + break; + // mpfr emax + case 'E': + flag=CP_FLAG_MPFR_EXP; + Eflag=1; + break; + // print value of variables + case 'v': + print_vars=1; + vflag=1; + break; + // number of threads + case 'C': + flag=CP_FLAG_THREADS; + Cflag=1; + break; + // use doubles instead of mpfr + case 'D': + *use_double=1; + Dflag=1; + // set prec to that of long doubles (for consistency when reading cli arguments with many digits) + mpfr_set_default_prec(64); + break; + // print version + case 'V': + Vflag=1; + break; + default: + fprintf(stderr, "unrecognized option '-%c'\n", *ptr); + print_usage(); + return(-1); + break; + } + } + } + // parameters + else if(flag==CP_FLAG_PARAMS){ + // read str later (after having set the MPFR precision and emax) + params_str=argv[i]; + flag=0; + } + // order of the integration + else if(flag==CP_FLAG_ORDER){ + ret=sscanf(argv[i],"%u",&tmp_uint); + if(ret!=1){ + fprintf(stderr, "error: '-O' should be followed by an unsigned int\n got '%s'\n",argv[i]); + return(-1); + } + *order=tmp_uint; + flag=0; + } + // tolerance + else if(flag==CP_FLAG_TOLERANCE){ + // read str later (after having set the MPFR precision and emax) + tolerance_str=argv[i]; + flag=0; + } + // maximal number of Newton steps + else if(flag==CP_FLAG_MAXITER){ + ret=sscanf(argv[i],"%u",maxiter); + if(ret!=1){ + fprintf(stderr, "error: '-N' should be followed by a positive int\n got '%s'\n",argv[i]); + return(-1); + } + flag=0; + } + // mpfr precision + else if(flag==CP_FLAG_MPFR_PREC){ + ret=sscanf(argv[i],"%ld",&tmp_lint); + if(ret!=1){ + fprintf(stderr, "error: '-P' should be followed by a long int\n got '%s'\n",argv[i]); + return(-1); + } + mpfr_set_default_prec(tmp_lint); + flag=0; + } + // mpfr emax + else if(flag==CP_FLAG_MPFR_EXP){ + ret=sscanf(argv[i],"%ld",&tmp_lint); + if(ret!=1){ + fprintf(stderr, "error: '-E' should be followed by a long int\n got '%s'\n",argv[i]); + return(-1); + } + mpfr_set_emax(tmp_lint); + flag=0; + } + // number of threads to use for the computation + else if(flag==CP_FLAG_THREADS){ + ret=sscanf(argv[i],"%u",threads); + if(ret!=1){ + fprintf(stderr, "error: '-C' should be followed by a positive int\n got '%s'\n",argv[i]); + return(-1); + } + flag=0; + } + // computation to run + else{ + if(str_cmp(argv[i], "phase")==1){ + *computation_nr=COMPUTATION_PHASE; + } + else if(str_cmp(argv[i], "z1-z2")==1){ + *computation_nr=COMPUTATION_ZZ; + } + else if(str_cmp(argv[i], "z1+z2")==1){ + *computation_nr=COMPUTATION_ZZZZ; + } + else{ + fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]); + print_usage(); + return(-1); + } + flag=0; + } + } + if(tolerance_str!=NULL){ + ret=mpfr_set_str(tolerance, tolerance_str, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "error: '-t' should be followed by an MPFR floating point number\n got '%s'\n", tolerance_str); + return(-1); + } + } + if(params_str!=NULL){ + ret=read_params(params, params_str); + if(ret<0){ + return(ret); + } + } + + // check for incompatible flags + if((Vflag==1 && (pflag!=0 || Oflag!=0 || tflag!=0 || Nflag!=0 || Pflag!=0 || Eflag!=0 || Cflag!=0 || Dflag!=0)) || \ + (*computation_nr!=COMPUTATION_ZZ && *computation_nr!=COMPUTATION_ZZZZ && Cflag==1) || \ + (Dflag==1 && (Pflag==1 || Eflag==1)) \ + ){ + print_usage(); + return(-1); + } + + // print version and exit + if(Vflag==1){ + printf("hhtop " VERSION "\n"); + printf("libinum " LIBINUM_VERSION "\n"); + if(vflag==1){ + // print datatype information + printf("\n\n"); + print_datatype_info(stdout); + } + return(1); + } + + // print variables + if(print_vars==1){ + fprintf(stderr, "t1="); + fprint_mpfr(stderr,params->t1); + fprintf(stderr, "\n"); + + fprintf(stderr, "t2="); + fprint_mpfr(stderr,params->t2); + fprintf(stderr, "\n"); + + fprintf(stderr, "lambda="); + fprint_mpfr(stderr,params->lambda); + fprintf(stderr, "\n"); + + fprintf(stderr, "phi="); + fprint_mpfr(stderr,params->phi); + fprintf(stderr, "\n"); + + fprintf(stderr, "sinphi="); + fprint_mpfr(stderr,params->sinphi); + fprintf(stderr, "\n"); + + fprintf(stderr, "W="); + fprint_mpfr(stderr,params->W); + fprintf(stderr, "\n"); + + fprintf(stderr, "omega=%d\n",params->omega); + + fprintf(stderr, "\ntolerance="); + fprint_mpfr(stderr,tolerance); + fprintf(stderr, "\n"); + + fprintf(stderr, "order of integration=%d\n", *order); + + fprintf(stderr, "\nMPFR precision=%ld\n", mpfr_get_default_prec()); + fprintf(stderr, "MPFR emax=%ld\n", mpfr_get_emax()); + + fprintf(stderr, "\n"); + } + + return(0); +} + +// compute the first loop correcion to the topological phase diagram +int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){ + mpfr_t val; + int ret; + args_integration args_int; + + // compute weights + ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &(args_int.abcissa), &(args_int.weights)); + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n try increasing the precision or the tolerance\n"); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n"); + return(ret); + } + + // set args + args_int.params=params; + mpfr_init(args_int.params.W); + + mpfr_init(val); + // initial value + mpfr_sqrt_ui(val, 3, MPFR_RNDN); + mpfr_mul_ui(val, val, 3, MPFR_RNDN); + mpfr_mul(val, val, params.sinphi, MPFR_RNDN); + mpfr_mul(val, val, params.t2, MPFR_RNDN); + if(params.omega==1){ + mpfr_neg(val, val, MPFR_RNDN); + } + + // compute root + ret=root_newton_inplace_mpfr(&val, &integration_wrapper, &d_integration_wrapper, tolerance, maxiter, &args_int); + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n try increasing the precision or the tolerance\n"); + mpfr_clear(val); + mpfr_clear(args_int.params.W); + array_mpfr_free(args_int.abcissa); + array_mpfr_free(args_int.weights); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n"); + mpfr_clear(val); + mpfr_clear(args_int.params.W); + array_mpfr_free(args_int.abcissa); + array_mpfr_free(args_int.weights); + return(ret); + } + + fprint_mpfr(stdout, val); + printf("\n"); + + mpfr_clear(val); + mpfr_clear(args_int.params.W); + array_mpfr_free(args_int.abcissa); + array_mpfr_free(args_int.weights); + + return(0); +} +// using double instead of mpfr +int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){ + long double tolerance_d; + hh_args_integration_double args_int; + long double val; + int ret; + + // convert mpfr to double + tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN); + hh_params_todouble(&(args_int.params), params); + + // compute weights + ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &(args_int.abcissa), &(args_int.weights)); + + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n"); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n"); + return(ret); + } + + + // initial value + val=-args_int.params.omega*3*sqrtl(3)*args_int.params.t2*args_int.params.sinphi; + + ret=root_newton_inplace_ldouble(&val, &hh_integration_wrapper_double, &hh_d_integration_wrapper_double, tolerance_d, maxiter, &args_int); + + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n"); + array_ldouble_free(args_int.abcissa); + array_ldouble_free(args_int.weights); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n"); + array_ldouble_free(args_int.abcissa); + array_ldouble_free(args_int.weights); + return(ret); + } + + printf("% .19Le\n",val); + + array_ldouble_free(args_int.abcissa); + array_ldouble_free(args_int.weights); + + return(0); +} + +// compute the sunrise diagram +int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){ + mpfr_t val; + int ret; + array_mpfr abcissa; + array_mpfr weights; + + // compute weights + ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &abcissa, &weights); + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n try increasing the precision, the tolerance, or the maximal number of Newton steps\n"); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n"); + return(ret); + } + + mpfr_init(val); + + // compute integral + ret=ss_integrate(&val, I, params, abcissa, weights, threads); + // return codes + if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered: the integrand is singular\n"); + mpfr_clear(val); + array_mpfr_free(abcissa); + array_mpfr_free(weights); + return(ret); + } + + fprint_mpfr(stdout, val); + printf("\n"); + + mpfr_clear(val); + array_mpfr_free(abcissa); + array_mpfr_free(weights); + + return(0); +} +// using doubles instead of mpfr +int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){ + long double tolerance_d; + hh_params_double params_d; + long double val; + int ret; + array_ldouble abcissa; + array_ldouble weights; + + // convert mpfr to double + tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN); + hh_params_todouble(¶ms_d, params); + + // compute weights + ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &abcissa, &weights); + // return codes + if(ret==LIBINUM_ERROR_MAXITER){ + fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n try increasing the tolerance, or the maximal number of Newton steps\n"); + return(ret); + } + else if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n"); + return(ret); + } + + // compute integral + ret=ss_integrate_double(&val, I, params_d, abcissa, weights, threads); + // return codes + if(ret==LIBINUM_ERROR_NAN){ + fprintf(stderr, "error: infinity encountered: the integrand is singular\n"); + array_ldouble_free(abcissa); + array_ldouble_free(weights); + return(ret); + } + + printf("% .19Le\n",val); + + array_ldouble_free(abcissa); + array_ldouble_free(weights); + + return(0); +} + diff --git a/src/parser.c b/src/parser.c new file mode 100644 index 0000000..9c4801a --- /dev/null +++ b/src/parser.c @@ -0,0 +1,274 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "parser.h" +#include <stdlib.h> +#include <stdio.h> +#include <stdarg.h> +// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears +#define MPFR_USE_VA_LIST +#include <mpfr.h> + +// read parameter string +#define P_VAR_T1 1 +#define P_VAR_T2 2 +#define P_VAR_LAMBDA 3 +#define P_VAR_OMEGA 4 +#define P_VAR_SINPHI 5 +#define P_VAR_W 6 +#define P_VAR_PHI 7 +int read_params(hh_params* params, char* str){ + char* ptr; + char* buffer; + char* buffer_ptr; + int ret; + int var_flag=0; + unsigned int reset_W=0; + unsigned int set_phi=0; + unsigned int set_sinphi=0; + mpfr_t tmp, tmp2; + + buffer=calloc(str_len(str), sizeof(char)); + buffer_ptr=buffer; + *buffer_ptr='\0'; + + + for(ptr=str; *ptr!='\0'; ptr++){ + switch(*ptr){ + // left side of equation + case '=': + if(str_cmp(buffer, "t1")==1){ + var_flag=P_VAR_T1; + } + else if(str_cmp(buffer, "t2")==1){ + var_flag=P_VAR_T2; + } + else if(str_cmp(buffer, "lambda")==1){ + var_flag=P_VAR_LAMBDA; + } + else if(str_cmp(buffer, "omega")==1){ + var_flag=P_VAR_OMEGA; + } + else if(str_cmp(buffer, "sinphi")==1){ + var_flag=P_VAR_SINPHI; + // reset W to -3*omega*sqrt(3)*sin(phi) provided W is not set explicitly + if(reset_W==0){ + reset_W=1; + } + // do not allow setting both sinphi and phi + set_sinphi=1; + } + else if(str_cmp(buffer, "phi")==1){ + var_flag=P_VAR_PHI; + // reset W to -3*omega*sqrt(3)*sin(phi) provided W is not set explicitly + if(reset_W==0){ + reset_W=1; + } + // do not allow setting both sinphi and phi + set_phi=1; + } + else if(str_cmp(buffer, "W")==1){ + var_flag=P_VAR_W; + // do not reset W + reset_W=2; + } + else{ + fprintf(stderr, "parsing error: unrecognized token '%s'\n", buffer); + free(buffer); + return(-1); + } + + // reset buffer + buffer_ptr=buffer; + *buffer_ptr='\0'; + break; + + // assign value + case ';': + if(var_flag==P_VAR_T1){ + ret=mpfr_set_str(params->t1, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: t1 must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + } + else if(var_flag==P_VAR_T2){ + ret=mpfr_set_str(params->t2, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: t2 must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + } + else if(var_flag==P_VAR_LAMBDA){ + ret=mpfr_set_str(params->lambda, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: lambda must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + } + else if(var_flag==P_VAR_SINPHI){ + if(set_phi==1){ + fprintf(stderr, "error: cannot set both phi and sinphi\n"); + free(buffer); + return(-1); + } + + ret=mpfr_set_str(params->sinphi, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: sinphi must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + // check value + if(mpfr_cmp_ui(params->sinphi,1)>0 || mpfr_cmp_si(params->sinphi,-1)<0){ + fprintf(stderr, "error: sinphi must be in [-1,1]\n"); + free(buffer); + return(-1); + } + // set phi + mpfr_asin(params->phi, params->sinphi, MPFR_RNDN); + } + else if(var_flag==P_VAR_PHI){ + if(set_sinphi==1){ + fprintf(stderr, "error: cannot set both phi and sinphi\n"); + free(buffer); + return(-1); + } + + ret=mpfr_set_str(params->phi, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: phi must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + // set sinphi + mpfr_sin(params->sinphi, params->phi, MPFR_RNDN); + } + else if(var_flag==P_VAR_W){ + ret=mpfr_set_str(params->W, buffer, 10, MPFR_RNDN); + if(ret<0){ + fprintf(stderr, "parsing error: W must be assigned to an MPFR floating point number\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + } + else if(var_flag==P_VAR_OMEGA){ + ret=sscanf(buffer, "%d", &(params->omega)); + if(ret!=1){ + fprintf(stderr, "parsing error: omega must be assigned to an integer\n got '%s'\n", buffer); + free(buffer); + return(-1); + } + // check value + if(params->omega!=1 && params->omega!=-1){ + fprintf(stderr, "error: omega must be either +1 or -1\n"); + return(-1); + } + } + else{ + fprintf(stderr, "parsing error: read right side of equation, but the matching token was not found\n"); + free(buffer); + return(-1); + } + + var_flag=0; + + // reset buffer + buffer_ptr=buffer; + *buffer_ptr='\0'; + break; + + // add to buffer + default: + buffer_ptr=str_addchar(buffer_ptr, *ptr); + break; + } + } + + // check that all variables were read + if(*buffer_ptr!='\0'){ + fprintf(stderr, "parsing error: trailing characters: '%s'\n", buffer); + free(buffer); + return(-1); + } + if(var_flag!=0){ + fprintf(stderr, "parsing error: empty assignment at the end of the string\n"); + free(buffer); + return(-1); + } + + // check that 3*abs(t2)<abs(t1) + mpfr_inits(tmp, tmp2, NULL); + mpfr_abs(tmp, params->t2, MPFR_RNDN); + mpfr_mul_ui(tmp, tmp, 3, MPFR_RNDN); + mpfr_abs(tmp2, params->t1, MPFR_RNDN); + if(mpfr_cmp(tmp, tmp2)>0){ + fprintf(stderr, "error: |t2| must be smaller than |t1|/3\n"); + return(-1); + } + + // if W was not set, reset its default to -3*omega*sqrt(3)*t2*sin(phi) + if(reset_W==1){ + mpfr_sqrt_ui(params->W, 3, MPFR_RNDN); + mpfr_mul_ui(params->W, params->W, 3, MPFR_RNDN); + mpfr_mul(params->W, params->W, params->sinphi, MPFR_RNDN); + mpfr_mul(params->W, params->W, params->t2, MPFR_RNDN); + if(params->omega==1){ + mpfr_neg(params->W, params->W, MPFR_RNDN); + } + } + + mpfr_clears(tmp, tmp2, NULL); + free(buffer); + return(0); +} + + +// length of a string +int str_len(char* str){ + char* ptr=str; + int ret=0; + while(*ptr!='\0'){ret++;ptr++;} + return(ret); +} + +// compare strings +int str_cmp(char* str1, char* str2){ + char* ptr1=str1; + char* ptr2=str2; + while(*ptr1==*ptr2 && *ptr1!='\0' && *ptr2!='\0'){ + ptr1++; + ptr2++; + } + if(*ptr1=='\0' && *ptr2=='\0'){ + return(1); + } + else{ + return(0); + } +} + +// append a character to the end of a string at the location pointed at by 'ptr' +char* str_addchar(char* ptr, const char c){ + *ptr=c; + ptr++; + *ptr='\0'; + return(ptr); +} + diff --git a/src/parser.h b/src/parser.h new file mode 100644 index 0000000..874889c --- /dev/null +++ b/src/parser.h @@ -0,0 +1,37 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + parser for hhtop +*/ + +#ifndef PARSER_H +#define PARSER_H + +#include "types.h" + +// read parameter string +int read_params(hh_params* params, char* str); + +// utilities +// length of a string +int str_len(char* str); +// compare strings +int str_cmp(char* str1, char* str2); +// append a character to the end of a string at the location pointed at by 'ptr' +char* str_addchar(char* ptr, const char c); + +#endif diff --git a/src/ss_integral.c b/src/ss_integral.c new file mode 100644 index 0000000..cb39ea8 --- /dev/null +++ b/src/ss_integral.c @@ -0,0 +1,958 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "ss_integral.h" +#include <pthread.h> + +// compute the integral +// takes a pointer to an integrand function 'I' as an argument. 'I' must have the following arguments +// int I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps, struct ss_cache cache) +int ss_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, unsigned int threads){ + mpfr_t val; + array_2_mpfr tmpss; + struct ss_cache cache; + unsigned int i; + + // compute pi and sqrt3 + mpfr_inits(cache.pi, cache.sqrt3, NULL); + mpfr_const_pi(cache.pi, MPFR_RNDN); + mpfr_sqrt_ui(cache.sqrt3, 3, MPFR_RNDN); + + mpfr_init(val); + + // init tmps + array_2_mpfr_init(&tmpss, 8*threads); + tmpss.length=8*threads; + for(i=0;i<threads;i++){ + // 8 array_mpfr's per thread + // do not allocate too much memory (I'm trying to keep a modicum of optimality) + array_mpfr_init(tmpss.values+8*i+0, 2); + array_mpfr_init(tmpss.values+8*i+1, 4); + array_mpfr_init(tmpss.values+8*i+2, 3); + array_mpfr_init(tmpss.values+8*i+3, 4); + array_mpfr_init(tmpss.values+8*i+4, 2); + array_mpfr_init(tmpss.values+8*i+5, 4); + array_mpfr_init(tmpss.values+8*i+6, 12); + array_mpfr_init(tmpss.values+8*i+7, 9); + } + + // compute A_FF, A_RF and A_RR + ssFF_integrate(out, I, params, abcissa, weights, cache, tmpss, threads); + + ssRF_integrate(&val, I, params, abcissa, weights, cache, tmpss, threads); + mpfr_mul_ui(val, val, 2, MPFR_RNDN); + mpfr_add(*out, *out, val, MPFR_RNDN); + + ssRR_integrate(&val, I, params, abcissa, weights, cache, tmpss, threads); + mpfr_add(*out, *out, val, MPFR_RNDN); + + mpfr_clears(cache.pi, cache.sqrt3, NULL); + mpfr_clear(val); + + array_2_mpfr_free(tmpss); + return(0); +} + + +// compute the integral A_FF +int ssFF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){ + ssFF_argsint_rho args; + int ret; + mpfr_t upper, lower; + array_pthread_t thread_ids; + + // init + thread_ids.length=0; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + args.cache=cache; + args.tmpss=tmpss; + args.thread_ids=&thread_ids; + + // bounds + mpfr_inits(upper, lower, NULL); + mpfr_set_ui(lower, 0, MPFR_RNDN); + mpfr_set_ui(upper, 1, MPFR_RNDN); + + if(threads==1){ + ret=integrate_gauss_mpfr(out, &ssFF_integrand_rho, lower, upper, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_mpfr(out, &ssFF_integrand_rho, lower, upper, abcissa, weights, &args, threads, args.thread_ids); + array_pthread_t_free(*(args.thread_ids)); + } + + // extra factor + mpfr_mul_ui(*out, *out, 6, MPFR_RNDN); + + mpfr_clears(upper, lower, NULL); + + return(ret); +} + +// integrand of the integral over rho in A_FF +int ssFF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){ + ssFF_argsint_theta nargs; + int ret; + ssFF_argsint_rho* argument=(ssFF_argsint_rho*)args; + array_mpfr* tmps1; + array_mpfr* tmps2; + int thread_id; + + // the number of the thread + if(argument->thread_ids->length>1){ + thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids)); + if(thread_id<0){ + fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n"); + return(-1); + } + } + else{ + thread_id=0; + } + + // tmps used in this function (not passed to integrand_theta) + tmps1=argument->tmpss.values+8*thread_id; + tmps2=argument->tmpss.values+8*thread_id+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.rho=rho; + + // copy tmps + // only pass tmps reserved for this thread + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+8*thread_id+2; + // 6 tmps in integrand_theta + nargs.tmpss.length=6; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN); + mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssFF_integrand_theta, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over theta in A_FF +int ssFF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){ + ssFF_argsint_psi nargs; + int ret; + ssFF_argsint_theta* argument=(ssFF_argsint_theta*)args; + unsigned int a; + // tmps used in this function (not passed to integrand_psi) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.rho=argument->rho; + nargs.theta=theta; + + // copy tmps + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 4 tmps in integrand_psi + nargs.tmpss.length=4; + + // alloc tmps if needed + array_mpfr_alloc_tmps(3, tmps1); + + // split the integral into four parts to improve the precision + mpfr_set_ui(*out, 0, MPFR_RNDN); + for(a=0;a<4;a++){ + // bounds + mpfr_set_ui(tmps1->values[0], a, MPFR_RNDN); + mpfr_sub_ui(tmps1->values[0], tmps1->values[0], 2, MPFR_RNDN); + mpfr_div_ui(tmps1->values[0], tmps1->values[0], 4, MPFR_RNDN); + mpfr_mul(tmps1->values[0], tmps1->values[0], nargs.cache.pi, MPFR_RNDN); + + mpfr_set_ui(tmps1->values[1], a, MPFR_RNDN); + mpfr_sub_ui(tmps1->values[1], tmps1->values[1], 1, MPFR_RNDN); + mpfr_div_ui(tmps1->values[1], tmps1->values[1], 4, MPFR_RNDN); + mpfr_mul(tmps1->values[1], tmps1->values[1], nargs.cache.pi, MPFR_RNDN); + + + ret=integrate_gauss_smarttmp_mpfr(tmps1->values+2, &ssFF_integrand_psi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + if(ret<0){ + break; + } + mpfr_add(*out, *out, tmps1->values[2], MPFR_RNDN); + } + + return(ret); +} + +// integrand of the integral over psi in A_FF +int ssFF_integrand_psi(mpfr_t* out, mpfr_t psi, void* args){ + ssFF_argsint_z nargs; + int ret; + ssFF_argsint_psi* argument=(ssFF_argsint_psi*)args; + // tmps used in this function (not passed to integrand_z) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.cache=argument->cache; + nargs.rho=argument->rho; + nargs.theta=argument->theta; + nargs.psi=psi; + + // copy tmps + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 2 tmps in integrand_z + nargs.tmpss.length=2; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN); + mpfr_set_ui(tmps1->values[1], 1, MPFR_RNDN); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssFF_integrand_z, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over z in A_FF +int ssFF_integrand_z(mpfr_t* out, mpfr_t z, void* args){ + ssFF_argsint_z* a=(ssFF_argsint_z*)args; + mpfr_t* tmps; + // tmps used in this function + array_mpfr* tmps1=a->tmpss.values; + array_mpfr* tmps2=a->tmpss.values+1; + + // alloc tmps if needed + array_mpfr_alloc_tmps(7, tmps1); + + tmps=tmps1->values; + + // r (store in tmps[0]) + mpfr_cos(tmps[0], a->psi, MPFR_RNDN); + mpfr_add_ui(tmps[1], tmps[0], 1, MPFR_RNDN); + mpfr_ui_sub(tmps[0], 1, tmps[0], MPFR_RNDN); + mpfr_div(tmps[0], tmps[0], tmps[1], MPFR_RNDN); + mpfr_mul(tmps[0], tmps[0], z, MPFR_RNDN); + mpfr_add_ui(tmps[0], tmps[0], 1, MPFR_RNDN); + mpfr_ui_sub(tmps[1], 1, z, MPFR_RNDN); + mpfr_mul(tmps[0], tmps[0], tmps[1], MPFR_RNDN); + + // varphi (store in tmps[1]) + mpfr_mul_ui(tmps[1], a->psi, 2, MPFR_RNDN); + mpfr_sin(tmps[2], tmps[1], MPFR_RNDN); + mpfr_cos(tmps[1], tmps[1], MPFR_RNDN); + mpfr_pow_ui(tmps[2], tmps[2], 2, MPFR_RNDN); + mpfr_mul(tmps[2], tmps[2], tmps[0], MPFR_RNDN); + mpfr_div_ui(tmps[2], tmps[2], 2, MPFR_RNDN); + mpfr_sub(tmps[1], tmps[1], tmps[2], MPFR_RNDN); + mpfr_acos(tmps[1], tmps[1], MPFR_RNDN); + if(mpfr_cmp_ui(a->psi, 0)<0){ + mpfr_neg(tmps[1], tmps[1], MPFR_RNDN); + } + + // jacobian + mpfr_mul_ui(*out, a->psi, 2, MPFR_RNDN); + mpfr_cos(*out, *out, MPFR_RNDN); + mpfr_mul(*out, *out, tmps[0], MPFR_RNDN); + mpfr_add_ui(*out, *out, 1, MPFR_RNDN); + mpfr_cos(tmps[2], a->psi, MPFR_RNDN); + mpfr_pow_ui(tmps[3], tmps[2], 2, MPFR_RNDN); + mpfr_mul(tmps[3], tmps[3], tmps[0], MPFR_RNDN); + mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN); + mpfr_sqrt(tmps[3], tmps[3], MPFR_RNDN); + mpfr_div(*out, *out, tmps[3], MPFR_RNDN); + mpfr_add_ui(tmps[2], tmps[2], 1, MPFR_RNDN); + mpfr_div(*out, *out, tmps[2], MPFR_RNDN); + mpfr_mul(*out, *out, tmps[0], MPFR_RNDN); + mpfr_mul(*out, *out, a->rho, MPFR_RNDN); + mpfr_mul(*out, *out, a->rho, MPFR_RNDN); + mpfr_mul(*out, *out, a->rho, MPFR_RNDN); + mpfr_mul_ui(*out, *out, 4, MPFR_RNDN); + + // cutoffs + ss_cutoff(tmps[2], a->rho, tmps[3], tmps[4]); + mpfr_mul(*out, *out, tmps[2], MPFR_RNDN); + mpfr_mul(tmps[2], a->rho, tmps[0], MPFR_RNDN); + ss_cutoff(tmps[2], tmps[2], tmps[3], tmps[4]); + mpfr_mul(*out, *out, tmps[2], MPFR_RNDN); + + // q1 (store in tmps[2]) + mpfr_add(tmps[2], a->theta, tmps[1], MPFR_RNDN); + mpfr_cos(tmps[2], tmps[2], MPFR_RNDN); + mpfr_mul(tmps[2], tmps[2], tmps[0], MPFR_RNDN); + mpfr_mul(tmps[2], tmps[2], a->rho, MPFR_RNDN); + mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN); + // q2 (store in tmps[3]) + mpfr_add(tmps[3], a->theta, tmps[1], MPFR_RNDN); + mpfr_sin(tmps[3], tmps[3], MPFR_RNDN); + mpfr_mul(tmps[3], tmps[3], tmps[0], MPFR_RNDN); + mpfr_mul(tmps[3], tmps[3], a->rho, MPFR_RNDN); + mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN); + if(a->params.omega==-1){ + mpfr_neg(tmps[3], tmps[3], MPFR_RNDN); + } + // p1 (store in tmps[0]) + mpfr_cos(tmps[0], a->theta, MPFR_RNDN); + mpfr_mul(tmps[0], tmps[0], a->rho, MPFR_RNDN); + mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN); + // p2 (store in tmps[1]) + mpfr_sin(tmps[1], a->theta, MPFR_RNDN); + mpfr_mul(tmps[1], tmps[1], a->rho, MPFR_RNDN); + mpfr_add_ui(tmps[1], tmps[1], 1, MPFR_RNDN); + if(a->params.omega==-1){ + mpfr_neg(tmps[1], tmps[1], MPFR_RNDN); + } + // F1 (store in tmps[4]) + mpfr_add(tmps[4], tmps[0], tmps[2], MPFR_RNDN); + mpfr_sub(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN); + // F2 (store in tmps[5]) + mpfr_add(tmps[5], tmps[1], tmps[3], MPFR_RNDN); + mpfr_sub_si(tmps[5], tmps[5], a->params.omega, MPFR_RNDN); + + // I + (*(a->I))(tmps[6], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache); + mpfr_mul(*out, *out, tmps[6], MPFR_RNDN); + + return(0); +} + + +// compute the integral A_RF +int ssRF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){ + ssRF_argsint_theta args; + int ret; + mpfr_t val; + mpfr_t lower, upper; + unsigned int a; + array_pthread_t thread_ids; + + // init + thread_ids.length=0; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + args.cache=cache; + args.tmpss=tmpss; + args.thread_ids=&thread_ids; + + mpfr_inits(val, lower, upper, NULL); + + // split the integral in three parts to improve the precision + mpfr_set_ui(*out, 0, MPFR_RNDN); + for(a=0;a<3;a++){ + // bounds + mpfr_set_ui(lower, a, MPFR_RNDN); + mpfr_mul_ui(lower, lower, 4, MPFR_RNDN); + mpfr_sub_ui(lower, lower, 1, MPFR_RNDN); + mpfr_div_ui(lower, lower, 6, MPFR_RNDN); + mpfr_mul(lower, lower, args.cache.pi, MPFR_RNDN); + + mpfr_set_ui(upper, a, MPFR_RNDN); + mpfr_mul_ui(upper, upper, 4, MPFR_RNDN); + mpfr_add_ui(upper, upper, 3, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + mpfr_mul(upper, upper, args.cache.pi, MPFR_RNDN); + + if(threads==1){ + ret=integrate_gauss_mpfr(&val, &ssRF_integrand_theta, lower, upper, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_mpfr(&val, &ssRF_integrand_theta, lower, upper, abcissa, weights, &args, threads, args.thread_ids); + array_pthread_t_free(*(args.thread_ids)); + } + if(ret<0){ + break; + } + + mpfr_add(*out, *out, val, MPFR_RNDN); + } + + // extra factor + mpfr_mul_ui(*out, *out, 3, MPFR_RNDN); + + mpfr_clears(val, lower, upper, NULL); + return(ret); +} + +// integrand of the integral over theta in A_RF +int ssRF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){ + ssRF_argsint_varphi nargs; + int ret; + ssRF_argsint_theta* argument=(ssRF_argsint_theta*)args; + // tmps used in this function (not passed to integrand_varphi) + array_mpfr* tmps1; + array_mpfr* tmps2; + int thread_id; + + // the number of the thread + if(argument->thread_ids->length>1){ + thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids)); + if(thread_id<0){ + fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n"); + return(-1); + } + } + else{ + thread_id=0; + } + + // tmps used in this function (not passed to integrand_varphi) + tmps1=argument->tmpss.values+8*thread_id; + tmps2=argument->tmpss.values+8*thread_id+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.theta=theta; + + // copy tmps + // only pass tmps reserved for this thread + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+8*thread_id+2; + // 6 tmps in integrand_varphi + nargs.tmpss.length=6; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN); + mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_varphi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RF +int ssRF_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args){ + ssRF_argsint_r nargs; + int ret; + ssRF_argsint_varphi* argument=(ssRF_argsint_varphi*)args; + // tmps used in this function (not passed to integrand_r) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + // copy tmps + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 4 tmps in integrand_r + nargs.tmpss.length=4; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN); + mpfr_set_ui(tmps1->values[1], 1, MPFR_RNDN); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_r, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RF +int ssRF_integrand_r(mpfr_t* out, mpfr_t r, void* args){ + ssRF_argsint_rho nargs; + int ret; + ssRF_argsint_r* argument=(ssRF_argsint_r*)args; + // tmps used in this function (not passed to integrand_rho) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.cache=argument->cache; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + // copy tmps + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 2 tmps in integrand_rho + nargs.tmpss.length=2; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN); + ss_R(tmps1->values[1], nargs.varphi, 0, nargs.cache); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRF_integrand_rho, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RF +int ssRF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){ + int nu; + ssRF_argsint_rho* a=(ssRF_argsint_rho*)args; + mpfr_t* tmps; + // tmps used in this function + array_mpfr* tmps1=a->tmpss.values; + array_mpfr* tmps2=a->tmpss.values+1; + + // alloc tmps if needed + array_mpfr_alloc_tmps(11, tmps1); + + tmps=tmps1->values; + + mpfr_set_ui(*out, 0, MPFR_RNDN); + for(nu=-1;nu<=1;nu=nu+2){ + // q1 (store in tmps[2]) + mpfr_cos(tmps[2], a->theta, MPFR_RNDN); + mpfr_mul(tmps[2], tmps[2], a->r, MPFR_RNDN); + mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN); + // q2 (store in tmps[3]) + mpfr_sin(tmps[3], a->theta, MPFR_RNDN); + mpfr_mul(tmps[3], tmps[3], a->r, MPFR_RNDN); + mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN); + if(a->params.omega==-1){ + mpfr_neg(tmps[3], tmps[3], MPFR_RNDN); + } + // F1 (store in tmps[4]) + mpfr_cos(tmps[4], a->varphi, MPFR_RNDN); + mpfr_mul(tmps[4], tmps[4], rho, MPFR_RNDN); + mpfr_add(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN); + // F2 (store in tmps[5]) + mpfr_sin(tmps[5], a->varphi, MPFR_RNDN); + mpfr_mul(tmps[5], tmps[5], rho, MPFR_RNDN); + mpfr_add_ui(tmps[5], tmps[5], 1, MPFR_RNDN); + if(nu==-1){ + mpfr_neg(tmps[5], tmps[5], MPFR_RNDN); + } + // p1 (store in tmps[0]) + mpfr_sub(tmps[0], tmps[4], tmps[2], MPFR_RNDN); + mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN); + // p2 (store in tmps[1]) + mpfr_sub(tmps[1], tmps[5], tmps[3], MPFR_RNDN); + mpfr_add_si(tmps[1], tmps[1], a->params.omega, MPFR_RNDN); + + // cutoff + // tmps[6]=(1-ss_cutoff(ss_norm(F1-q1,F2-q2)))*ss_cutoff(r) + mpfr_sub(tmps[6], tmps[4], tmps[2], MPFR_RNDN); + mpfr_sub(tmps[7], tmps[5], tmps[3], MPFR_RNDN); + ss_norm(tmps[8], tmps[6], tmps[7], a->cache, tmps[9], tmps[10]); + ss_cutoff(tmps[6], tmps[8], tmps[7], tmps[9]); + mpfr_ui_sub(tmps[6], 1, tmps[6], MPFR_RNDN); + ss_cutoff(tmps[7], a->r, tmps[8], tmps[9]); + mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + + (*(a->I))(tmps[7], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache); + + // r*rho*I*cutoffs + mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + mpfr_mul(tmps[6], tmps[6], a->r, MPFR_RNDN); + mpfr_mul(tmps[6], tmps[6], rho, MPFR_RNDN); + + // add to *out + mpfr_add(*out, *out, tmps[6], MPFR_RNDN); + } + + return(0); +} + + +// compute the integral A_RR +int ssRR_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads){ + ssRR_argsint_theta args; + int ret; + unsigned int a; + mpfr_t lower, upper; + mpfr_t val; + array_pthread_t thread_ids; + + // init + thread_ids.length=0; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + args.cache=cache; + args.tmpss=tmpss; + args.thread_ids=&thread_ids; + + mpfr_inits(val, lower, upper, NULL); + + // split the integral in three parts + mpfr_set_ui(*out, 0, MPFR_RNDN); + for(a=0;a<3;a++){ + // bounds + mpfr_set_ui(lower, a, MPFR_RNDN); + mpfr_mul_ui(lower, lower, 4, MPFR_RNDN); + mpfr_sub_ui(lower, lower, 1, MPFR_RNDN); + mpfr_div_ui(lower, lower, 6, MPFR_RNDN); + mpfr_mul(lower, lower, args.cache.pi, MPFR_RNDN); + + mpfr_set_ui(upper, a, MPFR_RNDN); + mpfr_mul_ui(upper, upper, 4, MPFR_RNDN); + mpfr_add_ui(upper, upper, 3, MPFR_RNDN); + mpfr_div_ui(upper, upper, 6, MPFR_RNDN); + mpfr_mul(upper, upper, args.cache.pi, MPFR_RNDN); + + // save sector + args.sector_theta=a; + + if(threads==1){ + ret=integrate_gauss_mpfr(&val, &ssRR_integrand_theta, lower, upper, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_mpfr(&val, &ssRR_integrand_theta, lower, upper, abcissa, weights, &args, threads, args.thread_ids); + array_pthread_t_free(*(args.thread_ids)); + } + if(ret<0){ + break; + } + + mpfr_add(*out, *out, val, MPFR_RNDN); + } + // extra factor + mpfr_mul_ui(*out, *out, 3, MPFR_RNDN); + + mpfr_clears(val, lower, upper, NULL); + + return(ret); +} + +// integrand of the integral over theta in A_RR +int ssRR_integrand_theta(mpfr_t* out, mpfr_t theta, void* args){ + ssRR_argsint_varphi nargs; + int ret; + ssRR_argsint_theta* argument=(ssRR_argsint_theta*)args; + array_mpfr* tmps1; + array_mpfr* tmps2; + int thread_id; + + // the number of the thread + if(argument->thread_ids->length>1){ + thread_id=array_pthread_t_find(pthread_self(), *(argument->thread_ids)); + if(thread_id<0){ + fprintf(stderr, "error: thread not found\nIf you see this message, then you've uncovered a serious bug...\n"); + return(-1); + } + } + else{ + thread_id=0; + } + + // tmps used in this function (not passed to integrand_varphi) + tmps1=argument->tmpss.values+8*thread_id; + tmps2=argument->tmpss.values+8*thread_id+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.sector_theta=argument->sector_theta; + nargs.theta=theta; + + // copy tmps + // only pass tmps reserved for this thread + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+8*thread_id+2; + // 6 tmps in integrand_varphi + nargs.tmpss.length=6; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_div_si(tmps1->values[0], nargs.cache.pi, -6, MPFR_RNDN); + mpfr_div_ui(tmps1->values[1], nargs.cache.pi, 2, MPFR_RNDN); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_varphi, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RR +int ssRR_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args){ + ssRR_argsint_r nargs; + int ret; + ssRR_argsint_varphi* argument=(ssRR_argsint_varphi*)args; + // tmps used in this function (not passed to integrand_r) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.cache=argument->cache; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + // copy tmps + // only pass tmps reserved for this thread + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 4 tmps in integrand_r + nargs.tmpss.length=4; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN); + ss_R(tmps1->values[1], nargs.theta, argument->sector_theta, nargs.cache); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_r, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RR +int ssRR_integrand_r(mpfr_t* out, mpfr_t r, void* args){ + ssRR_argsint_rho nargs; + int ret; + ssRR_argsint_r* argument=(ssRR_argsint_r*)args; + // tmps used in this function (not passed to integrand_rho) + array_mpfr* tmps1=argument->tmpss.values; + array_mpfr* tmps2=argument->tmpss.values+1; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.cache=argument->cache; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + // copy tmps + // only pass tmps reserved for this thread + // warning: 'nargs.tmpss' must not be freed or resized + nargs.tmpss.values=argument->tmpss.values+2; + // 2 tmps in integrand_rho + nargs.tmpss.length=2; + + // alloc tmps if needed + array_mpfr_alloc_tmps(2, tmps1); + + // bounds + mpfr_set_ui(tmps1->values[0], 0, MPFR_RNDN); + ss_R(tmps1->values[1], nargs.varphi, 0, nargs.cache); + + ret=integrate_gauss_smarttmp_mpfr(out, &ssRR_integrand_rho, tmps1->values[0], tmps1->values[1], argument->abcissa, argument->weights, tmps2, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RR +int ssRR_integrand_rho(mpfr_t* out, mpfr_t rho, void* args){ + int eta, nu; + ssRR_argsint_rho* a=(ssRR_argsint_rho*)args; + mpfr_t* tmps; + // tmps used in this function (not passed to integrand_rho) + array_mpfr* tmps1=a->tmpss.values; + array_mpfr* tmps2=a->tmpss.values+1; + + // alloc tmps if needed + array_mpfr_alloc_tmps(12, tmps1); + + tmps=tmps1->values; + + mpfr_set_ui(*out, 0, MPFR_RNDN); + for(eta=-1;eta<=1;eta=eta+2){ + for(nu=-1;nu<=1;nu=nu+2){ + // q1 (store in tmps[2]) + mpfr_cos(tmps[2], a->theta, MPFR_RNDN); + mpfr_mul(tmps[2], tmps[2], a->r, MPFR_RNDN); + mpfr_add(tmps[2], tmps[2], a->cache.sqrt3, MPFR_RNDN); + // q2 (store in tmps[3]) + mpfr_sin(tmps[3], a->theta, MPFR_RNDN); + mpfr_mul(tmps[3], tmps[3], a->r, MPFR_RNDN); + mpfr_add_ui(tmps[3], tmps[3], 1, MPFR_RNDN); + if(eta==-1){ + mpfr_neg(tmps[3], tmps[3], MPFR_RNDN); + } + // F1 (store in tmps[4]) + mpfr_cos(tmps[4], a->varphi, MPFR_RNDN); + mpfr_mul(tmps[4], tmps[4], rho, MPFR_RNDN); + mpfr_add(tmps[4], tmps[4], a->cache.sqrt3, MPFR_RNDN); + // F2 (store in tmps[5]) + mpfr_sin(tmps[5], a->varphi, MPFR_RNDN); + mpfr_mul(tmps[5], tmps[5], rho, MPFR_RNDN); + mpfr_add_ui(tmps[5], tmps[5], 1, MPFR_RNDN); + if(nu==-1){ + mpfr_neg(tmps[5], tmps[5], MPFR_RNDN); + } + // p1 (store in tmps[0]) + mpfr_sub(tmps[0], tmps[4], tmps[2], MPFR_RNDN); + mpfr_add(tmps[0], tmps[0], a->cache.sqrt3, MPFR_RNDN); + // p2 (store in tmps[1]) + mpfr_sub(tmps[1], tmps[5], tmps[3], MPFR_RNDN); + mpfr_add_si(tmps[1], tmps[1], a->params.omega, MPFR_RNDN); + + // cutoff + // tmps[6]=(1-ss_cutoff(ss_norm(F1-q1,F2-q2)))*(1-ss_cutoff(ss_norm(q1-SQRT3,q2-omega))) + mpfr_sub(tmps[6], tmps[4], tmps[2], MPFR_RNDN); + mpfr_sub(tmps[7], tmps[5], tmps[3], MPFR_RNDN); + ss_norm(tmps[8], tmps[6], tmps[7], a->cache, tmps[9], tmps[10]); + ss_cutoff(tmps[6], tmps[8], tmps[7], tmps[9]); + mpfr_ui_sub(tmps[6], 1, tmps[6], MPFR_RNDN); + mpfr_sub(tmps[7], tmps[2], a->cache.sqrt3, MPFR_RNDN); + mpfr_sub_si(tmps[8], tmps[3], a->params.omega, MPFR_RNDN); + ss_norm(tmps[9], tmps[7], tmps[8], a->cache, tmps[10], tmps[11]); + ss_cutoff(tmps[7], tmps[9], tmps[8], tmps[10]); + mpfr_ui_sub(tmps[7], 1, tmps[7], MPFR_RNDN); + mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + + (*(a->I))(tmps[7], tmps[0], tmps[1], tmps[2], tmps[3], tmps[4], tmps[5], a->params.t1, a->params.t2, a->params.phi, a->params.W, tmps2, a->cache); + + // r*rho*I*cutoffs + mpfr_mul(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + mpfr_mul(tmps[6], tmps[6], a->r, MPFR_RNDN); + mpfr_mul(tmps[6], tmps[6], rho, MPFR_RNDN); + + // add to *out + mpfr_add(*out, *out, tmps[6], MPFR_RNDN); + } + } + + return(0); +} + + +// R(theta) +int ss_R(mpfr_t out, mpfr_t theta, unsigned int sector, struct ss_cache cache){ + mpfr_set_ui(out, sector, MPFR_RNDN); + mpfr_mul_ui(out, out, 8, MPFR_RNDN); + mpfr_sub_ui(out, out, 1, MPFR_RNDN); + mpfr_mul(out, out, cache.pi, MPFR_RNDN); + mpfr_div_ui(out, out, 6, MPFR_RNDN); + mpfr_add(out, out, theta, MPFR_RNDN); + mpfr_cos(out, out, MPFR_RNDN); + mpfr_ui_div(out, 1, out, MPFR_RNDN); + return(0); +} + +// cutoff function +int ss_cutoff(mpfr_t out, mpfr_t x, mpfr_t tmp1, mpfr_t tmp2){ + + if(mpfr_cmp_d(x,0.5)<=0){ + mpfr_set_ui(out, 1, MPFR_RNDN); + return(0); + } + if(mpfr_cmp_ui(x,1)>=0){ + mpfr_set_ui(out, 0, MPFR_RNDN); + return(0); + } + + mpfr_ui_sub(tmp1, 1, x, MPFR_RNDN); + mpfr_d_div(tmp1, -0.5, tmp1, MPFR_RNDN); + mpfr_exp(tmp1, tmp1, MPFR_RNDN); + + mpfr_sub_d(tmp2, x, 0.5, MPFR_RNDN); + mpfr_d_div(tmp2, -0.5, tmp2, MPFR_RNDN); + mpfr_exp(tmp2, tmp2, MPFR_RNDN); + + mpfr_ui_sub(out, 1, tmp2, MPFR_RNDN); + mpfr_mul(out, out, tmp1, MPFR_RNDN); + mpfr_add(out, out, tmp2, MPFR_RNDN); + mpfr_div(out, tmp1, out, MPFR_RNDN); + + return(0); +} + +// periodic norm +int ss_norm(mpfr_t out, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1, mpfr_t tmp2){ + // split R^2 into equilateral triangles and find which triangle (k1,k2) is in + // tmp1=floor(k1/sqrt3)+1 + // tmp2=floor(k1/2/sqrt3-k2/2)+1 + // out=floor(k1/2/sqrt3+k2/2)+1 + mpfr_div(tmp1, k1, cache.sqrt3, MPFR_RNDN); + mpfr_sub(tmp2, tmp1, k2, MPFR_RNDN); + mpfr_add(out, tmp1, k2, MPFR_RNDN); + mpfr_floor(tmp1, tmp1); + mpfr_add_ui(tmp1, tmp1, 1, MPFR_RNDN); + mpfr_div_ui(tmp2, tmp2, 2, MPFR_RNDN); + mpfr_floor(tmp2, tmp2); + mpfr_add_ui(tmp2, tmp2, 1, MPFR_RNDN); + mpfr_div_ui(out, out, 2, MPFR_RNDN); + mpfr_floor(out, out); + mpfr_add_ui(out, out, 1, MPFR_RNDN); + + // translate (k1,k2) by -a*(SQRT3,3)-b(SQRT3,-3) + // a-b + // tmp2=b-a=floor((tmp2-out+1)/3); + mpfr_sub(tmp2, tmp2, out, MPFR_RNDN); + mpfr_add_ui(tmp2, tmp2, 1, MPFR_RNDN); + mpfr_div_ui(tmp2, tmp2, 3, MPFR_RNDN); + mpfr_floor(tmp2, tmp2); + // tmp1=b=floor((tmp1+b-a)/2) + mpfr_add(tmp1, tmp1, tmp2, MPFR_RNDN); + mpfr_div_ui(tmp1, tmp1, 2, MPFR_RNDN); + mpfr_floor(tmp1, tmp1); + // tmp2=a + mpfr_sub(tmp2, tmp1, tmp2, MPFR_RNDN); + + // out=(k1-sqrt3*(a+b))**2 + mpfr_add(out, tmp1, tmp2, MPFR_RNDN); + mpfr_mul(out, out, cache.sqrt3, MPFR_RNDN); + mpfr_sub(out, k1, out, MPFR_RNDN); + mpfr_pow_ui(out, out, 2, MPFR_RNDN); + // tmp1=(k2-3*(a-b))**2 + mpfr_sub(tmp1, tmp2, tmp1, MPFR_RNDN); + mpfr_mul_ui(tmp1, tmp1, 3, MPFR_RNDN); + mpfr_sub(tmp1, k2, tmp1, MPFR_RNDN); + mpfr_pow_ui(tmp1, tmp1, 2, MPFR_RNDN); + + // out=sqrt((k1-sqrt3*(a+b))**2+(k2-3*(a-b))**2) + mpfr_add(out, out, tmp1, MPFR_RNDN); + mpfr_sqrt(out, out, MPFR_RNDN); + + return(0); +} diff --git a/src/ss_integral.h b/src/ss_integral.h new file mode 100644 index 0000000..fb30ddd --- /dev/null +++ b/src/ss_integral.h @@ -0,0 +1,200 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + Computation of the sunrise diagram +*/ + +#ifndef SS_INTEGRAL_H +#define SS_INTEGRAL_H + +#include <libinum.h> +#include "types.h" + +// compute pi and sqrt3 ahead of time +struct ss_cache { + mpfr_t pi; + mpfr_t sqrt3; +}; + +// extra arguments for the integrations +// A_FF +typedef struct ssFF_argsint_rho { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + array_pthread_t* thread_ids; +} ssFF_argsint_rho; +typedef struct ssFF_argsint_theta { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr rho; +} ssFF_argsint_theta; +typedef struct ssFF_argsint_psi { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr rho; + mpfr_ptr theta; +} ssFF_argsint_psi; +typedef struct ssFF_argsint_z { + hh_params params; + TYPE_I; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr rho; + mpfr_ptr theta; + mpfr_ptr psi; +} ssFF_argsint_z; + +// A_RF +typedef struct ssRF_argsint_theta { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + array_pthread_t* thread_ids; +} ssRF_argsint_theta; +typedef struct ssRF_argsint_varphi { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr theta; +} ssRF_argsint_varphi; +typedef struct ssRF_argsint_r { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr theta; + mpfr_ptr varphi; +} ssRF_argsint_r; +typedef struct ssRF_argsint_rho { + hh_params params; + TYPE_I; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr theta; + mpfr_ptr varphi; + mpfr_ptr r; +} ssRF_argsint_rho; + +// A_RR +typedef struct ssRR_argsint_theta { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + array_pthread_t* thread_ids; + int sector_theta; +} ssRR_argsint_theta; +typedef struct ssRR_argsint_varphi { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + int sector_theta; + mpfr_ptr theta; +} ssRR_argsint_varphi; +typedef struct ssRR_argsint_r { + hh_params params; + TYPE_I; + array_mpfr abcissa; + array_mpfr weights; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr theta; + mpfr_ptr varphi; +} ssRR_argsint_r; +typedef struct ssRR_argsint_rho { + hh_params params; + TYPE_I; + struct ss_cache cache; + array_2_mpfr tmpss; + mpfr_ptr theta; + mpfr_ptr varphi; + mpfr_ptr r; +} ssRR_argsint_rho; + +// compute the integral +int ss_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, unsigned int threads); + +// compute the integral A_FF +int ssFF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss , unsigned int threads); +// integrand of the integral over rho in A_FF +int ssFF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args); +// integrand of the integral over theta in A_FF +int ssFF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args); +// integrand of the integral over psi in A_FF +int ssFF_integrand_psi(mpfr_t* out, mpfr_t psi, void* args); +// integrand of the integral over z in A_FF +int ssFF_integrand_z(mpfr_t* out, mpfr_t z, void* args); + +// compute the integral A_RF +int ssRF_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads); +// integrand of the integral over theta in A_RF +int ssRF_integrand_theta(mpfr_t* out, mpfr_t theta, void* args); +// integrand of the integral over varphi in A_RF +int ssRF_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args); +// integrand of the integral over r in A_RF +int ssRF_integrand_r(mpfr_t* out, mpfr_t r, void* args); +// integrand of the integral over rho in A_RF +int ssRF_integrand_rho(mpfr_t* out, mpfr_t rho, void* args); + +// compute the integral A_RR +int ssRR_integrate(mpfr_t* out, TYPE_I, hh_params params, array_mpfr abcissa, array_mpfr weights, struct ss_cache cache, array_2_mpfr tmpss, unsigned int threads); +// integrand of the integral over theta in A_RR +int ssRR_integrand_theta(mpfr_t* out, mpfr_t theta, void* args); +// integrand of the integral over varphi in A_RR +int ssRR_integrand_varphi(mpfr_t* out, mpfr_t varphi, void* args); +// integrand of the integral over r in A_RR +int ssRR_integrand_r(mpfr_t* out, mpfr_t r, void* args); +// integrand of the integral over rho in A_RR +int ssRR_integrand_rho(mpfr_t* out, mpfr_t rho, void* args); + +// R(theta) +int ss_R(mpfr_t out, mpfr_t theta, unsigned int sector, struct ss_cache cache); + +// cutoff function +int ss_cutoff(mpfr_t out, mpfr_t x, mpfr_t tmp1, mpfr_t tmp2); + +// periodic norm +int ss_norm(mpfr_t out, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1, mpfr_t tmp2); + +#endif + + diff --git a/src/ss_integral_double.c b/src/ss_integral_double.c new file mode 100644 index 0000000..00ab637 --- /dev/null +++ b/src/ss_integral_double.c @@ -0,0 +1,437 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "ss_integral_double.h" +#include <math.h> + +#define PI 3.1415926535897932385L +#define SQRT3 1.7320508075688772935L + +// compute the integral +int ss_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + long double val; + + // compute A_FF, A_RF and A_RR + ssFF_integrate_double(out, I, params, abcissa, weights, threads); + ssRF_integrate_double(&val, I, params, abcissa, weights, threads); + *out+=2*val; + ssRR_integrate_double(&val, I, params, abcissa, weights, threads); + *out+=val; + + return(0); +} + + +// compute the integral A_FF +int ssFF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssFF_argsint_rho_double args; + int ret; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + if(threads==1){ + ret=integrate_gauss_ldouble(out, &ssFF_integrand_rho_double, 0., 1., abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(out, &ssFF_integrand_rho_double, 0., 1., abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + + // extra factor + *out*=6; + + return(ret); +} + +// integrand of the integral over rho in A_FF +int ssFF_integrand_rho_double(long double* out, long double rho, void* args){ + ssFF_argsint_theta_double nargs; + int ret; + ssFF_argsint_rho_double* argument=(ssFF_argsint_rho_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.rho=rho; + + ret=integrate_gauss_ldouble(out, &ssFF_integrand_theta_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over theta in A_FF +int ssFF_integrand_theta_double(long double* out, long double theta, void* args){ + ssFF_argsint_psi_double nargs; + int ret; + ssFF_argsint_theta_double* argument=(ssFF_argsint_theta_double*)args; + long double val; + int a; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.rho=argument->rho; + nargs.theta=theta; + + // split the integral into four parts to improve the precision + *out=0; + for(a=0;a<4;a++){ + ret=integrate_gauss_ldouble(&val, &ssFF_integrand_psi_double, PI/4*(a-2), PI/4*(a-1), argument->abcissa, argument->weights, &nargs); + if(ret<0){ + break; + } + *out+=val; + } + + return(ret); +} + +// integrand of the integral over psi in A_FF +int ssFF_integrand_psi_double(long double* out, long double psi, void* args){ + ssFF_argsint_z_double nargs; + int ret; + ssFF_argsint_psi_double* argument=(ssFF_argsint_psi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.rho=argument->rho; + nargs.theta=argument->theta; + nargs.psi=psi; + + ret=integrate_gauss_ldouble(out, &ssFF_integrand_z_double, 0., 1., argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over z in A_FF +int ssFF_integrand_z_double(long double* out, long double z, void* args){ + ssFF_argsint_z_double* argument=(ssFF_argsint_z_double*)args; + long double p1,p2,q1,q2,F1,F2; + + // r and varphi (store in p1 and p2) + p1=(1-z)*(1+z*(1-cosl(argument->psi))/(1+cosl(argument->psi))); + p2=acosl(cos(2*argument->psi)-p1/2*sin(2*argument->psi)*sin(2*argument->psi)); + if(argument->psi<0.){ + p2=-p2; + } + + // jacobian + *out=4*argument->rho*argument->rho*argument->rho*p1*(1+p1*cosl(2*argument->psi))/(1+cosl(argument->psi))/sqrtl(1+p1*cosl(argument->psi)*cosl(argument->psi)); + + // cutoffs + *out*=ss_cutoff_double(argument->rho)*ss_cutoff_double(argument->rho*p1); + + q1=SQRT3+argument->rho*p1*cosl(argument->theta+p2); + q2=argument->params.omega*(1+argument->rho*p1*sinl(argument->theta+p2)); + p1=SQRT3+argument->rho*cosl(argument->theta); + p2=argument->params.omega*(1+argument->rho*sinl(argument->theta)); + F1=p1+q1-SQRT3; + F2=p2+q2-argument->params.omega; + + *out*=(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W); + + return(0); +} + + +// compute the integral A_RF +int ssRF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssRF_argsint_theta_double args; + int ret; + long double val; + int a; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + // split the integral in three parts to improve the precision + *out=0; + for(a=0;a<3;a++){ + if(threads==1){ + ret=integrate_gauss_ldouble(&val, &ssRF_integrand_theta_double, (4*a-1)*PI/6, (4*a+3)*PI/6, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(&val, &ssRF_integrand_theta_double, (4*a-1)*PI/6, (4*a+3)*PI/6, abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + if(ret<0){ + break; + } + + *out+=val; + } + + // extra factor + *out*=3; + + return(ret); +} + +// integrand of the integral over theta in A_RF +int ssRF_integrand_theta_double(long double* out, long double theta, void* args){ + ssRF_argsint_varphi_double nargs; + int ret; + ssRF_argsint_theta_double* argument=(ssRF_argsint_theta_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_varphi_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RF +int ssRF_integrand_varphi_double(long double* out, long double varphi, void* args){ + ssRF_argsint_r_double nargs; + int ret; + ssRF_argsint_varphi_double* argument=(ssRF_argsint_varphi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_r_double, 0., 1., argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RF +int ssRF_integrand_r_double(long double* out, long double r, void* args){ + ssRF_argsint_rho_double nargs; + int ret; + ssRF_argsint_r_double* argument=(ssRF_argsint_r_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + ret=integrate_gauss_ldouble(out, &ssRF_integrand_rho_double, 0., ss_R_double(nargs.varphi, 0), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RF +int ssRF_integrand_rho_double(long double* out, long double rho, void* args){ + int nu; + ssRF_argsint_rho_double* argument=(ssRF_argsint_rho_double*)args; + long double p1,p2,q1,q2,F1,F2; + + *out=0.; + for(nu=-1;nu<=1;nu=nu+2){ + q1=SQRT3+argument->r*cosl(argument->theta); + q2=argument->params.omega*(1+argument->r*sinl(argument->theta)); + F1=SQRT3+rho*cosl(argument->varphi); + F2=nu*(1+rho*sinl(argument->varphi)); + p1=SQRT3+F1-q1; + p2=argument->params.omega+F2-q2; + + *out+=argument->r*rho*(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W)*(1-ss_cutoff_double(ss_norm_double(F1-q1,F2-q2)))*ss_cutoff_double(argument->r); + } + + return(0); +} + + +// compute the integral A_RR +int ssRR_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads){ + ssRR_argsint_theta_double args; + int ret; + int a; + long double lower, upper; + long double val; + array_pthread_t thread_ids; + + args.params=params; + args.I=I; + args.abcissa=abcissa; + args.weights=weights; + + // split the integral in three parts + *out=0; + for(a=0;a<3;a++){ + lower=(4*a-1)*PI/6; + upper=(4*a+3)*PI/6; + + // save sector + args.sector_theta=a; + + if(threads==1){ + ret=integrate_gauss_ldouble(&val, &ssRR_integrand_theta_double, lower, upper, abcissa, weights, &args); + } + else{ + ret=integrate_gauss_multithread_ldouble(&val, &ssRR_integrand_theta_double, lower, upper, abcissa, weights, &args, threads, &thread_ids); + array_pthread_t_free(thread_ids); + } + if(ret<0){ + break; + } + + *out+=val; + } + // extra factor + *out*=3; + + return(ret); +} + +// integrand of the integral over theta in A_RR +int ssRR_integrand_theta_double(long double* out, long double theta, void* args){ + ssRR_argsint_varphi_double nargs; + int ret; + ssRR_argsint_theta_double* argument=(ssRR_argsint_theta_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.sector_theta=argument->sector_theta; + nargs.theta=theta; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_varphi_double, -PI/6, PI/2, argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over varphi in A_RR +int ssRR_integrand_varphi_double(long double* out, long double varphi, void* args){ + ssRR_argsint_r_double nargs; + int ret; + ssRR_argsint_varphi_double* argument=(ssRR_argsint_varphi_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.abcissa=argument->abcissa; + nargs.weights=argument->weights; + nargs.theta=argument->theta; + nargs.varphi=varphi; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_r_double, 0., ss_R_double(nargs.theta, argument->sector_theta), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over r in A_RR +int ssRR_integrand_r_double(long double* out, long double r, void* args){ + ssRR_argsint_rho_double nargs; + int ret; + ssRR_argsint_r_double* argument=(ssRR_argsint_r_double*)args; + + // new argument + nargs.params=argument->params; + nargs.I=argument->I; + nargs.theta=argument->theta; + nargs.varphi=argument->varphi; + nargs.r=r; + + ret=integrate_gauss_ldouble(out, &ssRR_integrand_rho_double, 0., ss_R_double(nargs.varphi, 0), argument->abcissa, argument->weights, &nargs); + + return(ret); +} + +// integrand of the integral over rho in A_RR +int ssRR_integrand_rho_double(long double* out, long double rho, void* args){ + int eta, nu; + ssRR_argsint_rho_double* argument=(ssRR_argsint_rho_double*)args; + long double p1,p2,q1,q2,F1,F2; + + *out=0.; + for(eta=-1;eta<=1;eta=eta+2){ + for(nu=-1;nu<=1;nu=nu+2){ + q1=SQRT3+argument->r*cosl(argument->theta); + q2=eta*(1+argument->r*sinl(argument->theta)); + F1=SQRT3+rho*cosl(argument->varphi); + F2=nu*(1+rho*sinl(argument->varphi)); + p1=SQRT3+F1-q1; + p2=argument->params.omega+F2-q2; + + *out+=argument->r*rho*(*(argument->I))(p1, p2, q1, q2, F1, F2, argument->params.t1, argument->params.t2, argument->params.phi, argument->params.W)*(1-ss_cutoff_double(ss_norm_double(F1-q1,F2-q2)))*(1-ss_cutoff_double(ss_norm_double(q1-SQRT3,q2-argument->params.omega))); + } + } + + return(0); +} + + +// R(theta) +long double ss_R_double(long double theta, int sector){ + return(1./cosl(theta+(8*sector-1)*PI/6)); +} + +// cutoff function +long double ss_cutoff_double(long double x){ + long double e1,e2; + + if(x<=0.5){ + return(1.); + } + if(x>=1){ + return(0.); + } + + e1=expl(-0.5/(1-x)); + e2=expl(-0.5/(x-0.5)); + return(e1/(e2+e1*(1-e2))); +} + +// periodic norm +long double ss_norm_double(long double k1, long double k2){ + long double n1, n2, n3; + long double tmp1, tmp2; + long double t1, t2; + + // split R^2 into equilateral triangles and find which triangle (k1,k2) is in + n1=floorl(k1/SQRT3)+1; + n2=floorl(k1/2/SQRT3-k2/2)+1; + n3=floorl(k1/2/SQRT3+k2/2)+1; + + // translate (k1,k2) by -a*(SQRT3,3)-b(SQRT3,-3) + // a-b + tmp1=-floorl((n2-n3+1)/3); + // b + tmp2=floorl((n1-tmp1)/2); + // a + tmp1=tmp1+tmp2; + t1=k1-SQRT3*(tmp1+tmp2); + t2=k2-3*(tmp1-tmp2); + + return(sqrtl(t1*t1+t2*t2)); +} diff --git a/src/ss_integral_double.h b/src/ss_integral_double.h new file mode 100644 index 0000000..67ec1de --- /dev/null +++ b/src/ss_integral_double.h @@ -0,0 +1,170 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + Computation of the sunrise diagram + + Takes a pointer to an integrand function 'I' as an argument. 'I' must have the following arguments + long double I(long double p1, long double p2, long double q1, long double q2, long double K1, long double K2, long double t1, long double t2, long double phi, long double W) +*/ + +#ifndef SS_INTEGRAL_DOUBLE_H +#define SS_INTEGRAL_DOUBLE_H + +#include <libinum.h> +#include "double_util.h" + +// extra arguments for the integrations +// A_FF +typedef struct ssFF_argsint_rho_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; +} ssFF_argsint_rho_double; +typedef struct ssFF_argsint_theta_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + long double rho; +} ssFF_argsint_theta_double; +typedef struct ssFF_argsint_psi_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + long double rho; + long double theta; +} ssFF_argsint_psi_double; +typedef struct ssFF_argsint_z_double { + hh_params_double params; + TYPE_I_DOUBLE; + long double rho; + long double theta; + long double psi; +} ssFF_argsint_z_double; + +// A_RF +typedef struct ssRF_argsint_theta_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; +} ssRF_argsint_theta_double; +typedef struct ssRF_argsint_varphi_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + long double theta; +} ssRF_argsint_varphi_double; +typedef struct ssRF_argsint_r_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + long double theta; + long double varphi; +} ssRF_argsint_r_double; +typedef struct ssRF_argsint_rho_double { + hh_params_double params; + TYPE_I_DOUBLE; + long double theta; + long double varphi; + long double r; +} ssRF_argsint_rho_double; + +// A_RR +typedef struct ssRR_argsint_theta_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + int sector_theta; +} ssRR_argsint_theta_double; +typedef struct ssRR_argsint_varphi_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + int sector_theta; + long double theta; +} ssRR_argsint_varphi_double; +typedef struct ssRR_argsint_r_double { + hh_params_double params; + TYPE_I_DOUBLE; + array_ldouble abcissa; + array_ldouble weights; + long double theta; + long double varphi; +} ssRR_argsint_r_double; +typedef struct ssRR_argsint_rho_double { + hh_params_double params; + TYPE_I_DOUBLE; + long double theta; + long double varphi; + long double r; +} ssRR_argsint_rho_double; + + +// compute the integral +int ss_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads); + +// compute the integral A_FF +int ssFF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads); +// integrand of the integral over rho in A_FF +int ssFF_integrand_rho_double(long double* out, long double rho, void* args); +// integrand of the integral over theta in A_FF +int ssFF_integrand_theta_double(long double* out, long double theta, void* args); +// integrand of the integral over psi in A_FF +int ssFF_integrand_psi_double(long double* out, long double psi, void* args); +// integrand of the integral over z in A_FF +int ssFF_integrand_z_double(long double* out, long double z, void* args); + +// compute the integral A_RF +int ssRF_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads); +// integrand of the integral over theta in A_RF +int ssRF_integrand_theta_double(long double* out, long double theta, void* args); +// integrand of the integral over varphi in A_RF +int ssRF_integrand_varphi_double(long double* out, long double varphi, void* args); +// integrand of the integral over r in A_RF +int ssRF_integrand_r_double(long double* out, long double r, void* args); +// integrand of the integral over rho in A_RF +int ssRF_integrand_rho_double(long double* out, long double rho, void* args); + +// compute the integral A_RR +int ssRR_integrate_double(long double* out, TYPE_I_DOUBLE, hh_params_double params, array_ldouble abcissa, array_ldouble weights, unsigned int threads); +// integrand of the integral over theta in A_RR +int ssRR_integrand_theta_double(long double* out, long double theta, void* args); +// integrand of the integral over varphi in A_RR +int ssRR_integrand_varphi_double(long double* out, long double varphi, void* args); +// integrand of the integral over r in A_RR +int ssRR_integrand_r_double(long double* out, long double r, void* args); +// integrand of the integral over rho in A_RR +int ssRR_integrand_rho_double(long double* out, long double rho, void* args); + +// R(theta; +long double ss_R_double(long double theta, int sector); + +// cutoff function +long double ss_cutoff_double(long double x); + +// periodic norm +long double ss_norm_double(long double k1, long double k2); + +#endif + diff --git a/src/types.h b/src/types.h new file mode 100644 index 0000000..0dec1fc --- /dev/null +++ b/src/types.h @@ -0,0 +1,35 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#ifndef TYPES_H +#define TYPES_H + +#include <mpfr.h> + +typedef struct hh_params { + int omega; + mpfr_t t1; + mpfr_t t2; + mpfr_t lambda; + mpfr_t sinphi; + mpfr_t phi; + mpfr_t W; +} hh_params; + +// format for I function (used to compute sunrise diagrams) +#define TYPE_I int (*I)(mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, mpfr_ptr, array_mpfr*, struct ss_cache) + +#endif diff --git a/src/zz_integral.c b/src/zz_integral.c new file mode 100644 index 0000000..ee42355 --- /dev/null +++ b/src/zz_integral.c @@ -0,0 +1,237 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "zz_integral.h" + + +// I for z1-z2 +int zz_I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps1, struct ss_cache cache){ + mpfr_t* tmps; + + // alloc tmps if needed + array_mpfr_alloc_tmps(9, tmps1); + + tmps=tmps1->values; + + // tmps[0]=alpha1(p1,p2) + zz_alpha1(tmps[0], p1, p2, cache, tmps[1]); + // tmps[1]=m(p1,p2) + zz_m(tmps[1], p1, p2, t2, phi, W, cache, tmps[2]); + // tmps[2]=alpha1(q1,q2) + zz_alpha1(tmps[2], q1, q2, cache, tmps[3]); + // tmps[3]=m(q1,q2) + zz_m(tmps[3], q1, q2, t2, phi, W, cache, tmps[4]); + // tmps[4]=alpha1(F1,F2) + zz_alpha1(tmps[4], F1, F2, cache, tmps[5]); + // tmps[5]=m(F1,F2) + zz_m(tmps[5], F1, F2, t2, phi, W, cache, tmps[6]); + + // tmps[6]=zeta(p) + zz_zeta(tmps[6], tmps[0], t2, phi); + // tmps[7]=zeta(q) + zz_zeta(tmps[7], tmps[2], t2, phi); + // tmps[8]=zeta(F) + zz_zeta(tmps[8], tmps[4], t2, phi); + // tmps[6]=Z=zeta(p)+zeta(q)-zeta(F) + mpfr_add(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + mpfr_sub(tmps[6], tmps[6], tmps[8], MPFR_RNDN); + + // tmps[0]=xi(p) + zz_xi(tmps[0], tmps[1], tmps[0], t1, tmps[7]); + // tmps[2]=xi(q) + zz_xi(tmps[2], tmps[3], tmps[2], t1, tmps[7]); + // tmps[4]=xi(F) + zz_xi(tmps[4], tmps[5], tmps[4], t1, tmps[7]); + + // tmps[1]=m(p)/xi(p) + zz_mxi(tmps[1], tmps[1], tmps[0]); + // tmps[3]=m(q)/xi(q) + zz_mxi(tmps[3], tmps[3], tmps[2]); + // tmps[5]=m(F)/xi(F) + zz_mxi(tmps[5], tmps[5], tmps[4]); + + // I=(xip+xiq+xiF)*(mp/xip+mq/xiq-mF/xiF-mp*mq*mF/xip/xiq/xiF)*Z/(Z**2-(xip+xiq+xiF)**2)**2/108 + mpfr_add(tmps[0], tmps[0], tmps[2], MPFR_RNDN); + mpfr_add(tmps[0], tmps[0], tmps[4], MPFR_RNDN); + mpfr_add(tmps[2], tmps[1], tmps[3], MPFR_RNDN); + mpfr_sub(tmps[2], tmps[2], tmps[5], MPFR_RNDN); + mpfr_mul(tmps[1], tmps[1], tmps[3], MPFR_RNDN); + mpfr_mul(tmps[1], tmps[1], tmps[5], MPFR_RNDN); + mpfr_sub(tmps[1], tmps[2], tmps[1], MPFR_RNDN); + mpfr_mul(tmps[1], tmps[0], tmps[1], MPFR_RNDN); + mpfr_mul(tmps[1], tmps[1], tmps[6], MPFR_RNDN); + mpfr_pow_ui(tmps[6], tmps[6], 2, MPFR_RNDN); + mpfr_pow_ui(tmps[0], tmps[0], 2, MPFR_RNDN); + mpfr_sub(tmps[0], tmps[6], tmps[0], MPFR_RNDN); + mpfr_pow_ui(tmps[0], tmps[0], 2, MPFR_RNDN); + mpfr_div(out, tmps[1], tmps[0], MPFR_RNDN); + mpfr_div_ui(out, out, 108, MPFR_RNDN); + + return(0); +} + + +// I for z1+z2 +int ZZ_I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps1, struct ss_cache cache){ + mpfr_t* tmps; + + // alloc tmps if needed + array_mpfr_alloc_tmps(9, tmps1); + + tmps=tmps1->values; + + // tmps[0]=alpha1(p1,p2) + zz_alpha1(tmps[0], p1, p2, cache, tmps[1]); + // tmps[1]=m(p1,p2) + zz_m(tmps[1], p1, p2, t2, phi, W, cache, tmps[2]); + // tmps[2]=alpha1(q1,q2) + zz_alpha1(tmps[2], q1, q2, cache, tmps[3]); + // tmps[3]=m(q1,q2) + zz_m(tmps[3], q1, q2, t2, phi, W, cache, tmps[4]); + // tmps[4]=alpha1(F1,F2) + zz_alpha1(tmps[4], F1, F2, cache, tmps[5]); + // tmps[5]=m(F1,F2) + zz_m(tmps[5], F1, F2, t2, phi, W, cache, tmps[6]); + + // tmps[6]=zeta(p) + zz_zeta(tmps[6], tmps[0], t2, phi); + // tmps[7]=zeta(q) + zz_zeta(tmps[7], tmps[2], t2, phi); + // tmps[8]=zeta(F) + zz_zeta(tmps[8], tmps[4], t2, phi); + // tmps[6]=Z**2=(zeta(p)+zeta(q)-zeta(F))**2 + mpfr_add(tmps[6], tmps[6], tmps[7], MPFR_RNDN); + mpfr_sub(tmps[6], tmps[6], tmps[8], MPFR_RNDN); + mpfr_pow_ui(tmps[6], tmps[6], 2, MPFR_RNDN); + + // tmps[0]=xi(p) + zz_xi(tmps[0], tmps[1], tmps[0], t1, tmps[7]); + // tmps[2]=xi(q) + zz_xi(tmps[2], tmps[3], tmps[2], t1, tmps[7]); + // tmps[4]=xi(F) + zz_xi(tmps[4], tmps[5], tmps[4], t1, tmps[7]); + + // tmps[1]=m(p)/xi(p) + zz_mxi(tmps[1], tmps[1], tmps[0]); + // tmps[3]=m(q)/xi(q) + zz_mxi(tmps[3], tmps[3], tmps[2]); + // tmps[5]=m(F)/xi(F) + zz_mxi(tmps[5], tmps[5], tmps[4]); + + // I=(1-mp*mF/xip/xiF-mq*mF/xiq/xiF+mp*mq/xip/xiq)*(Z**2+(xip+xiq+xiF)**2)/(Z**2-(xip+xiq+xiF)**2)**2/216 + mpfr_add(tmps[0], tmps[0], tmps[2], MPFR_RNDN); + mpfr_add(tmps[0], tmps[0], tmps[4], MPFR_RNDN); + mpfr_pow_ui(tmps[0], tmps[0], 2, MPFR_RNDN); + mpfr_mul(tmps[2], tmps[1], tmps[3], MPFR_RNDN); + mpfr_mul(tmps[4], tmps[3], tmps[5], MPFR_RNDN); + mpfr_sub(tmps[2], tmps[2], tmps[4], MPFR_RNDN); + mpfr_mul(tmps[4], tmps[1], tmps[5], MPFR_RNDN); + mpfr_sub(tmps[2], tmps[2], tmps[4], MPFR_RNDN); + mpfr_add_ui(out, tmps[2], 1, MPFR_RNDN); + mpfr_add(tmps[1], tmps[6], tmps[0], MPFR_RNDN); + mpfr_mul(out, out, tmps[1], MPFR_RNDN); + mpfr_sub(tmps[0], tmps[6], tmps[0], MPFR_RNDN); + mpfr_pow_ui(tmps[0], tmps[0], 2, MPFR_RNDN); + mpfr_div(out, out, tmps[0], MPFR_RNDN); + mpfr_div_ui(out, out, 216, MPFR_RNDN); + + return(0); +} + +// zeta(k1,k2) +// requires m and alpha1 to be computed beforehand +int zz_zeta(mpfr_t zeta, mpfr_t alpha1, mpfr_t t2, mpfr_t phi){ + mpfr_cos(zeta, phi, MPFR_RNDN); + mpfr_mul(zeta, zeta, t2, MPFR_RNDN); + mpfr_mul_ui(zeta, zeta, 2, MPFR_RNDN); + mpfr_mul(zeta, zeta, alpha1, MPFR_RNDN); + return(0); +} +// xi(k1,k2) +// requires m and alpha1 to be computed beforehand +// requires one initialized tmp number +int zz_xi(mpfr_t xi, mpfr_t m, mpfr_t alpha1, mpfr_t t1, mpfr_t tmp){ + mpfr_pow_ui(tmp, m, 2, MPFR_RNDN); + + mpfr_mul(xi, alpha1, t1, MPFR_RNDN); + mpfr_mul(xi, xi, t1, MPFR_RNDN); + mpfr_mul_ui(xi, xi, 2, MPFR_RNDN); + + mpfr_add(xi, xi, tmp, MPFR_RNDN); + mpfr_sqrt(xi, xi, MPFR_RNDN); + return(0); +} +// m(k1,k2) +// requires two initialized tmp numbers +int zz_m(mpfr_t m, mpfr_t k1, mpfr_t k2, mpfr_t t2, mpfr_t phi, mpfr_t W, struct ss_cache cache, mpfr_t tmp1){ + zz_alpha2(m, k1, k2, cache, tmp1); + mpfr_sin(tmp1, phi, MPFR_RNDN); + mpfr_mul(m, m, tmp1, MPFR_RNDN); + mpfr_mul(m, m, t2, MPFR_RNDN); + mpfr_mul_ui(m, m, 2, MPFR_RNDN); + mpfr_sub(m, W, m, MPFR_RNDN); + return(0); +} +// m(k1,k2)/xi(k1,k2) +int zz_mxi(mpfr_t mxi, mpfr_t m, mpfr_t xi){ + // if xi=0, then return 0 (m/xi->0 at pF) + if(mpfr_cmp_ui(xi,0)==0){ + mpfr_set_ui(mxi, 0, MPFR_RNDN); + return(0); + } + mpfr_div(mxi, m, xi, MPFR_RNDN); + + return(0); +} + +// alpha1(k1,k2) +// requires one initialized tmp number +int zz_alpha1(mpfr_t alpha1, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1){ + mpfr_mul(tmp1, k2, cache.pi, MPFR_RNDN); + mpfr_div_ui(tmp1, tmp1, 3, MPFR_RNDN); + mpfr_cos(tmp1, tmp1, MPFR_RNDN); + + mpfr_mul(alpha1, k1, cache.pi, MPFR_RNDN); + mpfr_div(alpha1, alpha1, cache.sqrt3, MPFR_RNDN); + mpfr_cos(alpha1, alpha1, MPFR_RNDN); + + mpfr_add(alpha1, alpha1, tmp1, MPFR_RNDN); + mpfr_mul(alpha1, alpha1, tmp1, MPFR_RNDN); + mpfr_mul_ui(alpha1, alpha1, 2, MPFR_RNDN); + mpfr_add_d(alpha1, alpha1, 0.5, MPFR_RNDN); + return(0); +} +// alpha2(k1,k2) +// requires one initialized tmp numbers +int zz_alpha2(mpfr_t alpha2, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1){ + mpfr_mul(tmp1, k2, cache.pi, MPFR_RNDN); + mpfr_div_ui(tmp1, tmp1, 3, MPFR_RNDN); + mpfr_cos(tmp1, tmp1, MPFR_RNDN); + + mpfr_mul(alpha2, k1, cache.pi, MPFR_RNDN); + mpfr_div(alpha2, alpha2, cache.sqrt3, MPFR_RNDN); + mpfr_cos(alpha2, alpha2, MPFR_RNDN); + + mpfr_sub(alpha2, alpha2, tmp1, MPFR_RNDN); + + mpfr_mul(tmp1, k2, cache.pi, MPFR_RNDN); + mpfr_div_ui(tmp1, tmp1, 3, MPFR_RNDN); + mpfr_sin(tmp1, tmp1, MPFR_RNDN); + + mpfr_mul(alpha2, alpha2, tmp1, MPFR_RNDN); + mpfr_mul_ui(alpha2, alpha2, 2, MPFR_RNDN); + return(0); +} diff --git a/src/zz_integral.h b/src/zz_integral.h new file mode 100644 index 0000000..b23c88a --- /dev/null +++ b/src/zz_integral.h @@ -0,0 +1,53 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + Computation of z1-z2 and z1+z2 +*/ + +#ifndef ZZ_INTEGRAL_H +#define ZZ_INTEGRAL_H + +#include "ss_integral.h" +#include <libinum.h> +#include <mpfr.h> + +// I for z1-z2 +int zz_I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps1, struct ss_cache cache); + +// I for z1+z2 +int ZZ_I(mpfr_t out, mpfr_t p1, mpfr_t p2, mpfr_t q1, mpfr_t q2, mpfr_t F1, mpfr_t F2, mpfr_t t1, mpfr_t t2, mpfr_t phi, mpfr_t W, array_mpfr* tmps1, struct ss_cache cache); + +// zeta(k1,k2) +int zz_zeta(mpfr_t zeta, mpfr_t alpha1, mpfr_t t2, mpfr_t phi); + +// xi(k1,k2) +int zz_xi(mpfr_t xi, mpfr_t m, mpfr_t alpha1, mpfr_t t1, mpfr_t tmp); + +// m(k1,k2) +int zz_m(mpfr_t m, mpfr_t k1, mpfr_t k2, mpfr_t t2, mpfr_t phi, mpfr_t W, struct ss_cache cache, mpfr_t tmp1); + +// m(k1,k2)/xi(k1,k2) +int zz_mxi(mpfr_t mxi, mpfr_t m, mpfr_t xi); + +// alpha1(k1,k2) +int zz_alpha1(mpfr_t alpha1, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1); + +// alpha2(k1,k2) +int zz_alpha2(mpfr_t alpha2, mpfr_t k1, mpfr_t k2, struct ss_cache cache, mpfr_t tmp1); + +#endif + diff --git a/src/zz_integral_double.c b/src/zz_integral_double.c new file mode 100644 index 0000000..9330ad1 --- /dev/null +++ b/src/zz_integral_double.c @@ -0,0 +1,89 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +#include "zz_integral_double.h" +#include <math.h> + +#define PI 3.1415926535897932385L +#define SQRT3 1.7320508075688772935L + +// I for z1-z2 +long double zz_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){ + long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF; + + // alpha1 (store in xi) + xip=zz_alpha1_double(p1,p2); + xiq=zz_alpha1_double(q1,q2); + xiF=zz_alpha1_double(F1,F2); + + // alpha2 (store in m) + mp=zz_alpha2_double(p1,p2); + mq=zz_alpha2_double(q1,q2); + mF=zz_alpha2_double(F1,F2); + + zetap=2*t2*cosl(phi)*xip; + zetaq=2*t2*cosl(phi)*xiq; + zetaF=2*t2*cosl(phi)*xiF; + + mp=W-2*t2*sinl(phi)*mp; + mq=W-2*t2*sinl(phi)*mq; + mF=W-2*t2*sinl(phi)*mF; + + xip=sqrtl(mp*mp+2*t1*t1*xip); + xiq=sqrtl(mq*mq+2*t1*t1*xiq); + xiF=sqrtl(mF*mF+2*t1*t1*xiF); + + return((xip+xiq+xiF)*(mp/xip+mq/xiq-mF/xiF-mp*mq*mF/xip/xiq/xiF)*(zetap+zetaq-zetaF)/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/108); +} + +// I for z1+z2 +long double ZZ_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){ + long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF; + + // alpha1 (store in xi) + xip=zz_alpha1_double(p1,p2); + xiq=zz_alpha1_double(q1,q2); + xiF=zz_alpha1_double(F1,F2); + + // alpha2 (store in m) + mp=zz_alpha2_double(p1,p2); + mq=zz_alpha2_double(q1,q2); + mF=zz_alpha2_double(F1,F2); + + zetap=2*t2*cosl(phi)*xip; + zetaq=2*t2*cosl(phi)*xiq; + zetaF=2*t2*cosl(phi)*xiF; + + mp=W-2*t2*sinl(phi)*mp; + mq=W-2*t2*sinl(phi)*mq; + mF=W-2*t2*sinl(phi)*mF; + + xip=sqrtl(mp*mp+2*t1*t1*xip); + xiq=sqrtl(mq*mq+2*t1*t1*xiq); + xiF=sqrtl(mF*mF+2*t1*t1*xiF); + + return((1-mp/xip*mF/xiF-mq/xiq*mF/xiF+mp/xip*mq/xiq)*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)+(xip+xiq+xiF)*(xip+xiq+xiF))/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/216); +} + +// bar alpha_1(k1,k2) +long double zz_alpha1_double(long double k1, long double k2){ + return(2*cosl(PI/3*k2)*(cosl(PI/SQRT3*k1)+cosl(PI/3*k2))+0.5); +} +// bar alpha_2(k1,k2) +long double zz_alpha2_double(long double k1, long double k2){ + return(2*sinl(PI/3*k2)*(cosl(PI/SQRT3*k1)-cosl(PI/3*k2))); +} + diff --git a/src/zz_integral_double.h b/src/zz_integral_double.h new file mode 100644 index 0000000..da21ef3 --- /dev/null +++ b/src/zz_integral_double.h @@ -0,0 +1,35 @@ +/* +Copyright 2016 Ian Jauslin + +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. +*/ + +/* + Computation of z1-z2 and z1+z2 (using doubles) +*/ + +#ifndef ZZ_INTEGRAL_DOUBLE_H +#define ZZ_INTEGRAL_DOUBLE_H + +// I for z1-z2 +long double zz_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W); + +// I for z1+z2 +long double ZZ_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W); + +// bar alpha_1(k1,k2) +long double zz_alpha1_double(long double k1, long double k2); +// bar alpha_2(k1,k2) +long double zz_alpha2_double(long double k1, long double k2); + +#endif |