Changeset 5243:9c6972d1055c

Timestamp:

07/01/07 14:15:09 (6 years ago)

Author:

William Stein <wstein@…>

Branch:

default

Message:

Reorganize and fix a bunch of basic arithmetic; make Li much much faster.

Location:

sage

Files:

• functions/all.py (modified) (1 diff)
• functions/constants.py (modified) (2 diffs)
• functions/transcendental.py (modified) (2 diffs)
• rings/arith.py (modified) (2 diffs)

sage/functions/all.py

@@ -6 +6 @@
 from transcendental import (exponential_integral_1,
                             gamma, gamma_inc, incomplete_gamma,
-                            zeta, zeta_symmetric)
+                            zeta, zeta_symmetric,
+                            Li,
+                            prime_pi, 
+                            number_of_primes)
 
 #from elementary import (cosine, sine, exponential,

sage/functions/constants.py

@@ -465 +465 @@
 pi = Pi()
 
+python_complex_i = complex(0,1)
 
 class I_class(Constant):
@@ -526 +527 @@
         return C.gen()
 
+    def __complex__(self):
+        return python_complex_i
+
     def _real_double_(self, R):
         raise TypeError        

sage/functions/transcendental.py

@@ -19 +19 @@
 
 import  sage.libs.pari.all
-pari = sage.libs.pari.all.pari
+from sage.libs.pari.all import pari, PariError
 import sage.rings.complex_field as complex_field
 import sage.rings.real_mpfr as real_field
+import sage.rings.real_double as real_double
 import sage.rings.complex_number
-
-from sage.rings.all import is_RealNumber, RealField, is_ComplexNumber, ComplexField
+from sage.gsl.integration import numerical_integral
+
+from sage.rings.all import (is_RealNumber, RealField,
+                            is_ComplexNumber, ComplexField,
+                            ZZ)
 
 def __prep_num(x):
@@ -219 +223 @@
 #    """
 #   return real_field.RealField(prec).pi()
+
+
+
+
+def Li(x, eps_rel=None, err_bound=False):
+    r"""
+    Return value of the function Li(x) as a real double field element.
+
+    This is the function
+    $$
+       \int_2^{x} dt / \log(t).
+    $$
+    
+    The function Li(x) is an approximation for the number
+    of primes up to $x$.  In fact, the famous Riemann
+    Hypothesis is equivalent to the statement that for
+    $x \geq 2.01$ we have
+    $$
+        |\pi(x) - Li(x)| \leq \sqrt{x} \log(x).
+    $$
+    For ``small'' $x$, $Li(x)$ is always slightly bigger than
+    $\pi(x)$.  However it is a theorem that there are (very large,
+    e.g., around $10^{316}$) values of $x$ so that $\pi(x) > Li(x)$.
+    See ``A new bound for the smallest x with $\pi(x) > li(x)$'',
+    Bays and Hudson, Mathematics of Computation, 69 (2000) 1285--1296.
+    
+    ALGORITHM: Computed numerically using GSL.
+
+    INPUT:
+        x -- a real number >= 2.
+
+    OUTPUT:
+        x -- a real double
+
+    EXAMPLES:
+        sage: Li(2)
+        0.0
+        sage: Li(5)
+        2.58942452992
+        sage: Li(1000)
+        176.56449421
+        sage: Li(10^5)
+        9628.76383727
+        sage: prime_pi(10^5)
+        9592
+        sage: Li(1)
+        Traceback (most recent call last):
+        ...
+        ValueError: Li only defined for x at least 2.    
+
+        sage: for n in range(1,7):
+        ...    print '%-10s%-10s%-20s'%(10^n, prime_pi(10^n), Li(10^n))
+        10        4         5.12043572467       
+        100       25        29.080977804        
+        1000      168       176.56449421        
+        10000     1229      1245.09205212       
+        100000    9592      9628.76383727       
+        1000000   78498     78626.5039957           
+    """
+    x = float(x)
+    if x < 2:
+        raise ValueError, "Li only defined for x at least 2."
+    if eps_rel:
+        ans = numerical_integral(_one_over_log, 2, float(x),
+                             eps_rel=eps_rel)
+    else:
+        ans = numerical_integral(_one_over_log, 2, float(x))
+    if err_bound:
+        return real_double.RDF(ans[0]), ans[1]
+    else:
+        return real_double.RDF(ans[0])
+    # Old PARI version -- much much slower
+    #x = RDF(x)
+    #return RDF(gp('intnum(t=2,%s,1/log(t))'%x))
+
+import math
+def _one_over_log(t):
+    return 1/math.log(t)    
+
+def prime_pi(x):
+    """
+    Return the number of primes $\leq x$.
+
+    EXAMPLES:
+        sage: prime_pi(7)
+        4
+        sage: prime_pi(100)
+        25
+        sage: prime_pi(1000)
+        168
+        sage: prime_pi(100000)
+        9592
+        sage: prime_pi(0.5)
+        0
+        sage: prime_pi(-10)
+        0
+    """
+    if x < 2:
+        return ZZ(0)
+    try:
+        return ZZ(pari(x).primepi())
+    except PariError:
+        pari.init_primes(pari(x)+1)
+        return ZZ(pari(x).primepi())
+
+number_of_primes = prime_pi
+
+

sage/rings/arith.py

@@ -207 +207 @@
     else:
         raise ValueError, "invalid choice of algorithm"
-
-def Li(x):
-    r"""
-    Return value of the function Li(x), which is by definition
-    $$
-       \int_2^{x} dt / \log(t).
-    $$
-    
-    The function Li(x) is an approximation for the number
-    of primes up to $x$.  In fact, the famous Riemann
-    Hypothesis is equivalent to the statement that for
-    $x \geq 2.01$ we have
-    $$
-        |\pi(x) - Li(x)| \leq \sqrt{x} \log(x).
-    $$
-    For ``small'' $x$, $Li(x)$ is always slightly bigger than
-    $\pi(x)$.  However it is a theorem that there are (very large,
-    e.g., around $10^{316}$) values of $x$ so that $\pi(x) > Li(x)$.
-    See ``A new bound for the smallest x with $\pi(x) > li(x)$'',
-    Bays and Hudson, Mathematics of Computation, 69 (2000) 1285--1296.
-    
-    ALGORITHM: Computed numerically using PARI.
-
-    INPUT:
-        x -- a real number >= 2.
-
-    OUTPUT:
-        x -- a real double
-
-    EXAMPLES:
-        sage: pari.init_primes(10^6)   # needed to compute prime_pi(10^6)
-        sage: for n in range(1,7):
-        ...    print '%-10s%-10s%-20s'%(10^n, prime_pi(10^n), Li(10^n))
-        10        4         5.12043572467       
-        100       25        29.080977804        
-        1000      168       176.56449421        
-        10000     1229      1245.09205212       
-        100000    9592      9628.76383727       
-        1000000   78498     78626.5039957           
-    """
-    from real_double import RDF
-    x = RDF(x)
-    return RDF(gp('intnum(t=2,%s,1/log(t))'%x))    
-
-def prime_pi(x):
-    """
-    Return the number of primes $\leq x$.
-
-    EXAMPLES:
-        sage: prime_pi(7)
-        4
-        sage: prime_pi(100)
-        25
-        sage: prime_pi(1000)
-        168
-        sage: prime_pi(100000)
-        9592
-        sage: prime_pi(0.5)
-        0
-        sage: prime_pi(-10)
-        0
-    """
-    if x < 2:
-        return integer_ring.ZZ(0)
-    try:
-        return integer_ring.ZZ(pari(x).primepi())
-    except PariError:
-        pari.init_primes(pari(x)+1)
-        return integer_ring.ZZ(pari(x).primepi())
-
-number_of_primes = prime_pi
 
 
@@ -2207 +2136 @@
                 return lower
             upper = mediant
-
-def number_of_partitions(n):
-    """
-    Return the number of partitions of the integer $n$.
-
-    To compute all the partitions of $n$ use \code{partitions(n)}.
-
-    EXAMPLES:
-        sage: number_of_partitions(3)
-        3
-        sage: number_of_partitions(10)
-        42
-        sage: number_of_partitions(40)
-        37338
-        sage: number_of_partitions(100)
-        190569292
-        sage: number_of_partitions(-5)
-        0
-        sage: number_of_partitions(0)
-        1
-    """
-    ZZ = integer_ring.ZZ
-    return ZZ(pari(ZZ(n)).numbpart())
-
-def partitions(n):
-    """
-    Generator of all the partitions of the integer $n$.
-
-    To compute the number of partitions of $n$ use
-    \code{number_of_partitions(n)}.
-
-    INPUT:
-        n -- int
-
-    EXAMPLES:
-        >> partitions(3)
-        <generator object at 0xab3b3eac>
-        sage: list(partitions(3))
-        [(1, 1, 1), (1, 2), (3,)]
-
-    AUTHOR: David Eppstein, Jan Van lent, George Yoshida; Python Cookbook 2, Recipe 19.16.
-    """
-    n == integer_ring.ZZ(n)
-    # base case of the recursion: zero is the sum of the empty tuple
-    if n == 0:
-        yield ( )
-        return
-    # modify the partitions of n-1 to form the partitions of n
-    for p in partitions(n-1):
-        yield (1,) + p
-        if p and (len(p) < 2 or p[1] > p[0]):
-            yield (p[0] + 1,) + p[1:]
-