# Category Archives: Expository

## 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

## A Jupyter Notebook from a SageMath tutorial

I gave an introduction to sage tutorial at the University of Warwick Computational Group seminar today, 2 November 2017. Below is a conversion of the sage/jupyter notebook I based the rest of the tutorial on. I said many things which are not included in the notebook, and during the seminar we added a few additional examples and took extra consideration to a few different calls. But for reference, the notebook is here.

The notebook itself (as a jupyter notebook) can be found and viewed on my github (link to jupyter notebook). When written, this notebook used a Sage 8.0.0.rc1 backend kernel and ran fine on the standard Sage 8.0 release , though I expect it to work fine with any recent official version of sage. The last cell requires an active notebook to be seen (or some way to export jupyter widgets to standalone javascript or something; this either doesn’t yet exist, or I am not aware of it).

I will also note that I converted the notebook for display on this website using jupyter’s nbconvert package. I have some CSS and syntax coloring set up that affects the display.

Good luck learning sage, and happy hacking.

# Sage¶

Sage (also known as SageMath) is a general purpose computer algebra system written on top of the python language. In Mathematica, Magma, and Maple, one writes code in the mathematica-language, the magma-language, or the maple-language. Sage is python.

But no python background is necessary for the rest of today’s guided tutorial. The purpose of today’s tutorial is to give an indication about how one really uses sage, and what might be available to you if you want to try it out.

I will spoil the surprise by telling you upfront the two main points I hope you’ll take away from this tutorial.

1. With tab-completion and documentation, you can do many things in sage without ever having done them before.
2. The ecosystem of libraries and functionality available in sage is tremendous, and (usually) pretty easy to use.

## Lightning Preview¶

Let’s first get a small feel for sage by seeing some standard operations and what typical use looks like through a series of trivial, mostly unconnected examples.

In [1]:
# Fundamental manipulations work as you hope

2+3

Out[1]:
5

You can also subtract, multiply, divide, exponentiate…

>>> 3-2
1
>>> 2*3
6
>>> 2^3
8
>>> 2**3 # (also exponentiation)
8



There is an order of operations, but these things work pretty much as you want them to work. You might try out several different operations.

Sage includes a lot of functionality, too. For instance,

In [2]:
factor(-1008)

Out[2]:
-1 * 2^4 * 3^2 * 7
In [3]:
list(factor(1008))

Out[3]:
[(2, 4), (3, 2), (7, 1)]

In the background, Sage is actually calling on pari/GP to do this factorization. Sage bundles lots of free and open source math software within it (which is why it’s so large), and provides a common access point. The great thing here is that you can often use sage without needing to know much pari/GP (or other software).

Sage knows many functions and constants, and these are accessible.

Posted in Expository, Mathematics, sage, sagemath | Tagged , , , , , | Leave a 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.

## Smooth Sums to Sharp Sums 1

In this note, I describe a combination of two smoothed integral transforms that has been very useful in my collaborations with Alex Walker, Chan Ieong Kuan, and Tom Hulse. I suspect that this particular technique was once very well-known. But we were not familiar with it, and so I describe it here.

In application, this is somewhat more complicated. But to show the technique, I apply it to reprove some classic bounds on $\text{GL}(2)$ $L$-functions.

This note is also available as a pdf. This was first written as a LaTeX document, and then modified to fit into wordpress through latex2jax.

## Introduction

Consider a Dirichlet series
$$D(s) = \sum_{n \geq 1} \frac{a(n)}{n^s}. \notag$$
Suppose that this Dirichlet series converges absolutely for $\Re s > 1$, has meromorphic continuation to the complex plane, and satisfies a functional equation of shape
$$\Lambda(s) := G(s) D(s) = \epsilon \Lambda(1-s), \notag$$
where $\lvert \epsilon \rvert = 1$ and $G(s)$ is a product of Gamma factors.

Dirichlet series are often used as a tool to study number theoretic functions with multiplicative properties. By studying the analytic properties of the Dirichlet series, one hopes to extract information about the coefficients $a(n)$. Some of the most common interesting information within Dirichlet series comes from partial sums
$$S(n) = \sum_{m \leq n} a(m). \notag$$
For example, the Gauss Circle and Dirichlet Divisor problems can both be stated as problems concerning sums of coefficients of Dirichlet series.

One can try to understand the partial sum directly by understanding the integral transform
$$S(n) = \frac{1}{2\pi i} \int_{(2)} D(s) \frac{X^s}{s} ds, \notag$$
a Perron integral. However, it is often challenging to understand this integral, as delicate properties concerning the convergence of the integral often come into play.

Instead, one often tries to understand a smoothed sum of the form
$$\sum_{m \geq 1} a(m) v(m) \notag$$
where $v(m)$ is a smooth function that vanishes or decays extremely quickly for values of $m$ larger than $n$. A large class of smoothed sums can be obtained by starting with a very nicely behaved weight function $v(m)$ and take its Mellin transform
$$V(s) = \int_0^\infty v(x) x^s \frac{dx}{x}. \notag$$
Then Mellin inversion gives that
$$\sum_{m \geq 1} a(m) v(m/X) = \frac{1}{2\pi i} \int_{(2)} D(s) X^s V(s) ds, \notag$$
as long as $v$ and $V$ are nice enough functions.

In this note, we will use two smoothing integral transforms and corresponding smoothed sums. We will use one smooth function $v_1$ (which depends on another parameter $Y$) with the property that
$$\sum_{m \geq 1} a(m) v_1(m/X) \approx \sum_{\lvert m – X \rvert < X/Y} a(m). \notag$$
And we will use another smooth function $v_2$ (which also depends on $Y$) with the property that
$$\sum_{m \geq 1} a(m) v_2(m/X) = \sum_{m \leq X} a(m) + \sum_{X < m < X + X/Y} a(m) v_2(m/X). \notag$$
Further, as long as the coefficients $a(m)$ are nonnegative, it will be true that
$$\sum_{X < m < X + X/Y} a(m) v_2(m/X) \ll \sum_{\lvert m – X \rvert < X/Y} a(m), \notag$$
which is exactly what $\sum a(m) v_1(m/X)$ estimates. Therefore
$$\label{eq:overall_plan} \sum_{m \leq X} a(m) = \sum_{m \geq 1} a(m) v_2(m/X) + O\Big(\sum_{m \geq 1} a(m) v_1(m/X) \Big).$$

Hence sufficient understanding of $\sum a(m) v_1(m/X)$ and $\sum a(m) v_2(m/X)$ allows one to understand the sharp sum
$$\sum_{m \leq X} a(m). \notag$$

Posted in Expository, Math.NT, Mathematics | | 3 Comments

## Computing $\pi$

This note was originally written in the context of my fall Math 100 class at Brown University. It is also available as a pdf note.

While investigating Taylor series, we proved that
\label{eq:base}
\frac{\pi}{4} = 1 – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} + \frac{1}{9} + \cdots

Let’s remind ourselves how. Begin with the geometric series

\frac{1}{1 + x^2} = 1 – x^2 + x^4 – x^6 + x^8 + \cdots = \sum_{n = 0}^\infty (-1)^n x^{2n}. \notag

(We showed that this has interval of convergence $\lvert x \rvert < 1$). Integrating this geometric series yields

\int_0^x \frac{1}{1 + t^2} dt = x – \frac{x^3}{3} + \frac{x^5}{5} – \frac{x^7}{7} + \cdots = \sum_{n = 0}^\infty (-1)^n \frac{x^{2n+1}}{2n+1}. \notag

Note that this has interval of convergence $-1 < x \leq 1$.

We also recognize this integral as

\int_0^x \frac{1}{1 + t^2} dt = \text{arctan}(x), \notag

one of the common integrals arising from trigonometric substitution. Putting these together, we find that

\text{arctan}(x) = x – \frac{x^3}{3} + \frac{x^5}{5} – \frac{x^7}{7} + \cdots = \sum_{n = 0}^\infty (-1)^n \frac{x^{2n+1}}{2n+1}. \notag

As $x = 1$ is within the interval of convergence, we can substitute $x = 1$ into the series to find the representation

\text{arctan}(1) = 1 – \frac{1}{3} + \frac{1}{5} – \frac{1}{7} + \cdots = \sum_{n = 0}^\infty (-1)^n \frac{1}{2n+1}. \notag

Since $\text{arctan}(1) = \frac{\pi}{4}$, this gives the representation for $\pi/4$ given in \eqref{eq:base}.

However, since $x=1$ was at the very edge of the interval of convergence, this series converges very, very slowly. For instance, using the first $50$ terms gives the approximation

\pi \approx 3.121594652591011. \notag

The expansion of $\pi$ is actually

\pi = 3.141592653589793238462\ldots \notag

So the first $50$ terms of \eqref{eq:base} gives two digits of accuracy. That’s not very good.

I think it is very natural to ask: can we do better? This series converges slowly — can we find one that converges more quickly?

| | 1 Comment

## “On Functions Whose Mean Value Abscissas are Midpoints, with Connections to Harmonic Functions” (with Paul Carter)

This is joint work with Paul Carter. Humorously, we completed this while on a cross-country drive as we moved the newly minted Dr. Carter from Brown to Arizona.

I’ve had a longtime fascination with the standard mean value theorem of calculus.

Mean Value Theorem
Suppose $f$ is a differentiable function. Then there is some $c \in (a,b)$ such that

\frac{f(b) – f(a)}{b-a} = f'(c).

The idea for this project started with a simple question: what happens when we interpret the mean value theorem as a differential equation and try to solve it? As stated, this is too broad. To narrow it down, we might specify some restriction on the $c$, which we refer to as the mean value abscissa, guaranteed by the Mean Value Theorem.

So I thought to try to find functions satisfying

\frac{f(b) – f(a)}{b-a} = f’ \left( \frac{a + b}{2} \right)

for all $a$ and $b$ as a differential equation. In other words, let’s try to find all functions whose mean value abscissas are midpoints.

This looks like a differential equation, which I only know some things about. But my friend and colleague Paul Carter knows a lot about them, so I thought it would be fun to ask him about it.

He very quickly told me that it’s essentially impossible to solve this from the perspective of differential equations. But like a proper mathematician with applied math leanings, he thought we should explore some potential solutions in terms of their Taylor expansions. Proceeding naively in this way very quickly leads to the answer that those (assumed smooth) solutions are precisely quadratic polynomials.

It turns out that was too simple. It was later pointed out to us that verifying that quadratic polynomials satisfy the midpoint mean value property is a common exercise in calculus textbooks, including the one we use to teach from at Brown. Digging around a bit reveals that this was even known (in geometric terms) to Archimedes.

So I thought we might try to go one step higher, and see what’s up with
\label{eq:original_midpoint}
\frac{f(b) – f(a)}{b-a} = f’ (\lambda a + (1-\lambda) b), \tag{1}

where $\lambda \in (0,1)$ is a weight. So let’s find all functions whose mean value abscissas are weighted averages. A quick analysis with Taylor expansions show that (assumed smooth) solutions are precisely linear polynomials, except when $\lambda = \frac{1}{2}$ (in which case we’re looking back at the original question).

That’s a bit odd. It turns out that the midpoint itself is distinguished in this way. Why might that be the case?

It is beneficial to look at the mean value property as an integral property instead of a differential property,

\frac{1}{b-a} \int_a^b f'(t) dt = f’\big(c(a,b)\big).

We are currently examining cases when $c = c_\lambda(a,b) = \lambda a + (1-\lambda b)$. We can see the right-hand side is differentiable by differentiating the left-hand side directly. Since any point can be a weighted midpoint, one sees that $f$ is at least twice-differentiable. One can actually iterate this argument to show that any $f$ satisfying one of the weighted mean value properties is actually smooth, justifying the Taylor expansion analysis indicated above.

An attentive eye might notice that the midpoint mean value theorem, written as the integral property

\frac{1}{b-a} \int_a^b f'(t) dt = f’ \left( \frac{a + b}{2} \right)

is exactly the one-dimensional case of the harmonic mean value property, usually written

\frac{1}{\lvert B_h \rvert} = \int_{B_h(x)} g(t) dV = g(x).

Here, $B_h(x)$ is the ball of radius $h$ and center $x$. Any harmonic function satisfies this mean value property, and any function satisfying this mean value property is harmonic.

From this viewpoint, functions satisfying our original midpoint mean value property~\eqref{eq:original_midpoint} have harmonic derivatives. But the only one-dimensional harmonic functions are affine functions $g(x) = cx + d$. This gives immediately that the set of solutions to~\eqref{eq:original_midpoint} are quadratic polynomials.

The weighted mean value property can also be written as an integral property. Trying to connect it similarly to harmonic functions led us to consider functions satisfying

\frac{1}{\lvert B_h \rvert} = \int_{B_h(x)} g(t) dV = g(c_\lambda(x,h)),

where $c_\lambda(x,h)$ should be thought of as some distinguished point in the ball $B_h(x)$ with a weight parameter $\lambda$. More specifically,

Are there weighted harmonic functions corresponding to a weighted harmonic mean value property?
In one dimension, the answer is no, as seen above. But there are many more multivariable harmonic functions [in fact, I’ve never thought of harmonic functions on $\mathbb{R}^1$ until this project, as they’re too trivial]. So maybe there are weighted harmonic functions in higher dimensions?

This ends up being the focus of the latter half of our paper. Unexpectedly (to us), an analogous methodology to our approach in the one-dimensional case works, with only a few differences.

It turns out that no, there are no weighted harmonic functions on $\mathbb{R}^n$ other than trivial extensions of harmonic functions from $\mathbb{R}^{n-1}$.

Harmonic functions are very special, and even more special than we had thought. The paper is a fun read, and can be found on the arxiv now. It has been accepted and will appear in American Mathematical Monthly.

## Math 420: Supplement on Gaussian Integers II

This is a secondary supplemental note on the Gaussian integers, written for my Spring 2016 Elementary Number Theory Class at Brown University. This note is also available as a pdf document.

In this note, we cover the following topics.

1. Assumed prerequisites from other lectures.
2. Which regular integer primes are sums of squares?
3. How can we classify all Gaussian primes?

1. Assumed Prerequisites

Although this note comes shortly after the previous note on the Gaussian integers, we covered some material from the book in the middle. In particular, we will assume use the results from chapters 20 and 21 from the textbook.

Most importantly, for ${p}$ a prime and ${a}$ an integer not divisible by ${p}$, recall the Legendre symbol ${\left(\frac{a}{p}\right)}$, which is defined to be ${1}$ if ${a}$ is a square mod ${p}$ and ${-1}$ if ${a}$ is not a square mod ${p}$. Then we have shown Euler’s Criterion, which states that

$$a^{\frac{p-1}{2}} \equiv \left(\frac{a}{p}\right) \pmod p, \tag{1}$$
and which gives a very efficient way of determining whether a given number ${a}$ is a square mod ${p}$.

We used Euler’s Criterion to find out exactly when ${-1}$ is a square mod ${p}$. In particular, we concluded that for each odd prime ${p}$, we have

$$\left(\frac{-1}{p}\right) = \begin{cases} 1 & \text{ if } p \equiv 1 \pmod 4 \ -1 & \text{ if } p \equiv 3 \pmod 4 \end{cases}. \tag{2}$$
Finally, we assume familiarity with the notation and ideas from the previous note on the Gaussian integers.

2. Understanding When ${p = a^2 + b^2}$.

Throughout this section, ${p}$ will be a normal odd prime. The case ${p = 2}$ is a bit different, and we will need to handle it separately. When used, the letters ${a}$ and ${b}$ will denote normal integers, and ${q_1,q_2}$ will denote Gaussian integers.

We will be looking at the following four statements.

1. ${p \equiv 1 \pmod 4}$
2. ${\left(\frac{-1}{p}\right) = 1}$
3. ${p}$ is not a Gaussian prime
4. ${p = a^2 + b^2}$

Our goal will be to show that each of these statements are equivalent. In order to show this, we will show that

$$(1) \implies (2) \implies (3) \implies (4) \implies (1). \tag{3}$$
Do you see why this means that they are all equivalent?

This naturally breaks down into four lemmas.

We have actually already shown one.

Lemma 1 ${(1) \implies (2)}$.

Proof: We have already proved this claim! This is exactly what we get from Euler’s Criterion applied to ${-1}$, as mentioned in the first section. $\Box$

There is one more that is somewhat straightfoward, and which does not rely on going up to the Gaussian integers.

Lemma 2 ${(4) \implies (1)}$.

Proof: We have an odd prime ${p}$ which is a sum of squares ${p = a^2 + b^2}$. If we look mod ${4}$, we are led to consider $$p = a^2 + b^2 \pmod 4. \tag{4}$$
What are the possible values of ${a^2 \pmod 4}$? A quick check shows that the only possibilites are ${a^2 \equiv 0, 1 \pmod 4}$.

So what are the possible values of ${a^2 + b^2 \pmod 4}$? We must have one of ${p \equiv 0, 1, 2 \pmod 4}$. Clearly, we cannot have ${p \equiv 0 \pmod 4}$, as then ${4 \mid p}$. Similarly, we cannot have ${p \equiv 2 \pmod 4}$, as then ${2 \mid p}$. So we necessarily have ${p \equiv 1 \pmod 4}$, which is what we were trying to prove. $\Box$

For the remaining two pieces, we will dive into the Gaussian integers.

Lemma 3 ${(2) \implies (3)}$.

Proof: As ${\left(\frac{-1}{p}\right) = 1}$, we know there is some ${a}$ so that ${a^2 \equiv -1 \pmod p}$. Rearranging, this becomes ${a^2 + 1 \equiv 0 \pmod p}$.

Over the normal integers, we are at an impasse, as all this tells us is that ${p \mid (a^2 + 1)}$. But if we suddenly view this within the Gaussian integers, then ${a^2 + 1}$ factors as ${a^2 + 1 = (a + i)(a – i)}$.

So we have that ${p \mid (a+i)(a-i)}$. If ${p}$ were a Gaussian prime, then we would necessarily have ${p \mid (a+i)}$ or ${p \mid (a-i)}$. (Do you see why?)

But is it true that ${p}$ divides ${a + i}$ or ${a – i}$? For instance, does ${p}$ divide ${a + i}$? No! If so, then ${\frac{a}{p} + \frac{i}{p}}$ would be a Gaussian integer, which is clearly not true.

So ${p}$ does not divide ${a + i}$ or ${a-i}$, and we must therefore conclude that ${p}$ is not a Gaussian prime. $\Box$

Lemma 4 ${(3) \implies (4)}$.

Proof: We now know that ${p}$ is not a Gaussian prime. In particular, this means that ${p}$ is not irreducible, and so it has a nontrivial factorization in the Gaussian integers. (For example, ${5}$ is a regular prime, but it is not a Gaussian prime. It factors as ${5 = (1 + 2i)(1 – 2i)}$ in the Gaussian integers.)

Let’s denote this nontrivial factorization as ${p = q_1 q_2}$. By nontrivial, we mean that neither ${q_1}$ nor ${q_2}$ are units, i.e. ${N(q_1), N(q_2) > 1}$. Taking norms, we see that ${N(p) = N(q_1) N(q_2)}$.

We can evaluate ${N(p) = p^2}$, so we have that ${p^2 = N(q_1) N(q_2)}$. Both ${N(q_1)}$ and ${N(q_2)}$ are integers, and their product is ${p^2}$. Yet ${p^2}$ has exactly two different factorizations: ${p^2 = 1 \cdot p^2 = p \cdot p}$. Since ${N(q_1), N(q_2) > 1}$, we must have the latter.

So we see that ${N(q_1) = N(q_2) = p}$. As ${q_1, q_2}$ are Gaussian integers, we can write ${q_1 = a + bi}$ for some ${a, b}$. Then since ${N(q_1) = p}$, we see that ${N(q_1) = a^2 + b^2}$. And so ${p}$ is a sum of squares, ending the proof. $\Box$

Notice that ${2 = 1 + 1}$ is also a sum of squares. Then all together, we can say the following theorem.

Theorem 5 A regular prime ${p}$ can be written as a sum of two squares, $$p = a^2 + b^2, \tag{5}$$
exactly when ${p = 2}$ or ${p \equiv 1 \pmod 4}$.

A remarkable aspect of this theorem is that it is entirely a statement about the behaviour of the regular integers. Yet in our proof, we used the Gaussian integers in a very fundamental way. Isn’t that strange?

You might notice that in the textbook, Dr. Silverman presents a proof that does not rely on the Gaussian integers. While interesting and clever, I find that the proof using the Gaussian integers better illustrates the deep connections between and around the structures we have been studying in this course so far. Everything connects!

Example 1 The prime ${5}$ is ${1 \pmod 4}$, and so ${5}$ is a sum of squares. In particular, ${5 = 1^2 + 2^2}$.

Example 2 The prime ${101}$ is ${1 \pmod 4}$, and so is a sum of squares. Our proof is not constructive, so a priori we do not know what squares sum to ${101}$. But in this case, we see that ${101 = 1^2 + 10^2}$.

Example 3 The prime ${97}$ is ${1 \pmod 4}$, and so it also a sum of squares. It’s less obvious what the squares are in this case. It turns out that ${97 = 4^2 + 9^2}$.

Example 4 The prime ${43}$ is ${3 \pmod 4}$, and so is not a sum of squares.

3. Classification of Gaussian Primes

In the previous section, we showed that each integer prime ${p \equiv 1 \pmod 4}$ actually splits into a product of two Gaussian numbers ${q_1}$ and ${q_2}$. In fact, since ${N(q_1) = p}$ is a regular prime, ${q_1}$ is a Gaussian irreducible and therefore a Gaussian prime (can you prove this? This is a nice midterm question.)

So in fact, ${p \equiv 1 \pmod 4}$ splits in to the product of two Gaussian primes ${q_1}$ and ${q_2}$.

In this way, we’ve found infinitely many Gaussian primes. Take a regular prime congruent to ${1 \pmod 4}$. Then we know that it splits into two Gaussian primes. Further, if we know how to write ${p = a^2 + b^2}$, then we know that ${q_1 = a + bi}$ and ${q_2 = a – bi}$ are those two Gaussian primes.

In general, we will find all Gaussian primes by determining their interaction with regular primes.

Suppose ${q}$ is a Gaussian prime. Then on the one hand, ${N(q) = q \overline{q}}$. On the other hand, ${N(q) = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}}$ is some regular integer. Since ${q}$ is a Gaussian prime (and so ${q \mid w_1 w_2}$ means that ${q \mid w_1}$ or ${q \mid w_2}$), we know that ${q \mid p_j}$ for some regular integer prime ${p_j}$.

So one way to classify Gaussian primes is to look at every regular integer prime and see which Gaussian primes divide it. We have figured this out for all primes ${p \equiv 1 \pmod 4}$. We can handle ${2}$ by noticing that ${2 = (1 + i) (1-i)}$. Both ${(1+i)}$ and ${(1-i)}$ are Gaussian primes.

The only primes left are those regular primes with ${p \equiv 3 \pmod 4}$. We actually already covered the key idea in the previous section.

Lemma 6 If ${p \equiv 3 \pmod 4}$ is a regular prime, then ${p}$ is also a Gaussian prime.

Proof: In the previous section, we showed that if ${p}$ is not a Gaussian prime, then ${p = a^2 + b^2}$ for some integers ${a,b}$, and then ${ p \equiv 1 \pmod 4}$. Since ${p \not \equiv 1 \pmod 4}$, we see that ${p}$ is a Gaussian prime. $\Box$

In total, we have classified all Gaussian primes.

Theorem 7 The Gaussian primes are given by

1. ${(1+i), (1-i)}$
2. Regular primes ${p \equiv 3 \pmod 4}$
3. The factors ${q_1 q_2}$ of a regular prime ${p \equiv 1 \pmod 4}$. Further, these primes are given by ${a \pm bi}$, where ${p = a^2 + b^2}$.

4. Concluding Remarks

I hope that it’s clear that the regular integers and the Gaussian integers are deeply connected and intertwined. Number theoretic questions in one constantly lead us to investigate the other. As one dives deeper into number theory, more and different integer-like rings appear, all deeply connected.

Each time I teach the Gaussian integers, I cannot help but feel the sense that this is a hint at a deep structural understanding of what is really going on. The interplay between the Gaussian integers and the regular integers is one of my favorite aspects of elementary number theory, which is one reason why I deviated so strongly from the textbook to include it. I hope you enjoyed it too.