text:
# Category Archives: Python

## Extra comparisons in Python’s timsort

Reading through the comp.lang.python mailing list, I saw an interesting question concerning the behavior of the default sorting algorithm in python. This led to this post.

Python uses timsort, a clever hybrid sorting algorithm with ideas borrowing from merge sort and (binary) insertion sort. A major idea in timsort is to use the structure of naturally occuring runs (consecutive elements in the list that are either monotone increasing or monotone decreasing) when sorting.

Let’s look at the following simple list.

10, 20, 5

A simple sorting algorithm is insertion sort, which just advances through the list and inserts each number into the correct spot. More explicitly, insertion sort would

1. Start with the first element, 10. As a list with one element, it is correctly sorted tautologically.
2. Now consider the second element, 20. We insert this into the correct position in the already-sorted list of previous elements. Here, that just means that we verify that 20 > 10, and now we have the sorted sublist consisting of 10, 20.
3. Now consider the third element, 5. We want to insert this into the correct position in the already-sorted list of previous elements. A naively easy way to do this is to either scan the list from the right or from the left, and insert into the correct place. For example, scanning from the right would mean that we compare 5 to the last element of the sublist, 20. As 5 < 20, we shift left and compare 5 to 10. As 5 < 10, we shift left again. As there is nothing left to compare against, we insert 5 at the beginning, yielding the sorted list 5, 10, 20.

How many comparisons did this take? This took 20 > 10, 5 < 20, and 5 < 10. This is three comparisons in total.

We can see this programmatically as well. Here is one implementation of insertion_sort, as described above.

def insertion_sort(lst):
'''
Sorts lst in-place. Note that this changes lst.
'''
for index in range(1, len(lst)):
current_value = lst[index]
position = index
while position > 0 and lst[position - 1] > current_value:
lst[position] = lst[position - 1]
position = position - 1
lst[position] = current_value

Let’s also create a simple Number class, which is just like a regular number, except that anytime a comparison is done it prints out the comparison. This will count the number of comparisons made for us.

class Number:
def __init__(self, value):
self.value = value

def __str__(self):
return str(self.value)

def __repr__(self):
return self.__str__()

def __lt__(self, other):
if self.value < other.value:
print("{} < {}".format(self, other))
return True
print("{} >= {}".format(self, other))
return False

def __eq__(self, other):
if self.value == other.value:
print("{} = {}".format(self, other))
return True
return False

def __gt__(self, other):
return not ((self == other) or (self < other))

def __le__(self, other):
return (self < other) or (self == other)

def __ge__(self, other):
return not (self < other)

def __nq__(self, other):
return not (self == other)

With this class and function, we can run

lst = [Number(10), Number(20), Number(5)]
insertion_sort(lst)
print(lst)

which will print

10 < 20
20 >= 5
10 >= 5
[5, 10, 20]

These are the three comparisons we were expecting to see.

Returning to python’s timsort — what happens if we call python’s default sorting method on this list? The code

lst = [Number(10), Number(20), Number(5)]
lst.sort()
print(lst)

prints

20 >= 10
5 < 20
5 < 20
5 < 10
[5, 10, 20]

There are four comparisons! And weirdly, the method checks that 5 < 20 twice in a row. What’s going on there?1

At its heart, this was at the core of the thread on comp.lang.python. Why are there extra comparisons in cases like this?

Poking around the implementation of timsort taught me a little bit more about timsort.2

Timsort approaches this sorting task in the following way.

1. First, timsort tries to identify how large the first run within the sequence is. So it keeps comparing terms until it finds one that is out of order. In this case, it compares 20 to 10 (finding that 20 > 10, and thus the run is increasing), and then compares 5 to 20 (finding that 5 < 20, and thus that 5 is not part of the same run as 10, 20). Now the run is identified, and there is one element left to incorporate.
2. Next timsort tries to insert 5 into the already-sorted run. It is more correct to say that timsort attempts to do a binary insertion, since one knows already that the run is sorted.3 In this binary insertion, timsort will compare 5 with the middle of the already-sorted run 10, 20. But this is a list of length 2, so what is its middle element? It turns out that timsort takes the latter element, 20, in this case. As 5 < 20, timsort concludes that 5 should be inserted somewhere in the first half of the run 10, 20, and not in the second half.
3. Of course, the first half consists entirely of 10. Thus the remaining comparison is to check that 5 < 10, and now the list is sorted.

We count4 all four of the comparisons. The doubled comparison is due to the two tasks of checking whether 5 is in the same run as 10, 20, and then of deciding through binary insertion where to place 5 in the smaller sublist of 10, 20.

Now that we’ve identified a doubled comparison, we might ask Why is it this way? Is this something that should change?

The short answer is it doesn’t really matter. A longer answer is that to apply this in general would cause additional comparisons to be made, since this applies in the case when the last element of the run agrees in value with the central value of the run (which may occur for longer lists if there are repeated values). Checking that this happens would probably either involve comparing the last element of the run with the central value (one extra comparison, so nothing is really saved anyway), or perhaps adding another data structure like a skip list (which seems sufficiently more complicated to not be worth the effort). Or it would only apply when sorting really short lists, in which case there isn’t much to worry about.

Learning a bit more about timsort made me realize that I could probably learn a lot by really understanding an implementation of timsort, or even a slightly simplified implementation. It’s a nice reminder that one can choose to optimize for certain situations or behaviors, and this might not cover all cases perfectly — and that’s ok.

## Cracking Codes with Python: A Book Review

How do you begin to learn a technical subject?

My first experience in “programming” was following a semi-tutorial on how to patch the Starcraft exe in order to make it understand replays from previous versions. I was about 10, and I cobbled together my understanding from internet mailing lists and chatrooms. The documentation was awful and the original description was flawed, and to make it worse, I didn’t know anything about any sort of programming yet. But I trawled these lists and chatroom logs and made it work, and learned a few things. Each time Starcraft was updated, the old replay system broke completely and it was necessary to make some changes, and I got pretty good at figuring out what changes were necessary and how to perform these patches.

On the other hand, my first formal experience in programming was taking a course at Georgia Tech many years later, in which a typical activity would revolve around an exciting topic like concatenating two strings or understanding object polymorphism. These were dry topics presented to us dryly, but I knew that I wanted to understand what was going on and so I suffered the straight-faced-ness of the class and used the course as an opportunity to build some technical depth.

Now I recognize that these two approaches cover most first experiences learning a technical subject: a motivated survey versus monographic study. At the heart of the distinction is a decision to view and alight on many topics (but not delving deeply in most) or to spend as much time as is necessary to understand completely each topic (and hence to not touch too many different topics). Each has their place, but each draws a very different crowd.

The book Cracking Codes with Python: An Introduction to Building and Breaking Ciphers by Al Sweigart1 is very much a motivated flight through various topics in programming and cryptography, and not at all a deep technical study of any individual topic. A more accurate (though admittedly less beckoning) title might be An Introduction to Programming Concepts Through Building and Breaking Ciphers in Python. The main goal is to promote programmatical thinking by exploring basic ciphers, and the medium happens to be python.

But ciphers are cool. And breaking them is cool. And if you think you might want to learn something about programming and you might want to learn something about ciphers, then this is a great book to read.

Sweigart has a knack for writing approachable descriptions of what’s going on without delving into too many details. In fact, in some sense Sweigart has already written this book before: his other books Automate the Boring Stuff with Python and Invent your own Computer Games with Python are similarly survey material using python as the medium, though with different motivating topics.

Each chapter of this book is centered around exploring a different aspect of a cipher, and introduces additional programming topics to do so. For example, one chapter introduces the classic Caesar cipher, as well as the “if”, “else”, and “elif” conditionals (and a few other python functions). Another chapter introduces brute-force breaking the Caesar cipher (as well as string formatting in python).

In each chapter, Sweigart begins by giving a high-level overview of the topics in that chapter, followed by python code which accomplishes the goal of the chapter, followed by a detailed description of what each block of code accomplishes. Readers get to see fully written code that does nontrivial things very quickly, but on the other hand the onus of code generation is entirely on the book and readers may have trouble adapting the concepts to other programming tasks. (But remember, this is more survey, less technical description). Further, Sweigart uses a number of good practices in his code that implicitly encourages good programming behaviors: modules are written with many well-named functions and well-named variables, and sufficiently modularly that earlier modules are imported and reused later.

But this book is not without faults. As with any survey material, one can always disagree on what topics are or are not included. The book covers five classical ciphers (Caesar, transposition, substitution, Vigenere, and affine) and one modern cipher (textbook-RSA), as well as the write-backwards cipher (to introduce python concepts) and the one-time-pad (presented oddly as a Vigenere cipher whose key is the same length as the message). For some unknown reason, Sweigart chooses to refer to RSA almost everywhere as “the public key cipher”, which I find both misleading (there are other public key ciphers) and giving insufficient attribution (the cipher is implemented in chapter 24, but “RSA” appears only once as a footnote in that chapter. Hopefully the reader was paying lots of attention, as otherwise it would be rather hard to find out more about it).

Further, the choice of python topics (and their order) is sometimes odd. In truth, this book is almost language agnostic and one could very easily adapt the presentation to any other scripting language, such as C.

In summary, this book is an excellent resource for the complete beginner who wants to learn something about programming and wants to learn something about ciphers. After reading this book, the reader will be a mid-beginner student of python (knee-deep is apt) and well-versed in classical ciphers. Should the reader feel inspired to learn more python, then he or she would probably feel comfortable diving into a tutorial or reference for their area of interest (like Full Stack Python if interested in web dev, or Python for Data Analysis if interested in data science). Or he or she might dive into a more complete monograph like Dive into Python or the monolithic Learn Python. Many fundamental topics (such as classes and objects, list comprehensions, data structures or algorithms) are not covered, and so “advanced” python resources would not be appropriate.

Further, should the reader feel inspired to learn more about cryptography, then I recommend that he or she consider Cryptanalysis by Gaines, which is a fun book aimed at diving deeper into classical pre-computer ciphers, or the slightly heavier but still fun resource would “Codebreakers” by Kahn. For much further cryptography, it’s necessary to develop a bit of mathematical maturity, which is its own hurdle.

This book is not appropriate for everyone. An experienced python programmer could read this book in an hour, skipping the various descriptions of how python works along the way. An experienced programmer who doesn’t know python could similarly read this book in a lazy afternoon. Both would probably do better reading either a more advanced overview of either cryptography or python, based on what originally drew them to the book.

## Hosting a Flask App on WebFaction on a Non-root Domain

Since I came to Warwick, I’ve been working extensively on the LMFDB, which uses python, sage, flask, and mongodb at its core. Thus I’ve become very familiar with flask. Writing a simple flask application is very quick and easy. So I thought it would be a good idea to figure out how to deploy a flask app on the server which runs this website, which is currently at WebFaction.

In short, it was not too hard, and now the app is set up for use. (It’s not a public tool, so I won’t link to it).

But there were a few things that I had to think figure out which I would quickly forget. Following the variety of information I found online, the only nontrivial aspect was configuring the site to run on a non-root domain (like davidlowryduda.com/subdomain instead of at davidlowryduda.com). I’m writing this so as to not need to figure this out when I write and hoost more flask apps. (Which I’ll almost certainly do, as it’s so straightforward).

There are some uninteresting things one must do on WebFaction.

2. Add a new application of type mod_wsgi (and the desired version of python, which is hopefully 3.6+).
3. Add this application to the desired website and subdomain in the WebFaction control panel.

After this, WebFaction will set up a skeleton “Hello World” mod_wsgi application with many reasonable server setting defaults. The remainder of the setup is done on the server itself.

In ~/webapps/application_name there will now appear

apache2/    # Apache config files and bin
htdocs/     # Default location where Apache looks for the app

We won’t change that structure. In htdocs1 there is a file index.py, which is where apache expects to find a python wsgi application called application. We will place the flask app along this structure and point to it in htdocs/index.py.

Usually I will use a virtualenv here. So in ~/webapps/application_name, I will run something like virtualenv flask_app_venv and virtualenv activate (or actually out of habit I frequently source the flask_app_venv/bin/activate file). Then pip install flask and whatever other python modules are necessary for the application to run. We will configure the server to use this virtual environment to run the app in a moment.

Copy the flask app, so that the resulting structure looks something like

~/webapps/application_name:

- apache2/
- htdocs/
- config.py
- libs/
- main/
- static/
- templates/
- __init__.py
- views.py
- models.my

I find it conceptually easiest if I have flask_app/main/__init__.py to directly contain the flask app to reference it by name in htdocs/index.py. It can be made elsewhere (for instance, perhaps in a file like flask_app/main/app.py, which appears to be a common structure), but I assume that it is at least imported in __init__.py.

For example, __init__.py might look something like

# application_name/flask_app/main/__init__.py

# ... other import statements from project if necessary

app.config.from_object('config')

# Importing the views for the rest of our site
# We do this here to avoid circular imports
# Note that I call it "main" where many call it "app"
from main import views

if __name__ == '__main__':
app.run()

The Flask constructor returns exactly the sort of wsgi application that apache expects. With this structure, we can edit the htdocs/index.py file to look like

# application_name/htdocs/index.py

import sys

# launching our app
from main import app as application

Now the server knows the correct wsgi_application to serve.

We must configure it to use our python virtual environment (and we’ll add a few additional convenience pieces). We edit /apache2/conf/httpd.conf as follows. Near the top of the file, certain modules are loaded. Add in the alias module, so that the modules look something like

#... other modules

This allows us to alias the root of the site. Since all site functionality is routed through htdocs/index.py, we want to think of the root / as beginning with /htdocs/index.py. At the end of the file

We now set the virtual environment to be used properly. There will be a set of lines containing names like WSGIDaemonProcess and WSGIProcessGroup. We edit these to refer to the correct python. WebFaction will have configured WSGIDaemonProcess to point to a local version of python by setting the python-path. Remove that, making that line look like

(or similar). We set the python path below, adding the line

I believe that this could also actually be done by setting puthon-path in WSGIDaemonProcess, but I find this more aesthetically pleasing.

We must also modify the “ section. Edit it to look like

<Directory>

It may very well be that I don't use the RewriteEngine at all, but if I do then this is where it's done. Script reloading is a nice convenience, especially while reloading and changing the app.

I note that it may be convenient to add an additional alias for static file hosting,

though I have not used this so far. (I get the same functionality through controlling the flask views appropriately).

The rest of this file has been setup by WebFaction for us upon creating the wsgi application.

## If the application is on a non-root domain...

If the application is to be run on a non-root domain, such as davidlowryduda.com/subdomain, then there is currently a problem. In flask, when using url getters like url_for, urls will be returned as though there is no subdomain. And thus all urls will be incorrect. It is necessary to alter provided urls in some way.

The way that worked for me was to insert a tiny bit of middleware in the wsgi_application. Alter htdocs/index.py to read

#application_name/htdocs/index.py

import sys

# subdomain url rerouting middleware
from webfaction_middleware import Middleware

from main import app

# set app through middleware
application = Middleware(app)

Now of course we need to write this middleware.

In application_name/flask_app, I create a file called webfaction_middleware.py, which reads

# application_name/flask_app/webfaction_middleware.py

class Middleware(object):  # python2 aware
def __init__(self, app):
self.app = app

def __call__(self, environ, start_response):
app_url = '/subdomain'
if app_url != '/':
environ['SCRIPT_NAME'] = app_url
return self.app(environ, start_response)

I now have a template file in which I keep app_url = '/' so that I can forget this and not worry, but that is where the subdomain url is prepended. Note that the leading slash is necessary. When I first tried using this, I omitted the leading slash. The application worked sometimes, and horribly failed in some other places. Some urls were correcty constructed, but most were not. I didn't try to figure out which ones were doomed to fail --- but it took me an embarassingly long time to realize that prepending a slash solved all problems.

The magical-names of environ and start_response are because the flask app is a wsgi_application, and this is the api of wsgi_applications generically.

Restart the apache server (/apache2/bin/restart) and go. Note that when incrementally making changes above, some changes can take a few minutes to fully propogate. It's only doing it the first time which takes some thought.

## Segregation, Gerrymandering, and Schelling’s Model

[This note is more about modeling some of the mathematics behind political events than politics themselves. And there are pretty pictures.]

Gerrymandering has become a recurring topic in the news. The Supreme Court of the US, as well as more state courts and supreme courts, is hearing multiple cases on partisan gerrymandering (all beginning with a case in Wisconsin).

Intuitively, it is clear that gerrymandering is bad. It allows politicians to choose their voters, instead of the other way around. And it allows the majority party to quash minority voices.

But how can one identify a gerrymandered map? To quote Justice Kennedy in his Concurrence the 2004 Supreme Court case Vieth v. Jubelirer:

When presented with a claim of injury from partisan gerrymandering, courts confront two obstacles. First is the lack of comprehensive and neutral principles for drawing electoral boundaries. No substantive definition of fairness in districting seems to command general assent. Second is the absence of rules to limit and confine judicial intervention. With uncertain limits, intervening courts–even when proceeding with best intentions–would risk assuming political, not legal, responsibility for a process that often produces ill will and distrust.

Later, he adds to the first obstacle, saying:

The object of districting is to establish “fair and effective representation for all citizens.” Reynolds v. Sims, 377 U.S. 533, 565—568 (1964). At first it might seem that courts could determine, by the exercise of their own judgment, whether political classifications are related to this object or instead burden representational rights. The lack, however, of any agreed upon model of fair and effective representation makes this analysis difficult to pursue.

From Justice Kennedy’s Concurrence emerges a theme — a “workable standard” of identifying gerrymandering would open up the possibility of limiting partisan gerrymandering through the courts. Indeed, at the core of the Wisconsin gerrymandering case is a proposed “workable standard”, based around the efficiency gap.

## Thomas Schelling and Segregation

In 1971, American economist Thomas Schelling (who later won the Nobel Prize in Economics in 2005) published Dynamic Models of Segregation (Journal of Mathematical Sociology, 1971, Vol 1, pp 143–186). He sought to understand why racial segregation in the United States seems so difficult to combat.

He introduced a simple model of segregation suggesting that even if each individual person doesn’t mind living with others of a different race, they might still choose to segregate themselves through mild preferences. As each individual makes these choices, overall segregation increases.

I write this post because I wondered what happens if we adapt Schelling’s model to instead model a state and its district voting map. In place of racial segregation, I consider political segregation. Supposing the district voting map does not change, I wondered how the efficiency gap will change over time as people further segregate themselves.

It seemed intuitive to me that political segregation (where people who had the same political beliefs stayed largely together and separated from those with different political beliefs) might correspond to more egregious cases of gerrymandering. But to my surprise, I was (mostly) wrong.

Let’s set up and see the model.

## Advent of Code: Day 4

This is a very short post in my collection working through this year’s Advent of Code challenges. Unlike the previous ones, this has no mathematical comments, as it was a very short exercise. This notebook is available in its original format on my github.

# Day 4: High Entropy Passphrases¶

Given a list of strings, determine how many strings have no duplicate words.

This is a classic problem, and it’s particularly easy to solve this in python. Some might use collections.Counter, but I think it’s more straightforward to use sets.

The key idea is that the set of words in a sentence will not include duplicates. So if taking the set of a sentence reduces its length, then there was a duplicate word.

In [1]:
with open("input.txt", "r") as f:

def count_lines_with_unique_words(lines):
num_pass = 0
for line in lines:
s = line.split()
if len(s) == len(set(s)):
num_pass += 1
return num_pass

count_lines_with_unique_words(lines)

Out[1]:
455

I think this is the first day where I would have had a shot at the leaderboard if I’d been gunning for it.

# Part 2¶

Let’s add in another constraint. Determine how many strings have no duplicate words, even after anagramming. Thus the string

abc bac

is not valid, since the second word is an anagram of the first. There are many ways to tackle this as well, but I will handle anagrams by sorting the letters in each word first, and then running the bit from part 1 to identify repeated words.

In [2]:
with open("input.txt", "r") as f:

sorted_lines = []
for line in lines:
sorted_line = ' '.join([''.join(l) for l in map(sorted, line.split())])
sorted_lines.append(sorted_line)

sorted_lines[:2]

Out[2]:
['bddjjow acimrv bcjjm anr flmmos fiosv',
'bcmnoxy dfinyzz dgmp dfgioy hinrrv eeklpuu adgpw kqv']
In [3]:
count_lines_with_unique_words(sorted_lines)

Out[3]:
186
Posted in Expository, Programming, Python | Tagged , , | 1 Comment

## Advent of Code: Day 3

This is the third notebook in my posts on the Advent of Code challenges. The notebook in its original format can be found on my github.

# Day 3: Spiral Memory¶

Numbers are arranged in a spiral

17  16  15  14  13
18   5   4   3  12
19   6   1   2  11
20   7   8   9  10
21  22  23---> ...

Given an integer n, what is its Manhattan Distance from the center (1) of the spiral? For instance, the distance of 3 is $2 = 1 + 1$, since it’s one space to the right and one space up from the center.

Here’s my idea. The bottom right corner of the $k$th layer is the integer $(2k+1)^2$, since that’s how many integers are contained within that square. The other three corners in that layer are $(2k+1)^2 – 2k, (2k+1)^2 – 4k$, and $(2k+1)^2 – 6k$. Finally, the closest spot on the $k$th layer to the origin is at distance $k$: these are the four “axis locations” halfway between the corners, at $(2k+1)^2 – k, (2k+1)^2 – 3k, (2k+1)^2 – 5k$, and $(2k+1)^2 – 7k$.

For instance when $k = 1$, the bottom right is $(2 + 1)^2 = 9$, and the four “axis locations” are $9 – 1, 9 – 3, 9-5$, and $9-7$. The “axis locations” are $k$ away, and the corners are $2k$ away.

So I will first find which layer the number is on. Then I’ll figure out which side it’s on, and then how far away it is from the nearest “axis location” or “corner”.

My given number happens to be 289326.

In [1]:
import math

def find_lowest_larger_odd_square(n):
upper = math.ceil(n**.5)
if upper %2 == 0:
upper += 1
return upper

In [2]:
assert find_lowest_larger_odd_square(39) == 7
assert find_lowest_larger_odd_square(26) == 7
assert find_lowest_larger_odd_square(25) == 5

In [3]:
find_lowest_larger_odd_square(289326)

Out[3]:
539
In [4]:
539**2 - 289326

Out[4]:
1195

It happens to be that our integer is very close to an odd square.
The square is $539^2$, and the distance to that square is $538$ from the center.

Note that $539 = 2(269) + 1$, so this is the $269$th layer of the square.
The previous corner to $539^2$ is $539^2 – 538$, and the previous corner to that is $539^2 – 2\cdot538 = 539^2 – 1076$.
This is the nearest corner.
How far away from the square is this corner?

## Advent of Code: Day 2

This is the second notebook in my posts on the Advent of Code challenges. This notebook in its original format can be found on my github.

# Day 2: Corruption Checksum, part I¶

You are given a table of integers. Find the difference between the maximum and minimum of each row, and add these differences together.

There is not a lot to say about this challenge. The plan is to read the file linewise, compute the difference on each line, and sum them up.

In [1]:
with open("input.txt", "r") as f:
lines[0]

Out[1]:
'5048\t177\t5280\t5058\t4504\t3805\t5735\t220\t4362\t1809\t1521\t230\t772\t1088\t178\t1794\n'
In [2]:
l = lines[0]
l = l.split()
l

Out[2]:
['5048',
'177',
'5280',
'5058',
'4504',
'3805',
'5735',
'220',
'4362',
'1809',
'1521',
'230',
'772',
'1088',
'178',
'1794']
In [3]:
def max_minus_min(line):
'''Compute the difference between the largest and smallest integer in a line'''
line = list(map(int, line.split()))
return max(line) - min(line)

def sum_differences(lines):
'''Sum the value of max_minus_min for each line in lines'''
return sum(max_minus_min(line) for line in lines)

In [4]:
testcase = ['5 1 9 5','7 5 3', '2 4 6 8']
assert sum_differences(testcase) == 18

In [5]:
sum_differences(lines)

Out[5]:
58975

## Mathematical Interlude¶

In line with the first day’s challenge, I’m inclined to ask what we should “expect.” But what we should expect is not well-defined in this case. Let us rephrase the problem in a randomized sense.

Suppose we are given a table, $n$ lines long, where each line consists of $m$ elements, that are each uniformly randomly chosen integers from $1$ to $10$. We might ask what is the expected value of this operation, of summing the differences between the maxima and minima of each row, on this table. What should we expect?

As each line is independent of the others, we are really asking what is the expected value across a single row. So given $m$ integers uniformly randomly chosen from $1$ to $10$, what is the expected value of the maximum, and what is the expected value of the minimum?

### Expected Minimum¶

Let’s begin with the minimum. The minimum is $1$ unless all the integers are greater than $2$. This has probability
$$1 – \left( \frac{9}{10} \right)^m = \frac{10^m – 9^m}{10^m}$$
of occurring. We rewrite it as the version on the right for reasons that will soon be clear.
The minimum is $2$ if all the integers are at least $2$ (which can occur in $9$ different ways for each integer), but not all the integers are at least $3$ (each integer has $8$ different ways of being at least $3$). Thus this has probability
$$\frac{9^m – 8^m}{10^m}.$$
Continuing to do one more for posterity, the minimum is $3$ if all the integers are at least $3$ (each integer has $8$ different ways of being at least $3$), but not all integers are at least $4$ (each integer has $7$ different ways of being at least $4$). Thus this has probability

$$\frac{8^m – 7^m}{10^m}.$$

And so on.

Recall that the expected value of a random variable is

$$E[X] = \sum x_i P(X = x_i),$$

so the expected value of the minimum is

$$\frac{1}{10^m} \big( 1(10^m – 9^m) + 2(9^m – 8^m) + 3(8^m – 7^m) + \cdots + 9(2^m – 1^m) + 10(1^m – 0^m)\big).$$

This simplifies nicely to

$$\sum_ {k = 1}^{10} \frac{k^m}{10^m}.$$

### Expected Maximum¶

The same style of thinking shows that the expected value of the maximum is

$$\frac{1}{10^m} \big( 10(10^m – 9^m) + 9(9^m – 8^m) + 8(8^m – 7^m) + \cdots + 2(2^m – 1^m) + 1(1^m – 0^m)\big).$$

This simplifies to

$$\frac{1}{10^m} \big( 10 \cdot 10^m – 9^m – 8^m – \cdots – 2^m – 1^m \big) = 10 – \sum_ {k = 1}^{9} \frac{k^m}{10^m}.$$

### Expected Difference¶

Subtracting, we find that the expected difference is

$$9 – 2\sum_ {k=1}^{9} \frac{k^m}{10^m}.$$

From this we can compute this for each list-length $m$. It is good to note that as $m \to \infty$, the expected value is $9$. Does this make sense? Yes, as when there are lots of values we should expect one to be a $10$ and one to be a $1$. It’s also pretty straightforward to see how to extend this to values of integers from $1$ to $N$.

Looking at the data, it does not appear that the integers were randomly chosen. Instead, there are very many relatively small integers and some relatively large integers. So we shouldn’t expect this toy analysis to accurately model this problem — the distribution is definitely not uniform random.
But we can try it out anyway.

## Advent of Code: Day 1

I thoroughly enjoyed reading through Peter Norvig’s extraordinarily clean and nice solutions to the Advent of Code challenge last year. Inspired by his clean, literate programming style and the convenience of jupyter notebook demonstrations, I will look at several of these challenges in my own jupyter notebooks.

My background and intentions aren’t the same as Peter Norvig’s: his expertise dwarfs mine. And timezones are not kind to those of us in the UK, and thus I won’t be competing for a position on the leaderboards. These are to be fun. And sometimes there are tidbits of math that want to come out of the challenges.

Enough of that. Let’s dive into the first day.

# Day 1: Inverse Captcha, Part 1¶

Given a sequence of digits, find the sum of those digits which match the following digit. The sequence is presumed circular, so the first digit may match the last digit.

This would probably be done the fastest by looping through the sequence.

In [1]:
with open('input.txt', 'r') as f:
seq = seq.strip()
seq[:10]

Out[1]:
'1118313623'
In [2]:
def sum_matched_digits(s):
"Sum of digits which match following digit, and first digit if it matches last digit"
total = 0
for a,b in zip(s, s[1:]+s[0]):
if a == b:
total += int(a)

They provide a few test cases which we use to test our method against.

In [3]:
assert sum_matched_digits('1122') == 3
assert sum_matched_digits('1111') == 4
assert sum_matched_digits('1234') == 0
assert sum_matched_digits('91212129') == 9

For fun, this is a oneline version.

Posted in Expository, Programming, Python | Tagged , , | 1 Comment

# Interfacing sage and the LMFDB — a prototype¶

The lmfdb and sagemath are both great things, but they don’t currently talk to each other. Much of the lmfdb calls sage, but the lmfdb also includes vast amounts of data on $L$-functions and modular forms (hence the name) that is not accessible from within sage.
This is an example prototype of an interface to the lmfdb from sage. Keep in mind that this is a prototype and every aspect can change. But we hope to show what may be possible in the future. If you have requests, comments, or questions, please request/comment/ask either now, or at my email: david@lowryduda.com.

Note that this notebook is available on http://davidlowryduda.com or https://gist.github.com/davidlowryduda/deb1f88cc60b6e1243df8dd8f4601cde, and the code is available at https://github.com/davidlowryduda/sage2lmfdb

Let’s dive into an example.

In [1]:
# These names will change
from sage.all import *
import LMFDB2sage.elliptic_curves as lmfdb_ecurve

In [2]:
lmfdb_ecurve.search(rank=1)

Out[2]:
[Elliptic Curve defined by y^2 + x*y = x^3 - 887688*x - 321987008 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 + 10795*x - 97828 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 - 2294115305*x - 42292668425178 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 - 3170*x - 49318 over Rational Field,
Elliptic Curve defined by y^2 + y = x^3 + 1050*x - 26469 over Rational Field,
Elliptic Curve defined by y^2 + x*y = x^3 - x^2 - 1240542*x - 531472509 over Rational Field,
Elliptic Curve defined by y^2 + y = x^3 - x^2 + 8100*x - 263219 over Rational Field,
Elliptic Curve defined by y^2 + x*y = x^3 + 637*x - 68783 over Rational Field,
Elliptic Curve defined by y^2 + y = x^3 + x^2 + 36*x - 380 over Rational Field,
Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2535*x - 49982 over Rational Field]
This returns 10 elliptic curves of rank 1. But these are a bit different than sage’s elliptic curves.

In [3]:
Es = lmfdb_ecurve.search(rank=1)
E = Es[0]
print(type(E))

<class 'LMFDB2sage.ell_lmfdb.EllipticCurve_rational_field_lmfdb_with_category'>

Note that the class of an elliptic curve is an lmfdb ElliptcCurve. But don’t worry, this is a subclass of a normal elliptic curve. So we can call the normal things one might call on an elliptic curve.

th

In [4]:
# Try autocompleting the following. It has all the things!
print(dir(E))

['CPS_height_bound', 'CartesianProduct',
'Chow_form', 'Hom',
'Jacobian', 'Jacobian_matrix',
'Lambda', 'Np',
'S_integral_points', '_AlgebraicScheme__A',
'_AlgebraicScheme__divisor_group', '_AlgebraicScheme_subscheme__polys',
'_EllipticCurve_generic__ainvs', '_EllipticCurve_generic__b_invariants',
'_EllipticCurve_generic__base_ring', '_EllipticCurve_generic__discriminant',
'_EllipticCurve_generic__is_over_RationalField', '_EllipticCurve_generic__multiple_x_denominator',
'_EllipticCurve_generic__multiple_x_numerator', '_EllipticCurve_rational_field__conductor_pari',
'_EllipticCurve_rational_field__generalized_congruence_number', '_EllipticCurve_rational_field__generalized_modular_degree',
'_EllipticCurve_rational_field__gens', '_EllipticCurve_rational_field__modular_degree',
'_EllipticCurve_rational_field__np', '_EllipticCurve_rational_field__rank',
'_EllipticCurve_rational_field__regulator', '_EllipticCurve_rational_field__torsion_order',
'__class__', '__cmp__', '__contains__', '__delattr__',
'__dict__', '__dir__', '__div__', '__doc__',
'__eq__', '__format__', '__ge__', '__getattribute__',
'__getitem__', '__getstate__', '__gt__', '__hash__',
'__init__', '__le__', '__lt__', '__make_element_class__',
'__module__', '__mul__', '__ne__', '__new__',
'__nonzero__', '__pari__', '__pow__', '__pyx_vtable__',
'__rdiv__', '__reduce__', '__reduce_ex__', '__repr__',
'__rmul__', '__setattr__', '__setstate__', '__sizeof__',
'__str__', '__subclasshook__', '__temporarily_change_names', '__truediv__',
'_ascii_art_', '_assign_names', '_axiom_', '_axiom_init_',
'_base', '_base_ring', '_base_scheme', '_best_affine_patch',
'_cache__point_homset', '_cache_an_element', '_cache_key', '_check_satisfies_equations',
'_cmp_', '_coerce_map_from_', '_coerce_map_via', '_coercions_used',
'_compute_gens', '_convert_map_from_', '_convert_method_name', '_defining_names',
'_defining_params_', '_doccls', '_element_constructor', '_element_constructor_',
'_element_constructor_from_element_class', '_element_init_pass_parent', '_factory_data', '_first_ngens',
'_forward_image', '_fricas_', '_fricas_init_', '_gap_',
'_gap_init_', '_generalized_congmod_numbers', '_generic_coerce_map', '_generic_convert_map',
'_get_action_', '_get_local_data', '_giac_', '_giac_init_',
'_gp_', '_gp_init_', '_heegner_best_tau', '_heegner_forms_list',
'_heegner_index_in_EK', '_homset', '_init_category_', '_initial_action_list',
'_initial_coerce_list', '_initial_convert_list', '_interface_', '_interface_init_',
'_interface_is_cached_', '_internal_coerce_map_from', '_internal_convert_map_from', '_introspect_coerce',
'_is_category_initialized', '_is_valid_homomorphism_', '_isoclass', '_json',
'_kash_', '_kash_init_', '_known_points', '_latex_',
'_lmfdb_label', '_lmfdb_regulator', '_macaulay2_', '_macaulay2_init_',
'_magma_init_', '_maple_', '_maple_init_', '_mathematica_',
'_mathematica_init_', '_maxima_', '_maxima_init_', '_maxima_lib_',
'_maxima_lib_init_', '_modsym', '_modular_symbol_normalize', '_morphism',
'_multiple_of_degree_of_isogeny_to_optimal_curve', '_multiple_x_denominator', '_multiple_x_numerator', '_names',
'_pari_', '_pari_init_', '_point', '_point_homset',
'_polymake_', '_polymake_init_', '_populate_coercion_lists_', '_r_init_',
'_reduce_model', '_reduce_point', '_reduction', '_refine_category_',
'_repr_', '_repr_option', '_repr_type', '_sage_',
'_scale_by_units', '_set_conductor', '_set_cremona_label', '_set_element_constructor',
'_set_gens', '_set_modular_degree', '_set_rank', '_set_torsion_order',
'_shortest_paths', '_singular_', '_singular_init_', '_symbolic_',
'_test_an_element', '_test_cardinality', '_test_category', '_test_elements',
'_test_elements_eq_reflexive', '_test_elements_eq_symmetric', '_test_elements_eq_transitive', '_test_elements_neq',
'_test_eq', '_test_new', '_test_not_implemented_methods', '_test_pickling',
'_test_some_elements', '_tester', '_torsion_bound', '_unicode_art_',
'_unset_category', '_unset_coercions_used', '_unset_embedding', 'a1',
'a2', 'a3', 'a4', 'a6',
'a_invariants', 'abelian_variety', 'affine_patch', 'ainvs',
'algebra', 'ambient_space', 'an', 'an_element',
'analytic_rank', 'analytic_rank_upper_bound', 'anlist', 'antilogarithm',
'ap', 'aplist', 'arithmetic_genus', 'automorphisms',
'b2', 'b4', 'b6', 'b8',
'b_invariants', 'base', 'base_extend', 'base_field',
'base_morphism', 'base_ring', 'base_scheme', 'c4',
'c6', 'c_invariants', 'cartesian_product', 'categories',
'category', 'change_ring', 'change_weierstrass_model', 'cm_discriminant',
'codimension', 'coerce', 'coerce_embedding', 'coerce_map_from',
'complement', 'conductor', 'congruence_number', 'construction',
'convert_map_from', 'coordinate_ring', 'count_points', 'cremona_label',
'database_attributes', 'database_curve', 'db', 'defining_ideal',
'defining_polynomial', 'defining_polynomials', 'degree', 'descend_to',
'dimension', 'dimension_absolute', 'dimension_relative', 'discriminant',
'division_field', 'division_polynomial', 'division_polynomial_0', 'divisor',
'divisor_group', 'divisor_of_function', 'dual', 'dump',
'dumps', 'element_class', 'elliptic_exponential', 'embedding_center',
'embedding_morphism', 'eval_modular_form', 'excellent_position', 'formal',
'formal_group', 'fundamental_group', 'galois_representation', 'gen',
'gens', 'gens_certain', 'gens_dict', 'gens_dict_recursive',
'genus', 'geometric_genus', 'get_action', 'global_integral_model',
'has_base', 'has_cm', 'has_coerce_map_from', 'has_global_minimal_model',
'has_good_reduction', 'has_good_reduction_outside_S', 'has_multiplicative_reduction', 'has_nonsplit_multiplicative_reduction',
'has_rational_cm', 'has_split_multiplicative_reduction', 'hasse_invariant', 'heegner_discriminants',
'heegner_discriminants_list', 'heegner_index', 'heegner_index_bound', 'heegner_point',
'heegner_point_height', 'heegner_sha_an', 'height', 'height_function',
'height_pairing_matrix', 'hom', 'hyperelliptic_polynomials', 'identity_morphism',
'inject_variables', 'integral_model', 'integral_points', 'integral_short_weierstrass_model',
'integral_weierstrass_model', 'integral_x_coords_in_interval', 'intersection', 'intersection_multiplicity',
'intersection_points', 'intersects_at', 'irreducible_components', 'is_atomic_repr',
'is_coercion_cached', 'is_complete_intersection', 'is_conversion_cached', 'is_exact',
'is_global_integral_model', 'is_global_minimal_model', 'is_good', 'is_integral',
'is_irreducible', 'is_isogenous', 'is_isomorphic', 'is_local_integral_model',
'is_minimal', 'is_on_curve', 'is_ordinary', 'is_ordinary_singularity',
'is_p_integral', 'is_p_minimal', 'is_parent_of', 'is_projective',
'is_singular', 'is_smooth', 'is_supersingular', 'is_transverse',
'is_x_coord', 'isogenies_prime_degree', 'isogeny', 'isogeny_class',
'isogeny_codomain', 'isogeny_degree', 'isogeny_graph', 'isomorphism_to',
'isomorphisms', 'j_invariant', 'kodaira_symbol', 'kodaira_type',
'kodaira_type_old', 'kolyvagin_point', 'label', 'latex_name',
'latex_variable_names', 'lift_x', 'lll_reduce', 'lmfdb_page',
'local_coordinates', 'local_data', 'local_integral_model', 'local_minimal_model',
'lseries', 'lseries_gross_zagier', 'manin_constant', 'matrix_of_frobenius',
'modular_degree', 'modular_form', 'modular_parametrization', 'modular_symbol',
'modular_symbol_numerical', 'modular_symbol_space', 'multiplication_by_m', 'multiplication_by_m_isogeny',
'multiplicity', 'mwrank', 'mwrank_curve', 'neighborhood',
'newform', 'ngens', 'non_minimal_primes', 'nth_iterate',
'objgen', 'objgens', 'optimal_curve', 'orbit',
'pari_mincurve', 'period_lattice', 'plane_projection', 'plot',
'point', 'point_homset', 'point_search', 'point_set',
'pollack_stevens_modular_symbol', 'preimage', 'projection', 'prove_BSD',
'quartic_twist', 'rank', 'rank_bound', 'rank_bounds',
'rational_parameterization', 'rational_points', 'real_components', 'reduce',
'reduction', 'register_action', 'register_coercion', 'register_conversion',
'register_embedding', 'regulator', 'regulator_of_points', 'rename',
'reset_name', 'root_number', 'rst_transform', 'satisfies_heegner_hypothesis',
'saturation', 'save', 'scale_curve', 'selmer_rank',
'sextic_twist', 'sha', 'short_weierstrass_model', 'silverman_height_bound',
'simon_two_descent', 'singular_points', 'singular_subscheme', 'some_elements',
'specialization', 'structure_morphism', 'supersingular_primes', 'tamagawa_exponent',
'tamagawa_number', 'tamagawa_number_old', 'tamagawa_numbers', 'tamagawa_product',
'tamagawa_product_bsd', 'tangents', 'tate_curve', 'three_selmer_rank',
'torsion_order', 'torsion_points', 'torsion_polynomial', 'torsion_subgroup',
'two_descent', 'two_descent_simon', 'two_division_polynomial', 'two_torsion_rank',
'union', 'variable_name', 'variable_names', 'weierstrass_p',
'weil_restriction', 'zeta_series']

This gives quick access to some data that is not stored within the LMFDB, but which is relatively quickly computable. For example,

In [5]:
E.defining_ideal()

Out[5]:
Ideal (-x^3 + x*y*z + y^2*z + 887688*x*z^2 + 321987008*z^3) of Multivariate Polynomial Ring in x, y, z over Rational Field
But one of the great powers is that there are some things which are computed and stored in the LMFDB, and not in sage. We can now immediately give many examples of rank 3 elliptic curves with:

In [6]:
Es = lmfdb_ecurve.search(conductor=11050, torsion_order=2)
print("There are {} curves returned.".format(len(Es)))
E = Es[0]
print(E)

There are 10 curves returned.
Elliptic Curve defined by y^2 + x*y + y = x^3 - 3476*x - 79152 over Rational Field

And for these curves, the lmfdb contains data on its rank, generators, regulator, and so on.

In [7]:
print(E.gens())
print(E.rank())
print(E.regulator())

[(-34 : 17 : 1)]
1
1.63852610029

In [8]:
res = []
%time for E in Es: res.append(E.gens()); res.append(E.rank()); res.append(E.regulator())

CPU times: user 971 ms, sys: 6.82 ms, total: 978 ms
Wall time: 978 ms

That’s pretty fast, and this is because all of this was pulled from the LMFDB when the curves were returned by the search() function.
In this case, elliptic curves over the rationals are only an okay example, as they’re really well studied and sage can compute much of the data very quickly. On the other hand, through the LMFDB there are millions of examples and corresponding data at one’s fingertips.

### This is where we’re really looking for input.¶

Think of what you might want to have easy access to through an interface from sage to the LMFDB, and tell us. We’re actively seeking comments, suggestions, and requests. Elliptic curves over the rationals are a prototype, and the LMFDB has lots of (much more challenging to compute) data. There is data on the LMFDB that is simply not accessible from within sage.
email: david@lowryduda.com, or post an issue on https://github.com/LMFDB/lmfdb/issues

## Now let’s describe what’s going on under the hood a little bit¶

There is an API for the LMFDB at http://beta.lmfdb.org/api/. This API is a bit green, and we will change certain aspects of it to behave better in the future. A call to the API looks like

http://beta.lmfdb.org/api/elliptic_curves/curves/?rank=i1&conductor=i11050

The result is a large mess of data, which can be exported as json and parsed.
But that’s hard, and the resulting data are not sage objects. They are just strings or ints, and these require time and thought to parse.
So we created a module in sage that writes the API call and parses the output back into sage objects. The 22 curves given by the above API call are the same 22 curves returned by this call:

In [9]:
Es = lmfdb_ecurve.search(rank=1, conductor=11050, max_items=25)
print(len(Es))
E = Es[0]

22

The total functionality of this search function is visible from its current documentation.

In [10]:
# Execute this cell for the documentation
print(lmfdb_ecurve.search.__doc__)

    Search the LMFDB for an elliptic curve.

Note that all inputs are optional, but at least one input is necessary.

INPUT:

-  label=l -- a string l representing a label in the LMFDB.

-  degree=d -- an int d giving the minimum degree of a
parameterization of the modular curve

-  conductor=c -- an int c giving the conductor of the curve

-  min_conductor=mc -- an int mc giving a lower bound on the
conductor for desired curves

-  max_conductor=mc -- an int mc giving an upper bound on the
conductor for desired curves

-  torsion_order=t -- an int t giving the order of the torsion
subgroup of the curve

-  rank=r -- an int r giving the rank of the curve

-  regulator=f -- a float f giving the regulator of the curve

-  max_items=m -- an int m (default: 10, max: 100) indicating the
maximum number of results to return

-  base_item=b -- an int b (default: 0) specifying where to start
returning values from. The search will begin by returning the bth
curve. Combined with max_items to return data in chunks.

-  sort=s -- a string s specifying what database field to sort the

EXAMPLES::

sage: Es = search(conductor=11050, rank=2)
[Elliptic Curve defined by y^2 + x*y = x^3 - x^2 - 442*x + 1716 over Rational Field, Elliptic Curve defined by y^2 + x*y = x^3 - x^2 + 1558*x + 11716 over Rational Field]
sage: E = E[0]
sage: E.conductor()
11050

In [11]:
# So, for instance, one could perform the following search, finding a unique elliptic curve
lmfdb_ecurve.search(rank=2, torsion_order=3, degree=4608)

Out[11]:
[Elliptic Curve defined by y^2 + y = x^3 + x^2 - 5155*x + 140756 over Rational Field]

### What if there are no curves?¶

If there are no curves satisfying the search criteria, then a message is displayed and that’s that. These searches may take a couple of seconds to complete.
For example, no elliptic curve in the database has rank 5.

In [12]:
lmfdb_ecurve.search(rank=5)

No fields were found satisfying input criteria.

### How does one step through the data?¶

Right now, at most 100 curves are returned in a single API call. This is the limit even from directly querying the API. But one can pass in the argument base_item (the name will probably change… to skip? or perhaps to offset?) to start returning at the base_itemth element.

In [13]:
from pprint import pprint
pprint(lmfdb_ecurve.search(rank=1, max_items=3))              # The last item in this list
print('')
pprint(lmfdb_ecurve.search(rank=1, max_items=3, base_item=2)) # should be the first item in this list

[Elliptic Curve defined by y^2 + x*y = x^3 - 887688*x - 321987008 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 + 10795*x - 97828 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 - 2294115305*x - 42292668425178 over Rational Field]

[Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 - 2294115305*x - 42292668425178 over Rational Field,
Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 - 3170*x - 49318 over Rational Field,
Elliptic Curve defined by y^2 + y = x^3 + 1050*x - 26469 over Rational Field]

Included in the documentation is also a bit of hopefulness. Right now, the LMFDB API does not actually accept max_conductor or min_conductor (or arguments of that type). But it will sometime. (This introduces a few extra difficulties on the server side, and so it will take some extra time to decide how to do this).

In [14]:
lmfdb_ecurve.search(rank=1, min_conductor=500, max_conductor=10000)  # Not implemented

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-14-3d98f2cf7a13> in <module>()
----> 1 lmfdb_ecurve.search(rank=Integer(1), min_conductor=Integer(500), max_conductor=Integer(10000))  # Not implemented

/home/djlowry/Dropbox/EllipticCurve_LMFDB/LMFDB2sage/elliptic_curves.py in search(**kwargs)
76             kwargs[item]
77             raise NotImplementedError("This would be a great thing to have, " +
---> 78                 "but the LMFDB api does not yet provide this functionality.")
79         except KeyError:
80             pass

NotImplementedError: This would be a great thing to have, but the LMFDB api does not yet provide this functionality.
Our EllipticCurve_rational_field_lmfdb class constructs a sage elliptic curve from the json and overrides (somem of the) the default methods in sage if there is quicker data available on the LMFDB. In principle, this new object is just a sage object with some slightly different methods.
Generically, documentation and introspection on objects from this class should work. Much of sage’s documentation carries through directly.

In [15]:
print(E.gens.__doc__)

        Return generators for the Mordell-Weil group E(Q) *modulo*
torsion.

.. warning::

If the program fails to give a provably correct result, it
prints a warning message, but does not raise an
exception. Use :meth:~gens_certain to find out if this
warning message was printed.

INPUT:

- proof -- bool or None (default None), see
proof.elliptic_curve or sage.structure.proof

- verbose - (default: None), if specified changes the
verbosity of mwrank computations

- rank1_search - (default: 10), if the curve has analytic
rank 1, try to find a generator by a direct search up to
this logarithmic height.  If this fails, the usual mwrank
procedure is called.

- algorithm -- one of the following:

- 'mwrank_shell' (default) -- call mwrank shell command

- 'mwrank_lib' -- call mwrank C library

- only_use_mwrank -- bool (default True) if False, first
attempts to use more naive, natively implemented methods

- use_database -- bool (default True) if True, attempts to
find curve and gens in the (optional) database

- descent_second_limit -- (default: 12) used in 2-descent

- sat_bound -- (default: 1000) bound on primes used in
saturation.  If the computed bound on the index of the
points found by two-descent in the Mordell-Weil group is
greater than this, a warning message will be displayed.

OUTPUT:

- generators - list of generators for the Mordell-Weil
group modulo torsion

IMPLEMENTATION: Uses Cremona's mwrank C library.

EXAMPLES::

sage: E = EllipticCurve('389a')
sage: E.gens()                 # random output
[(-1 : 1 : 1), (0 : 0 : 1)]

A non-integral example::

sage: E = EllipticCurve([-3/8,-2/3])
sage: E.gens() # random (up to sign)
[(10/9 : 29/54 : 1)]

A non-minimal example::

sage: E = EllipticCurve('389a1')
sage: E1 = E.change_weierstrass_model([1/20,0,0,0]); E1
Elliptic Curve defined by y^2 + 8000*y = x^3 + 400*x^2 - 320000*x over Rational Field
sage: E1.gens() # random (if database not used)
[(-400 : 8000 : 1), (0 : -8000 : 1)]

Modified methods should have a note indicating that the data comes from the LMFDB, and then give sage’s documentation. This is not yet implemented. (So if you examine the current version, you can see some incomplete docstrings like regulator().)

In [16]:
print(E.regulator.__doc__)

        Return the regulator of the curve. This is taken from the lmfdb if available.

NOTE:
In later implementations, this docstring will probably include the
docstring from sage's regular implementation. But that's not
currently the case.

## This concludes our demo of an interface between sage and the LMFDB.¶

Thank you, and if you have any questions, comments, or concerns, please find me/email me/raise an issue on LMFDB’s github.

## A Brief Notebook on Cryptography

This is a post written for my Spring 2016 Number Theory class at Brown University. It was written and originally presented from a Jupyter Notebook, and changed for presentation on this site. There is a version of the notebook available at github.

# Cryptography¶

Recall the basic setup of cryptography. We have two people, Anabel and Bartolo. Anabel wants to send Bartolo a secure message. What do we mean by “secure?” We mean that even though that dastardly Eve might intercept and read the transmitted message, Eve won’t learn anything about the actual message Anabel wants to send to Bartolo.

Before the 1970s, the only way for Anabel to send Bartolo a secure message required Anabel and Bartolo to get together beforehand and agree on a secret method of communicating. To communicate, Anabel first decides on a message. The original message is sometimes called the plaintext. She then encrypts the message, producing a ciphertext.

She then sends the ciphertext. If Eve gets hold of hte ciphertext, she shouldn’t be able to decode it (unless it’s a poor encryption scheme).

When Bartolo receives the ciphertext, he can decrypt it to recover the plaintext message, since he agreed on the encryption scheme beforehand.

## Caesar Shift¶

The first known instance of cryptography came from Julius Caesar. It was not a very good method. To encrypt a message, he shifted each letter over by 2, so for instance “A” becomes “C”, and “B” becomes “D”, and so on.

Let’s see what sort of message comes out.

alpha = "abcdefghijklmnopqrstuvwxyz".upper()
punct = ",.?:;'\n "

from functools import partial

def shift(l, s=2):
l = l.upper()
return alpha[(alpha.index(l) + s) % 26]

def caesar_shift_encrypt(m, s=2):
m = m.upper()
c = "".join(map(partial(shift, s=s), m))
return c

def caesar_shift_decrypt(c, s=-2):
c = c.upper()
m = "".join(map(partial(shift, s=s), c))
return m

print "Original Message: HI"
print "Ciphertext:", caesar_shift_encrypt("hi")

Original Message: HI
Ciphertext: JK

m = """To be, or not to be, that is the question:
Whether 'tis Nobler in the mind to suffer
The Slings and Arrows of outrageous Fortune,
Or to take Arms against a Sea of troubles,
And by opposing end them."""

m = "".join([l for l in m if not l in punct])

print "Original Message:"
print m

print
print "Ciphertext:"
tobe_ciphertext = caesar_shift_encrypt(m)
print tobe_ciphertext

Original Message:
TobeornottobethatisthequestionWhethertisNoblerinthemindtosuffer
TheSlingsandArrowsofoutrageousFortuneOrtotakeArmsagainstaSeaof
troublesAndbyopposingendthem

Ciphertext:
VQDGQTPQVVQDGVJCVKUVJGSWGUVKQPYJGVJGTVKUPQDNGTKPVJGOKPFVQUWHHGT
VJGUNKPIUCPFCTTQYUQHQWVTCIGQWUHQTVWPGQTVQVCMGCTOUCICKPUVCUGCQH
VTQWDNGUCPFDAQRRQUKPIGPFVJGO

print "Decrypted first message:", caesar_shift_decrypt("JK")

Decrypted first message: HI

print "Decrypted second message:"
print caesar_shift_decrypt(tobe_ciphertext)

Decrypted second message:
TOBEORNOTTOBETHATISTHEQUESTIONWHETHERTISNOBLERINTHEMINDTOSUFFER
THESLINGSANDARROWSOFOUTRAGEOUSFORTUNEORTOTAKEARMSAGAINSTASEAOF
TROUBLESANDBYOPPOSINGENDTHEM

Is this a good encryption scheme? No, not really. There are only 26 different things to try. This can be decrypted very quickly and easily, even if entirely by hand.

## Substitution Cipher¶

A slightly different scheme is to choose a different letter for each letter. For instance, maybe “A” actually means “G” while “B” actually means “E”. We call this a substitution cipher as each letter is substituted for another.

import random
permutation = list(alpha)
random.shuffle(permutation)

# Display the new alphabet
print alpha
subs = "".join(permutation)
print subs

ABCDEFGHIJKLMNOPQRSTUVWXYZ
EMJSLZKAYDGTWCHBXORVPNUQIF

def subs_cipher_encrypt(m):
m = "".join([l.upper() for l in m if not l in punct])
return "".join([subs[alpha.find(l)] for l in m])

def subs_cipher_decrypt(c):
c = "".join([l.upper() for l in c if not l in punct])
return "".join([alpha[subs.find(l)] for l in c])

m1 = "this is a test"

print "Original message:", m1
c1 = subs_cipher_encrypt(m1)

print
print "Encrypted Message:", c1

print
print "Decrypted Message:", subs_cipher_decrypt(c1)

Original message: this is a test

Encrypted Message: VAYRYREVLRV

Decrypted Message: THISISATEST

print "Original message:"
print m

print
c2 = subs_cipher_encrypt(m)
print "Encrypted Message:"
print c2

print
print "Decrypted Message:"
print subs_cipher_decrypt(c2)

Original message:
TobeornottobethatisthequestionWhethertisNoblerinthemindtosuffer
TheSlingsandArrowsofoutrageousFortuneOrtotakeArmsagainstaSeaof
troublesAndbyopposingendthem

Encrypted Message:
VHMLHOCHVVHMLVAEVYRVALXPLRVYHCUALVALOVYRCHMTLOYCVALWYCSVHRPZZLO
VALRTYCKRECSEOOHURHZHPVOEKLHPRZHOVPCLHOVHVEGLEOWREKEYCRVERLEHZ
VOHPMTLRECSMIHBBHRYCKLCSVALW

Decrypted Message:
TOBEORNOTTOBETHATISTHEQUESTIONWHETHERTISNOBLERINTHEMINDTOSUFFER
THESLINGSANDARROWSOFOUTRAGEOUSFORTUNEORTOTAKEARMSAGAINSTASEAOF
TROUBLESANDBYOPPOSINGENDTHEM

Is this a good encryption scheme? Still no. In fact, these are routinely used as puzzles in newspapers or puzzle books, because given a reasonably long message, it’s pretty easy to figure it out using things like the frequency of letters.

For instance, in English, the letters RSTLNEAO are very common, and much more common than other letters. So one might start to guess that the most common letter in the ciphertext is one of these. More powerfully, one might try to see which pairs of letters (called bigrams) are common, and look for those in the ciphertext, and so on.

From this sort of thought process, encryption schemes that ultimately rely on a single secret alphabet (even if it’s not our typical alphabet) fall pretty quickly. So… what about polyalphabetic ciphers? For instance, what if each group of 5 letters uses a different set of alphabets?

This is a great avenue for exploration, and there are lots of such encryption schemes that we won’t discuss in this class. But a class on cryptography (or a book on cryptography) would certainly go into some of these. It might also be a reasonable avenue of exploration for a final project.

## The German Enigma¶

One very well-known polyalphabetic encryption scheme is the German Enigma used before and during World War II. This was by far the most complicated cryptosystem in use up to that point, and the story of how it was broken is a long and tricky one. Intial successes towards breaking the Enigma came through the work of Polish mathematicians, fearful (and rightfully so) of the Germans across the border. By 1937, they had built replicas and understood many details of the Enigma system. But in 1938, the Germans shifted to a more secure and complicated cryptosystem. Just weeks before the German invasion of Poland and the beginning of World War II, Polish mathematicians sent their work and notes to mathematicians in France and Britain, who carried out this work.

The second major progress towards breaking the Enigma occurred largely in Bletchley Park in Britain, a communication center with an enormous dedicated effort to breaking the Enigma. This is where the tragic tale of Alan Turing, recently popularized through the movie The Imitation Game, begins. This is also the origin tale for modern computers, as Alan Turing developed electromechanical computers to help break the Enigma.

The Enigma worked by having a series of cogs or rotors whose positions determined a substitution cipher. After each letter, the positions were changed through a mechanical process. An Enigma machine is a very impressive machine to look at [and the “computer” Alan Turing used to help break them was also very impressive].

Below, I have implemented an Enigma, by default set to 4 rotors. I don’t expect one to understand the implementation. The interesting part is how meaningless the output message looks. Note that I’ve kept the spacing and punctuation from the original message for easier comparison. Really, you wouldn’t do this.

The plaintext used for demonstration is from Wikipedia’s article on the Enigma.

from random import shuffle,randint,choice
from copy import copy
num_alphabet = range(26)

def en_shift(l, n):                         # Rotate cogs and arrays
return l[n:] + l[:n]

class Cog:                                  # Each cog has a substitution cipher
def __init__(self):
self.shuf = copy(num_alphabet)
shuffle(self.shuf)                  # Create the individual substition cipher
return                              # Really, these were not random

def subs_in(self,i):                    # Perform a substition
return self.shuf[i]

def subs_out(self,i):                   # Perform a reverse substition
return self.shuf.index(i)

def rotate(self):                       # Rotate the cog by 1.
self.shuf = en_shift(self.shuf, 1)

def setcog(self,a):                     # Set up a particular substitution
self.shuf = a

class Enigma:
self.numcogs = numcogs
self.cogs = []
self.oCogs = []                     # "Original Cog positions"

for i in range(0,self.numcogs):     # Create the cogs
self.cogs.append(Cog())
self.oCogs.append(self.cogs[i].shuf)

refabet = copy(num_alphabet)
self.reflector = copy(num_alphabet)
while len(refabet) > 0:             # Pair letters in the reflector
a = choice(refabet)
refabet.remove(a)

b = choice(refabet)
refabet.remove(b)

self.reflector[a] = b
self.reflector[b] = a

def print_setup(self): # Print out substituion setup.
print "Enigma Setup:\nCogs: ",self.numcogs,"\nCog arrangement:"
for i in range(0,self.numcogs):
print self.cogs[i].shuf
print "Reflector arrangement:\n",self.reflector,"\n"

def reset(self):
for i in range(0,self.numcogs):
self.cogs[i].setcog(self.oCogs[i])

def encode(self,text):
t = 0     # Ticker counter
ciphertext=""
for l in text.lower():
num = ord(l) % 97
# Handle special characters for readability
if (num>25 or num<0):
ciphertext += l
else:
pass

else:
# Pass through cogs, reflect, then return through cogs
t += 1
for i in range(self.numcogs):
num = self.cogs[i].subs_in(num)

num = self.reflector[num]

for i in range(self.numcogs):
num = self.cogs[self.numcogs-i-1].subs_out(num)
ciphertext += "" + chr(97+num)

# Rotate cogs
for i in range(self.numcogs):
if ( t % ((i*6)+1) == 0 ):
self.cogs[i].rotate()
return ciphertext.upper()

plaintext="""When placed in an Enigma, each rotor can be set to one of 26 possible positions.
When inserted, it can be turned by hand using the grooved finger-wheel, which protrudes from
the internal Enigma cover when closed. So that the operator can know the rotor's position,
each had an alphabet tyre (or letter ring) attached to the outside of the rotor disk, with
26 characters (typically letters); one of these could be seen through the window, thus indicating
the rotational position of the rotor. In early models, the alphabet ring was fixed to the rotor
disk. A later improvement was the ability to adjust the alphabet ring relative to the rotor disk.
The position of the ring was known as the Ringstellung ("ring setting"), and was a part of the
initial setting prior to an operating session. In modern terms it was a part of the
initialization vector."""

# Remove newlines for encryption
pt = "".join([l.upper() for l in plaintext if not l == "\n"])
# pt = "".join([l.upper() for l in plaintext if not l in punct])

x=enigma(4)
#x.print_setup()

print "Original Message:"
print pt

print
ciphertext = x.encode(pt)
print "Encrypted Message"
print ciphertext

print
# Decryption and encryption are symmetric, so to decode we reset and re-encrypt.
x.reset()
decipheredtext = x.encode(ciphertext)
print "Decrypted Message:"
print decipheredtext

Original Message:
WHEN PLACED IN AN ENIGMA, EACH ROTOR CAN BE SET TO ONE OF 26 POSSIBLE POSITIONS.
WHEN INSERTED, IT CAN BE TURNED BY HAND USING THE GROOVED FINGER-WHEEL, WHICH
PROTRUDES FROM THE INTERNAL ENIGMA COVER WHEN CLOSED. SO THAT THE OPERATOR CAN
KNOW THE ROTOR'S POSITION, EACH HAD AN ALPHABET TYRE (OR LETTER RING) ATTACHED
TO THE OUTSIDE OF THE ROTOR DISK, WITH 26 CHARACTERS (TYPICALLY LETTERS); ONE
OF THESE COULD BE SEEN THROUGH THE WINDOW, THUS INDICATING THE ROTATIONAL
POSITION OF THE ROTOR. IN EARLY MODELS, THE ALPHABET RING WAS FIXED TO THE
ROTOR DISK. A LATER IMPROVEMENT WAS THE ABILITY TO ADJUST THE ALPHABET RING
RELATIVE TO THE ROTOR DISK. THE POSITION OF THE RING WAS KNOWN AS THE
RINGSTELLUNG ("RING SETTING"), AND WAS A PART OF THE INITIAL SETTING PRIOR TO
AN OPERATING SESSION. IN MODERN TERMS IT WAS A PART OF THE INITIALIZATION VECTOR.

Encrypted Message
RYQY LWMVIH OH EK GGWFBO, PDZN FNALL ZEP AP IUD JM FEV UG 26 HGRKODYT XWOLIYTSA.
TJQK KBCNVMHG, VS YBB XI BYZTVU XP BGWP IFNIY UQZ IQUPTKJ QXPVYE-EAOVT, PYMNN
RXFIOWYVK LWNN OJL JXZBTRVP YMQCXX STJQF KXNZ LXSWDC. LK KGRH ZKT RDZUMCWA GKY
LYKC LLV ZNFAD'Z BAXLDFTH, QSHA NBY RZ SXEXDAMP FXUS (WC RWHZOX ARWP) VKYNDJQP
WL OLW ICBALPI RV ISD GXSDQ ZKMJ, OHNN 26 IBGUHNYIQE (GDBNTEBWP QGECUNA); QPG
SQ JYUQX NDBWB NM AAUF RUNKYDL OLD LMPZQV, ZMKP SQZFDEHKCI KLJ AKPZLRAZZX
QXPJEXLV HS AFG IIMGK. NT PVAYP RFOKYV, XUY CHZAQEXG DUSS CKN YENSR FX PLS CYMZK
YHTB. J RLZLB DNIDWNWAIOO HOK PGM HVXXHIJ VV SCTNIC QGM TFKPQYPQ XFQU PSAECHTX
RA QUW CUNRV JCUW. LCP GFKZLLXW LN YWF OAQM CMF YVTQH EI JAU LZTXEYDGDFLQ
("CDBS TQHUEIR"), BLN QIG O UHJJ DV DQY KQGVFKI EBEFRCT WEZWG LJ WC NIAUKIYQJ
EHZLZLS. TY XPZSBL IPGWN ZS IUY E YQMD BW NVF QYWSDMTRSNSHYO GJGNUL.

Decrypted Message:
WHEN PLACED IN AN ENIGMA, EACH ROTOR CAN BE SET TO ONE OF 26 POSSIBLE POSITIONS.
WHEN INSERTED, IT CAN BE TURNED BY HAND USING THE GROOVED FINGER-WHEEL, WHICH
PROTRUDES FROM THE INTERNAL ENIGMA COVER WHEN CLOSED. SO THAT THE OPERATOR
CAN KNOW THE ROTOR'S POSITION, EACH HAD AN ALPHABET TYRE (OR LETTER RING)
ATTACHED TO THE OUTSIDE OF THE ROTOR DISK, WITH 26 CHARACTERS (TYPICALLY
LETTERS); ONE OF THESE COULD BE SEEN THROUGH THE WINDOW, THUS INDICATING THE
ROTATIONAL POSITION OF THE ROTOR. IN EARLY MODELS, THE ALPHABET RING WAS FIXED
TO THE ROTOR DISK. A LATER IMPROVEMENT WAS THE ABILITY TO ADJUST THE ALPHABET
RING RELATIVE TO THE ROTOR DISK. THE POSITION OF THE RING WAS KNOWN AS THE
RINGSTELLUNG ("RING SETTING"), AND WAS A PART OF THE INITIAL SETTING PRIOR TO
AN OPERATING SESSION. IN MODERN TERMS IT WAS A PART OF THE INITIALIZATION VECTOR.

The advent of computers brought in a paradigm shift in approaches towards cryptography. Prior to computers, one of the ways of maintaining security was to come up with a hidden key and a hidden cryptosystem, and keep it safe merely by not letting anyone know anything about how it actually worked at all. This has the short cute name security through obscurity. As the number of type of attacks on cryptosystems are much, much higher with computers, a different model of security and safety became necessary.

It is interesting to note that it is not obvious that security through obscurity is always bad, as long as it’s really really well-hidden. This is relevant to some problems concerning current events and cryptography.

# Public Key Cryptography¶

The new model begins with a slightly different setup. We should think of Anabel and Bartolo as sitting on opposite sides of a classroom, trying to communicate securely even though there are lots of people in the middle of the classroom who might be listening in. In particular, Anabel has something she wants to tell Bartolo.

Instead of keeping the cryptosystem secret, Bartolo tells everyone (in our metaphor, he shouts to the entire classroom) a public key K, and explains how to use it to send him a message. Anabel uses this key to encrypt her message. She then sends this message to Bartolo.

If the system is well-designed, no one will be able to understand the ciphertext even though they all know how the cryptosystem works. This is why the system is called Public Key.

Bartolo receives this message and (using something only he knows) he decrypts the message.

We will learn one such cryptosystem here: RSA, named after Rivest, Shamir, and Addleman — the first major public key cryptosystem.

## RSA¶

Bartolo takes two primes such as $p = 12553$ and $q = 13007$. He notes their product
$$m = pq = 163276871$$
and computes $\varphi(m)$,
$$\varphi(m) = (p-1)(q-1) = 163251312.$$
Finally, he chooses some integer $k$ relatively prime to $\varphi(m)$, like say
$$k = 79921.$$
Then the public key he distributes is
$$(m, k) = (163276871, 79921).$$

In order to send Bartolo a message using this key, Anabel must convert her message to numbers. She might use the identification A = 11, B = 12, C = 13, … and concatenate her numbers. To send the word CAB, for instance, she would send 131112. Let’s say that Anabel wants to send the message

NUMBER THEORY IS THE QUEEN OF THE SCIENCES

Then she needs to convert this to numbers.

conversion_dict = dict()
alpha = "abcdefghijklmnopqrstuvwxyz".upper()
curnum = 11
for l in alpha:
conversion_dict[l] = curnum
curnum += 1

print "Original Message:"
msg = "NUMBERTHEORYISTHEQUEENOFTHESCIENCES"
print msg
print

def letters_to_numbers(m):
return "".join([str(conversion_dict[l]) for l in m.upper()])

print "Numerical Message:"
msg_num = letters_to_numbers(msg)
print msg_num

Original Message:
NUMBERTHEORYISTHEQUEENOFTHESCIENCES

Numerical Message:
2431231215283018152528351929301815273115152425163018152913191524131529

So Anabel’s message is the number

$$2431231215283018152528351929301815273115152425163018152913191524131529$$

which she wants to encrypt and send to Bartolo. To make this manageable, she cuts the message into 8-digit numbers,

$$24312312, 15283018, 15252835, 19293018, 15273115, 15242516, 30181529, 13191524, 131529.$$

To send her message, she takes one of the 8-digit blocks and raises it to the power of $k$ modulo $m$. That is, to transmit the first block, she computes

$$24312312^{79921} \equiv 13851252 \pmod{163276871}.$$

# Secret information
p = 12553
q = 13007
phi = (p-1)*(q-1) # varphi(pq)

# Public information
m = p*q # 163276871
k = 79921

print pow(24312312, k, m)

13851252

She sends this number
$$13851252$$
to Bartolo (maybe by shouting. Even though everyone can hear, they cannot decrypt it). How does Bartolo decrypt this message?

He computes $\varphi(m) = (p-1)(q-1)$ (which he can do because he knows $p$ and $q$ separately), and then finds a solution $u$ to
$$uk = 1 + \varphi(m) v.$$
This can be done quickly through the Euclidean Algorithm.

def extended_euclidean(a,b):
if b == 0:
return (1,0,a)
else :
x, y, gcd = extended_euclidean(b, a % b) # Aside: Python has no tail recursion
return y, x - y * (a // b),gcd           # But it does have meaningful stack traces

# This version comes from Exercise 6.3 in the book, but without recursion
def extended_euclidean2(a,b):
x = 1
g = a
v = 0
w = b
while w != 0:
q = g // w
t = g - q*w
s = x - q*v
x,g = v,w
v,w = s,t
y = (g - a*x) / b
return (x,y,g)

def modular_inverse(a,m) :
x,y,gcd = extended_euclidean(a,m)
if gcd == 1 :
return x % m
else :
return None

print "k, p, q:", k, p, q
print
u = modular_inverse(k,(p-1)*(q-1))
print u

k, p, q: 79921 12553 13007

145604785

In particular, Bartolo computes that his $u = 145604785$. To recover the message, he takes the number $13851252$ sent to him by Anabel and raises it to the $u$ power. He computes
$$13851252^{u} \equiv 24312312 \pmod{pq},$$
which we can see must be true as
$$13851252^{u} \equiv (24312312^{k})^u \equiv 24312312^{1 + \varphi(pq)v} \equiv 24312312 \pmod{pq}.$$
In this last step, we have used Euler’s Theorem to see that
$$24312312^{\varphi(pq)v} \equiv 1 \pmod{pq}.$$

# Checking this power explicitly.
print pow(13851252, 145604785, m)

24312312

Now Bartolo needs to perform this process for each 8-digit chunk that Anabel sent over. Note that the work is very easy, as he computes the integer $u$ only once. Each other time, he simply computes $c^u \pmod m$ for each ciphertext $c$, and this is very fast with repeated-squaring.

We do this below, in an automated fashion, step by step.

First, we split the message into 8-digit chunks.

# Break into chunks of 8 digits in length.
def chunk8(message_number):
cp = str(message_number)
ret_list = []
while len(cp) > 7:
ret_list.append(cp[:8])
cp = cp[8:]
if cp:
ret_list.append(cp)
return ret_list

msg_list = chunk8(msg_num)
print msg_list

['24312312', '15283018', '15252835', '19293018', '15273115',
'15242516', '30181529', '13191524', '131529']

This is a numeric representation of the message Anabel wants to send Bartolo. So she encrypts each chunk. This is done below

# Compute ciphertexts separately on each 8-digit chunk.
def encrypt_chunks(chunked_list):
ret_list = []
for chunk in chunked_list:
#print chunk
#print int(chunk)
ret_list.append(pow(int(chunk), k, m))
return ret_list

cipher_list = encrypt_chunks(msg_list)
print cipher_list

[13851252, 14944940, 158577269, 117640431, 139757098, 25099917,
88562046, 6640362, 10543199]

This is the encrypted message. Having computed this, Anabel sends this message to Bartolo.

To decrypt the message, Bartolo uses his knowledge of $u$, which comes from his ability to compute $\varphi(pq)$, and decrypts each part of the message. This looks like below.

# Decipher the ciphertexts all in the same way
def decrypt_chunks(chunked_list):
ret_list = []
for chunk in chunked_list:
ret_list.append(pow(int(chunk), u, m))
return ret_list

decipher_list = decrypt_chunks(cipher_list)
print decipher_list

[24312312, 15283018, 15252835, 19293018, 15273115, 15242516,
30181529, 13191524, 131529]

Finally, Bartolo concatenates these numbers together and translates them back into letters. Will he get the right message back?

alpha = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

# Collect deciphered texts into a single list, and translate back into letters.
def chunks_to_letters(chunked_list):
s = "".join([str(chunk) for chunk in chunked_list])
ret_str = ""
while s:
ret_str += alpha[int(s[:2])-11].upper()
s = s[2:]
return ret_str

print chunks_to_letters(decipher_list)

NUMBERTHEORYISTHEQUEENOFTHESCIENCES

Yes! Bartolo successfully decrypts the message and sees that Anabel thinks that Number Theory is the Queen of the Sciences. This is a quote from Gauss, the famous mathematician who has been popping up again and again in this course.

## Why is this Secure?¶

Let’s pause to think about why this is secure.

What if someone catches the message in transit? Suppose Eve is eavesdropping and hears Anabel’s first chunk, $13851252$. How would she decrypt it?

Eve knows that she wants to solve
$$x^k \equiv 13851252 \pmod {pq}$$
or rather
$$x^{79921} \equiv 13851252 \pmod {163276871}.$$
How can she do that? We can do that because we know to raise this to a particular power depending on $\varphi(163276871)$. But Eve doesn’t know what $\varphi(163276871)$ is since she can’t factor $163276871$. In fact, knowing $\varphi(163276871)$ is exactly as hard as factoring $163276871$.

But if Eve could somehow find $79921$st roots modulo $163276871$, or if Eve could factor $163276871$, then she would be able to decrypt the message. These are both really hard problems! And it’s these problems that give the method its security.

More generally, one might use primes $p$ and $q$ that are each about $200$ digits long, and a fairly large $k$. Then the resulting $m$ would be about $400$ digits, which is far larger than we know how to factor effectively. The reason for choosing a somewhat large $k$ is for security reasons that are beyond the scope of this segment. The relevant idea here is that since this is a publicly known encryption scheme, many people have strenghtened it over the years by making it more secure against every clever attack known. This is another, sometimes overlooked benefit of public-key cryptography: since the methodology isn’t secret, everyone can contribute to its security — indeed, it is in the interest of anyone desiring such security to contribute. This is sort of the exact opposite of Security through Obscurity.

The nature of code being open, public, private, or secret is also very topical in current events. Recently, Volkswagen cheated in its car emissions-software and had it report false outputs, leading to a large scandal. Their software is proprietary and secret, and the deliberate bug went unnoticed for years. Some argue that this means that more software, and especially software that either serves the public or that is under strong jurisdiction, should be publicly viewable for analysis.

Another relevant current case is that the code for most voting machines in the United States is proprietary and secret. Hopefully they aren’t cheating!

On the other side, many say that it is necessary for companies to be able to have secret software for at least some time period to recuperate the expenses that go into developing the software. This is similar to how drug companies have patents on new drugs for a number of years. This way, a new successful drug pays for its development since the company can charge above the otherwise-market-rate.

Further, many say that opening some code would open it up to attacks from malicious users who otherwise wouldn’t be able to see the code. Of course, this sounds a lot like trying for security through obscurity.

This is a very relevant and big topic, and the shape it takes over the next few years may very well have long-lasting impacts.

## Condensed Version¶

Now that we’ve gone through it all once, we have a condensed RSA system set up with our $p, q,$ and $k$ from above. To show that this can be done quickly with a computer, let’s do another right now.

Let us encrypt, transmit, and decrypt the message

“I have never done anything useful. No discovery of mine has made, or is likely to make, directly or indirectly, for good or ill, the least difference to the amenity of the world”.

First, we prepare the message for conversion to numbers.

message = """I have never done anything useful. No discovery of mine has made,
or is likely to make, directly or indirectly, for good or ill, the least
difference to the amenity of the world"""

message = "".join([l.upper() for l in message if not l in "\n .,"])
print "Our message:\n"+message

Our message:
TOMAKEDIRECTLYORINDIRECTLYFORGOODORILLTHELEASTDIFFERENCETOTHE
AMENITYOFTHEWORLD

Now we convert the message to numbers, and transform those numbers into 8-digit chunks.

numerical_message = letters_to_numbers(message)
print "Our message, converted to numbers:"
print numerical_message
print

plaintext_chunks = chunk8(numerical_message)
print "Separated into 8-digit chunks:"
print plaintext_chunks

Our message, converted to numbers:
1918113215241532152814252415112435301819241731291516312224251419
2913253215283525162319241518112923111415252819292219211522353025
2311211514192815133022352528192414192815133022351625281725251425
2819222230181522151129301419161615281524131530253018151123152419
303525163018153325282214

Separated into 8-digit chunks:
['19181132', '15241532', '15281425', '24151124', '35301819',
'24173129', '15163122', '24251419', '29132532', '15283525',
'16231924', '15181129', '23111415', '25281929', '22192115',
'22353025', '23112115', '14192815', '13302235', '25281924',
'14192815', '13302235', '16252817', '25251425', '28192222',
'30181522', '15112930', '14191616', '15281524', '13153025',
'30181511', '23152419', '30352516', '30181533', '25282214']

We encrypt each chunk by computing $P^k \bmod {m}$ for each plaintext chunk $P$.

ciphertext_chunks = encrypt_chunks(plaintext_chunks)
print ciphertext_chunks

[99080958, 142898385, 80369161, 11935375, 108220081, 82708130,
158605094, 96274325, 154177847, 121856444, 91409978, 47916550,
155466420, 92033719, 95710042, 86490776, 15468891, 139085799,
68027514, 53153945, 139085799, 68027514, 9216376, 155619290,
83776861, 132272900, 57738842, 119368739, 88984801, 83144549,
136916742, 13608445, 92485089, 89508242, 25375188]

This is the message that Anabel can sent Bartolo.

To decrypt it, he computes $c^u \bmod m$ for each ciphertext chunk $c$.

deciphered_chunks = decrypt_chunks(ciphertext_chunks)
print "Deciphered chunks:"
print deciphered_chunks

Deciphered chunks:
[19181132, 15241532, 15281425, 24151124, 35301819, 24173129,
15163122, 24251419, 29132532, 15283525, 16231924, 15181129,
23111415, 25281929, 22192115, 22353025, 23112115, 14192815,
13302235, 25281924, 14192815, 13302235, 16252817, 25251425,
28192222, 30181522, 15112930, 14191616, 15281524, 13153025,
30181511, 23152419, 30352516, 30181533, 25282214]

Finally, he translates the chunks back into letters.

decoded_message = chunks_to_letters(deciphered_chunks)
print "Decoded Message:"
print decoded_message

Decoded Message:
TOMAKEDIRECTLYORINDIRECTLYFORGOODORILLTHELEASTDIFFERENCETOTHE
AMENITYOFTHEWORLD

Even with large numbers, RSA is pretty fast. But one of the key things that one can do with RSA is securely transmit secret keys for other types of faster-encryption that don’t work in a public-key sense. There is a lot of material in this subject, and it’s very important.

The study of sending and receiving secret messages is called Cryptography. There are lots of interesting related topics, some of which we’ll touch on in class.

The study of analyzing and breaking cryptosystems is called Cryptanalysis, and is something I find quite fun. But it’s also quite intense.

I should mention that in practice, RSA is performed a little bit differently to make it more secure.

Aside: this is the first time I’ve converted a Jupyter Notebook, and it turns out it’s very easy. However, there are a few formatting annoyances that I didn’t consider, but which are too minor to spend time fixing. But now I know!
| Tagged , , , , | 4 Comments

format: plain


text:
# (Plain)TeX

This is not a simple (Plain)TeX Reference page or a quick guide to TeX or friends: this is my own reference sheet in wich you can find what I don't find in internet or The TeXBook.

TeX (not LaTeX) is here not the very sistem, but the code of the text to compile with (Plain)TeX (not LaTeX): that is, how do the difficult things simple in (Plain)TeX.

I quote from http://www.topology.org/soft/tex.html (20-IX-07):

"TeX is the venerable and still revolutionary typesetting software from Donald Knuth.
My documents written in TeX in 1985 still work exactly the same now as then.
And TeX is still by far the best typesetting system.
No versionitis. No bugs. And it's free and open.
TeX is what started the open source revolution, not Richard Stallman's FSF."

You can freely copy and paste or download the examplain.tex.

## 1. Definition of \alph to do enumerations with letters always in (Plain)TeX:

\def\alph #1{ %
\ifcase #1 %
\or a\or b\or c\or d\or e\or f\or g\or h\or i %
\or j\or k\or l\or m\or n\or o\or p\or q\or r %
\or s\or t\or u\or v\or w\or x\or y\or z %
\fi %
}

## 2. To do automatic numbering of footnotes define this:

\newcount\nnnota
\nnnota= 0
\def\note{\global\advance\nnnota
by 1\footnote{$^{\the\nnnota}$}}

And then write so:

This is the first note\note{first}, and this is the second\note{second}, and this is the third\note{third}.

## 3. Now do the same but with letters:

\newcount\nnota
\nnota= 0
\def\note{\global\advance\nnota
by 1\footnote{$^{\alph{\the\nnota}}$}}

pay attention to the different number of "n" in the definitions (\nnnota and \nnota).

## 4. Do you think that the two series could coexist?

Yes. Here you have:

\newcount\nnota \nnota= 0
\def\anote{\global\advance\nnota
by 1\footnote{$^{\alph{\the\nnota}}$}}


text:
# Scalar multiplication Rule of matrices

## Formula

$k \times \big[\,e_{ij}\,\big] \,=\, \big[k \times e_{ij}\big]$

The multiplication of a matrix by a scalar is equal to the multiplication of each entry (or element) by the scalar in the matrix, is called the scalar multiplication rule of the matrices.

## Introduction

Let $k$ be a constant and the matrix $M$ represents a matrix of the order $m \times n$.

$M$ $\,=\,$ $\begin{bmatrix} e_{11} & e_{12} & e_{13} & \cdots & e_{1n}\\ e_{21} & e_{22} & e_{23} & \cdots & e_{2n}\\ e_{31} & e_{32} & e_{33} & \cdots & e_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ e_{m1} & e_{m2} & e_{m3} & \cdots & e_{mn} \end{bmatrix}$

Let’s assume that the matrix $M$ is multiplied by the scalar $k$.

$\implies$ $k \times M$ $\,=\,$ $k \times \begin{bmatrix} e_{11} & e_{12} & e_{13} & \cdots & e_{1n}\\ e_{21} & e_{22} & e_{23} & \cdots & e_{2n}\\ e_{31} & e_{32} & e_{33} & \cdots & e_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ e_{m1} & e_{m2} & e_{m3} & \cdots & e_{mn} \end{bmatrix}$

The product of a scalar and a matrix is equal to the product of each element and scalar in a matrix.

$\implies$ $k \times M$ $\,=\,$ $\begin{bmatrix} k \times e_{11} & k \times e_{12} & k \times e_{13} & \cdots & k \times e_{1n}\\ k \times e_{21} & k \times e_{22} & k \times e_{23} & \cdots & k \times e_{2n}\\ k \times e_{31} & k \times e_{32} & k \times e_{33} & \cdots & k \times e_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ k \times e_{m1} & k \times e_{m2} & k \times e_{m3} & \cdots & k \times e_{mn} \end{bmatrix}$

This mathematical property of matrices is called the scalar multiplication of the matrices. It is used as a formula to multiply any matrix by any scalar.

If the matrix $M$ is simply expressed as $\big[e_{ij}\big]$, then the scalar multiplication rule can be simply written as follows.

$\implies$ $k \big[e_{ij}\big] \,=\, \big[k e_{ij}\big]$

Here, the letters $i$ and $j$ represent “the number of the row” and “the number of the column”.

### Examples

Let’s learn how to use the scalar multiplication rule of the matrices from some understandable examples.

$(1).\,\,\,$ $4$ $\times$ $\begin{bmatrix} 1 & 7 \\ -2 & 6 \end{bmatrix}$

In this example, the matrix of the order $2$ is multiplied by a scalar $4$. It can be evaluated by multiplying each entry in the matrix by the scalar $4$.

$\implies$ $4$ $\times$ $\begin{bmatrix} 1 & 7 \\ -2 & 6 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} 4 \times 1 & 4 \times 7 \\ 4 \times (-2) & 4 \times 6 \end{bmatrix}$

$\,\,\,\therefore\,\,\,\,\,\,$ $4$ $\times$ $\begin{bmatrix} 1 & 7 \\ -2 & 6 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} 4 & 28 \\ -8 & 24 \end{bmatrix}$

$(2).\,\,\,$ $-1$ $\times$ $\begin{bmatrix} 2 & 0 & 3 & 5 \\ 5 & 1 & 8 & 4 \\ 7 & 3 & 2 & 6 \end{bmatrix}$

In this case, a matrix of the order $3 \times 4$ is multiplied by the scalar $-1$. So, multiply each entry in this matrix by the number -1.

$\implies$ $-1$ $\times$ $\begin{bmatrix} 2 & 0 & 3 & 5 \\ 5 & 1 & 8 & 4 \\ 7 & 3 & 2 & 6 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} (-1) \times 2 & (-1) \times 0 & (-1) \times 3 & (-1) \times 5 \\ (-1) \times 5 & (-1) \times 1 & (-1) \times 8 & (-1) \times 4 \\ (-1) \times 7 & (-1) \times 3 & (-1) \times 2 & (-1) \times 6 \end{bmatrix}$

$\,\,\,\therefore\,\,\,\,\,\,$ $-1$ $\times$ $\begin{bmatrix} 2 & 0 & 3 & 5 \\ 5 & 1 & 8 & 4 \\ 7 & 3 & 2 & 6 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} -2 & 0 & -3 & -5 \\ -5 & -1 & -8 & -4 \\ -7 & -3 & -2 & -6 \end{bmatrix}$

$(3).\,\,\,$ $6$ $\times$ $\begin{bmatrix} -9 & 5 & -2 \\ 5 & 8 & -3 \\ 6 & -1 & 0 \end{bmatrix}$

Use the scalar multiplication of matrices formula for multiplying the square matrix of the order $3$ by the scalar $6$.

$\implies$ $6$ $\times$ $\begin{bmatrix} -9 & 5 & -2 \\ 5 & 8 & -3 \\ 6 & -1 & 0 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} 6 \times (-9) & 6 \times 5 & 6 \times (-2) \\ 6 \times 5 & 6 \times 8 & 6 \times (-3) \\ 6 \times 6 & 6 \times (-1) & 6 \times 0 \end{bmatrix}$

$\,\,\,\therefore\,\,\,\,\,\,$ $6$ $\times$ $\begin{bmatrix} -9 & 5 & -2 \\ 5 & 8 & -3 \\ 6 & -1 & 0 \end{bmatrix}$ $\,=\,$ $\begin{bmatrix} -54 & 30 & -12 \\ 30 & 48 & -18 \\ 36 & -6 & 0 \end{bmatrix}$

Latest Math Topics
Email subscription
Math Doubts is a best place to learn mathematics and from basics to advanced scientific level for students, teachers and researchers. Know more


text:
Does Mathematica implement Risch algorithm? If it does, in which cases? [duplicate]

I am asking this questions because when trying to evaluate the integrals:

\begin{eqnarray} (1) \qquad & \int \frac{x}{\sqrt{x^4+10x^2-96x-71}}\ \mathrm{d}x \\ (2) \qquad & \int \frac{x^2+2x+1+(3x+1)\sqrt{x+\ln(x)}}{x\left(x+\sqrt{x+\ln(x)}\right)\sqrt{x+\ln(x)}}\ \mathrm{d}x \end{eqnarray}

In Mathematica 11 using

Integrate[x/Sqrt[x^4 + 10 x^2 - 96 x - 71], x]

for (1) I obtain an evaluation involving non-elementary functions, while for (2) the expression

Integrate[(x^2 + 2 x + 1 + (3 x + 1) Sqrt[x + Log[x]])/(x (x + Sqrt[x + Log[x]]) Sqrt[x + Log[x]]), x]

is not even evaluated.

In Maple 2016 I get the same results as in Mathematica. In Rubi 4.11 neither are evaluated. The two indefinite integrals can be expressed with elementary functions.

ADDENDUM: The antiderivative of the integrand (1) is correctly expressed with elementary functions using Axiom. Why can't the two powerful Computer Algebra Systems listed above, obtain the same result as Axiom?

• Provide the formatted code of the integration. – Edmund Mar 15 '17 at 1:22
• @Nasser He did not go into detail with his answer, I am looking for something more elaborate. – user372003 Mar 15 '17 at 1:52
• I remember I read from somewhere "no CAS has completely implemented Risch Algorithm" – happy fish Mar 15 '17 at 2:06
• Although the documentation stated "For indefinite integrals, an extended version of the Risch algorithm is used whenever both the integrand and integral can be expressed in terms of elementary functions" from long time ago, I have found some elementary examples Mathematica cannot solve. That may imply Mathematica did not fully implement the algorithm. – happy fish Mar 15 '17 at 3:19
• (2) is a duplicate of (109690), (127488) – Michael E2 Mar 15 '17 at 12:02

For (1) DSolve returns an answer quickly, so this "powerful Computer Algebra System" can do it:

DSolve[y'[x] == x/Sqrt[x^4 + 10 x^2 - 96 x - 71], y, x]
(*  {{y -> Function[{x}, x - Log[x/Sqrt[-71 - 96 x + 10 x^2 + x^4]]]}}  *)

(2) is a duplicate of Why can't Mathematica integrate this? as well as How to do this integral?.

All that seems to be available to the user can be found here:

http://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html#23196

For more details, one might need help from Wolfram Support.

• The OP wasn't saying MMA's Integrate[ ] doesn't return an answer for (1), he was saying it doesn't return an answer in terms of elementary functions (by which I believe he means one without Root objects). As s/he notes, MMA does return answer with Integrate[ ]. I checked and, except for the lack of an explicit constant, Integrate's solution is identical to the one given by DSolve[ ]. – theorist Apr 14 '18 at 0:23

This is not a answer ,but you can simplify yours integral.

For (2):

func = (1 + 2 x + x^2 + (1 + 3 x) Sqrt[x + Log[x]])/(x Sqrt[x + Log[x]] (x + Sqrt[x + Log[x]]));

func1 = Apart[func] // ExpandAll

Plus @@ Table[Integrate[func1[[n]], x], {n, 1, Length[func1]}]

$$\int \frac{1}{x^2 \log (x)-\log ^2(x)-x \log (x)} \, dx-\int \frac{x}{x^2 \log (x)-\log ^2(x)-x \log (x)} \, dx-\int \frac{\log (x)}{x^2 \log (x)-\log ^2(x)-x \log (x)} \, dx+2 \int \frac{x \log (x)}{x^2 \log (x)-\log ^2(x)-x \log (x)} \, dx+\int \frac{x \sqrt{x+\log (x)}}{x^2 \log (x)-\log ^2(x)-x \log (x)} \, dx+\int \frac{\sqrt{x+\log (x)}}{x^2 \log ^2(x)-\log ^3(x)-x \log ^2(x)} \, dx-\int \frac{x \sqrt{x+\log (x)}}{x^2 \log ^2(x)-\log ^3(x)-x \log ^2(x)} \, dx-2 \int \frac{\sqrt{x+\log (x)}}{x^2-x-\log (x)} \, dx+\log (\log (x))+\int \frac{1}{\log ^2(x) \sqrt{x+\log (x)}} \, dx+\int \frac{1}{\sqrt{x+\log (x)}} \, dx-2 \int \frac{1}{\log (x) \sqrt{x+\log (x)}} \, dx+\int \frac{\sqrt{x+\log (x)}}{x \log (x)} \, dx$$

• Is the simplification just the Log[Log[x]] term? It seems to me the integrand has been broken apart into terms which cannot all be integrated in finite terms, unlike the original problem. Or am I missing something? – Michael E2 Mar 16 '17 at 12:46