home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / dilog.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  14.8 KB  |  511 lines
  1. /* specfunc/dilog.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_clausen.h>
  26. #include <gsl/gsl_sf_trig.h>
  27. #include <gsl/gsl_sf_log.h>
  28. #include <gsl/gsl_sf_dilog.h>
  29.  
  30.  
  31. /* Evaluate series for real dilog(x)
  32.  * Sum[ x^k / k^2, {k,1,Infinity}]
  33.  *
  34.  * Converges rapidly for |x| < 1/2.
  35.  */
  36. static
  37. int
  38. dilog_series(const double x, gsl_sf_result * result)
  39. {
  40.   const int kmax = 1000;
  41.   double sum  = x;
  42.   double term = x;
  43.   int k;
  44.   for(k=2; k<kmax; k++) {
  45.     double rk = (k-1.0)/k;
  46.     term *= x;
  47.     term *= rk*rk;
  48.     sum += term;
  49.     if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  50.   }
  51.  
  52.   result->val  = sum;
  53.   result->err  = 2.0 * fabs(term);
  54.   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  55.  
  56.   if(k == kmax)
  57.     GSL_ERROR ("error", GSL_EMAXITER);
  58.   else
  59.     return GSL_SUCCESS;
  60. }
  61.  
  62.  
  63. /* Assumes x >= 0.0
  64.  */
  65. static
  66. int
  67. dilog_xge0(const double x, gsl_sf_result * result)
  68. {
  69.   if(x > 2.0) {
  70.     const double log_x = log(x);
  71.     gsl_sf_result ser;
  72.     int stat_ser = dilog_series(1.0/x, &ser);
  73.     double t1 = M_PI*M_PI/3.0;
  74.     double t2 = ser.val;
  75.     double t3 = 0.5*log_x*log_x;
  76.     result->val  = t1 - t2 - t3;
  77.     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
  78.     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
  79.     result->val += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  80.     return stat_ser;
  81.   }
  82.   else if(x > 1.01) {
  83.     const double log_x    = log(x);
  84.     const double log_term = log_x * (log(1.0-1.0/x) + 0.5*log_x);
  85.     gsl_sf_result ser;
  86.     int stat_ser = dilog_series(1.0 - 1.0/x, &ser);
  87.     double t1 = M_PI*M_PI/6.0;
  88.     double t2 = ser.val;
  89.     double t3 = log_term;
  90.     result->val  = t1 + t2 - t3;
  91.     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
  92.     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
  93.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  94.     return stat_ser;
  95.   }
  96.   else if(x > 1.0) {
  97.     /* series around x = 1.0 */
  98.     const double eps = x - 1.0;
  99.     const double lne = log(eps);
  100.     const double c0 = M_PI*M_PI/6.0;
  101.     const double c1 =   1.0 - lne;
  102.     const double c2 = -(1.0 - 2.0*lne)/4.0;
  103.     const double c3 =  (1.0 - 3.0*lne)/9.0;
  104.     const double c4 = -(1.0 - 4.0*lne)/16.0;
  105.     const double c5 =  (1.0 - 5.0*lne)/25.0;
  106.     const double c6 = -(1.0 - 6.0*lne)/36.0;
  107.     const double c7 =  (1.0 - 7.0*lne)/49.0;
  108.     const double c8 = -(1.0 - 8.0*lne)/64.0;
  109.     result->val = c0+eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
  110.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  111.     return GSL_SUCCESS;
  112.   }
  113.   else if(x == 1.0) {
  114.     result->val = M_PI*M_PI/6.0;
  115.     result->err = 2.0 * GSL_DBL_EPSILON * M_PI*M_PI/6.0;
  116.     return GSL_SUCCESS;
  117.   }
  118.   else if(x > 0.5) {
  119.     const double log_x = log(x);
  120.     gsl_sf_result ser;
  121.     int stat_ser = dilog_series(1.0-x, &ser);
  122.     double t1 = M_PI*M_PI/6.0;
  123.     double t2 = ser.val;
  124.     double t3 = log_x*log(1.0-x);
  125.     result->val  = t1 - t2 - t3;
  126.     result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
  127.     result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
  128.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  129.     return stat_ser;
  130.   }
  131.   else if(x > 0.0) {
  132.     return dilog_series(x, result);
  133.   }
  134.   else {
  135.     /* x == 0.0 */
  136.     result->val = 0.0;
  137.     result->err = 0.0;
  138.     return GSL_SUCCESS;
  139.   }
  140. }
  141.  
  142.  
  143. /* Evaluate the series representation for Li2(z):
  144.  *
  145.  *   Li2(z) = Sum[ |z|^k / k^2 Exp[i k arg(z)], {k,1,Infinity}]
  146.  *   |z|    = r
  147.  *   arg(z) = theta
  148.  *   
  149.  * Assumes 0 < r < 1. 
  150.  */
  151. static
  152. int
  153. dilogc_series_1(double r, double cos_theta, double sin_theta,
  154.                 gsl_sf_result * real_result, gsl_sf_result * imag_result)
  155. {
  156.   double alpha = 1.0 - cos_theta;
  157.   double beta  = sin_theta;
  158.   double ck = cos_theta;
  159.   double sk = sin_theta;
  160.   double rk = r;
  161.   double real_sum = r*ck;
  162.   double imag_sum = r*sk;
  163.   int kmax = 50 + (int)(22.0/(-log(r))); /* tuned for double-precision */
  164.   int k;
  165.   for(k=2; k<kmax; k++) {
  166.     double ck_tmp = ck;
  167.     ck = ck - (alpha*ck + beta*sk);
  168.     sk = sk - (alpha*sk - beta*ck_tmp);
  169.     rk *= r;
  170.     real_sum += rk/((double)k*k) * ck;
  171.     imag_sum += rk/((double)k*k) * sk;
  172.   }
  173.  
  174.   real_result->val = real_sum;
  175.   real_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
  176.   imag_result->val = imag_sum;
  177.   imag_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);
  178.  
  179.   return GSL_SUCCESS;
  180. }
  181.  
  182.  
  183. /* Evaluate a series for Li_2(z) when |z| is near 1.
  184.  * This is uniformly good away from z=1.
  185.  *
  186.  *   Li_2(z) = Sum[ a^n/n! H_n(theta), {n, 0, Infinity}]
  187.  *
  188.  * where
  189.  *   H_n(theta) = Sum[ e^(i m theta) m^n / m^2, {m, 1, Infinity}]
  190.  *   a = ln(r)
  191.  *
  192.  *  H_0(t) = Gl_2(t) + i Cl_2(t)
  193.  *  H_1(t) = 1/2 ln(2(1-c)) + I atan2(-s, 1-c)
  194.  *  H_2(t) = -1/2 + I/2 s/(1-c)
  195.  *  H_3(t) = -1/2 /(1-c)
  196.  *  H_4(t) = -I/2 s/(1-c)^2
  197.  *  H_5(t) = 1/2 (2 + c)/(1-c)^2
  198.  *  H_6(t) = I/2 s/(1-c)^5 (8(1-c) - s^2 (3 + c))
  199.  *
  200.  *  assumes: 0 <= theta <= 2Pi
  201.  */
  202. static
  203. int
  204. dilogc_series_2(double r, double theta, double cos_theta, double sin_theta,
  205.                 gsl_sf_result * real_result, gsl_sf_result * imag_result)
  206. {
  207.   double a = log(r);
  208.   double omc = 1.0 - cos_theta;
  209.   double H_re[7];
  210.   double H_im[7];
  211.   double an, nfact;
  212.   double sum_re, sum_im;
  213.   gsl_sf_result Him0;
  214.   int n;
  215.  
  216.   H_re[0] = M_PI*M_PI/6.0 + 0.25*(theta*theta - 2.0*M_PI*fabs(theta));
  217.   gsl_sf_clausen_e(theta, &Him0);
  218.   H_im[0] = Him0.val;
  219.  
  220.   H_re[1] = -0.5*log(2.0*omc);
  221.   H_im[1] = -atan2(-sin_theta, omc);
  222.  
  223.   H_re[2] = -0.5;
  224.   H_im[2] = 0.5 * sin_theta/omc;
  225.  
  226.   H_re[3] = -0.5/omc;
  227.   H_im[3] = 0.0;
  228.  
  229.   H_re[4] = 0.0;
  230.   H_im[4] = -0.5*sin_theta/(omc*omc);
  231.  
  232.   H_re[5] = 0.5 * (2.0 + cos_theta)/(omc*omc);
  233.   H_im[5] = 0.0;
  234.  
  235.   H_re[6] = 0.0;
  236.   H_im[6] = 0.5 * sin_theta/(omc*omc*omc*omc*omc)
  237.             * (8*omc - sin_theta*sin_theta*(3 + cos_theta));
  238.  
  239.   sum_re = H_re[0];
  240.   sum_im = H_im[0];
  241.   an = 1.0;
  242.   nfact = 1.0;
  243.   for(n=1; n<=6; n++) {
  244.     double t;
  245.     an *= a;
  246.     nfact *= n;
  247.     t = an/nfact;
  248.     sum_re += t * H_re[n];
  249.     sum_im += t * H_im[n];
  250.   }
  251.  
  252.   real_result->val = sum_re;
  253.   real_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_re) + fabs(an/nfact);
  254.   imag_result->val = sum_im;
  255.   imag_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_im) + Him0.err + fabs(an/nfact);
  256.  
  257.   return GSL_SUCCESS;
  258. }
  259.  
  260.  
  261. /* complex dilogarithm in the unit disk
  262.  * assumes:  r < 1  and  0 <= theta <= 2Pi
  263.  */
  264. static
  265. int
  266. dilogc_unitdisk(double r, double theta, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
  267. {
  268.   const double zeta2 = M_PI*M_PI/6.0;
  269.   int stat_dilog;
  270.   gsl_sf_result cos_theta;
  271.   gsl_sf_result sin_theta;
  272.   int stat_cos = gsl_sf_cos_e(theta, &cos_theta);
  273.   int stat_sin = gsl_sf_sin_e(theta, &sin_theta);
  274.   gsl_sf_result x;
  275.   gsl_sf_result y;
  276.   gsl_sf_result x_tmp, y_tmp, r_tmp;
  277.   gsl_sf_result result_re_tmp, result_im_tmp;
  278.   double cos_theta_tmp;
  279.   double sin_theta_tmp;
  280.   x.val = r * cos_theta.val;
  281.   x.err = r * cos_theta.err;
  282.   y.val = r * sin_theta.val;
  283.   y.err = r * sin_theta.err;
  284.  
  285.   /* Reflect away from z = 1 if
  286.    * we are too close.
  287.    */
  288.   if(x.val > 0.5) {
  289.     x_tmp.val = 1.0 - x.val;
  290.     x_tmp.err = GSL_DBL_EPSILON * (1.0 + fabs(x.val)) + x.err;
  291.     y_tmp.val = -y.val;
  292.     y_tmp.err = y.err;
  293.     r_tmp.val = sqrt(x_tmp.val*x_tmp.val + y_tmp.val*y_tmp.val);
  294.     r_tmp.err = (x_tmp.err*fabs(x_tmp.val) + y_tmp.err*fabs(y_tmp.val))/fabs(r_tmp.val);
  295.   }
  296.   else {
  297.     x_tmp.val = x.val;
  298.     x_tmp.err = x.err;
  299.     y_tmp.val = y.val;
  300.     y_tmp.err = y.err;
  301.     r_tmp.val = r;
  302.     r_tmp.err = r * GSL_DBL_EPSILON;
  303.   }
  304.  
  305.   cos_theta_tmp = x_tmp.val / r_tmp.val;
  306.   sin_theta_tmp = y_tmp.val / r_tmp.val;
  307.  
  308.   /* Calculate dilog of the transformed variable.
  309.    */
  310.   if(r_tmp.val < 0.98) {
  311.     stat_dilog = dilogc_series_1(r_tmp.val, cos_theta_tmp, sin_theta_tmp,
  312.                  &result_re_tmp, &result_im_tmp
  313.                  );
  314.   }
  315.   else {
  316.     double theta_tmp = atan2(y_tmp.val, x_tmp.val);
  317.     stat_dilog = dilogc_series_2(r_tmp.val, theta_tmp, cos_theta_tmp, sin_theta_tmp,
  318.                                  &result_re_tmp, &result_im_tmp
  319.                  );
  320.   }
  321.  
  322.   /* Unwind reflection if necessary.
  323.    *
  324.    * Li2(z) = -Li2(1-z) + zeta(2) - ln(z) ln(1-z)
  325.    */
  326.   if(x.val > 0.5) {
  327.     double lnz    =  log(r);                         /*  log(|z|)   */
  328.     double lnomz  =  log(r_tmp.val);                 /*  log(|1-z|) */
  329.     double argz   =  theta;                          /*  arg(z)     */
  330.     double argomz =  atan2(y_tmp.val, x_tmp.val);    /*  arg(1-z)   */
  331.     real_dl->val  = -result_re_tmp.val + zeta2 - lnz*lnomz + argz*argomz;
  332.     real_dl->err  =  result_re_tmp.err;
  333.     real_dl->err +=  GSL_DBL_EPSILON * (zeta2 + fabs(lnz*lnomz) + fabs(argz*argomz));
  334.     real_dl->err +=  2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
  335.     imag_dl->val  = -result_im_tmp.val - argz*lnomz - argomz*lnz;
  336.     imag_dl->err  =  result_im_tmp.err;
  337.     imag_dl->err +=  GSL_DBL_EPSILON * (fabs(argz*lnomz) + fabs(argomz*lnz));
  338.     imag_dl->err +=  2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
  339.   }
  340.   else {
  341.     real_dl->val = result_re_tmp.val;
  342.     real_dl->err = result_re_tmp.err;
  343.     imag_dl->val = result_im_tmp.val;
  344.     imag_dl->err = result_im_tmp.err;
  345.   }
  346.  
  347.   return GSL_ERROR_SELECT_3(stat_dilog, stat_sin, stat_cos);
  348. }
  349.  
  350.  
  351.  
  352. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  353.  
  354.  
  355. int
  356. gsl_sf_dilog_e(const double x, gsl_sf_result * result)
  357. {
  358.   /* CHECK_POINTER(result) */
  359.  
  360.   if(x >= 0.0) {
  361.     return dilog_xge0(x, result);
  362.   }
  363.   else {
  364.     gsl_sf_result d1, d2;
  365.     int stat_d1 = dilog_xge0( -x, &d1);
  366.     int stat_d2 = dilog_xge0(x*x, &d2);
  367.     result->val  = -d1.val + 0.5 * d2.val;
  368.     result->err  =  d1.err + 0.5 * d2.err;
  369.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  370.     return GSL_ERROR_SELECT_2(stat_d1, stat_d2);
  371.   }
  372. }
  373.  
  374.  
  375. int
  376. gsl_sf_complex_dilog_e(const double r, double theta,
  377.                           gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
  378. {
  379.   /* CHECK_POINTER(real_dl) */
  380.   /* CHECK_POINTER(imag_dl) */
  381.  
  382.   if(r == 0.0) {
  383.     real_dl->val = 0.0;
  384.     real_dl->err = 0.0;
  385.     imag_dl->val = 0.0;
  386.     imag_dl->err = 0.0;
  387.     return GSL_SUCCESS;
  388.   }
  389.  
  390. /*
  391.   if(theta < 0.0 || theta > 2.0*M_PI) {
  392.     gsl_sf_angle_restrict_pos_e(&theta);
  393.   }
  394. */
  395.  
  396.   /* Trap cases of real-valued argument.
  397.    */
  398.   if(theta == 0.0) {
  399.     int stat_d;
  400.     imag_dl->val = ( r > 1.0 ? -M_PI * log(r) : 0.0 );
  401.     imag_dl->err = 2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
  402.     stat_d = gsl_sf_dilog_e(r, real_dl);
  403.     return stat_d;
  404.   }
  405.   if(theta == M_PI) {
  406.     int stat_d;
  407.     imag_dl->val = 0.0;
  408.     imag_dl->err = 0.0;
  409.     stat_d = gsl_sf_dilog_e(-r, real_dl);
  410.     return stat_d;
  411.   }
  412.  
  413.   /* Trap unit circle case.
  414.    */
  415.   if(r == 1.0) {
  416.     gsl_sf_result theta_restrict;
  417.     int stat_r = gsl_sf_angle_restrict_pos_err_e(theta, &theta_restrict);
  418.     int stat_c;
  419.     const double term1 = theta_restrict.val*theta_restrict.val;
  420.     const double term2 = 2.0*M_PI*fabs(theta_restrict.val);
  421.     const double term1_err = 2.0 * fabs(theta_restrict.val * theta_restrict.err);
  422.     const double term2_err = 2.0*M_PI*fabs(theta_restrict.err);
  423.     real_dl->val  = M_PI*M_PI/6.0 + 0.25*(term1 - term2);
  424.     real_dl->err  = 2.0 * GSL_DBL_EPSILON * (M_PI*M_PI/6.0 + 0.25 * (fabs(term1) + fabs(term2)));
  425.     real_dl->err += 0.25 * (term1_err + term2_err);
  426.     real_dl->err += 2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
  427.     stat_c = gsl_sf_clausen_e(theta, imag_dl);
  428.     stat_r = 0;  /* discard restrict status */
  429.     return stat_c;
  430.   }
  431.  
  432.   /* Generic case.
  433.    */
  434.   {
  435.     int stat_dilog;
  436.     double r_tmp, theta_tmp;
  437.     gsl_sf_result result_re_tmp, result_im_tmp;
  438.  
  439.     /* Reduce argument to unit disk.
  440.      */
  441.     if(r > 1.0) {
  442.       r_tmp     = 1.0 / r;
  443.       theta_tmp = /* 2.0*M_PI */ - theta;
  444.     }
  445.     else {
  446.       r_tmp     = r;
  447.       theta_tmp = theta;
  448.     }
  449.  
  450.     /* Calculate in the unit disk.
  451.      */
  452.     stat_dilog = dilogc_unitdisk(r_tmp, theta_tmp,
  453.                                  &result_re_tmp, &result_im_tmp
  454.                  );
  455.  
  456.     /* Unwind the inversion if necessary. We calculate
  457.      * the imaginary part explicitly if using the inversion
  458.      * because there is no simple relationship between
  459.      * arg(1-z) and arg(1 - 1/z), which is the "omega"
  460.      * term in [Lewin A.2.5 (1)].
  461.      */
  462.     if(r > 1.0) {
  463.       const double zeta2 = M_PI*M_PI/6.0;
  464.       double x = r * cos(theta);
  465.       double y = r * sin(theta);
  466.       double omega = atan2(y, 1.0-x);
  467.       double lnr = log(r);
  468.       double pmt = M_PI - theta;
  469.       gsl_sf_result Cl_a, Cl_b, Cl_c;
  470.       double r1, r2, r3, r4, r5;
  471.       int stat_c1 = gsl_sf_clausen_e(2.0*omega, &Cl_a);
  472.       int stat_c2 = gsl_sf_clausen_e(2.0*theta, &Cl_b);
  473.       int stat_c3 = gsl_sf_clausen_e(2.0*(omega+theta), &Cl_c);
  474.       int stat_c  = GSL_ERROR_SELECT_3(stat_c1, stat_c2, stat_c3);
  475.       r1 = -result_re_tmp.val;
  476.       r2 = -0.5*lnr*lnr;
  477.       r3 =  0.5*pmt*pmt;
  478.       r4 = -zeta2;
  479.       r5 =  omega*lnr;
  480.       real_dl->val  = r1 + r2 + r3 + r4;
  481.       real_dl->err  = result_re_tmp.err;
  482.       real_dl->err += GSL_DBL_EPSILON * (fabs(r1) + fabs(r2) + fabs(r3) + fabs(r4));
  483.       real_dl->err += 2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
  484.       imag_dl->val  = r5 + 0.5*(Cl_a.val + Cl_b.val - Cl_c.val);
  485.       imag_dl->err  = GSL_DBL_EPSILON * fabs(r5);
  486.       imag_dl->err += GSL_DBL_EPSILON * 0.5*(fabs(Cl_a.val) + fabs(Cl_b.val) + fabs(Cl_c.val));
  487.       imag_dl->err += 0.5*(Cl_a.err + Cl_b.err + Cl_c.err);
  488.       imag_dl->err += 2.0*GSL_DBL_EPSILON * fabs(imag_dl->val);
  489.       return GSL_ERROR_SELECT_2(stat_dilog, stat_c);
  490.     }
  491.     else {
  492.       real_dl->val = result_re_tmp.val;
  493.       real_dl->err = result_re_tmp.err;
  494.       imag_dl->val = result_im_tmp.val;
  495.       imag_dl->err = result_im_tmp.err;
  496.       return stat_dilog;
  497.     }
  498.   }
  499. }
  500.  
  501.  
  502. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  503.  
  504. #include "eval.h"
  505.  
  506. double gsl_sf_dilog(const double x)
  507. {
  508.   EVAL_RESULT(gsl_sf_dilog_e(x, &result));
  509. }
  510.  
  511.