Category Archives: Programming

Programming Masthead

I maintain the following programming projects:

HNRSS: (source), a HackerNews RSS generator written in python. HNRSS periodically updates RSS feeds from the HN frontpage and best list. It also attempts to automatically summarize the link (if there is a link) and includes the top five comments, all to make it easier to determine whether it’s worth checking out.

LaTeX2Jax: (source), a tool to convert LaTeX documents to HTML with MathJax. This is a modification of the earlier MSE2WP, which converts Math.StackExchange flavored markdown to WordPress+MathJax compatible html. In particular, this is more general, and allows better control of the resulting html by exposing more CSS elements (that generically aren’t available on free WordPress setups). This is what is used for all math posts on this site.

MSE2WP: (source), a tool to convert Math.Stackexchange flavored markdown to WordPress+MathJax compatible html. This was once written for the Math.Stackexchange Community Blog. But as that blog is shutting down, there is much less of a purpose for this script. Note that this began as a modified version of latex2wp.

 

I actively contribute to:

python-markdown2: (source),  a fast and complete python implementation of markdown, with a few additional features.

 

And I generally support or have contributed to:

SageMath: (main site), a free and open source system of tools for mathematics. Some think of it as a free alternative to the “Big M’s” — Maple, Mathematica, Magma.

Matplotlib: (main site), a plotting library in python. Most of the static plots on this site were creating using matplotlib.

crouton: (source), a tool for making Chromebooks, which by default are very limited in capability, into hackable linux laptops. This lets you directly run Linux on the device at the same time as having ChromeOS installed. The only cost is that there is absolutely no physical security at all (and every once in a while a ChromeOS update comes around and breaks lots of things). It’s great!

 

Below, you can find my most recent posts tagged “Programming” on this site.

I will note the following posts which have received lots of positive feedback.

  1. A Notebook Preparing for a Talk at Quebec-Maine
  2. A Brief Notebook on Cryptography
  3. Computing pi with Tools from Calculus (which includes computational tidbits, though no actual programming).
Posted in Programming | Leave a 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

Experimenting with latex2html5: PSTricks to HTML interactivity

I recently learned about about latex2html5, a javascript library which allows one to write LaTeX and PSTricks to produce interactive objects on websites.At its core, it functions in a similar way to MathJax, which is what I use to generate mathematics on this (and my other) sites. As an example of MathJax, I can write the following.

$$ \int_0^1 f(x) dx = F(1) – F(0). $$

The dream of latex2html5 is to be able to describe a diagram using the language of PSTricks inside LaTeX, throw in a bit of sugar to describe how interactivity should work on the web, and then render this to a beautiful svg using javascript.

Unfortunately, I did not try to make this work on WordPress (as WordPress is a bit finicky about how it interacts with javascript). So instead, I wrote a more detailed description about latex2html5, including some examples and some criticisms, on my non-Wordpress website david.lowryduda.com.

 

Posted in Programming | Leave a comment

Revealing zero in fully homomorphic encryption is a Bad Thing

When I was first learning number theory, cryptography seemed really fun and really practical. I thought elementary number theory was elegant, and that cryptography was an elegant application. As I continued to learn more about mathematics, and in particular modern mathematics, I began to realize that decades of instruction and improvement (and perhaps of more useful points of view) have simplified the presentation of elementary number theory, and that modern mathematics is less elegant in presentation.

Similarly, as I learned more about cryptography, I learned that though the basic ideas are very simple, their application is often very inelegant. For example, the basis of RSA follows immediately from Euler’s Theorem as learned while studying elementary number theory, or alternately from Lagrange’s Theorem as learned while studying group theory or abstract algebra. And further, these are very early topics in these two areas of study!

But a naive implementation of RSA is doomed (For that matter, many professional implementations have their flaws too). Every now and then, a very clever expert comes up with a new attack on popular cryptosystems, generating new guidelines and recommendations. Some guidelines make intuitive sense [e.g. don’t use too small of an exponent for either the public or secret keys in RSA], but many are more complicated or designed to prevent more sophisticated attacks [especially side-channel attacks].

In the summer of 2013, I participated in the ICERM IdeaLab working towards more efficient homomorphic encryption. We were playing with existing homomorphic encryption schemes and trying to come up with new methods. One guideline that we followed is that an attacker should not be able to recognize an encryption of zero. This seems like a reasonable guideline, but I didn’t really understand why, until I was chatting with others at the 2017 Joint Mathematics Meetings in Atlanta.

It turns out that revealing zero isn’t just against generally sound advice. Revealing zero is a capital B capital T Bad Thing.

Basic Setup

For the rest of this note, I’ll try to identify some of this reasoning.

In a typical cryptosystem, the basic setup is as follows. Andrew has a message that he wants to send to Beatrice. So Andrew converts the message into a list of numbers $M$, and uses some sort of encryption function $E(\cdot)$ to encrypt $M$, forming a ciphertext $C$. We can represent this as $C = E(M)$. Andrew transmits $C$ to Beatrice. If an eavesdropper Eve happens to intercept $C$, it should be very hard for Eve to recover any information about the original message from $C$. But when Beatrice receives $C$, she uses a corresponding decryption function $D(\cdot)$ to decrypt $C$, $M = d(C)$.

Often, the encryption and decryption techniques are based on number theoretic or combinatorial primitives. Some of these have extra structure (or at least they do with basic implementation). For instance, the RSA cryptosystem involves a public exponent $e$, a public mod $N$, and a private exponent $d$. Andrew encrypts the message $M$ by computing $C = E(M) \equiv M^e \bmod N$. Beatrice decrypts the message by computing $M = C^d \equiv M^{ed} \bmod N$.

Notice that in the RSA system, given two messages $M_1, M_2$ and corresponding ciphertexts $C_1, C_2$, we have that
\begin{equation}
E(M_1 M_2) \equiv (M_1 M_2)^e \equiv M_1^e M_2^e \equiv E(M_1) E(M_2) \pmod N. \notag
\end{equation}
The encryption function $E(\cdot)$ is a group homomorphism. This is an example of extra structure.

A fully homomorphic cryptosystem has an encryption function $E(\cdot)$ satisfying both $E(M_1 + M_2) = E(M_1) + E(M_2)$ and $E(M_1M_2) = E(M_1)E(M_2)$ (or more generally an analogous pair of operations). That is, $E(\cdot)$ is a ring homomorphism.

This extra structure allows for (a lot of) extra utility. A fully homomorphic $E(\cdot)$ would allow one to perform meaningful operations on encrypted data, even though you can’t read the data itself. For example, a clinic could store (encrypted) medical information on an external server. A doctor or nurse could pull out a cellphone or tablet with relatively little computing power or memory and securely query the medical data. Fully homomorphic encryption would allow one to securely outsource data infrastructure.

A different usage model suggests that we use a different mental model. So suppose Alice has sensitive data that she wants to store for use on EveCorp’s servers. Alice knows an encryption method $E(\cdot)$ and a decryption method $D(\cdot)$, while EveCorp only ever has mountains of ciphertexts, and cannot read the data [even though they have it].

Why revealing zero is a Bad Thing

Let us now consider some basic cryptographic attacks. We should assume that EveCorp has access to a long list of plaintext messages $M_i$ and their corresponding ciphertexts $C_i$. Not everything, but perhaps from small leaks or other avenues. Among the messages $M_i$ it is very likely that there are two messages $M_1, M_2$ which are relatively prime. Then an application of the Euclidean Algorithm gives a linear combination of $M_1$ and $M_2$ such that
\begin{equation}
M_1 x + M_2 y = 1 \notag
\end{equation}
for some integers $x,y$. Even though EveCorp doesn’t know the encryption method $E(\cdot)$, since we are assuming that they have access to the corresponding ciphertexts $C_1$ and $C_2$, EveCorp has access to an encryption of $1$ using the ring homomorphism properties:
\begin{equation}\label{eq:encryption_of_one}
E(1) = E(M_1 x + M_2 y) = x E(M_1) + y E(M_2) = x C_1 + y C_2.
\end{equation}
By multiplying $E(1)$ by $m$, EveCorp has access to a plaintext and encryption of $m$ for any message $m$.

Now suppose that EveCorp can always recognize an encryption of $0$. Then EveCorp can mount a variety of attacks exposing information about the data it holds.

For example, EveCorp can test whether a particular message $m$ is contained in the encrypted dataset. First, EveCorp generates a ciphertext $C_m$ for $m$ by multiplying $E(1)$ by $m$, as in \eqref{eq:encryption_of_one}. Then for each ciphertext $C$ in the dataset, EveCorp computes $C – C_m$. If $m$ is contained in the dataset, then $C – C_m$ will be an encryption of $0$ for the $C$ corresponding to $m$. EveCorp recognizes this, and now knows that $m$ is in the data. To be more specific, perhaps a list of encrypted names of medical patients appears in the data, and EveCorp wants to see if JohnDoe is in that list. If they can recognize encryptions of $0$, then EveCorp can access this information.

And thus it is unacceptable for external entities to be able to consistently recognize encryptions of $0$.

Up to now, I’ve been a bit loose by saying “an encryption of zero” or “an encryption of $m$”. The reason for this is that to protect against recognition of encryptions of $0$, some entropy is added to the encryption function $E(\cdot)$, making it multivalued. So if we have a message $M$ and we encrypt it once to get $E(M)$, and we encrypt $M$ later and get $E'(M)$, it is often not true that $E(M) = E'(M)$, even though they are both encryptions of the same message. But these systems are designed so that it is true that $C(E(M)) = C(E'(M)) = M$, so that the entropy doesn’t matter.

This is a separate matter, and something that I will probably return to later.

Posted in Crypto, Math.NT, Mathematics, Programming | Tagged , | Leave a comment

A Notebook Preparing for a Talk at Quebec-Maine

This is a notebook containing a representative sample of the code I used to  generate the results and pictures presented at the Quebec-Maine Number Theory Conference on 9 October 2016. It was written in a Jupyter Notebook using Sage 7.3, and later converted for presentation on this site.
There is a version of the notebook available on github. Alternately, a static html version without WordPress formatting is available here. Finally, this notebook is also available in pdf form.
The slides for my talk are available here.

Testing for a Generalized Conjecture on Iterated Sums of Coefficients of Cusp Forms

Let $f$ be a weight $k$ cusp form with Fourier expansion

$$ f(z) = \sum_{n \geq 1} a(n) e(nz). $$

Deligne has shown that $a(n) \ll n^{\frac{k-1}{2} + \epsilon}$. It is conjectured that

$$ S_f^1(n) := \sum_{m \leq X} a(m) \ll X^{\frac{k-1}{2} + \frac{1}{4} + \epsilon}. $$

It is known that this holds on average, and we recently showed that this holds on average in short intervals.
(See HKLDW1, HKLDW2, and HKLDW3 for details and an overview of work in this area).
This is particularly notable, as the resulting exponent is only 1/4 higher than that of a single coefficient.
This indicates extreme cancellation, far more than what is implied merely by the signs of $a(n)$ being random.

It seems that we also have

$$ \sum_{m \leq X} S_f^1(m) \ll X^{\frac{k-1}{2} + \frac{2}{4} + \epsilon}. $$

That is, the sum of sums seems to add in only an additional 1/4 exponent.
This is unexpected and a bit mysterious.

The purpose of this notebook is to explore this and higher conjectures.
Define the $j$th iterated sum as

$$ S_f^j(X) := \sum_{m \leq X} S_f^{j-1} (m).$$

Then we numerically estimate bounds on the exponent $\delta(j)$ such that

$$ S_f^j(X) \ll X^{\frac{k-1}{2} + \delta(j) + \epsilon}. $$

In [1]:
# This was written in SageMath 7.3 through a Jupyter Notebook.
# Jupyter interfaces to sage by loading it as an extension
%load_ext sage

# sage plays strangely with ipython. This re-allows inline plotting
from IPython.display import display, Image

We first need a list of coefficients of one (or more) cusp forms.
For initial investigation, we begin with a list of 50,000 coefficients of the weight $12$ cusp form on $\text{SL}(2, \mathbb{Z})$, $\Delta(z)$, i.e. Ramanujan’s delta function.
We will use the data associated to the 50,000 coefficients for pictoral investigation as well.

We will be performing some numerical investigation as well.
For this, we will use the first 2.5 million coefficients of $\Delta(z)$

In [2]:
# Gather 10 coefficients for simple checking
check_10 = delta_qexp(11).coefficients()
print check_10

fiftyk_coeffs = delta_qexp(50000).coefficients()
print fiftyk_coeffs[:10] # these match expected

twomil_coeffs = delta_qexp(2500000).coefficients()
print twomil_coeffs[:10] # these also match expected
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
In [3]:
# Function which iterates partial sums from a list of coefficients

def partial_sum(baselist):
    ret_list = [baselist[0]]
    for b in baselist[1:]:
        ret_list.append(ret_list[-1] + b)
    return ret_list

print check_10
print partial_sum(check_10) # Should be the partial sums
[1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
[1, -23, 229, -1243, 3587, -2461, -19205, 65275, -48368, -164288]
In [4]:
# Calculate the first 10 iterated partial sums
# We store them in a single list list, `sums_list`
# the zeroth elelemnt of the list is the array of initial coefficients
# the first element is the array of first partial sums, S_f(n)
# the second element is the array of second iterated partial sums, S_f^2(n)

fiftyk_sums_list = []
fiftyk_sums_list.append(fiftyk_coeffs) # zeroth index contains coefficients
for j in range(10):                    # jth index contains jth iterate
    fiftyk_sums_list.append(partial_sum(fiftyk_sums_list[-1]))
    
print partial_sum(check_10)
print fiftyk_sums_list[1][:10]         # should match above
    
twomil_sums_list = []
twomil_sums_list.append(twomil_coeffs) # zeroth index contains coefficients
for j in range(10):                    # jth index contains jth iterate
    twomil_sums_list.append(partial_sum(twomil_sums_list[-1]))
    
print twomil_sums_list[1][:10]         # should match above
[1, -23, 229, -1243, 3587, -2461, -19205, 65275, -48368, -164288]
[1, -23, 229, -1243, 3587, -2461, -19205, 65275, -48368, -164288]
[1, -23, 229, -1243, 3587, -2461, -19205, 65275, -48368, -164288]

As is easily visible, the sums alternate in sign very rapidly.
For instance, we believe tha the first partial sums should change sign about once every $X^{1/4}$ terms in the interval $[X, 2X]$.
In this exploration, we are interested in the sizes of the coefficients.
But in HKLDW3, we investigated some of the sign changes of the partial sums.

Now seems like a nice time to briefly look at the data we currently have.
What do the first 50 thousand coefficients look like?
So we normalize them, getting $A(n) = a(n)/n^{5.5}$ and plot these coefficients.

In [5]:
norm_list = []
for n,e in enumerate(fiftyk_coeffs, 1):
    normalized_element = 1.0 * e / (1.0 * n**(5.5))
    norm_list.append(normalized_element)
print norm_list[:10]
1
In [6]:
# Make a quick display
normed_coeffs_plot = scatter_plot(zip(range(1,60000), norm_list), markersize=.02)
normed_coeffs_plot.save("normed_coeffs_plot.png")
display(Image("normed_coeffs_plot.png"))

Since some figures will be featuring prominently in the talk I’m giving at Quebec-Maine, let us make high-quality figures now.

 

(more…)

  1. 00000000000000, -0.530330085889911, 0.598733612492945, -0.718750000000000, 0.691213333204735, -0.317526448138560, -0.376547696558964, 0.911504835123284, -0.641518061271148, -0.366571226366719
Posted in Math.NT, Mathematics, Open, Programming, sagemath | Tagged , , , | 1 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

Programming

I announce the beginning of my (likely intermittent) programming posts.

I’m slowly getting better at python and java, though I spend most of my programming time on python. [And in particular, in sage or using matplotlib].

This category will grow with time. But for now, the most important thing I can link to is my github.

 

Posted in Programming | Tagged , | 1 Comment