Square-cube-free-integers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
   ],
   "source": [
    "#Functions used below\n",
    "def rup(x): #Rounds a real number x up with 4 decimal digits of precision.\n",
    "    return ceil(10000*x)/10000\n",
    "\n",
    "def rdn(x): #Rounds a real number x down with 4 decimal digits of precision.\n",
    "    return floor(10000*x)/10000\n",
    "\n",
    "def rupstr(x): #Returns a string containing the rounded up digits of x (5 digits)\n",
    "    d=floor(log(abs(x),10))\n",
    "    return str(N(rup(x),digits=max(d+5,4)))\n",
    "\n",
    "def rdnstr(x): #Returns a string containing the rounded up digits of x (5 digits)\n",
    "    d=floor(log(abs(x),10))\n",
    "    return str(N(rdn(x),digits=max(d+5,4)))\n",
    "\n",
    "def ispos(x):\n",
    "    if x>0:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0\n",
    "    \n",
    "sigma2(x,y) = y/log(x*y)*(1+1.5/log(x*y))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7934 1.2604 1.3860 1.4430\n"
     ]
    }
   ],
   "source": [
    "# Section 3\n",
    "# Lemma 5\n",
    "cL5 = rdn((0.4995)^(1/3))\n",
    "# Corollary 2\n",
    "cC2 = rup((0.4995)^(-1/3))\n",
    "# Lemma 6\n",
    "cL6 = rdn((1.997*4/3)^(1/3))\n",
    "# Corollary 3\n",
    "cC3 = 2*rup((1.997*4/3)^(-1/3))\n",
    "# Lemma 7\n",
    "cL7(x) = (2/5) - 1/(20*x)\n",
    "# Corollary 4\n",
    "cC4(x) = 2/cL7(x)^(1/5)\n",
    "print(rdnstr(cL5), rupstr(cC2), rdnstr(cL6), rupstr(cC3))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lemma 8 constant: 0.6600\n",
      "Lemma 9 constants: 1.3401 0.2900\n",
      "Equations (4.30) and (4.31) are not mutually exclusive.\n",
      "Constants appearing in Lemmas 10-12: 0.09990 0.06640 2.8467 0.1170\n"
     ]
    }
   ],
   "source": [
    "# Section 4\n",
    "# Lemma 8\n",
    "m = 5            # lower bound on m\n",
    "l = 1.05         # upper bound on lambda\n",
    "if m < 3.005*l^5/4:\n",
    "    print(\"Failed contradiction!\")\n",
    "cL8 = rdn((2 - 2*(l-1)/m)/3)\n",
    "print(\"Lemma 8 constant:\",rdnstr(cL8))\n",
    "# Lemma 9\n",
    "cL9i = rup(l^6)\n",
    "cL9ii = rdn(0.5-l/m)\n",
    "# Check that (4.5) and (4.6) are mutually exclusive as claimed\n",
    "if m*cL9ii < cL9i:\n",
    "        print(\"Equations (4.5) and (4.6) are not mutually exclusive.\")\n",
    "print(\"Lemma 9 constants:\", rupstr(cL9i), rupstr(cL9ii))\n",
    "# Lemma 10\n",
    "m = 3\n",
    "l = 1.04\n",
    "if 9.99*m < (3*l-1)*l^6:\n",
    "    print(\"Failed contradiction!\")\n",
    "if (3*l-1)/m^2 < 1:\n",
    "    cL10 = 0.0999\n",
    "# Lemma 11\n",
    "m = 8\n",
    "if rdn(24/l^7)*cL10*2 < (6*l+0.11)/m:\n",
    "    print(\"Failed contradiction!\")\n",
    "cL11(x,y) = rdn((2-(6*x+0.11)/(100*y))/30)\n",
    "# Lemma 12\n",
    "m = 10\n",
    "cL12i(x) = rup(2*(x)^9)\n",
    "cL12ii(x,y) = rdn(1/3-2*x^2/y)\n",
    "# Check whether (4.30) and (4.31) are mutually exclusive\n",
    "# The choice in the paper should output a message\n",
    "if m*cL12ii(l,m) < cL12i(l):\n",
    "    print(\"Equations (4.30) and (4.31) are not mutually exclusive.\")\n",
    "print(\"Constants appearing in Lemmas 10-12:\", rdnstr(cL10), rdnstr(cL11(l,m)), rupstr(cL12i(l)), rdnstr(cL12ii(l,m)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Constant in Eq. 5.1: 1.3860\n",
      "Lower bound for b: 0.8706\n",
      "Eq 5.12: 1.4959\n",
      "Eq 5.13: 0.1452\n",
      "Following display: 0.2173\n",
      "Eq 5.14: 2.9918 0.2904 0.4346\n",
      "Lower bound for B: 10.2850 1193\n",
      "Eq 5.20: 0.01540 0.0001000 0.2649\n",
      "Lower Bound for A: 11.8501\n",
      "Eq 5.22: 0.4083 0.2370\n",
      "Eq 5.4: 0.5298 0.3400 0.03080\n",
      "Eq 5.5: 1.2105 1.4182\n"
     ]
    }
   ],
   "source": [
    "# Section 5.1\n",
    "C2 = 11           # Main constant in the definition of h\n",
    "# Proposition 1\n",
    "l = 1045/1000         # Value of lambda\n",
    "m = 55/10           # Value of m\n",
    "k0 = 116\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "MC1 = 5           # Coefficient of split point in (5.5)\n",
    "d1 = 17/100         # Value of delta in (5.11)\n",
    "# From Sect. 3\n",
    "cL6 = rdn((1.997*4/3)^(1/3))\n",
    "# From Sect. 4, updated\n",
    "cL8a = rdn(cL8^(1/3))\n",
    "cL9ii = rdn(1/2-l/m)\n",
    "print(\"Constant in Eq. 5.1:\", rdnstr(cL6))\n",
    "print(\"Lower bound for b:\",rdnstr(cL8a))\n",
    "#print(rupstr(cL9i), rupstr(cL9ii))\n",
    "\n",
    "# Computed constants appearing in the proof of Proposition 1:\n",
    "cP1512 = rup(l^6/cL8a)             # Eq 5.12\n",
    "cP1513 = rup((l-1)/cL9ii)          # Eq 5.13\n",
    "cP151x = rup(cP1512*cP1513)        \n",
    "print(\"Eq 5.12:\", rupstr(cP1512))\n",
    "print(\"Eq 5.13:\", rupstr(cP1513))\n",
    "print(\"Following display:\", rupstr(cP151x)) \n",
    "print(\"Eq 5.14:\", rupstr(2*cP1512), rupstr(2*cP1513), rupstr(2*cP151x))             # Eq 5.14\n",
    "bB = floor(m*d1*C2*k0)                          # Lower bound on B\n",
    "print(\"Lower bound for B:\",rdnstr(m*d1*C2), bB)\n",
    "\n",
    "cP1517 = rup(2*cP151x*(1002/1000*d1)^(4/3)/(8/3))    # Eq 5.17\n",
    "cP1518 = rup(2*cP1513*(1002/1000*d1)^4/8)                 # Eq 5.18\n",
    "cP1520 = rup((l-1)/d1+cP1518+1.5*x0^(-1/5))\n",
    "print(\"Eq 5.20:\", rupstr(cP1517), rupstr(cP1518), rupstr(cP1520))                                 # Eq 5.20 \n",
    "\n",
    "cP1521 = 2023/10000                                               # Eq 5.21\n",
    "bA = rdn(cL6*MC1^(4/3))                                       # Lower bound on A\n",
    "cP1522a = rup(0.6*2^(-8/3)*(1/2-1/bA)^(-5/3))\n",
    "cP1522b = rup(cP1522a*1.386^(-5/3))\n",
    "print(\"Lower Bound for A:\",rdnstr(bA))\n",
    "print(\"Eq 5.22:\",rupstr(cP1522a), rupstr(cP1522b))  \n",
    "                                                              # Eq 5.22\n",
    "\n",
    "print(\"Eq 5.4:\",rupstr(2*cP1520), rupstr(2*d1), rupstr(2*cP1517))                  # Eq 5.4\n",
    "cP155a = rup(4*cP1512*cP1521)\n",
    "cP155b = rup(4*cP1512*cP1522b)\n",
    "print(\"Eq 5.5:\", rupstr(cP155a), rupstr(cP155b))              # Eq 5.5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower Bound on log(x) for small H: 150\n",
      "Small H, Small M: 0.1034 0.6019 -89.7886 0.1582\n",
      "Small H, Large M: 0.1148 1.8055 -10.9463 0.2790\n",
      "Large H: \t  0.2378 2.4073 -98.1700 0.3976\n",
      "Eq 6.3: 0.09010 0.0009000\n",
      "Sigma_0: -0.05950\n",
      "Sigma_2 (Large H): 0.1797\n",
      "Sigma_3 (Large H): 0.3985\n",
      "Total (Large H): 0.9710\n",
      "Sigma_2 (Small H): 0.1461\n",
      "Sigma_3 (Small H): 0.4381\n",
      "Total (Small H): 0.9770\n"
     ]
    }
   ],
   "source": [
    "# Section 6.1\n",
    "C3 = 11           # Main constant in the definition of h\n",
    "# Proof of Theorem 1\n",
    "l = 1045/1000         # Value of lambda\n",
    "m = 55/10           # Value of m\n",
    "MC1 = 5           # Coefficient of split point in (5.5)\n",
    "d1 = 17/100         # Value of delta in (5.11)\n",
    "  \n",
    "#k1 = 150\n",
    "k0 = 116\n",
    "x0 = exp(k0)      # Regular lower bound on x\n",
    "\n",
    "\n",
    "f1(x)=x^(1/20)-(m*C3/MC1)*log(x)\n",
    "k1=floor(find_root(f1(e^x),100,300))\n",
    "print(\"Lower Bound on log(x) for small H:\",k1)\n",
    "x1 = exp(k1)      # Enhanced lower bound on x\n",
    "k1u = k1+1          # Upper bound on k1 for small H\n",
    "\n",
    "cT161a = cP155a*MC1^(7/3)*x1^(-1/12)/(1-l^(-7/3))\n",
    "cT161b = 2*cP1517*(m*C3*k1)^(-1/3)/(1-l^(-1/3))\n",
    "cT161c = rup(2*cP1520/(20*ln(l)))\n",
    "cT161d = (2*MC1*d1*x1^(-3/20)/(1-l^(-1)))\n",
    "cT161e = ((2*cP1520)*(1-log(m*C3*k1/MC1)/log(l)))\n",
    "cT161I = rup(cT161a+cT161b)+rup(cT161c)/C3\n",
    "print(\"Small H, Small M:\",rupstr(cT161a+cT161b), rupstr(cT161c), rupstr(cT161d+cT161e),rupstr(cT161I))\n",
    "if cT161d + cT161e >0:\n",
    "    print(\"Error!\")\n",
    "\n",
    "cT161f = cP155b*x1^(-1/15)/(1-l^(-1/9))\n",
    "cT161g = 2*cP1517*MC1^(-1/3)*x1^(-1/60)/(1-l^(-1/3))\n",
    "cT161h = rup(3*2*cP1520/(20*ln(l)))\n",
    "cT161j = 2*d1/(1-l^(-1))\n",
    "cT161k = ((2*cP1520)*(1-log(MC1)/log(l)))\n",
    "cT161II = rup(rup(cT161f+cT161g)+cT161h/C3)\n",
    "print(\"Small H, Large M:\", rupstr(cT161f+cT161g), rupstr(cT161h), rupstr(cT161j+cT161k),rupstr(cT161II))\n",
    "if cT161j + cT161k >0:\n",
    "    print(\"Error!\")\n",
    "\n",
    "\n",
    "cT162a = cP155b*x0^(-1/15)/(1-l^(-1/9))\n",
    "cT162b = 2*cP1517*(m*C3*k0)^(-1/3)/(1-l^(-1/3))\n",
    "cT162c = rup(2*cP1520/(5*ln(l)))\n",
    "cT162d = 2*d1/(1-l^(-1))\n",
    "cT162e = (2*cP1520)*(1-log(m*C3*log(x0))/log(l))\n",
    "cT162f = 10*rup((cT162d+cT162e)/10)\n",
    "cT162 = rup(rup(cT162a+cT162b)+cT162c/C3) + rup(cT162f/k1u/C3)\n",
    "print(\"Large H: \\t \",rupstr(cT162a+cT162b), rupstr(cT162c), rupstr(cT162f), rupstr(cT162))\n",
    "\n",
    "# From Sect. 3\n",
    "cL4 = 4995/10000\n",
    "cC1 = rup((l-1)/cL4)\n",
    "cT163a = cC1/(1-l^(-2))\n",
    "cT163b = x0^(-1/5)/log(l)*(1/10+(1/2*log(2)+log(l))/k0)    #Double check?\n",
    "cT163 = rup((cT163a/k0+cT163b)/C3)\n",
    "print(\"Eq 6.3:\",rupstr(cC1),rupstr(cT163))\n",
    "\n",
    "#From Sect. 2\n",
    "bJ = 120\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C3*k0*x0^(1/5)\n",
    "sigma0 = rup(1 - prod([1 - 1/p^2 for p in prList]) - sum([1/p^2 for p in prList]) + 2^bK/h0)\n",
    "sigma1 = 4523/10000\n",
    "\n",
    "print(\"Sigma_0:\",rupstr(sigma0))\n",
    "print(\"Sigma_2 (Large H):\",rupstr(sigma2(h0,m)))\n",
    "print(\"Sigma_3 (Large H):\",rupstr(cT162 + cT163))\n",
    "print(\"Total (Large H):\",rupstr(sigma0+sigma1+rup(sigma2(h0,m)) + rup(cT162 + cT163)))\n",
    "\n",
    "h1 = C3*k1*x1^(1/5)\n",
    "\n",
    "print(\"Sigma_2 (Small H):\",rupstr(sigma2(h1,m)))\n",
    "print(\"Sigma_3 (Small H):\",rupstr(cT161I + cT161II + cT163))\n",
    "print(\"Total (Small H):\",rupstr(sigma0+sigma1+rup(sigma2(h1,m)) + rup(cT161I + cT161II + cT163)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S2 Bound: 0.4272\n",
      "Sigma_3: 0.4331\n",
      "Sigma_0: -0.05430\n",
      "Sigma_2: 0.1580\n",
      "Total: 0.9891\n",
      "This constant proves Theorem 1 from e^41 until e^109.7200\n"
     ]
    }
   ],
   "source": [
    "# Section 6.2\n",
    "# Proposition 3\n",
    "l = 1025/1000         # Value of lambda\n",
    "m = 175/100          # Value of m\n",
    "k0 = 41\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "C2 = 5            \n",
    "C3 = 11\n",
    "bJ = 20\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C2*x0^(1/4)\n",
    "\n",
    "# From Sect. 3\n",
    "cL6 = rdn((1997/1000*4/3)^(1/3))\n",
    "cC3 = 2*rup((1997/1000*4/3)^(-1/3))\n",
    "cP31 = rup(cC3*(l-1)/(1-l^(-1/3))/(m*C2^4)^(1/3))\n",
    "cP32 = rup(k0/(2*log(l)*h0))\n",
    "print(\"S2 Bound:\",rupstr(cP31))\n",
    "print(\"Sigma_3:\",rupstr(cP31+cP32))\n",
    "\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^2 for p in prList]) - sum([1/p^2 for p in prList]) + 2^bK/h0))\n",
    "print(\"Sigma_0:\",rupstr(sigma0))\n",
    "print(\"Sigma_2:\",rupstr(sigma2(h0,m)))\n",
    "Tot=rup(sigma0)+sigma1+rup(sigma2(h0,m) + cP31+cP32)\n",
    "\n",
    "print(\"Total:\",rupstr(Tot))\n",
    "km1p2=floor(100*(find_root((e^x)^(1/20)*x^(-1)-C3/C2,40,120)))/100\n",
    "if Tot<1:\n",
    "    print(\"This constant proves Theorem 1 from e^\"+str(k0)+\" until e^\"+rdnstr(km1p2))\n",
    "else:\n",
    "    print(\"Error!\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S2 Bound: 0.4423\n",
      "Sigma_3: 0.4424\n",
      "Sigma_0: -0.05940\n",
      "Sigma_2: 0.1571\n",
      "Total: 0.9924\n",
      "This constant proves Theorem 1 from e^109 until e^116.3900\n"
     ]
    }
   ],
   "source": [
    "# Proposition 4\n",
    "l = 10001/10000        # Value of lambda\n",
    "m = 45/10           # Value of m\n",
    "k0 = 109\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "C2 = 38/10\n",
    "C3 = 11\n",
    "bJ = 100\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C2*x0^(1/4)\n",
    "\n",
    "# From Sect. 3\n",
    "cL6 = rdn((1997/1000*4/3)^(1/3))\n",
    "cC3 = 2*rup((1997/1000*4/3)^(-1/3))\n",
    "cP31 = rup(cC3*(l-1)/(1-l^(-1/3))/(m*C2^4)^(1/3))\n",
    "cP32 = rup(k0/(2*log(l)*h0))\n",
    "print(\"S2 Bound:\",rupstr(cP31))\n",
    "print(\"Sigma_3:\",rupstr(cP31+cP32))\n",
    "\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^2 for p in prList]) - sum([1/p^2 for p in prList]) + 2^bK/h0))\n",
    "print(\"Sigma_0:\",rupstr(sigma0))\n",
    "print(\"Sigma_2:\",rupstr(sigma2(h0,m)))\n",
    "Tot=rup(sigma0)+sigma1+rup(sigma2(h0,m) + cP31+cP32)\n",
    "\n",
    "print(\"Total:\",rupstr(Tot))\n",
    "km1p2=floor(100*(find_root((e^x)^(1/20)*x^(-1)-C3/C2,40,120)))/100\n",
    "if Tot<1:\n",
    "    print(\"This constant proves Theorem 1 from e^\"+str(k0)+\" until e^\"+rdnstr(km1p2))\n",
    "else:\n",
    "    print(\"Error!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T_3 \\cap I Bound: 0.06640 7.0298\n",
      "Num Intervals Bound: 0.2929\n",
      "Eq 5.24 2.0591\n",
      "B lower bound: 2090\n",
      "Eq 5.26: 0.01870\n",
      "Eq 5.27: 1.15521866233743e-6\n",
      "Eq 5.29: 0.2106\n",
      "A lower bound: 0.8306 26.6513\n",
      "Eq 5.31: 0.2107 0.3677\n",
      "Eq 5.7: 0.4212 0.1900 0.03740\n",
      "Eq 5.8: 0.9519 5.1698\n"
     ]
    }
   ],
   "source": [
    "# Section 5.2\n",
    "C3 = 5            # Main constant in the definition of h\n",
    "# Proposition 2\n",
    "l = 104/100          # Value of lambda\n",
    "m = 11            # Value of m\n",
    "k0 = 200\n",
    "x0 = exp(k0)      # Enhanced lower bound on x\n",
    "MC2 = 18          # Coefficient of split point in (5.8)\n",
    "d2 = 19/100         # Value of delta in (5.23)\n",
    "# From Sect. 3\n",
    "cL7 = rdn((0.4 - 0.05/m)^(1/5))\n",
    "# From Sect. 4, updated\n",
    "# cL11 = 0.0665\n",
    "cL11a = cL11(l,m)^(1/3)\n",
    "\n",
    "# Computed constants appearing in the proof of Proposition 2:\n",
    "cP2524i = rup(2*l^9/cL11a)\n",
    "cP2524ii = rup((l-1)/cL12ii(l,m))\n",
    "cP2524x = rup(cP2524i*cP2524ii)\n",
    "print(\"T_3 \\\\cap I Bound:\",rdnstr(cL11(l,m)),rupstr(cP2524i))\n",
    "print(\"Num Intervals Bound:\", rupstr(cP2524ii))                             # Eq 5.24\n",
    "print(\"Eq 5.24\", rupstr(cP2524x))\n",
    "bB = floor(m*d2*C3*k0)                                        # Lower bound on B\n",
    "print(\"B lower bound:\",bB)\n",
    "\n",
    "cP2526 = rup(cP2524x*(1.001*d2)^2/4)             # Eq 5.26\n",
    "cP2527 = cP2524ii*(1.001*d2)^6/12                # Eq 5.27\n",
    "cP2529 = rup((l-1)/d2+cP2527+1.5*x0^(-1/7))\n",
    "print(\"Eq 5.26:\",rupstr(cP2526))\n",
    "print(\"Eq 5.27:\",n(cP2527))\n",
    "print(\"Eq 5.29:\",rupstr(cP2529))                                 # Eq 5.29 \n",
    "\n",
    "cP2530 = 0.0677                                               # Eq 5.30\n",
    "bA = rdn(cL7*MC2^(6/5))                                       # Lower bound on A\n",
    "cP2531a = rup((1/48)*(1/2-1/bA)^(-3))\n",
    "cP2531b = rup(cP2531a*cL7^(-3))\n",
    "print(\"A lower bound:\", rdnstr(cL7),rdnstr(bA))\n",
    "print(\"Eq 5.31:\",rupstr(cP2531a), rupstr(cP2531b))  \n",
    "                                                              # Eq 5.31\n",
    "\n",
    "print(\"Eq 5.7:\", rupstr(2*cP2529), rupstr(d2), rupstr(2*cP2526))          # Eq 5.7 \n",
    "cP258a = rup(2*cP2524i*cP2530)\n",
    "cP258b = rup(2*cP2524i*cP2531b)\n",
    "print(\"Eq 5.8:\",rupstr(cP258a), rupstr(cP258b))\n",
    "                                                              # Eq 5.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "K1: 284\n",
      "Small H, Small M: 0.1552 0.2557 -72.2396 0.2064\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Small H, Large M: 0.1180 1.2785 -25.6791 0.3737\n",
      "Large H:\t  0.2742 1.5342 -94.5742 0.5148\n",
      "Eq 7.3: 0.1202 0.001100\n",
      "Sigma_3 (Large H): 0.5159\n",
      "Sigma_3 (Small H): 0.5812\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2 (Large H): 0.3020\n",
      "Sigma_2 (Small H): 0.2256\n",
      "Total (Large H): 0.9861\n",
      "Total (Small H): 0.9750\n"
     ]
    }
   ],
   "source": [
    "# Section 7.1 \n",
    "C3 = 5           # Main constant in the definition of h\n",
    "# Proof of Theorem 2\n",
    "l = 104/100          # Value of lambda\n",
    "m = 11            # Value of m\n",
    "k0 = 200\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "MC2 = 18          # Coefficient of split point in (5.8)\n",
    "d2 = 19/100         # Value of delta in (5.23)\n",
    "\n",
    "f3(x)=x^(1/42)-(m*C3/MC2)*log(x)\n",
    "k1=floor(find_root(f3(e^x),100,600))\n",
    "#k1 = 287\n",
    "print(\"K1:\",k1)\n",
    "x1 = exp(k1)      # Lower bound on x\n",
    "k1u = k1+1          # Upper bound on x in this case\n",
    "\n",
    "cT271a = cP258a*MC2^(11/3)*x1^(-1/18)/(1-l^(-11/3))\n",
    "cT271b = 2*cP2526*(m*C3*k1)^(-1/3)/(1-l^(-1/3))\n",
    "cT271c = rup(2*cP2529/(42*ln(l)))\n",
    "cT271d = (2*MC2*d2*x1^(-5/42)/(1-l^(-1)))\n",
    "cT271e = ((2*cP2529)*(1-log(m*C3*k1/MC2)/log(l)))\n",
    "cT271I = rup(rup(cT271a+cT271b)+cT271c/C3)\n",
    "print(\"Small H, Small M:\",rupstr(cT271a+cT271b), rupstr(cT271c), rupstr(cT271d+cT271e), rupstr(cT271I))\n",
    "if cT271d+cT271e >0:\n",
    "    print(\"Error!\")\n",
    "\n",
    "cT271f = cP258b*x1^(-1/21)/(1-l^(-1/15))\n",
    "cT271g = 2*cP2526*MC2^(-1/3)*x1^(-1/126)/(1-l^(-1/3))\n",
    "cT271h = rup(5*2*cP2529/(42*ln(l)))\n",
    "cT271j = d2/(1-l^(-1))\n",
    "cT271k = (2*cP2529)*(1-log(MC2)/log(l))\n",
    "cT271II = rup(rup(cT271f+cT271g)+cT271h/C3)\n",
    "cT271 = cT271I + cT271II\n",
    "print(\"Small H, Large M:\",rupstr(cT271f+cT271g), rupstr(cT271h), rupstr(cT271j+cT271k), rupstr(cT271II))\n",
    "if cT271j+cT271k >0:\n",
    "    print(\"Error!\")\n",
    "\n",
    "\n",
    "cT272a = cP258b*x0^(-1/21)/(1-l^(-1/15))\n",
    "cT272b = 2*cP2526*(m*C3*k0)^(-1/3)/(1-l^(-1/3))\n",
    "cT272c = rup(2*cP2529/(7*ln(l)))\n",
    "cT272d = d2/(1-l^(-1))\n",
    "cT272e = (2*cP2529)*(1-log(m*C3*k0)/log(l))\n",
    "cT272f = cT272d+cT272e\n",
    "cT272 = rup(rup(cT272a+cT272b)+cT272c/C3) + rup(cT272f/k1u/C3)\n",
    "print(\"Large H:\\t \",rupstr(cT272a+cT272b), rupstr(cT272c), rupstr(cT272f), rupstr(cT272))\n",
    "\n",
    "# From Sect. 3\n",
    "cL4 = 333/1000\n",
    "cC1 = rup((l-1)/cL4)\n",
    "cT273a = cC1/(1-l^(-3))\n",
    "cT273b = x0^(-1/7)/(3*ln(l))*(1/7+log(2*l^3)/k0)\n",
    "cT273 = rup((cT273a/k0+cT273b)/C3)\n",
    "print(\"Eq 7.3:\",rupstr(cC1), rupstr(cT273))\n",
    "\n",
    "print(\"Sigma_3 (Large H):\",rupstr(cT272 + cT273))\n",
    "print(\"Sigma_3 (Small H):\",rupstr(cT271 + cT273))\n",
    "\n",
    "#From Sect. 2\n",
    "bJ = 100\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C3*k0*x0^(1/7)\n",
    "h1 = C3*k1*x1^(1/7)\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^3 for p in prList]) - sum([1/p^3 for p in prList]) + 2^bK/h0))\n",
    "sigma1 = 1748/10000\n",
    "print(\"Sigma_0:\",rupstr(sigma0))\n",
    "print(\"Sigma_2 (Large H):\",rupstr(sigma2(h0,m)))\n",
    "print(\"Sigma_2 (Small H):\", rupstr(sigma2(h1,m)))\n",
    "print(\"Total (Large H):\",rupstr(sigma0+sigma1+ rup(sigma2(h0,m)) + rup(cT272 + cT273)))\n",
    "print(\"Total (Small H):\",rupstr(sigma0+sigma1+ rup(sigma2(h1,m)) + rup(cT271 + cT273)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S Bound: 0.3225\n",
      "Sigma_3: 0.3354\n",
      "Sigma_0: -0.004700\n",
      "Sigma_2: 0.3146\n",
      "Total: 0.8201\n",
      "This constant proves the theorem from e^41 until e^95.8900\n"
     ]
    }
   ],
   "source": [
    "Cmt2=5\n",
    "# Proposition 5.1\n",
    "l = 106/100         # Value of lambda\n",
    "m = 27/10          # Value of m\n",
    "k0 = 41\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "C3 = 2            \n",
    "bJ = 6\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C3*x0^(1/5)\n",
    "\n",
    "# From Sect. 3\n",
    "cL5 = rdn((4995/10000)^(1/3))\n",
    "cC2 = rup((4995/10000)^(-1/3))\n",
    "cP51 = rup(cC2*(l-1)/(1-l^(-2/3))/(m*C3)^(2/3)/C3)\n",
    "cP52 = 2*k0/(15*log(l)*h0)\n",
    "print(\"S Bound:\",rupstr(cP51)) \n",
    "print(\"Sigma_3:\",rupstr(cP51+cP52))\n",
    "\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^3 for p in prList]) - sum([1/p^3 for p in prList]) + 2^bK/h0))\n",
    "print(\"Sigma_0:\",rupstr(sigma0))\n",
    "print(\"Sigma_2:\", rupstr(sigma2(h0,m)))\n",
    "Tot=rup(sigma0)+sigma1+rup(sigma2(h0,m)) + cP51+cP52\n",
    "print(\"Total:\",rupstr(Tot))\n",
    "km0p2=floor(100*(find_root((e^x)^(2/35)*x^(-1)-Cmt2/C3,1,500)))/100\n",
    "if Tot<1:\n",
    "    print(\"This constant proves the theorem from e^\"+str(k0)+\" until e^\"+rdnstr(km0p2))\n",
    "else:\n",
    "    print(\"Error!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficent of (xH^-1)^(1/5): 0.07260\n",
      "Sigma_3: 0.5891\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2: 0.2207\n",
      "Total: 0.9780\n",
      "This constant proves the theorem from e^95 until e^191.6100\n"
     ]
    }
   ],
   "source": [
    "# Proposition 6 part 1\n",
    "l = 103/100       # Value of lambda\n",
    "m = 4           # Value of m\n",
    "k0 = 95\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "Cmt2 = 5           # Constant used in Theorem 2\n",
    "C3 = 10            \n",
    "bJ = 40\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C3*x0^(1/6)\n",
    "\n",
    "# From Sect. 3\n",
    "cL7(x) = (2/5) - 1/(20*x)\n",
    "cP60 = rup(2*cL7(m)^(-1/5)*(l-1))\n",
    "cP61 = rup(cP60/(1-l^(-1/5))/(m*C3)^(1/5)/C3)\n",
    "cP62 = rup(2*k0/(6*log(l)*h0))\n",
    "print(\"Coefficent of (xH^-1)^(1/5):\",rupstr(cP60))\n",
    "print(\"Sigma_3:\", rupstr(cP61+cP62))\n",
    "\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^3 for p in prList]) - sum([1/p^3 for p in prList]) + 2^bK/h0))\n",
    "print(\"Sigma_0:\", rupstr(sigma0))\n",
    "print(\"Sigma_2:\", rupstr(sigma2(h0,m)))\n",
    "Tot=sigma0+sigma1+sigma2(h0,m) + cP61+cP62\n",
    "print(\"Total:\",rupstr(Tot))\n",
    "km1p2=floor(100*(find_root((e^x)^(1/42)*x^(-1)-Cmt2/C3,100,700)))/100\n",
    "if Tot<1:\n",
    "    print(\"This constant proves the theorem from e^\"+str(k0)+\" until e^\"+rdnstr(km1p2))\n",
    "else:\n",
    "    print(\"Error!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coefficent of (xH^-1)^(1/5): 0.02420\n",
      "Sigma_3: 0.6767\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2: 0.1465\n",
      "Total: 0.9914\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This constant proves the theorem from e^191 until e^200.3000\n"
     ]
    }
   ],
   "source": [
    "# Proposition 6 part 2\n",
    "l = 101/100       # Value of lambda\n",
    "m = 5           # Value of m\n",
    "k0 = 191\n",
    "x0 = exp(k0)      # Lower bound on x\n",
    "Cmt2 = 5           # Constant used in Theorem 2\n",
    "C3 = 85/10            \n",
    "bJ = 40\n",
    "bK = prime_pi(bJ)\n",
    "prList = prime_range(1,bJ)\n",
    "h0 = C3*x0^(1/6)\n",
    "\n",
    "# From Sect. 3\n",
    "cL7(x) = (2/5) - 1/(20*x)\n",
    "cP60 = rup(2*cL7(m)^(-1/5)*(l-1))\n",
    "cP61 = rup(cP60/(1-l^(-1/5))/(m*C3)^(1/5)/C3)\n",
    "cP62 = rup(2*k0/(6*log(l)*h0))\n",
    "print(\"Coefficent of (xH^-1)^(1/5):\",rupstr(cP60))\n",
    "print(\"Sigma_3:\", rupstr(cP61+cP62))\n",
    "\n",
    "sigma0 = n(rup(1 - prod([1 - 1/p^3 for p in prList]) - sum([1/p^3 for p in prList]) + 2^bK/h0))\n",
    "print(\"Sigma_0:\", rupstr(sigma0))\n",
    "print(\"Sigma_2:\", rupstr(sigma2(h0,m)))\n",
    "Tot=sigma0+sigma1+sigma2(h0,m) + cP61+cP62\n",
    "print(\"Total:\",rupstr(Tot))\n",
    "km1p2=floor(100*(find_root((e^x)^(1/42)*x^(-1)-Cmt2/C3,100,700)))/100\n",
    "if Tot<1:\n",
    "    print(\"This constant proves the theorem from e^\"+str(k0)+\" until e^\"+rdnstr(km1p2))\n",
    "else:\n",
    "    print(\"Error!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m: 20.0000000000000\n",
      "Constant in Eq. 5.1: 1.3860\n",
      "Lower bound for b: 0.8706\n",
      "Eq 5.12: 1.2936\n",
      "Eq 5.13: 0.04460\n",
      "Following display: 0.05770\n",
      "Eq 5.14: 2.5872 0.08920 0.1154\n",
      "Lower bound for B: 30.0000 12000\n",
      "Eq 5.20: 0.008800 0.0001000 0.06680\n",
      "Lower Bound for A: 11.8501\n",
      "Eq 5.22: 0.4083 0.2370\n",
      "Eq 5.4: 0.1336 0.6000 0.01760\n",
      "Eq 5.5: 1.0468 1.2264\n",
      "Small H, Small M: 0.07830 0.3374 -60.4992 0.1458\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Small H, Large M: 0.002000 1.0120 19.8754 0.2175\n",
      "Eq 6.3: 0.04010 0.0006000\n",
      "1.10812447687870e38 20.0000000000000\n",
      "Sigma_0: -0.05960\n",
      "Sigma_2 (Small H): 0.2245\n",
      "Sigma_3 (Small H): 0.3639\n",
      "Total (C = 5, k_0 = 400 ): 0.9811\n",
      "\n",
      "\n",
      "\n",
      "m: 42.4264068711929\n",
      "Constant in Eq. 5.1: 1.3860\n",
      "Lower bound for b: 0.8706\n",
      "Eq 5.12: 1.2936\n",
      "Eq 5.13: 0.04210\n",
      "Following display: 0.05450\n",
      "Eq 5.14: 2.5872 0.08420 0.1090\n",
      "Lower bound for B: 50.9116 91641\n",
      "Eq 5.20: 0.02080 0.001400 0.03480\n",
      "Lower Bound for A: 11.8501\n",
      "Eq 5.22: 0.4083 0.2370\n",
      "Eq 5.4: 0.06960 1.2000 0.04160\n",
      "Eq 5.5: 1.0468 1.2264\n",
      "Small H, Small M: 0.1183 0.1758 -36.2266 0.2062\n",
      "Small H, Large M: 0.0001000 0.5273 55.6130 0.2719\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eq 6.3: 0.04010 0.0003000\n",
      "7.98575507113880e159 42.4264068711929\n",
      "Sigma_0: -0.05960\n",
      "Sigma_2 (Small H): 0.1146\n",
      "Sigma_3 (Small H): 0.4784\n",
      "Total (C = 2, k_0 = 1800 ): 0.9857\n",
      "\n",
      "\n",
      "\n",
      "m: 707.106781186548\n",
      "Constant in Eq. 5.1: 1.3860\n",
      "Lower bound for b: 0.8706\n",
      "Eq 5.12: 1.2936\n",
      "Eq 5.13: 0.04020\n",
      "Following display: 0.05210\n",
      "Eq 5.14: 2.5872 0.08040 0.1042\n",
      "Lower bound for B: 615.1828 307591449\n",
      "Eq 5.20: 0.03260 0.005900 0.02890\n",
      "Lower Bound for A: 11.8501\n",
      "Eq 5.22: 0.4083 0.2370\n",
      "Eq 5.4: 0.05780 1.7400 0.06520\n",
      "Eq 5.5: 1.0468 1.2264\n",
      "Small H, Small M: 0.01410 0.1460"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " -52.6969 0.1601\n",
      "Small H, Large M: 0.0001000 0.4379 84.1002 0.4381\n",
      "Eq 6.3: 0.04010 0.0001000\n",
      "1.40333168021306e43435 707.106781186548\n",
      "Sigma_0: -0.05960\n",
      "Sigma_2 (Small H): 0.007100\n",
      "Sigma_3 (Small H): 0.5983\n",
      "Total (C = 1, k_0 = 500000 ): 0.9981\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Theorem 3:\n",
    "# Section 5.1\n",
    "for (C3,k0,d1) in [(5,400,30/100),(2,1800,60/100),(1,500000,87/100)]:\n",
    "\n",
    "    # Proposition 1\n",
    "    l = 102/100         # Value of lambda\n",
    "    m = sqrt(k0)\n",
    "    x0 = exp(k0)      # Lower bound on x\n",
    "    MC1 = 5           # Coefficient of split point in (5.5)\n",
    "    # From Sect. 3\n",
    "    print(\"m:\",n(m))\n",
    "    cL6 = rdn((1.997*4/3)^(1/3))\n",
    "    # From Sect. 4, updated\n",
    "    cL8a = rdn(cL8^(1/3))\n",
    "    cL9ii = rdn(1/2-l/m)\n",
    "    print(\"Constant in Eq. 5.1:\", rdnstr(cL6))\n",
    "    print(\"Lower bound for b:\",rdnstr(cL8a))\n",
    "    #print(rupstr(cL9i), rupstr(cL9ii))\n",
    "\n",
    "    # Computed constants appearing in the proof of Proposition 1:\n",
    "    cP1512 = rup(l^6/cL8a)             # Eq 5.12\n",
    "    cP1513 = rup((l-1)/cL9ii)          # Eq 5.13\n",
    "    cP151x = rup(cP1512*cP1513)        \n",
    "    print(\"Eq 5.12:\", rupstr(cP1512))\n",
    "    print(\"Eq 5.13:\", rupstr(cP1513))\n",
    "    print(\"Following display:\", rupstr(cP151x)) \n",
    "    print(\"Eq 5.14:\", rupstr(2*cP1512), rupstr(2*cP1513), rupstr(2*cP151x))             # Eq 5.14\n",
    "    bB = floor(m*d1*C3*k0)                          # Lower bound on B\n",
    "    print(\"Lower bound for B:\",rdnstr(m*d1*C3), bB)\n",
    "\n",
    "    cP1517 = rup(2*cP151x*(1002/1000*d1)^(4/3)/(8/3))    # Eq 5.17\n",
    "    cP1518 = rup(2*cP1513*(1002/1000*d1)^4/8)                 # Eq 5.18\n",
    "    cP1520 = rup((l-1)/d1+cP1518+1.5*x0^(-1/5))\n",
    "    print(\"Eq 5.20:\", rupstr(cP1517), rupstr(cP1518), rupstr(cP1520))                                 # Eq 5.20 \n",
    "\n",
    "    cP1521 = 2023/10000                                               # Eq 5.21\n",
    "    bA = rdn(cL6*MC1^(4/3))                                       # Lower bound on A\n",
    "    cP1522a = rup(0.6*2^(-8/3)*(1/2-1/bA)^(-5/3))\n",
    "    cP1522b = rup(cP1522a*1.386^(-5/3))\n",
    "    print(\"Lower Bound for A:\",rdnstr(bA))\n",
    "    print(\"Eq 5.22:\",rupstr(cP1522a), rupstr(cP1522b))  \n",
    "                                                                  # Eq 5.22\n",
    "\n",
    "    print(\"Eq 5.4:\",rupstr(2*cP1520), rupstr(2*d1), rupstr(2*cP1517))                  # Eq 5.4\n",
    "    cP155a = rup(4*cP1512*cP1521)\n",
    "    cP155b = rup(4*cP1512*cP1522b)\n",
    "    print(\"Eq 5.5:\", rupstr(cP155a), rupstr(cP155b))              # Eq 5.5\n",
    "\n",
    "\n",
    "\n",
    "    k1=k0 # \n",
    "    x1=exp(k1)\n",
    "\n",
    "    cT161a = cP155a*MC1^(7/3)*x1^(-1/12)/(1-l^(-7/3))\n",
    "    cT161b = 2*cP1517*(m*C3*k1)^(-1/3)/(1-l^(-1/3))\n",
    "    cT161c = rup(2*cP1520/(20*ln(l)))\n",
    "    cT161d = (2*MC1*d1*x1^(-3/20)/(1-l^(-1)))\n",
    "    cT161e = ((2*cP1520)*(1-log(m*C3*k1/MC1)/log(l)))\n",
    "    cT161I = rup(cT161a+cT161b)+rup(cT161c)/C3\n",
    "    print(\"Small H, Small M:\",rupstr(cT161a+cT161b), rupstr(cT161c), rupstr(cT161d+cT161e),rupstr(cT161I))\n",
    "    if cT161d + cT161e >0:\n",
    "        print(\"Error!!\")\n",
    "\n",
    "    cT161f = cP155b*x1^(-1/15)/(1-l^(-1/9))\n",
    "    cT161g = 2*cP1517*MC1^(-1/3)*x1^(-1/60)/(1-l^(-1/3))\n",
    "    cT161h = rup(3*2*cP1520/(20*ln(l)))\n",
    "    cT161j = 2*d1/(1-l^(-1))\n",
    "    cT161k = ((2*cP1520)*(1-log(MC1)/log(l)))\n",
    "    cT161II = rup(rup(cT161f+cT161g)+cT161h/C3+max(0,rup(cT161j+cT161k)/(C2*k0)))\n",
    "    print(\"Small H, Large M:\", rupstr(cT161f+cT161g), rupstr(cT161h), rupstr(cT161j+cT161k),rupstr(cT161II))\n",
    "\n",
    "\n",
    "\n",
    "    # From Sect. 3\n",
    "    cL4 = 4995/10000\n",
    "    cC1 = rup((l-1)/cL4)\n",
    "    cT163a = cC1/(1-l^(-2))\n",
    "    cT163b = x0^(-1/5)/log(l)*(1/10+(1/2*log(2)+log(l))/k0)    \n",
    "    cT163 = rup((cT163a/k0+cT163b)/C3)\n",
    "    print(\"Eq 6.3:\",rupstr(cC1),rupstr(cT163))\n",
    "\n",
    "    #From Sect. 2\n",
    "    bJ = 120\n",
    "    bK = prime_pi(bJ)\n",
    "    prList = prime_range(1,bJ)\n",
    "    h0 = C3*k0*x0^(1/5)\n",
    "    sigma0 = rup(1 - prod([1 - 1/p^2 for p in prList]) - sum([1/p^2 for p in prList]) + 2^bK/h0)\n",
    "    sigma1 = 4523/10000\n",
    "\n",
    "\n",
    "    h1 = C3*k1*x1^(1/5)\n",
    "    print(n(h1),n(m))\n",
    "\n",
    "    print(\"Sigma_0:\",rupstr(sigma0))\n",
    "    print(\"Sigma_2 (Small H):\",rupstr(sigma2(h1,m)))\n",
    "    print(\"Sigma_3 (Small H):\",rupstr(cT161I + cT161II + cT163))\n",
    "    print(\"Total (C = \"+str(C3)+\", k_0 =\",k0,\"):\",rupstr(sigma0+sigma1+rup(sigma2(h1,m)) + rup(cT161I + cT161II + cT163)))\n",
    "    print(\"\\n\\n\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m: 23.4520787991171\n",
      "T_3 \\cap I Bound: 0.06650 5.0262\n",
      "Num Intervals Bound: 0.008100\n",
      "Eq 5.24 0.04080\n",
      "B lower bound: 9802\n",
      "Eq 5.26: 0.001500\n",
      "Eq 5.27: 2.04460687796423e-6\n",
      "Eq 5.29: 0.005300\n",
      "A lower bound: 0.8316 26.6834\n",
      "Eq 5.31: 0.2106 0.3662\n",
      "Eq 5.7: 0.01060 0.3800 0.003000\n",
      "Eq 5.8: 0.6806 3.6812\n",
      "Small H, Small M: 0.1525 0.1264 -38.5464 0.2157\n",
      "Small H, Large M: 0.02190 0.6316 175.0564 0.3377\n",
      "Large H:\t  0.1525 0.7579 136.4993 0.7710\n",
      "Eq 7.3: 0.006100 0.001000\n",
      "Sigma_3 (Large H): 0.7720\n",
      "Sigma_3 (Small H): 0.5544\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2 (Large H): 0.2688\n",
      "Sigma_2 (Small H): 0.2688\n",
      "Total (Large H): 1.2090\n",
      "Total (C = 2, k_0 = 550 ): 0.9914\n",
      "\n",
      "\n",
      "\n",
      "m: 47.9583152331272\n",
      "T_3 \\cap I Bound: 0.06660 5.0236\n",
      "Num Intervals Bound: 0.006900\n",
      "Eq 5.24 0.03470\n",
      "B lower bound: 72800\n",
      "Eq 5.26: 0.003800\n",
      "Eq 5.27: 0.0000478118912283077\n",
      "Eq 5.29: 0.003100\n",
      "A lower bound: 0.8321 26.6994\n",
      "Eq 5.31: 0.2106 0.3656\n",
      "Eq 5.7: 0.006200 0.6600 0.007600\n",
      "Eq 5.8: 0.6802 3.6733\n",
      "Small H, Small M: 0.2381"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 0.07390 -27.0547 0.3120\n",
      "Small H, Large M: 0.0001000 0.3695 321.6971 0.3696\n",
      "Large H:\t  0.2381 0.4433 294.6362 1.7153\n",
      "Eq 7.3: 0.006100 0.0005000\n",
      "Sigma_3 (Large H): 1.7158\n",
      "Sigma_3 (Small H): 0.6821"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2 (Large H): 0.1416\n",
      "Sigma_2 (Small H): 0.1416\n",
      "Total (Large H): 2.0256\n",
      "Total (C = 1, k_0 = 2300 ): 0.9919\n",
      "\n",
      "\n",
      "\n",
      "m: 273.861278752583\n",
      "T_3 \\cap I Bound: 0.06660 5.0236\n",
      "Num Intervals Bound: 0.006200\n",
      "Eq 5.24 0.03120\n",
      "B lower bound: 9242818\n",
      "Eq 5.26: 0.006400\n",
      "Eq 5.27: 0.000276229441263427\n",
      "Eq 5.29: 0.002500\n",
      "A lower bound: 0.8324 26.7090\n",
      "Eq 5.31: 0.2106 0.3652\n",
      "Eq 5.7: 0.005000 0.9000 0.01280\n",
      "Eq 5.8: 0.6802 3.6693\n",
      "Small H, Small M: 0.08850 0.05960 -33.1639 0.2077\n",
      "Small H, Large M: 0.0001000 0.2980 443.6719 0.5961\n",
      "Large H:\t  0.08850 0.3575 410.5029 3.6843\n",
      "Eq 7.3: 0.006100 0.0001000\n",
      "Sigma_3 (Large H): 3.6844\n",
      "Sigma_3 (Small H): 0.8039\n",
      "Sigma_0: -0.006600\n",
      "Sigma_2 (Large H): 0.02560\n",
      "Sigma_2 (Small H): 0.02560\n",
      "Total (Large H): 3.8782\n",
      "Total (C = 1/2, k_0 = 75000 ): 0.9977\n",
      "\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Theorem 4\n",
    "# Section 5.2\n",
    "for (C3,k0,d2) in [(2,550,38/100),(1,2300,66/100),(1/2,75000,90/100)]:\n",
    "    # Proposition 2\n",
    "    l = 1002/1000          # Value of lambda\n",
    "    m = sqrt(k0)            # Value of m\n",
    "    print(\"m:\",n(m))\n",
    "    x0 = exp(k0)      # Enhanced lower bound on x\n",
    "    MC2 = 18          # Coefficient of split point in (5.8)\n",
    "    # From Sect. 3\n",
    "    cL7 = rdn((0.4 - 0.05/m)^(1/5))\n",
    "    # From Sect. 4, updated\n",
    "    # cL11 = 0.0665\n",
    "    cL11a = cL11(l,m)^(1/3)\n",
    "\n",
    "    # Computed constants appearing in the proof of Proposition 2:\n",
    "    cP2524i = rup(2*l^9/cL11a)\n",
    "    cP2524ii = rup((l-1)/cL12ii(l,m))\n",
    "    cP2524x = rup(cP2524i*cP2524ii)\n",
    "    print(\"T_3 \\\\cap I Bound:\",rdnstr(cL11(l,m)),rupstr(cP2524i))\n",
    "    print(\"Num Intervals Bound:\", rupstr(cP2524ii))                             # Eq 5.24\n",
    "    print(\"Eq 5.24\", rupstr(cP2524x))\n",
    "    bB = floor(m*d2*C3*k0)                                        # Lower bound on B\n",
    "    print(\"B lower bound:\",bB)\n",
    "\n",
    "    cP2526 = rup(cP2524x*(1.001*d2)^2/4)             # Eq 5.26\n",
    "    cP2527 = cP2524ii*(1.001*d2)^6/12                # Eq 5.27\n",
    "    cP2529 = rup((l-1)/d2+cP2527+1.5*x0^(-1/7))\n",
    "    print(\"Eq 5.26:\",rupstr(cP2526))\n",
    "    print(\"Eq 5.27:\",n(cP2527))\n",
    "    print(\"Eq 5.29:\",rupstr(cP2529))                                 # Eq 5.29 \n",
    "\n",
    "    cP2530 = 0.0677                                               # Eq 5.30\n",
    "    bA = rdn(cL7*MC2^(6/5))                                       # Lower bound on A\n",
    "    cP2531a = rup((1/48)*(1/2-1/bA)^(-3))\n",
    "    cP2531b = rup(cP2531a*cL7^(-3))\n",
    "    print(\"A lower bound:\", rdnstr(cL7),rdnstr(bA))\n",
    "    print(\"Eq 5.31:\",rupstr(cP2531a), rupstr(cP2531b))  \n",
    "                                                                  # Eq 5.31\n",
    "\n",
    "    print(\"Eq 5.7:\", rupstr(2*cP2529), rupstr(d2), rupstr(2*cP2526))          # Eq 5.7 \n",
    "    cP258a = rup(2*cP2524i*cP2530)\n",
    "    cP258b = rup(2*cP2524i*cP2531b)\n",
    "    print(\"Eq 5.8:\",rupstr(cP258a), rupstr(cP258b))\n",
    "                                                                  # Eq 5.8\n",
    "\n",
    "    # Section 7.1 \n",
    "    k1=k0\n",
    "    x1 = exp(k1)      # Lower bound on x\n",
    "\n",
    "\n",
    "    cT271a = cP258a*MC2^(11/3)*x1^(-1/18)/(1-l^(-11/3))\n",
    "    cT271b = 2*cP2526*(m*C3*k1)^(-1/3)/(1-l^(-1/3))\n",
    "    cT271c = rup(2*cP2529/(42*ln(l)))\n",
    "    cT271d = (2*MC2*d2*x1^(-5/42)/(1-l^(-1)))\n",
    "    cT271e = ((2*cP2529)*(1-log(m*C3*k1/MC2)/log(l)))\n",
    "    cT271I = rup(rup(cT271a+cT271b)+cT271c/C3)\n",
    "    print(\"Small H, Small M:\",rupstr(cT271a+cT271b), rupstr(cT271c), rupstr(cT271d+cT271e), rupstr(cT271I))\n",
    "    if cT271d+cT271e >0:\n",
    "        print(\"Error!\")\n",
    "\n",
    "    cT271f = cP258b*x1^(-1/21)/(1-l^(-1/15))\n",
    "    cT271g = 2*cP2526*MC2^(-1/3)*x1^(-1/126)/(1-l^(-1/3))\n",
    "    cT271h = rup(5*2*cP2529/(42*ln(l)))\n",
    "    cT271j = d2/(1-l^(-1))\n",
    "    cT271k = (2*cP2529)*(1-log(MC2)/log(l))\n",
    "    cT271II = rup(rup(cT271f+cT271g)+cT271h/C3)\n",
    "    cT271 = cT271I + cT271II\n",
    "    print(\"Small H, Large M:\",rupstr(cT271f+cT271g), rupstr(cT271h), rupstr(cT271j+cT271k), rupstr(cT271II))\n",
    "\n",
    "\n",
    "    cT272a = cP258b*x0^(-1/21)/(1-l^(-1/15))\n",
    "    cT272b = 2*cP2526*(m*C3*k0)^(-1/3)/(1-l^(-1/3))\n",
    "    cT272c = rup(2*cP2529/(7*ln(l)))\n",
    "    cT272d = d2/(1-l^(-1))\n",
    "    cT272e = (2*cP2529)*(1-log(m*C3*k0)/log(l))\n",
    "    cT272f = cT272d+cT272e\n",
    "    cT272 = rup(rup(cT272a+cT272b)+cT272c/C3) + rup(cT272f/k1u/C3)\n",
    "    print(\"Large H:\\t \",rupstr(cT272a+cT272b), rupstr(cT272c), rupstr(cT272f), rupstr(cT272))\n",
    "\n",
    "    # From Sect. 3\n",
    "    cL4 = 333/1000\n",
    "    cC1 = rup((l-1)/cL4)\n",
    "    cT273a = cC1/(1-l^(-3))\n",
    "    cT273b = x0^(-1/7)/(3*ln(l))*(1/7+log(2*l^3)/k0)\n",
    "    cT273 = rup((cT273a/k0+cT273b)/C3)\n",
    "    print(\"Eq 7.3:\",rupstr(cC1), rupstr(cT273))\n",
    "\n",
    "    print(\"Sigma_3 (Large H):\",rupstr(cT272 + cT273))\n",
    "    print(\"Sigma_3 (Small H):\",rupstr(cT271 + cT273))\n",
    "\n",
    "    #From Sect. 2\n",
    "    bJ = 20\n",
    "    bK = prime_pi(bJ)\n",
    "    prList = prime_range(1,bJ)\n",
    "    h0 = C3*k0*x0^(1/7)\n",
    "    h1 = C3*k1*x1^(1/7)\n",
    "    sigma0 = n(rup(1 - prod([1 - 1/p^3 for p in prList]) - sum([1/p^3 for p in prList]) + 2^bK/h0))\n",
    "    sigma1 = 1748/10000\n",
    "    print(\"Sigma_0:\",rupstr(sigma0))\n",
    "    print(\"Sigma_2 (Large H):\",rupstr(sigma2(h0,m)))\n",
    "    print(\"Sigma_2 (Small H):\", rupstr(sigma2(h1,m)))\n",
    "    print(\"Total (Large H):\",rupstr(sigma0+sigma1+ rup(sigma2(h0,m)) + rup(cT272 + cT273)))\n",
    "    print(\"Total (C = \"+str(C3)+\", k_0 =\",k0,\"):\",rupstr(sigma0+sigma1+ rup(sigma2(h1,m)) + rup(cT271 + cT273)))\n",
    "    print(\"\\n\\n\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.7",
   "language": "sagemath",
   "metadata": {
    "cocalc": {
     "description": "Open-source mathematical software system",
     "priority": 10,
     "url": "https://www.sagemath.org/"
    }
   },
   "name": "sage-9.7",
   "resource_dir": "/ext/jupyter/kernels/sage-9.7"
  },
  "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}