a-semismooth.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Semismooth\n",
    "\n",
    "We use [Banks and Shparlinski](https://arxiv.org/pdf/math/0601460.pdf) to approximate the density of semismooth numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "from scipy.integrate import quad\n",
    "import matplotlib.pyplot as plt\n",
    "PRECISION=255"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best way we found to evaluate the [Buchstab function](https://en.wikipedia.org/wiki/Buchstab_function) was to tabulate it. The errors do not seem to have a noticeable impact on the final estimates.\n",
    "\n",
    "The table can be recomputed using the Julia script [`buchstab.jl`](buchstab.jl)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "buchstabtab = []\n",
    "with open('buchstab.tsv') as f:\n",
    "    for x, y in csv.reader(f, dialect=csv.excel_tab):\n",
    "        buchstabtab.append((float(x),float(y)))\n",
    "\n",
    "gamma = euler_gamma.n(PRECISION)\n",
    "buchstab_lim = e**(-euler_gamma).n(PRECISION)\n",
    "log2 = log(2).n(PRECISION)\n",
    "\n",
    "def buchstab(u):\n",
    "    '''Buchstab function'''\n",
    "    if u < 1:\n",
    "        return 0\n",
    "    elif u <= 2:\n",
    "        return 1/u\n",
    "    elif u < 7:\n",
    "        i, x1, y1 = next((i,x,y) for i, (x, y) in enumerate(buchstabtab) if u < x)\n",
    "        x0, y0 = buchstabtab[i-1]\n",
    "        return y0 + u*(y1 - y0)/(x1-x0)\n",
    "    else:\n",
    "        return buchstab_lim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"partial convolutions\" of the Buchstab function ω and the [Dickman function](https://en.wikipedia.org/wiki/Dickman_function) ρ, as defined by Banks & Sparlinski, evaluated by numerical integration\n",
    "\n",
    "$$C(u,v) = \\int_v^∞ ω(u-s)ρ(s)ds$$\n",
    "$$C'(u,v) = \\int_v^∞ ω(u-s)ρ'(s)ds$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def C(u,v):\n",
    "    integral, err = quad(lambda s: buchstab(u-s)*dickman_rho(s), v, u-1)\n",
    "    return integral\n",
    "\n",
    "def C1(u, v):\n",
    "    integral, err = quad(lambda s: -buchstab(u-s)*dickman_rho(s-1)/s, v, u-1)\n",
    "    return integral"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The smooth number density functions.\n",
    "\n",
    "$$Ψ(α,β) = ρ(α/β)$$\n",
    "\n",
    "is the classical density of $x^β$-smooth numbers in the interval $[1,x^α]$.\n",
    "\n",
    "$$Θ(α,β,γ) = ρ(α/β) + C(α/β, γ/β) - γ \\frac{C'(α/β, γ/β)}{β}$$\n",
    "\n",
    "is the function introduced by Banks & Sparlinksy, estimating the density of numbers in the interval $[1,e^α]$ that have a $e^β$-smooth divisor larger than $e^γ$.\n",
    "\n",
    "For convenience, the latter is modified to accept exponents in base $2$ instead of base $e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Psi(a, b):\n",
    "    return dickman_rho(a/b)\n",
    "\n",
    "def Theta(a, b, c):\n",
    "    return dickman_rho(a/b) + C(a/b, c/b) - gamma/log2 * C1(a/b, c/b) / b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function $Θ(a,b,c)$ above is the same as the semi-smooth counting function $ψ(a, a-c, b)$.\n",
    "\n",
    "We use it to approximate the limit $G(1/u,2/u)$ in Drew's thesis: `u` just what you'd expect; `x` is a large constant: the larger, the more accurate the approximation... you may want to play a bit with it, as I'm not sure how fast the limit converges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def G(u, x = 10^20):\n",
    "    return Theta(x, x/u, x-2*x/u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check that the formula matches the table in Drew's thesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(2.10000000000000, 0.948866710938980),\n",
       " (2.20000000000000, 0.895875321922736),\n",
       " (2.30000000000000, 0.841529473096127),\n",
       " (2.40000000000000, 0.786218441980182),\n",
       " (2.50000000000000, 0.730246154979384),\n",
       " (2.60000000000000, 0.673851847528524),\n",
       " (2.70000000000000, 0.617225068815002),\n",
       " (2.80000000000000, 0.560516765435068),\n",
       " (2.90000000000000, 0.503847590588720)]"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[(i/10., G(i/10)) for i in range(21, 30)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.447314214932594,\n",
       " 0.223819493954199,\n",
       " 0.0963990059356955,\n",
       " 0.0365730650779422,\n",
       " 0.0124134827486074,\n",
       " 0.00109226674235895,\n",
       " 3.66265164509229e-6,\n",
       " 5.38286135968268e-9,\n",
       " 4.25491179975736e-12,\n",
       " 6.53407244918750e-19,\n",
       " 2.41550417023604e-26]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[G(i) for i in [3, 3.5, 4, 4.5, 5, 6, 8, 10, 12, 16, 20]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we finally check the values you extrapolated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-104.671702604978"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log(G(23), 2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-101.353787485191"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log(G(22.5), 2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-125.002622333603"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log(G(26), 2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-128.457523309872"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log(G(26.5), 2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So your extrapolation looks accurate, unsurprisingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.1",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

buchstab.jl

using DifferentialEquations, ArgParse

function buchstab_model(dω,ω,h,p,u)
    dω[1] = ( h(p, u-1)[1] - ω[1] ) / u
end

function buchstabtab(up, res=1, prec=40)
    setprecision(prec) do
        buchstab_history(p, u) = 1/u
        buchstab_initial = [1/BigFloat(2)]
        buchstab_problem = DDEProblem(buchstab_model, buchstab_initial, buchstab_history,
                                      (BigFloat(2), BigFloat(up)+3); constant_lags=[1])
        alg = MethodOfSteps(Tsit5())
        buchstab = solve(buchstab_problem, alg)
        return [(x, buchstab(x)[1]) for x in (BigFloat(2):BigFloat(res):BigFloat(up))]
    end
end

s = ArgParseSettings()
@add_arg_table s begin
    "--resolution", "-r"
        help = "Evaluate at points spaced by this resolution"
        arg_type = BigFloat
        default = 1/BigFloat(10)
    "--precision", "-p"
        help = "Use this much floating point precision"
        arg_type = Int
        default = 40
    "bound"
        help = "Evaluate up to bound"
        arg_type = Int
        required = true
end

args = parse_args(ARGS, s)
prec = args["precision"]
table = buchstabtab(args["bound"], args["resolution"], prec)
for (x,bx) in table
    print(x, "\t", bx, "\n")
end

buchstab.tsv