OKlibrary  0.2.1.6
Asymptotics.mac
Go to the documentation of this file.
00001 /* Oliver Kullmann, 8.5.2010 (Swansea) */
00002 /* Copyright 2010 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/NumberTheory/Lisp/PrimeNumbers.mac")$
00023 
00024 /* *********************************************
00025    * The formula from [Grosswald, Hagis, 1979] *
00026    *********************************************
00027 */
00028 
00029 /* The task is to approximate
00030      nhyp_arithprog_primes_hg(k,n) = n_arithprog_primes_c(k,n),
00031    the number of arithmetic progressions of length k in the first n prime
00032    numbers.
00033 */
00034 
00035 /* First computing the coefficients C_k = C_gh(k,inf), defined as
00036      C_k = (prod_{p <= k} 1/p * (p/(p-1))^(k-1)) *
00037            (prod_{p > k} (p/(p-1))^(k-1) * (p-k+1)/p)
00038    where p ranges over the prime numbers.
00039 
00040    The first factor is computed by C_gh_fin(k), the second by
00041    C_gh_inf(k,max_p), where only p <= max_p is considered.
00042  */
00043 
00044 /* The "finite factor": */
00045 factor_C_gh_fin(k) := buildq([k],lambda([p], p^(k-2)/(p-1)^(k-1)))$
00046 
00047 C_gh_fin(k) := prod_l(map(factor_C_gh_fin(k), primes_first(count_primes(k))))$
00048 
00049 /* The "infinite factor": */
00050 
00051 factor_C_gh_inf(k) := buildq([k],lambda([p], p^(k-2)*(p-k+1)/(p-1)^(k-1)))$
00052 factorbf_C_gh_inf(k) := buildq([k],lambda([p], bfloat(p^(k-2)*(p-k+1)/(p-1)^(k-1))))$
00053 
00054 /* Prerequisite: max_p >= k, k >= 2. */
00055 C_gh_inf(k,max_p) := block(
00056  [L : rest(primes_first(count_primes(max_p)), count_primes(k))],
00057   prod_l(map(factor_C_gh_inf(k), L)))$
00058 /* Using floating-point arithmetic (with the current value of fpprec): */
00059 C_gh_inf_hp(k,max_p) := block(
00060  [L : rest(primes_first(count_primes(max_p)), count_primes(k))],
00061   if emptyp(L) then 1.0b0 else
00062   prod_l(map(factorbf_C_gh_inf(k), L)))$
00063 /* Upper bound on the absolute value of the relative error of 
00064    C_gh_inf_hp(k,max_p) (as compared to C_gh_inf(k,max_p)), where num_p
00065    is the number of primes considered.
00066    Prerequisite: eps = 10^(-fpprec+1), the "machine epsilon", fulfils
00067      eps <= min(2/3, 1/(8*num_p - 6)).
00068    (Since fpprec >= 16, for all feasible computations at Maxima-level this
00069    is automatically fulfilled.)
00070 */
00071 relerror_C_gh_inf_hp(k,num_p) := block([eps : 10^(-fpprec+1)],
00072   bfloat(eps * (1 + 4*(num_p-1))))$
00073 /* Returning a pair [a,b], representing an open interval containing the true
00074    value:
00075 */
00076 C_gh_inf_hpint(k,max_p) := block(
00077  [L : rest(primes_first(count_primes(max_p)), count_primes(k)), res, eps],
00078   res : if emptyp(L) then 1.0b0 else
00079     prod_l(map(factorbf_C_gh_inf(k), L)),
00080   eps : relerror_C_gh_inf_hp(k,length(L)),
00081   res*[1-eps,1+eps])$
00082 
00083 
00084 /* The final computations: */
00085 
00086 C_gh(k,max_p) := C_gh_fin(k) * C_gh_inf(k,max_p)$
00087 /* Using floating-point arithmetic: */
00088 C_gh_hp(k,max_p,decimal_digits) := block([fpprec : decimal_digits],
00089   C_gh_fin(k) * C_gh_inf_hp(k,max_p))$
00090 /* Returning a pair [a,b], representing an open interval containing the true
00091    value: */
00092 C_gh_hpint(k,max_p,decimal_digits) := block(
00093  [fpprec : decimal_digits,
00094   eps : 10^(-decimal_digits+1)],
00095   4*[-eps,+eps] + C_gh_fin(k) * C_gh_inf_hpint(k,max_p))$
00096 
00097 
00098 /* The current approximations of the C_k for k >= 2: */
00099 C_gh_values : [
00100  1,
00101  1.3203236316942413054,
00102  2.8582485957224816583,
00103  4.1511808632468886479,
00104  10.131794950034614014,
00105  17.298612311683576106,
00106  53.971948300560721615,
00107  148.55162866536732961,
00108  336.03432675383279702,
00109  511.42228206774875224,
00110  1312.3197113260945073,
00111  2364.5989633652830187,
00112  7820.6000304765722324,
00113  22938.908633119341404,
00114  55651.462555721561157,
00115  91555.111230322724999,
00116  256474.85986744035674,
00117  510992.01033894385383,
00118  1900972.5849978144616
00119 ]$
00120 
00121 /* The resulting estimation of the number of arithmetic progressions of
00122    length k in the first n prime numbers:
00123 */
00124 gh_nhyp_arithprog_primes_hg(k,n) := if k-1 > length(C_gh_values) then unknown
00125  else float(C_gh_values[k-1]/2/(k-1) * n^2/log(n)^(k-2))$
00126 
00127 /* Using the original formula: */
00128 ghorig_nhyp_arithprog_primes_hg(k,n) :=
00129  if k-1 > length(C_gh_values) then unknown
00130  else block([p : unrank_primes(n)],
00131  float(C_gh_values[k-1]/2/(k-1) * p^2/log(p)^k))$
00132 
00133 
00134 /* ********************************************
00135    * The approximation from [Granville, 2006] *
00136    ********************************************
00137 */
00138 
00139 /* The approximation of unrank_primes(greentao([k]) from
00140    [Granville, 2006] ("ur" stands for "unranked"): */
00141 approxgv_grt1ur(k) := (k/2*exp(1-%gamma))^(k/2)$
00142 /* The derived approximation of greentao([k]), using the approximation
00143    of count_primes via the logarithmic integral:
00144 */
00145 approxgv_grt1Li_hp(k,decimal_digits) := block([fpprec : decimal_digits],
00146  expand(bfloat(Li(approxgv_grt1ur(k)))))$
00147 
00148