- Introduced derivative of Li2.
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Thu, 6 Jul 2000 23:34:42 +0000 (23:34 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Thu, 6 Jul 2000 23:34:42 +0000 (23:34 +0000)
- Introduced nontrivial part of series expansion of Li2 around 0.  (Read the
  comments in the source to see why it's nontrivial!)  Missing: the branch
  point and cut.

ginac/inifcns.cpp
ginac/pseries.cpp
ginac/utils.cpp
ginac/utils.h

index 85a7365fa980ece03a596b9b9cff70724c8da529..877c2b852710a49acc2ea2e61bc141385afbd423 100644 (file)
@@ -132,16 +132,75 @@ REGISTER_FUNCTION(csgn, eval_func(csgn_eval).
 
 static ex Li2_eval(const ex & x)
 {
+    // Li2(0) -> 0
     if (x.is_zero())
         return x;
+    // Li2(1) -> Pi^2/6
     if (x.is_equal(_ex1()))
-        return power(Pi, _ex2()) / _ex6();
+        return power(Pi,_ex2())/_ex6();
+    // Li2(1/2) -> Pi^2/12 - log(2)^2/2
+    if (x.is_equal(_ex1_2()))
+        return power(Pi,_ex2())/_ex12() + power(log(_ex2()),_ex2())*_ex_1_2();
+    // Li2(-1) -> -Pi^2/12
     if (x.is_equal(_ex_1()))
-        return -power(Pi, _ex2()) / _ex12();
+        return -power(Pi,_ex2())/_ex12();
+    // Li2(I) -> -Pi^2/48+Catalan*I
+    if (x.is_equal(I))
+        return power(Pi,_ex2())/_ex_48() + Catalan*I;
+    // Li2(-I) -> -Pi^2/48-Catalan*I
+    if (x.is_equal(-I))
+        return power(Pi,_ex2())/_ex_48() - Catalan*I;
     return Li2(x).hold();
 }
 
-REGISTER_FUNCTION(Li2, eval_func(Li2_eval));
+static ex Li2_deriv(const ex & x, unsigned deriv_param)
+{
+    GINAC_ASSERT(deriv_param==0);
+    
+    // d/dx Li2(x) -> -log(1-x)/x
+    return -log(1-x)/x;
+}
+
+static ex Li2_series(const ex &x, const relational &rel, int order)
+{
+    const ex x_pt = x.subs(rel);
+    if (!x_pt.is_zero() && !x_pt.is_equal(_ex1()))
+        throw do_taylor();  // caught by function::series()
+    // First case: x==0 (derivatives have poles)
+    if (x_pt.is_zero()) {
+        // method:
+        // The problem is that in d/dx Li2(x==0) == -log(1-x)/x we cannot 
+        // simply substitute x==0.  The limit, however, exists: it is 1.  We
+        // also know all higher derivatives' limits: (d/dx)^n Li2(x) == n!/n^2.
+        // So the primitive series expansion is Li2(x==0) == x + x^2/4 + x^3/9
+        // and so on.
+        // We first construct such a primitive series expansion manually in
+        // a dummy symbol s and then insert the argument's series expansion
+        // for s.  Reexpanding the resulting series returns the desired result.
+        const symbol s;
+        ex ser;
+        // construct manually the primitive expansion
+        for (int i=1; i<order; ++i)
+            ser += pow(s,i)/pow(numeric(i),numeric(2));
+        // substitute the argument's series expansion
+        ser = ser.subs(s==x.series(rel,order));
+        // maybe that was terminanting, so add a proper order term
+        epvector nseq;
+        nseq.push_back(expair(Order(_ex1()), numeric(order)));
+        ser += pseries(rel, nseq);
+        // reexpand will collapse the series again
+        ser = ser.series(rel,order);
+        return ser;
+    }
+    // second problematic case: x real, >=1 (branch cut)
+    return pseries();
+    // TODO: Li2_series should do something around branch point?
+    // Careful: may involve logs!
+}
+
+REGISTER_FUNCTION(Li2, eval_func(Li2_eval).
+                       derivative_func(Li2_deriv).
+                       series_func(Li2_series));
 
 //////////
 // trilogarithm
index 0782be07ebd34aaeb108e2358b6fd8328422360f..d12907a2cbb43d6c1efc20ad06ca6744fc1d1ee8 100644 (file)
@@ -454,7 +454,7 @@ ex basic::series(const relational & r, int order) const
         seq.push_back(expair(coeff, numeric(0)));
     
     int n;
-    for (n=1; n<order; n++) {
+    for (n=1; n<order; ++n) {
         fac = fac.mul(numeric(n));
         deriv = deriv.diff(*s).expand();
         if (deriv.is_zero()) {
index f82addff4328eceffab6fe3d9d522096f7f80cac..5a92a31c17b4b879e4ea741f4d5f0ece33ae3726 100644 (file)
@@ -53,6 +53,20 @@ int compare_pointers(const void * a, const void * b)
 // `construct on first use' chest of numbers
 //////////
 
+// numeric -120
+const numeric & _num_120(void)
+{
+    const static ex e = ex(numeric(-120));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+const ex & _ex_120(void)
+{
+    static ex * e = new ex(_num_120());
+    return *e;
+}
+
 // numeric -60
 const numeric & _num_60(void)
 {
@@ -67,17 +81,17 @@ const ex & _ex_60(void)
     return *e;
 }
 
-// numeric -120
-const numeric & _num_120(void)
+// numeric -48
+const numeric & _num_48(void)
 {
-    const static ex e = ex(numeric(-120));
+    const static ex e = ex(numeric(-48));
     const static numeric * n = static_cast<const numeric *>(e.bp);
     return *n;
 }
 
-const ex & _ex_120(void)
+const ex & _ex_48(void)
 {
-    static ex * e = new ex(_num_120());
+    static ex * e = new ex(_num_48());
     return *e;
 }
 
@@ -683,6 +697,20 @@ const ex & _ex30(void)
     return *e;
 }
 
+// numeric  48
+const numeric & _num48(void)
+{
+    const static ex e = ex(numeric(48));
+    const static numeric * n = static_cast<const numeric *>(e.bp);
+    return *n;
+}
+
+const ex & _ex48(void)
+{
+    static ex * e = new ex(_num48());
+    return *e;
+}
+
 // numeric  60
 const numeric & _num60(void)
 {
index 13fa2c8b956bafc17d8ca33640b8575467b9ba44..107a6f5282d6eb15aa500acb7bd893c47e7c6f57 100644 (file)
@@ -174,6 +174,8 @@ const numeric & _num_120(void);   // -120
 const ex & _ex_120(void);
 const numeric & _num_60(void);    // -60
 const ex & _ex_60(void);
+const numeric & _num_48(void);    // -48
+const ex & _ex_48(void);
 const numeric & _num_30(void);    // -30
 const ex & _ex_30(void);
 const numeric & _num_25(void);    // -25
@@ -260,6 +262,8 @@ const numeric & _num25(void);     //  25
 const ex & _ex25(void);
 const numeric & _num30(void);     //  30
 const ex & _ex30(void);
+const numeric & _num48(void);     //  48
+const ex & _ex48(void);
 const numeric & _num60(void);     //  60
 const ex & _ex60(void);
 const numeric & _num120(void);    //  120