]> www.ginac.de Git - ginac.git/blobdiff - ginac/inifcns_nstdsums.cpp
function_options::nparams is now always correctly initialized, even if no
[ginac.git] / ginac / inifcns_nstdsums.cpp
index 6217d8b5b2901ace0eb1cd547781cdaf03d0f465..7b47e4910005d0ff88e5bda6a2872bfaefd40c44 100644 (file)
@@ -46,7 +46,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -189,7 +189,7 @@ cln::cl_N Li2_do_sum(const cln::cl_N& x)
 {
        cln::cl_N res = x;
        cln::cl_N resbuf;
-       cln::cl_N num = x;
+       cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits));
        cln::cl_I den = 1; // n^2 = 1
        unsigned i = 3;
        do {
@@ -208,7 +208,7 @@ cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
 {
        std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
        cln::cl_N u = -cln::log(1-x);
-       cln::cl_N factor = u;
+       cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
        cln::cl_N res = u - u*u/4;
        cln::cl_N resbuf;
        unsigned i = 1;
@@ -226,7 +226,7 @@ cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
 // calculates Li(n,x), n>2 without Xn
 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
 {
-       cln::cl_N factor = x;
+       cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits));
        cln::cl_N res = x;
        cln::cl_N resbuf;
        int i=2;
@@ -245,7 +245,7 @@ cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 {
        std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
        cln::cl_N u = -cln::log(1-x);
-       cln::cl_N factor = u;
+       cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
        cln::cl_N res = u;
        cln::cl_N resbuf;
        unsigned i=2;
@@ -461,7 +461,7 @@ static ex Li_evalf(const ex& x1, const ex& x2)
                                return Li(x1, x2).hold();
                        }
                        conv *= x2.op(i);
-                       if (conv >= 1) {
+                       if (abs(conv) >= 1) {
                                return Li(x1, x2).hold();
                        }
                }
@@ -846,7 +846,9 @@ cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t&
        }
 
        // should be done otherwise
-       cln::cl_N xf = x * cln::cl_float(1, prec);
+       cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+       cln::cl_N xf = x * one;
+       //cln::cl_N xf = x * cln::cl_float(1, prec);
 
        cln::cl_N res;
        cln::cl_N resbuf;
@@ -937,9 +939,8 @@ numeric S_num(int n, int p, const numeric& x)
        else if (!x.imag().is_rational())
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 
-
        // [Kol] (5.3)
-       if (cln::realpart(value) < -0.5) {
+       if ((cln::realpart(value) < -0.5) || (n == 0)) {
 
                cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
                                   * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
@@ -1027,6 +1028,10 @@ static ex S_eval(const ex& n, const ex& p, const ex& x)
                        return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
                }
        }
+       if (n.is_zero()) {
+               // [Kol] (5.3)
+               return pow(-log(1-x), p) / factorial(p);
+       }
        return S(n, p, x).hold();
 }
 
@@ -1086,6 +1091,10 @@ REGISTER_FUNCTION(S,
 // anonymous namespace for helper functions
 namespace {
 
+       
+// regulates the pole (used by 1/x-transformation)
+symbol H_polesign("IMSIGN");
+
 
 // convert parameters from H to Li representation
 // parameters are expected to be in expanded form, i.e. only 0, 1 and -1
@@ -1141,10 +1150,6 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
                        }
                }
        }
-       for (; acc > 1; acc--) {
-               throw std::runtime_error("ERROR!");
-               m.append(0);
-       }
        
        return has_negative_parameters;
 }
@@ -1603,7 +1608,7 @@ struct map_trafo_H_1overx : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops())
+                                               return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                }
@@ -1930,6 +1935,11 @@ static ex H_evalf(const ex& x1, const ex& x2)
                        // x -> 1/x
                        map_trafo_H_1overx trafo;
                        res *= trafo(H(m, xtemp));
+                       if (cln::imagpart(x) <= 0) {
+                               res = res.subs(H_polesign == -I*Pi);
+                       } else {
+                               res = res.subs(H_polesign == I*Pi);
+                       }
                }
 
                // simplify result
@@ -2141,9 +2151,9 @@ ex convert_H_to_Li(const ex& m, const ex& x)
        map_trafo_H_reduce_trailing_zeros filter;
        map_trafo_H_convert_to_Li filter2;
        if (is_a<lst>(m)) {
-               return filter2(filter(H(m, x).hold())).eval();
+               return filter2(filter(H(m, x).hold()));
        } else {
-               return filter2(filter(H(lst(m), x).hold())).eval();
+               return filter2(filter(H(lst(m), x).hold()));
        }
 }
 
@@ -2625,7 +2635,7 @@ static void zeta1_print_latex(const ex& m_, const print_context& c)
 }
 
 
-unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta").
+unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1).
                                 evalf_func(zeta1_evalf).
                                 eval_func(zeta1_eval).
                                 derivative_func(zeta1_deriv).
@@ -2762,7 +2772,7 @@ static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c
 }
 
 
-unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta").
+unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2).
                                 evalf_func(zeta2_evalf).
                                 eval_func(zeta2_eval).
                                 derivative_func(zeta2_deriv).