```python
1+1
```




    2




```python
import re
from functools import partial
```

## Preparing the List of Prime Gaps

The file `raw_prime_gaps_list.txt` contains the list of known lead prime gaps for all gaps up to $9998$. The data is uniformly stored. A typical row looks like the following, and the meaning of each element of the row is described below.

    6* CFC Glaisher 1877   1.91     2  23
    || |   |        |      |        |  |
    || |   |        |      |        |  ---> first prime in the prime gap
    || |   |        |      |        ------> number of digits of the first prime
    || |   |        |      ---------------> "merit", which is bigger if the gap is bigger wrt prime size
    || |   |        ----------------------> year of discovery
    || |   -------------------------------> discoverer
    || -----------------------------------> C if "conventional gap" (all here are)
    ||                                      F if known first occurrence of gap size, N if not, ? if unknown
    ||                                      C if certified primes, P if probabilistic primes
    |-------------------------------------> * if a "maximal gap", larger than all previous gaps
    --------------------------------------> size of gap
    
Sometimes the first prime is stored not as a number, but as an expression in the form $a \cdot b\# / c - d$, where $b\#$ denotes a *primorial*, or as an expression in the form $2^a + b$. We first read in the list. We extract those rows which are maximal (i.e. has * following the gapsize). Then we keep the gap size and the first prime.


```python
with open("raw_prime_gaps_list.txt", "r") as f:
    lines = f.readlines()
print(lines[3])
```

         6* CFC Glaisher 1877   1.91     2  23
    



```python
def is_maximal_gap(splitlist):
    if "*" in splitlist[0]:
        return True
    return False

def extract_gap_and_prime(line):
    splitlist = line.split(None, 6) # Split on whitespace (None), and make at most 6 splits
                                    # Splitting at most 6 times keeps the expression for the first prime together
    gap_str = splitlist[0]
    first_prime_str = splitlist[6]
    first_prime_str = first_prime_str.rstrip()
    
    #assert gap_str[-1] == "*"
    if gap_str[-1] == "*":
        gap = gap_str[:-1]
    else:
        gap = gap_str
    gap = int(gap)
    
    if "*" in first_prime_str:
        first_prime = extract_primorial(first_prime_str)
    elif "^" in first_prime_str:
        first_prime = extract_two_power_prime(first_prime_str)
    else:
        first_prime = int(first_prime_str)
        
    return gap, first_prime


def memoize(f):
    memo = {}
    def helper(x):
        if x not in memo:            
            memo[x] = f(x)
        return memo[x]
    return helper


primes_to_541 = primes_first_n(100)
@memoize
def primorial(n):
    return prod([p for p in primes_to_541 if p<= n])


def extract_primorial(expression):
    pattern = re.compile(r"(\w+)\W*\*\W*(\w+)\W*\#\W*/\W*(\w+)\W*-\W*(\w+)")
    m = pattern.match(expression)
    a,b,c,d = map(int, m.groups())
    return a*primorial(b)/c - d


def extract_two_power_prime(expression):
    pattern = re.compile(r"2\^(\w+)\W*\+\W*(\w+)")
    m = pattern.match(expression)
    exp, addend = map(int, m.groups())
    return 2**exp + addend
```


```python
# First we test the method on just the first 10 prime gaps
gaps = []
gaps_dict = dict()
for line in lines[:10]:
    gap, first_prime = extract_gap_and_prime(line)
    gaps.append(gap)
    gaps_dict[gap] = [first_prime, first_prime + gap]
```


```python
# These are the gaps
gaps
```




    [1, 2, 4, 6, 8, 10, 12, 14, 16, 18]




```python
# Clear prime pair storage
gaps_dict
```




    {1: [2, 3],
     2: [3, 5],
     4: [7, 11],
     6: [23, 29],
     8: [89, 97],
     10: [139, 149],
     12: [199, 211],
     14: [113, 127],
     16: [1831, 1847],
     18: [523, 541]}




```python
# First we test the method on just the first 10 prime gaps
gaps = []
gaps_dict = dict()
for line in lines:
    gap, first_prime = extract_gap_and_prime(line)
    gaps.append(gap)
    gaps_dict[gap] = [first_prime, first_prime + gap]
```


```python
# Choose random representative to check
gaps_dict[116]
```




    [5845193, 5845309]



## Setting up the model functions

We will model the distribution of primes in a few different ways. The nature of the distribution of primes is uncertain. By the prime number theorem, the average prime gap looks like
$$ X, X + \log X.$$
We expect instead that there are many gaps that are much larger. Cramer thought that prime gaps were never larger than
$$ X, X + \log^2 X.$$
It turns out that he wasn't correct, but it may be true that there is a constant $B$ such that gaps are never larger than
$$ X, X + B \log^2 X.$$
It is known that this constant, if it exists, is at least $2e^{-\gamma} \approx 1.1$.
It is also conjectured that no such constant exists, but that for any $A > 2$ there exists a corresponding constant (which I also denote by $B$) such that gaps are never larger than
$$ X, X + B \log^A X.$$

Recall in the background that we conjecture that for each $\alpha$ there exists a constant $C(\alpha)$ such that
$$ p_{n+1}^\alpha - p_n^\alpha \leq C(\alpha).$$

We will model how gaps behave according to these models and compare them to the actual behavior.
To do this, we define the function
$$ f(A, B, \alpha, x) := (X + B \log^A X)^\alpha - X^\alpha.$$
For each fixed $A, B$, and $\alpha$, maximizing $f_{A, B, \alpha}(x)$ gives a largest-possible estimate of the value of $C(\alpha)$. We say largest-possible because this presupposes that there are primes at exactly the maximizing value of $x$, which corresponds to $x$ and $x + B \log^A x$ being primes themselves.

Note that if Cramer's Conjecture were true (for instance), then it may be that the $C(\alpha)$ predictions we get from the Cramer model might be slightly larger than reality, and in general each model errs on the side of overestimating rather than underestimating.

We will set up model functions with values of $(A, B)$ being $(1,1)$ (PNT), $(2,1)$ (Naive Cramer), $(2,2e^{-\gamma})$ (Current Cramer), $(2, 10)$ (Cramer), and $(2.5, 1)$ (Adleman-McCurley).

### Maximizing the model functions

Maximizing $f_{A,B,\alpha}(x)$ is done naively. Compute the derivative and find when it's zero. There is a unique zero greater than $1$, which is our maximum point. A first-order estimation shows that $f_{A,B,\alpha}$ is maximized near $m(A,B,\alpha) = e^{A/(1-\alpha)}$. Note that $B$ does not affect this first-order estimation, so I do not expect the three variants of Cramer to behave significantly differently.

We have that
$$ \begin{align}
f(A, B, \alpha, x) &= (x + B \log^A x)^\alpha - x^\alpha, \\
f'(A, B, \alpha, x) &= \alpha(x + B \log^A x)^{\alpha - 1} \left( 1 + AB \frac{\log^{A-1} x}{x} \right) - \alpha x^{\alpha - 1}.
\end{align}$$

Note that $f' = 0$ (for $x > 1$) is the same as $f'/(\alpha x^{\alpha - 1}) = 0$, which can be rearranged into
$$ g_{A, B, \alpha}(x) := \left( 1 + AB \frac{\log^{A-1} x}{x} \right) - \left( 1 + B \frac{\log^A x}{x} \right)^{1 - \alpha} = 0.$$

This is a pretty well-behaved function, and numerical root solvers should be able to handle it. This is what we implement now.


```python
def F(A, B, alpha, x):
    return (x + B*log(x)**A)**alpha - x**alpha
```


```python
def g(A, B, alpha, x):
    return (1 + A * B * log(x)**(A-1)/x) - (1 + B*log(x)**A / x)**(1-alpha)
```


```python
find_root(partial(g,2,1,.931), 2, 10**17)
```




    3875382310049.5977



Actually, this caps out at 17 digit precision and won't suffice for $\alpha > 0.93$. But for really large $x$, the binomial expansion should become more and more accurate, leading to $e^{A/(1-\alpha)}$ being a better and better estimate of the maximizing $x$. For $\alpha > 0.93$, I need to use a different technique. This runs into 64-bit precision arithmetic problems. To solve this, for $\alpha > 0.9$, I will use a homegrown bisection rootfinder in a high precision arithmetic sage environment.


```python
RF = RealField(1000)
def RF_F(A, B, alpha, x):
    return (RF(x) + B*RF(log(x))**A)**alpha - RF(x)**alpha
```


```python
start = exp(2/.005)
```


```python
RF_F(2, 1, .995, start)
```




    21545.3770912687021748108986403427084075846916931792212525384307640395688240237337048227691696804791326950035399102608956974569724349727991073964722581240384122762528927904117476919639364452675336235413841573393011675588938510358312010132614449997916430297316353906271329030935986082731144694523430359




```python
def RF_g(A, B, alpha, x):
    return (RF(1) + A * B * RF(log(x))**(A-1)/RF(x)) - (RF(1) + B*RF(log(x))**A / RF(x))**(1-alpha)
```


```python
def RF_gsim(x):
    return RF_g(2,1,0.8,x)
```


```python
start = exp(2/.2)
```


```python
RF_gsim(start)
```




    1.64444546432196045785420540155893726063200274572267336535350345517820986134145980935118162918424805451405780491393962735635632068351784816308968439342270391699018760425327496611550914706761553155070977740717893810316453496034348416368377098158251454957552350910050505616367388235938291760692342296640e-6




```python
find_root(RF_gsim, 2, 10**17)
```




    22423.973931885273




```python
# Note that this requires that func(2*start) < 0 to work properly
# In our case, this means that the `start` guess can't be too bad
def bisection_rootfinding(func, start, iterations=1, goal=10**(-7)):
    if start*2**(-iterations) < goal:
        return start
    if func(start) > 0:
        return bisection_rootfinding(func, start + start * 2**(-iterations), iterations=iterations+1, goal=goal)
    return bisection_rootfinding(func, start - start * 2**(-iterations), iterations=iterations+1, goal=goal)
```


```python
bisection_rootfinding(RF_gsim, start)
```




    22423.9739319633



Notice in particular that this is essentially the same as scipy's find_root result, so my method seems to work pretty well.

Now let's try it on something larger.


```python
def RF_gsim2(x):
    return RF_g(2,1,0.95,x)
```


```python
start = RF(exp(2/.05))
```


```python
find_root(RF_gsim2, 2, 10**16)
```




    9999999710581220.0




```python
bisection_rootfinding(RF_gsim2, start)
```




    2.35385266837042022843327073581445209550189717806195566505913760542398087081961498860906821406821821384500066808118082097402013219336066587283757576744855862678903728952660301215698355426390611932122470916392998635276106218830903304888702478640727370013260828376486514621410432957225533987162136976562e17




```python
fr = find_root(RF_gsim2, 2, 10**16)
br = bisection_rootfinding(RF_gsim2, start)
```


```python
br/fr
```




    23.538527364953389




```python
RF_F(2,1,0.95,br)
```




    205.709630519650881691924521792765193253506704662129143200009919811578164796694847391279127980866813238127336707532494092877547455355110417039290138055484169795377093169784536561532277570158710114624712809121313730522096004967342309066003863302396648990201591070015928728914879882408923972678419871216




```python
RF_F(2,1,0.95,fr)
```




    204.359554823890368403363980346311818555686040610577329335810851723660578683578255065736205720488148424074134614592141667900060298707050938806333994911863866233402912888249520944989763989818618341662786563588774804722056915716107185579765431262102230328894524070859044620196971152617523013227880281789



This is interesting --- there is a region in which find_root fails due to precision reasons, but not in an obviously detectable way. But as this clearly shows, my bisection method suffices.

## Data Creation

The plan is to loop through $\alpha$ from $0.005$ to $0.995$, compute everything, and compare everything. But for many initial $\alpha$, everything is boring. Let's first find exactly when $\alpha$ is boring.


```python
def estimate_c_alpha(alpha):
    MAX = 0
    MAXIMIZING_PRIME = 0
    MAXIMIZING_GAP = 0
    for gap in gaps:  # Note: It is possible to optimize this selection based on vals for other alpha
        a,b = gaps_dict[gap]
        diff = RF(b)**alpha - RF(a)**alpha
        if diff > MAX:
            MAX = diff
            MAXIMIZING_PRIME = a
            MAXIMIZING_GAP = gap
    return MAX, MAXIMIZING_GAP, MAXIMIZING_PRIME
```


```python
for a in range(20):
    alpha = 0.05 + 0.04 * a
    print alpha, estimate_c_alpha(alpha)
```

    0.0500000000000000 (0.0273310781848302800658590880398522815347590177971289917771832998742533538405129347101901809408899390154774514006487716615848197146224848388750750143388527967917247425289525849885984701315952214994461518869467334113477417254612430110708981668255496105519249618121865605658797754277991125986546424926508, 2, 3)
    0.0900000000000000 (0.0519370809347252520685872649136154113889575896783660909311486774384967763294410100201853785718930335741505701554636256358169356949220163741754562328432777159041140235202911046252875832381116056987398716939257761619006205913299310880723598291675100978615964192750372577866405938486631637439508364320858, 2, 3)
    0.130000000000000 (0.0792030225749460145408332846210646211349709553727122695589775991379642815865105358693876664136993146934495792581811221651860043234387972830839584009007088431244195480057541946715422880580136690530288636745487275691029720040890241179152279971140525502557180885918952288021915731111539843731655078009650, 2, 3)
    0.170000000000000 (0.111181263499968631320136517388785040808714392234047453641809201804708639290279254938652323470002276346079710485543325649046090275258826662034012458243570622205524629138788775066951458339047395559511537749810240846304395803420508311453056782745024686529577255723133918595809482942348237543967846160003, 4, 7)
    0.210000000000000 (0.149826273160475994962243083701533895900842934406639563544922783813024602459916508959343188251692967568113386862740598596871970918514700800403775752958239520903512136772589541047142530312391410824945291551508017637758651939264125623653517866914330922645399426878943623423072451500026730493789321506636, 4, 7)
    0.250000000000000 (0.194583725140086155235410867344262555198847337413419396460571687675616808191366806740479594064082037831924757984080937834304785500452520906853755534860850295188536908645060113128466605462778750838482301927791913256899345519152556683960954609815854565747156297795965378511888898199097339590595739963184, 4, 7)
    0.290000000000000 (0.246248813819103170253557250449355251524182358470425963859211259303798615289095890248478723442933301746121823090879497185382980418432530720592316569672647132718753876297567208875018396019594127419469968573820605380516809172280817095923124075568117875266835223189694717618200541375491919546679627098467, 4, 7)
    0.330000000000000 (0.305711343206456535065347320087438603464924475308724063528038459208099421322271293811067461182623734565140057698384504321768168205679408867081213453787162538160335263784242371458734357856213201079531522786702433012164371829170292653473361351573652392532724606923212418835000408796683348784457264111052, 4, 7)
    0.370000000000000 (0.373966429543365074817695190057943133748849296316293016605715723578219057277796879857315638543939944511494607661966046966136173460500005119867241883996572229561550012814088131390300062064484716155148271213476199228912371725988542658561793549101843469396474674653770477560665905219304663796828624873265, 4, 7)
    0.410000000000000 (0.452126376701729070189230306846564166572718811750208095294640063601049113790103772310658620779822483542466716249347768517447112764806799712472932886358915375650933223964845094374376513000612507694391310907365957562029322741530071576739195666770017192298200796425941866160240244584347982547800294338410, 4, 7)
    0.450000000000000 (0.541433850160674856517448019977648670845387468609120824704039687062082365058862967523552961325586337836348173605227514038998558396631077325859615480296792512772133696891677065764305911703916867104804927875994085360984285663433909652187556181498138595496488010425855281207978804879931061852800799756917, 4, 7)
    0.490000000000000 (0.643276488631938665905352690029550515258841518006625633245321825184230438817001007398166180674071026761014493558595088627772662382841185013505094743174300363527368853334206213217354029111657745924317354222327336246191131256221501182273919522940639882386609873185598418596322399256320888450121570394212, 4, 7)
    0.530000000000000 (0.782274375222570659961116752409324914028545518530568244320104431746104885735461769616383026384078842399462241244581615655826225425861047007355692369799948454253659291398201594313709988207531393383294538319815566762354200807118161016244669292398154414432780957805186418730929458670499355583639366807496, 14, 113)
    0.570000000000000 (1.01884142574031551509021338668081500059651084095766586110657227080467560613268941623508491511136916158752831640137800592349753586674611323376787825530400982051690086491024022896463721852472889392683539307104501345623645811741468826667843990073735076008410738714094795624533929203128446699606735615606, 14, 113)
    0.610000000000000 (1.32041627446091249684010282433791853388208830739656554504138875808096002722343829956112453368302449982309552879586646808802152034816717620334622823023810628771492706007214793761849531553602264771273615622992631932925771633549693623574091063218430747286446606981048690064876572566265298842449766093359, 14, 113)
    0.650000000000000 (1.77606150694455767148184879493035887223003417395280279720283805453747271273934846538180951981092679591381731044243677198786743664042472880933235210451979387660865985682252547374928040543948333428791944282726895450499622924451542431454744301895332923950015913854758146588145372911146101685185764793123, 34, 1327)
    0.690000000000000 (2.51494758180812497524004824914538864215457135214415764156856021955949292849607775801439209422484553476893581354288954398402859547490062743970891503842286272825447055396110162612007004241446273913524824849127987779554312507146909153549001997795834074088281582457832963082371713882142574398014399748687, 34, 1327)
    0.730000000000000 (3.54926107221611583152257524364752677913475564845984515579046097207485619605424351240800845901580940996317144955978790853632298886397526616816851581340281228865107002973080947704010636374478688987465249159849949548386146640449395044703904284644139242341066041035459498716689367467428185462776529563374, 34, 1327)
    0.770000000000000 (5.12182230348182730576342401351876551208525168860138306791801121718534370070661507654496625833636423338279349496907637710728402502004783185640546392968476514453708497086637272112232436492441948446521882798746778382139248643241899662976156695772601498121414397776266212838948354134427604190418086056576, 72, 31397)
    0.810000000000000 (8.15293966895682191988905194405319109615825254285457443436102887983737009629372030499598020896742552803941325831592633996341564615822777249900088377876688336896985842163801399806797855932027735326506916167614518727288169134555219040681812889496872163795527757013603506051256509476820385589025707472966, 72, 31397)


It appears that starting at $\alpha = 0.005$ will be fine, though not a lot happens until $\alpha \approx 0.5$. The semiimportant data there is really just the value of $C(\alpha)$.

Let us now set up the data structures to store the data. These will be a set of dictionaries, contained in a "model" dictionary.


```python
experimental_model = dict()
PNT_model = dict()
naive_cramer_model = dict()
current_cramer_model = dict()
weak_cramer_model = dict()
AM_model = dict()

models = dict()
models["experimental"] = experimental_model
models["PNT"] = PNT_model
models["naive"] = naive_cramer_model
models["current"] = current_cramer_model
models["weak"] = weak_cramer_model
models["AM"] = AM_model
```

In each model dictionary, the keys will be values of $\alpha$ from $0.005$ to $0.995$, and the values will be `(C(alpha), gap, lower_prime)` tuples. That's that.

I'll store the data as pickled dictionaries, and as json dumps (for external portability).


```python
RF = RealField(1000)
def RF_F(A, B, alpha, x):
    return (RF(x) + B*RF(log(x))**A)**alpha - RF(x)**alpha

def RF_g(A, B, alpha, x):
    return (RF(1) + A * B * RF(log(x))**(A-1)/RF(x)) - (RF(1) + B*RF(log(x))**A / RF(x))**(1-alpha)


def maximize_model(A, B, alpha):
    fun = lambda x: RF_F(A, B, alpha, x)
    deriv = lambda x: RF_g(A, B, alpha, x)
    estimate_lower_prime = 2*RF(exp(A/(1 - alpha)))
    lower_prime = bisection_rootfinding(deriv, start)
    Calpha = fun(lower_prime)
    gap = B*log(lower_prime)**A
    return (Calpha, gap, lower_prime)
```


```python
def make_experimental_data(alpha):
    Calpha, gap, lower_prime = estimate_c_alpha(alpha)
    return (Calpha, gap, lower_prime)

def make_model_data(A, B, alpha):
    return maximize_model(A,B,alpha)

def make_PNT_model_data(alpha):
    return make_model_data(1,1,alpha)

def make_naive_model_data(alpha):
    return make_model_data(2,1,alpha)

def make_current_model_data(alpha):
    return make_model_data(2,1.1229189671337703396, alpha)

def make_weak_model_data(alpha):
    return make_model_data(2, 10, alpha)

def make_AM_model_data(alpha):
    return make_model_data(2.1, 1, alpha)
```


```python
for a in range(1000):
    alpha = 0.000 + a/1000
    models["experimental"][alpha] = make_experimental_data(alpha)
    models["PNT"][alpha] = make_PNT_model_data(alpha)
    models["naive"][alpha] = make_naive_model_data(alpha)
    models["current"][alpha] = make_current_model_data(alpha)
    models["weak"][alpha] = make_weak_model_data(alpha)
    models["AM"][alpha] = make_AM_model_data(alpha)
    if a%50 == 0:
        print a
```

    0
    50
    100
    150
    200
    250
    300
    350
    400
    450
    500
    550
    600
    650
    700
    750
    800
    850
    900
    950



```python
#t = models.keys()[0]
#for k in models[t].keys():
#    print k
#    for key in models.keys():
#        print key, "{}".format(float(models[key][k][0]))
#    print
```


```python
alist = models['experimental'].keys()
```


```python
alist = sorted(alist)
alist[:10]
```




    [0.000000000000000,
     0.00100000000000000,
     0.00200000000000000,
     0.00300000000000000,
     0.00400000000000000,
     0.00500000000000000,
     0.00600000000000000,
     0.00700000000000000,
     0.00800000000000000,
     0.00900000000000000]




```python
models['experimental'][alist[-50]]
```




    (186.266539364655933626490915855618094964571041849497438999997826389079814853291919149871340027107619054174485385512758967203039410422848728780559014112645466156493093420011269068962647210601068265065442552646499480648773915312419366944230198770879264447877675770996717230251871662312084943372625852009,
     1132,
     1693182318746371)




```python
exp_plotys = [models['experimental'][a][0] for a in alist]
cur_plotys = [models['current'][a][0] for a in alist]
weak_plotys = [models['weak'][a][0] for a in alist]
naive_plotys = [models['naive'][a][0] for a in alist]
pnt_plotys = [models['PNT'][a][0] for a in alist]
am_plotys = [models['AM'][a][0] for a in alist]
```


```python
exp_plotys[:10]
```




    [0,
     0.000511517768525288490431299564019530566521268528293360224090648225671358406177645208427960618154537414232201941218446425306818783533501147565303119802041928147556730502807612639983069318331337966366989636147947059894660635945989823381565534871988347392491415512712789530278446363938930989574701960971573,
     0.00102442172449282884332508157534822491563469058472756719379918434497150349164072269796738895556658217479601021664620351202402928247792955884906721383928526283258088661094103784976857592343906755172954042356731975418539998880830592392230883140406784021001485475079790186985430084990582641020616554607448,
     0.00153871471873605261609594567254478753444623540399582289121122633476589317565823494066369723602547767528506869164326737509390842483499917786372877323718544377009408577779985644905224201012604546334079646010266208402550801858445699818597717680152286230410860527428679105037060870377855630342605362083545,
     0.00205439960735970041704645998947656177812210389930506653845904193767730286099795243878248319312669702013552227025385042847607239098634986006808941199232325363397106263902487262040597462454010379778333332955140235088061999175180377117872404786751318950384094489121442812986482906623907565586558415282809,
     0.00257147925174905936109663014134097336834442393609311620079119771972767194500766471612512274531782092337622866174440543746300567206334808347961436348689967356555975464533599100596026029910740394817019762901927392208024221858943093784678802392647510480682890612867455763120157088180379124648103287551767,
     0.00308995651857921622546796559200596882746424951676754671589696002353232888475187180207023465138402128824739465016627718117073708463068539941891599637570574693047123871730515453760687505233732008693440140525073857761664831298381288803400485841363604104199122966657955365153196391365677979590315607030346,
     0.00360983427982432633151290676273828177065418823632371050144958102026278975990287331278230039244296151938098161366750880745556329436894367278846820410099015999510215259570924904252430315259755884537232301410584943365318499431178635391607582265002963417075775796112480314734163131768771943838626823279178,
     0.00413111541276689817892355442974284642043215842675225341315076577505076312445243624385504765217391148512516335857815162425101960944675790086942557310550369600595084936465873338520976161756438902253630154479581694290320640417906763915305024844757572628498939719652771807118350659671983418068743491981217,
     0.00465380280000709340462715100328122372167500192805071772269677766578089916306974130382072879733511866037963297784029049269126452696333254096221820907436597556468574696079353308055141116896900016154860397150953509265025577577698312210423534852671473702537036156003845806675868441170706672772677146830196]




```python
P = scatter_plot(zip(alist, exp_plotys), xmin=0, xmax=0.995, ymin=0, ymax=1000, markersize=2, facecolor='black', edgecolor='')
P += scatter_plot(zip(alist, cur_plotys), markersize=2, facecolor='red', edgecolor='')
P += scatter_plot(zip(alist, weak_plotys), markersize=2, facecolor='green', edgecolor='')
P += scatter_plot(zip(alist, pnt_plotys), markersize=2, facecolor='blue', edgecolor='')
P += scatter_plot(zip(alist, am_plotys), markersize=2, facecolor='orange', edgecolor='')
P += scatter_plot(zip(alist, naive_plotys), markersize=2, facecolor='purple', edgecolor='')
P
```


```python
alpha_in_8_to_10 = []
for a in range(200):
    alpha = 0.8 + 0.001 * a
    print "\r" + str(alpha),
    alpha_in_8_to_10.append((alpha, estimate_c_alpha(alpha)[0]))
```

    0.999000000000000    



```python
alpha_in_8_to_10[:10]
```




    [(0.800000000000000,
      7.26014575499710583007752218476427663828926179611157953282726566695049107265377270133560692997036040530500477768867195065477012290917719308708782018245669078582044688607873551418902354119005826714341041787179613542691696812731368111577302778025140783938715803646904732042037141424376103751080011998710),
     (0.801000000000000,
      7.34488929791419664639472074261502192281957152607915337406435419916315200585249998380204802464158172422854707558003676849509272095842610708260861187394075385302864085424482425636215460071693157771352179388596036196527196233733066745050675507537225353771965938807629237803255305635045814031962805504380),
     (0.802000000000000,
      7.43061042252953118360499889585345215018025636058791613024779180241959349552833645769717749958099242955401006554977231966370929595157074587453403662539834509389164431402822851597631998341664334362076793184715469597355362017466344732729542069893132901288343947966821105638284381796933618371826846838619),
     (0.803000000000000,
      7.51732029848634200607433174510508482594122453398141176966268473529751351184646813784800405777531754550192381100544366450785598014002895673006395936770365346627588379366606543940945237679157583436394422090354325265492733676752398517946425200692808124042667052044810697692986390060537707486364599143525),
     (0.804000000000000,
      7.60503022204000734801611994175910058752179636206651996147193043335584596511723929814175288234159427739125200776567201269508854261086012056659928873030403103232887807707163451358405790873269189649627354113518436152672779483523449395753734423912537899862066995252477954706024219583312729407926179911563),
     (0.805000000000000,
      7.69375161748367094826638210122414507148859404247720099871848487689468214713202082463994969124267991718825786287112188912678100536191574433336629908447710876306600467603662495511710583014759510914202253333621545659830430946106714917813648256271009955928561611769164808920368003266509206523186762233538),
     (0.806000000000000,
      7.78349603858982244721779300896340370719651892262359702987062149173543437496559797298960185871544210722387598361933245542587426523378352635215594880061012627819890365697867503048798428438187426312900797894815467141758355394778799680454552418211658510893648886609371746014810249794042746393852202824007),
     (0.807000000000000,
      7.87427517006801615180146452499849061927596528455360826721739657562644794302766771087638979639602449116574992077357141394603976062924856101812872569383559596854294631178413183904768796450844864789738475061388541959034575255714643195898000598548584414383123038816412410548231650309595727849497961986431),
     (0.808000000000000,
      7.96610082903890794669462114539672129478305173791675241547842421561455113734040619186986942287577193152223869834403630008542757211046233366303597315991416245196747104402192182372699339898663234679669488422469870090098190800507409854987061382181736808634586083222536568831129544352980099058547400185061),
     (0.809000000000000,
      8.05898496652479212401601208637250723827102374629869062269780036497368626661197110233759851953267391453351914055670274814265889392222687741155299567683025316275204849480434304438911646122416249013640276960552932240889258471885514419724654940837698763355241134565033813625655292442352013103404866524715)]




```python
alphaprime = []
yprime = []
for tup in alpha_in_8_to_10:
    alpha, y = tup
    alphaprime.append(log(1-alpha))
    yprime.append(log(y))
```


```python
alphaprime[:5]
```




    [-1.60943791243410,
     -1.61445045425764,
     -1.61948824828760,
     -1.62455155024415,
     -1.62964061975162]




```python
yprime[:5]
```




    [1.98239990508105792432914206584034425361291820164801420571371447279489000805005017210551940329748997949158276953758402139559706347895794635989738106647433102557324635006924478567083315201989122813678874733027136858749688442475437260290841954181034549727611229670749280412218485884946091336845146816553,
     1.99400473775108917242607592672379876282962602758486599176148576772489434639025636624056896545690094280224195809733448464124202508093422016580455567072652267572463293430462118712221270594276058468607769040072042105653811588064721872113903111728451820492918543610398905495968245410433690439953389342567,
     2.00560801181928100231139813531008896336733947147889255572229770276025059061492979832016543643104245292920992912477217573183288125389017314482489589303875686728543697919670636242183802728761179540256596711814791281490102327649928484622757922109228990643321882068823404558914950808859897896899546791065,
     2.01720973117000196837249208594601178602418005180730883439142747157620936053011217274044740680743596917763962963449407257327502987553886810995997876898715267399577354310049293809033854648554959973418779247307890041118469773630291760102404488143229295135303350509732825944937643716389433950777297807287,
     2.02880989967311768574223509850438224081778861713460381994857172328879726399734062796210317666818748089330102932762383113852809191013114520866613652369138698179850038119705032754770196695094060071695316927452740911833479460919282968072815444250088745055725905599645333532839537328901190948844738968695]




```python
zip(alphaprime[:5], yprime[:5])
```




    [(-1.60943791243410,
      1.98239990508105792432914206584034425361291820164801420571371447279489000805005017210551940329748997949158276953758402139559706347895794635989738106647433102557324635006924478567083315201989122813678874733027136858749688442475437260290841954181034549727611229670749280412218485884946091336845146816553),
     (-1.61445045425764,
      1.99400473775108917242607592672379876282962602758486599176148576772489434639025636624056896545690094280224195809733448464124202508093422016580455567072652267572463293430462118712221270594276058468607769040072042105653811588064721872113903111728451820492918543610398905495968245410433690439953389342567),
     (-1.61948824828760,
      2.00560801181928100231139813531008896336733947147889255572229770276025059061492979832016543643104245292920992912477217573183288125389017314482489589303875686728543697919670636242183802728761179540256596711814791281490102327649928484622757922109228990643321882068823404558914950808859897896899546791065),
     (-1.62455155024415,
      2.01720973117000196837249208594601178602418005180730883439142747157620936053011217274044740680743596917763962963449407257327502987553886810995997876898715267399577354310049293809033854648554959973418779247307890041118469773630291760102404488143229295135303350509732825944937643716389433950777297807287),
     (-1.62964061975162,
      2.02880989967311768574223509850438224081778861713460381994857172328879726399734062796210317666818748089330102932762383113852809191013114520866613652369138698179850038119705032754770196695094060071695316927452740911833479460919282968072815444250088745055725905599645333532839537328901190948844738968695)]




```python
var('A')
var('B')
model(x) = B - A*x
```


```python
find_fit(zip(alphaprime, yprime), model)
```




    [A == 1.6670411794480962, B == -0.28990653254233467]



This indicates that the exponent from the linear fit is 1.66, which is not so good. I suspect what's happening is the lack of sufficiently large prime gaps means that the experimental data for alpha very near $1$ is simply incorrect. So let's cut out the highest portion of the range.


```python
# From 0.8 to 0.95
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.38952780953096, B == -1.8741114946012887]




```python
# From 0.85 to 0.95
LOWER = 50
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.3838016290728183, B == -1.8595471688135894]



Wow, look at the stability of that result.


```python
# From 0.8 to 0.9
LOWER = 0
UPPER = 100
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.381841660510989, B == -1.8594294921692536]




```python
# From 0.8 to 0.85
LOWER = 0
UPPER = 50
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.3446674868637274, B == -1.7982001515169852]



The sheer stability of this in these ranges indicates to me that maybe there is something here. What if we eek just a little bit more out of the upper range?


```python
# From 0.8 to 0.96
LOWER = 0
UPPER = 160
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.3517182353898307, B == -1.7986473277851818]




```python
# From 0.8 to 0.97
LOWER = 0
UPPER = 170
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
find_fit(zip(alphaprime_9, yprime_9), model)
```




    [A == 2.2790797480846114, B == -1.6504264401594178]



The guess really plummets when we increase from 0.96 to 0.97, so I would guess that this is around where the prime gaps peter off.

Let's use the model data on the range $\alpha \leq 0.8 < 0.95$. On this range, we have

    A' = 2.3895
    B' = -1.874
    
I use $A'$ and $B'$ to denote that this is what the model estimates. The estimate for $A$ is exactly $A'$, so we guess that $A \approx 2.3895$. The estimate for $B$ is a bit more complicated, and work in my notebook RN5.28 shows that $B \approx \exp(B' - A' \log(A'/e))$.


```python
# From 0.8 to 0.95
var('A')
var('B')
model(x) = B - A*x
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
res = find_fit(zip(alphaprime_9, yprime_9), model, solution_dict=True)

```


```python
Aprime = res[A]
Bprime = res[B]
print Aprime, Bprime
```

    2.38952780953 -1.8741114946



```python
exper_A = Aprime
exper_B = exp(Bprime - Aprime * log(Aprime / e))
exper_B = exper_B.n()
print exper_A, exper_B
```

    2.38952780953 0.208859966269026


Hmm, interesting. The growth estimate appears to be overcompensating and tossing more into the exponent and taking away from the coefficient.

As an additional experiment, shat if we force the leading coefficient to be $-1$? I'm not quite sure what this would mean, but sometimes one must play a trombone to petunias.


```python
# From 0.8 to 0.95
var('A')
var('B')
model(x) = -1 - A*x
LOWER = 0
UPPER = 150
alphaprime_9 = alphaprime[LOWER:UPPER]
yprime_9 = yprime[LOWER:UPPER]
res = find_fit(zip(alphaprime_9, yprime_9), model, solution_dict=True)
Aprime = res[A]
print Aprime
```

    1.99402560961



```python
exper_A = Aprime
exper_B = exp(-1 - Aprime * log(Aprime / e))
exper_B = exper_B.n()
print exper_A, exper_B
```

    1.99402560961 0.682384387756181


It's amusingly close to causing the exponent to be $2$, but otherwise it is not particularly meaningful.


```python

```
