1 files changed, 237 insertions, 0 deletions
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);
+}