Category Archives: Python

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:
    lines = f.readlines()
    
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:
    lines = f.readlines()
    
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 , , | Leave a 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?

(more…)

Posted in Expository, Programming, Python | Tagged , , | Leave a comment

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 = f.readlines()
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.

(more…)

Posted in Expository, Programming, Python | Tagged , , | Leave a comment

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 = f.read()
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)
    return total

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.

(more…)

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

Sage Days 87 Demo: Interfacing between sage and the LMFDB

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',
'_Hom_', '__add__', '__cached_methods', '__call__',
'__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__',
'__weakref__', '_abstract_element_class', '_adjust_heegner_index', '_an_element_',
'_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',
'_normalize_padic_lseries', '_octave_', '_octave_init_', '_p_primary_torsion_basis',
'_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',
'global_minimal_model', 'global_minimality_class', 'has_additive_reduction', 'has_bad_reduction',
'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_quadratic_twist', 'is_quartic_twist', 'is_semistable', 'is_sextic_twist',
'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',
'minimal_discriminant_ideal', 'minimal_model', 'minimal_quadratic_twist', 'mod5family',
'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',
'ordinary_model', 'ordinary_primes', 'padic_E2', 'padic_height',
'padic_height_pairing_matrix', 'padic_height_via_multiply', 'padic_lseries', 'padic_regulator',
'padic_sigma', 'padic_sigma_truncated', 'parent', 'pari_curve',
'pari_mincurve', 'period_lattice', 'plane_projection', 'plot',
'point', 'point_homset', 'point_search', 'point_set',
'pollack_stevens_modular_symbol', 'preimage', 'projection', 'prove_BSD',
'q_eigenform', 'q_expansion', 'quadratic_transform', 'quadratic_twist',
'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']
All the things
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 ``b``th
       curve. Combined with ``max_items`` to return data in chunks.

    -  ``sort=s`` -- a string ``s`` specifying what database field to sort the
       results on. See the LMFDB api for more info.

    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.
XKCD's automation

Posted in Expository, LMFDB, Math.NT, Mathematics, Programming, Python, sagemath | Tagged , , , , , , , | Leave a comment

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:  
    def __init__(self, numcogs,readability=True):  
        self.readability = readability  
        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):  
                if (self.readability):
                    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:
IHAVENEVERDONEANYTHINGUSEFULNODISCOVERYOFMINEHASMADEORISLIKELY
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:
IHAVENEVERDONEANYTHINGUSEFULNODISCOVERYOFMINEHASMADEORISLIKELY
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!
Posted in Brown University, Expository, Math 420, Math.NT, Mathematics, Programming, Python, Teaching | Tagged , , , , | 4 Comments

Comments vs documentation in python

In my first programming class we learned python. It went fine (I thought), I got the idea, and I moved on (although I do now see that much of what we did was not ‘pythonic’).

But now that I’m returning to programming (mostly in python), I see that I did much of it all wrong. One of my biggest surprises was how wrong I was about comments. Too many comments are terrible. Redundant comments make code harder to maintain. If the code is too complex to understand without comments, it’s probably just bad code.

That’s not so hard, right? You read some others’ code, see their comment conventions, and move on. You sort of get into zen moments where the code becomes very clear and commentless, and there is bliss. But then we were at a puzzling competition, and someone wanted to use a piece of code I’d written. Sure, no problem! And the source was so beautifully clear that it almost didn’t even need a single comment to understand.

But they didn’t want to read the code. They wanted to use the code. Comments are very different than documentation. The realization struck me, and again I had it all wrong. In hindsight, it seems so obvious! I’ve programmed java, I know about javadoc. But no one had ever actually wanted to use my code before (and wouldn’t have been able to easily if they had)!

Enter pydoc and sphinx. These are tools that allow html API to be generated from comment docstrings in the code itself. There is a cost – some comments with specific formatting are below each method or class. But it’s reasonable, I think.

Pydoc ships with python, and is fine for single modules or small projects. But for large projects, you’ll need something more. In fact, even the documentation linked above for pydoc was generated with sphinx:

The pydoc documentation is  ‘Created using Sphinx 1.0.7.’

This isn’t to say that pydoc is bad. But I didn’t want to limit myself. Python uses sphinx, so I’ll give it a try too.

And thus I (slightly excessively, to get the hang of it) comment on my solutions to Project Euler on my github. The current documentation looks alright, and will hopefully look better as I get the hang of it.

Full disclosure – this was originally going to be a post on setting up sphinx-generated documentation on github’s pages automatically. I got sidetracked – that will be the next post.

Posted in Programming, Python | Tagged , , , , , | Leave a comment