OKlibrary  0.2.1.6
PrimeNumbers.mac
Go to the documentation of this file.
00001 /* Oliver Kullmann, 20.9.2008 (Swansea) */
00002 /* Copyright 2008, 2009, 2010, 2011 Oliver Kullmann
00003 This file is part of the OKlibrary. OKlibrary is free software; you can redistribute
00004 it and/or modify it under the terms of the GNU General Public License as published by
00005 the Free Software Foundation and included in this library; either version 3 of the
00006 License, or any later version. */
00007 
00022 oklib_include("OKlib/ComputerAlgebra/Hypergraphs/Lisp/Generators/GreenTao.mac")$
00023 oklib_include("OKlib/ComputerAlgebra/NumberTheory/Lisp/Auxiliary.mac")$
00024 
00025 
00026 /* **********************************
00027    * Finding and enumerating primes *
00028    **********************************
00029 */
00030 
00031 /* Maxima has: next_prime, prev_prime, primep (controlled by
00032    primep_number_of_tests )
00033 */
00034 
00035 /* The memoised version of primep: */
00036 mprimep[n] := primep(n)$
00037 
00038 /* The list of primes in the closed interval from a to b: */
00039 primes_interval(a,b) := block([L : [], p : if a <= 1 then 1 else a-1],
00040   p : next_prime(p),
00041   while p <= b do (
00042     L : cons(p,L),
00043     p : next_prime(p)
00044   ),
00045   reverse(L))$
00046 
00047 /* The list of the first n prime numbers (n an integer): */
00048 primes_first(n) := if n <= 0 then [] else
00049   block([p : 2, L : [2]],
00050     n : n - 1,
00051     while n > 0 do (
00052       p : next_prime(p),
00053       L : cons(p,L),
00054       n : n - 1
00055     ),
00056     reverse(L))$
00057 
00058 /* Computing the i-th prime number (i a natural number >= 1): */
00059 unrank_primes(i) := block([p : 2],
00060   while i > 1 do (
00061     p : next_prime(p),
00062     i : i - 1
00063   ),
00064   p)$
00065 
00066 /* For a natural number n >= 1 compute the smallest index i >= 1
00067    such that unrank_primes(i) >= n:
00068 */
00069 rank_primes(n) := block([i : 1],
00070   while n > 2 do (
00071     n : prev_prime(n),
00072     i : i + 1
00073   ),
00074   i)$
00075 /* Ranking a set or list P of natural numbers: */
00076 rank_primes_s(P) := map(rank_primes, P)$
00077 /* Remark: If n is not a prime number, then we have
00078    rank_primes(n) = count_primes(n) + 1, while otherwise we have
00079    equality.
00080 */
00081 
00082 /* For a natural number n >= 0 or inf compute the numbers of primes p
00083    with p <= n (in mathematics this is denoted by pi(n)):
00084 */
00085 count_primes(n) := if n <= 1 then 0
00086  elseif n=inf then inf
00087  elseif primep(n) then rank_primes(n)
00088  else rank_primes(prev_prime(n))$
00089 
00090 
00091 /* The product of the prime numbers less than or equal to n >= 0
00092    (in mathematics this is often denoted by Pi(n)): */
00093 product_primes(n) := block([pr : 1, p : 2],
00094  while p <= n do (pr : pr * p, p : next_prime(p)),
00095  return(pr))$
00096 
00097 
00098 /* **********************************
00099    * The Prime Number Theorem (PNT) *
00100    **********************************
00101 */
00102 
00103 /* The simplest approximation of count_primes(x) via the PNT: */
00104 approx_count_primes_0(x) := x / log(x)$
00105 approx_count_primes_0_hp(x,decimal_digits) := block([fpprec : decimal_digits],
00106  bfloat(x / log(x)))$
00107 approx_count_primes_1(x) := x / (log(x)-1)$
00108 approx_count_primes_1_hp(x,decimal_digits) := block([fpprec : decimal_digits],
00109  bfloat(x / (log(x)-1)))$
00110 
00111 /* The simplest approximation of unrank_primes(x) via the PNT: */
00112 approx_unrank_primes_0(x) := x * log(x)$
00113 
00114 /* The approximation of count_primes(x) via the logarithmic integral: */
00115 Li(x) := expintegral_li(x)$
00116 /* The floating point version, using bfloats with decimal digits: */
00117 Li_hp(x,decimal_digits) := block([fpprec : decimal_digits],
00118  expand(bfloat(Li(x))))$
00119 /* The higher logarithmic integrals (Li(x) = Lih(x,1)): */
00120 /* TODO: the integration should start with 0. */
00121 Lih(x,m) := integrate(log(t)^(-m),t,2,x)$
00122 Lih_hp(x,m,decimal_digits) := block([fpprec : decimal_digits],
00123  expand(bfloat(Lih(x,m))))$
00124 
00125 /* Accordingly we provide abbreviations for the exponential integral: */
00126 Ei(x) := expintegral_ei(x)$
00127 
00128 /* Intervals containing count_primes(x) */
00129 
00130 /* Using the most basic bounds around x/log(x) (where x can be real): */
00131 count_primes_int_0(x) :=
00132  if x<2 then [0,0] elseif x<3 then [1,1] elseif x<5 then [2,2]
00133  else [x/ld(x), 2*x/log(x)]$
00134 
00135 /* From the Riemann-conjecture it follows that the quotient 
00136      abs(c - Li(x)) / (sqrt(x) * log(x)
00137    is bounded. The following function computes the maximal quotient
00138    in the interval from a >= 2 to b (a, b must be integers), together with
00139    the x where this maximum is assumed:
00140 */
00141 max_quotient_riemann_error(a,b) := block(
00142  [c : count_primes(a-1), maxq : minf, maxx : undef],
00143   for x : a thru b do (
00144     if primep(x) then c:c+1,
00145     q : float(abs(c - Li(x)) / (sqrt(x) * log(x))),
00146     if q > maxq then (maxq : q, maxx : x)
00147   ),
00148   return([maxq, maxx]))$
00149 /* The Riemann conjecture is equivalent to the assertion that the result of
00150      max_quotient_riemann_error(2657,inf)
00151    is less then 1/(8 pi) ~ .03978873577297384.
00152    The corresponding interval (x here is an integer):
00153 */
00154 count_primes_int_riemann(x) :=
00155  if x < 2657 then block([c : count_primes(x)], [c,c])
00156  else block([d : 1/(8*%pi) * sqrt(x) * log(x)],
00157    Li(x) + [-d,+d])$
00158 /* In other words, the Riemann conjecture is equivalent to count_primes(x)
00159    for all x lying in the (open) intervall computed by 
00160    count_primes_int_riemann(x).
00161 */
00162 
00163 
00164 /* **************************
00165    * Additive number theory *
00166    **************************
00167 */
00168 
00169 /* In ComputerAlgebra/Hypergraphs/Lisp/Generators/GreenTao.mac we have
00170    arithprog_primes_finish[k,n].
00171 */
00172 
00173 /* The number of arithmetic progressions of length k in the set of
00174    prime numbers with the n-th prime number as last element ("nc"
00175    for "non-cumulative").
00176 */
00177 /* Prerequisite: n, k natural numbers, n >= 1, k >= 0. */
00178 n_arithprog_primes_nc[k,n] := 
00179  block([prp : product_primes(k), prk : primep(k), p : unrank_primes(n)],
00180   n_arithprog_primes_nc_[k,n])$
00181 
00182 /* Implementation functions (better not to be called directly, since wrong
00183    usage can introduce false results).
00184 */
00185 /* Inheriting prp, prk and p: */
00186 n_arithprog_primes_nc_[k,n] := n_arithprog_primes_nc_nm(k,n)$
00187 /* Inheriting again prp, prk and p: */
00188 n_arithprog_primes_nc_nm(k,n) :=
00189  if k=0 then 0 elseif k=1 then 1 elseif k=2 then n-1
00190  elseif n <= k then 0 else
00191   length(arithprog_primes_finish_nm(k,p))$
00192 
00193 /* The number of arithmetic progressions of length k in the set of
00194    prime numbers ("c" for "cumulative"):
00195 */
00196 n_arithprog_primes_c(k,n) :=
00197  if k <= 2 then binomial(n,k) else
00198  block([prp : product_primes(k), prk : primep(k), p : 1, sum : 0, count],
00199    for i : 1 thru n do (
00200      p : next_prime(p),
00201      count : if p <= k then 0 else n_arithprog_primes_nc_[k,i],
00202      sum : sum + count
00203    ),
00204    sum)$
00205 /* See ComputerAlgebra/RamseyTheory/Lisp/GreenTao/Asymptotics.mac
00206    for approximations of n_arithprog_primes_c(k,n).
00207 */
00208 
00209 
00210 /* The list of the numbers of arithmetic progressions of length k in 
00211    the set of the first i prime numbers for i from 1 to n ("c" for
00212    "cumulative"):
00213 */
00214 ln_arithprog_primes_c(k,n) :=
00215  if k <= 2 then create_list(binomial(i,k),i,1,n) else
00216  block([prp : product_primes(k), prk : primep(k), p : 1, res : [], sum : 0, count],
00217    for i : 1 thru n do (
00218      p : next_prime(p),
00219      count : if p <= k then 0 else n_arithprog_primes_nc_[k,i],
00220      sum : sum + count, res : cons(sum,res)
00221    ),
00222    reverse(res))$
00223 
00224 /* The smallest natural number n >= 0 such that amongst the first n
00225    prime numbers there exists an arithmetic progression of size k:
00226 */
00227 first_arithprog_primes(k) :=
00228  if k <= 2 then k else
00229   block([prp : product_primes(k), prk : primep(k), 
00230    n : rank_primes(k), p : if primep(k) then k else next_prime(k), 
00231    found : false],
00232    while not found do (
00233      p : next_prime(p), n : n + 1,
00234      found : not emptyp(arithprog_primes_finish[k,p])
00235    ),
00236    n)$
00237 
00238 
00239 /* *****************
00240    * Prime factors *
00241    *****************
00242 */
00243 
00244 /* The exponent of integer p <> 0,1 for integer n <> 0: */
00245 expfact(n,p) := if n=0 then 0 else block([e:0],
00246  while mod(n,p)=0 do (n:n/p, e:e+1), e)$
00247 
00248