Select the correct newform via traces in CMF GP and Magma code snippets

lmfdb/classical_modular_forms/code-form.yaml

@@ -36,11 +36,6 @@ initialize-newspace-common: &initialize-newspace-common
     [N,k,chi] = [{N},{k},Mod({conrey_index},{N})]
     mf = mfinit([N,k,chi],0)
     lf = mfeigenbasis(mf)
-  magma: |
-    // Please install CHIMP (https://github.com/edgarcosta/CHIMP) if you want to run this code
-    chi := DirichletCharacter("{N}.{conrey_index}");
-    S:= CuspForms(chi, {k});
-    N := Newforms(S);
 
 initialize-newspace-weight-1:
   <<: *initialize-newspace-common
@@ -57,12 +52,14 @@ initialize-newspace-weight-not-1:
     from sage.modular.dirichlet import DirichletCharacter
     H = DirichletGroup({N}, base_ring=CyclotomicField({sage_zeta_order}))
     chi = DirichletCharacter(H, H._module({sage_genvalues}))
-    N = Newforms(chi, {k}, names="a") 
-
-coeff-field:
-  comment: Coefficient field, relative polynomial
-  pari: |
-    f.mod \\ as an extension of the cyclotomic field Q(t)/Phi
+    N = Newforms(chi, {k}, names="a")
+  magma: |
+    // Please install CHIMP (https://github.com/edgarcosta/CHIMP) if you want to run this code
+    chi := DirichletCharacter("{N}.{conrey_index}");
+    // Magma's newform constructor for spaces of cusp forms can silently return
+    // incorrect eigenvalues, so we decompose the space of modular symbols instead
+    M := NewSubspace(CuspidalSubspace(ModularSymbols(chi, {k}, -1)));
+    D := NewformDecomposition(M);
 
 relative_polynomial:
   sage: K.relative_polynomial()
@@ -73,7 +70,18 @@ newform:
     {sage_traces}
     f = next(g for g in N if [g.coefficient(i+1).trace() for i in range({sage_trace_bound})] == traces)
   pari: |
-    f = lf[1] \\ Warning: the index may be different
+    {gp_traces}
+    abstrace(a,d) = my(s=1); while(type(a) == "t_POLMOD", s *= poldegree(a.mod); a = trace(a)); a*d/s
+    tracevec(g,n) = my(p=mfparams(g), v=mfcoefs(g,n)); vector(n, i, abstrace(v[i+1], poldegree(p[4])*poldegree(p[5])))
+    f = [g | g <- lf, tracevec(g,#traces) == traces][1]
+  magma: |
+    {magma_traces}
+    f := [d : d in D | [Integers()|Trace(Trace(Coefficient(e,j))) : j in [1..#traces]] eq traces where e := qEigenform(d, #traces+1)][1];
+
+coeff-field:
+  comment: Coefficient field, relative polynomial
+  pari: |
+    f.mod \\ as an extension of the cyclotomic field Q(t)/Phi
 
 coefficient_field:
   sage: |
@@ -90,6 +98,8 @@ qexp:
     f.q_expansion() # note that sage often uses an isomorphic number field
   pari: |
     mfcoefs(f, 20)
+  magma: |
+    qEigenform(f, 20);
 
 embeddings:
   comment: embeddings in the coefficient field
@@ -104,3 +114,25 @@ l-function:
   pari: |
     L = lfunmf(mf,f);
     lfun(L,1)
+
+# Code snippet tests for newform pages; all of these labels are in spaces with
+# several newform orbits, so they exercise the trace-based selection
+# (for 37.2.a the Pari/GP eigenbasis order differs from the LMFDB order,
+# and for 80.5.h the Magma order does, see LMFDB/lmfdb#6059)
+snippet_test:
+  test37_2_a_a:
+    label: 37.2.a.a
+    langs:
+      - sage
+      - gp
+    url: ModularForm/GL2/Q/holomorphic/download_code_newform/37.2.a.a/{lang}
+  test95_1_d_a:
+    label: 95.1.d.a
+    langs:
+      - gp
+    url: ModularForm/GL2/Q/holomorphic/download_code_newform/95.1.d.a/{lang}
+  test80_5_h_b:
+    label: 80.5.h.b
+    langs:
+      - magma
+    url: ModularForm/GL2/Q/holomorphic/download_code_newform/80.5.h.b/{lang}

lmfdb/classical_modular_forms/download.py

@@ -313,6 +313,10 @@ def download_code(self, label, lang):
         if data is None:
             return abort(404, "Label not found: %s" % label)
         form = WebNewform(data)
+        if lang not in form.code_langs:
+            # e.g. Magma for a weight 1 newform, where no snippet can be
+            # generated; emitting a header-only script would be misleading
+            return abort(404, "%s code is not available for %s" % (Fullname[lang], label))
         code = form.code
         # 'initialize-newspace-common' is only a YAML merge anchor base (see
         # code-form.yaml): its snippet is inlined into both weight-specific

lmfdb/classical_modular_forms/test_cmf.py

@@ -1,5 +1,6 @@
 
 from lmfdb.tests import LmfdbTest
+from ast import literal_eval
 import unittest
 
 
@@ -699,17 +700,124 @@ def get(label, lang):
             return self.tc.get('%s/%s/%s' % (base, label, lang)).get_data(as_text=True)
 
         # weight 2 form: the init commands must each appear once...
-        assert get('417.2.d.a', 'magma').count('CuspForms(chi, 2)') == 1
+        assert get('417.2.d.a', 'magma').count('NewformDecomposition(') == 1
         pari = get('417.2.d.a', 'pari')
         assert pari.count('mf = mfinit(') == 1 and pari.count('mfeigenbasis(') == 1, pari
         sage = get('417.2.d.a', 'sage')
         assert sage.count('Newforms(chi, 2, names="a")') == 1, sage
         assert 'cuspidal_submodule' not in sage  # ...and the weight 1 sage variant must not leak in
 
         # weight 1 form: the cuspidal-submodule sage variant is used instead
-        assert get('23.1.b.a', 'magma').count('CuspForms(chi, 1)') == 1
         pari = get('23.1.b.a', 'pari')
         assert pari.count('mf = mfinit(') == 1 and pari.count('mfeigenbasis(') == 1, pari
         sage = get('23.1.b.a', 'sage')
         assert sage.count('cuspidal_submodule().basis()') == 1, sage
         assert 'Newforms(chi' not in sage
+
+    def test_code_download_magma_uses_modular_symbols(self):
+        """The Magma snippet must select the newform from a modular symbols
+        decomposition. Newforms of a space of cusp forms can silently return
+        incorrect eigenvalues in current Magma (e.g. for 32.2.g, 110.2.f and
+        102.2.h), so that computation must not appear at all."""
+        base = '/ModularForm/GL2/Q/holomorphic/download_code_newform'
+        for label in ['417.2.d.a', '80.5.h.b']:
+            magma = self.tc.get('%s/%s/magma' % (base, label)).get_data(as_text=True)
+            assert magma.count('NewformDecomposition(') == 1, magma
+            assert magma.count('ModularSymbols(chi,') == 1, magma
+            # the unreliable computation must not be present, even unused
+            assert 'CuspForms' not in magma, magma
+            assert 'Newforms(S)' not in magma, magma
+            # the form is selected by matching traces, and everything
+            # downstream uses the selected form
+            assert 'f := [d : d in D |' in magma, magma
+            assert 'eq traces' in magma, magma
+            assert 'qEigenform(f, 20);' in magma, magma
+
+        # 80.5.h.b is the regression case from LMFDB/lmfdb#6059: Magma's
+        # decomposition lists this orbit before 80.5.h.a, so an unselected
+        # snippet returns the wrong form. The two orbits first differ at a_5.
+        magma = self.tc.get('%s/80.5.h.b/magma' % base).get_data(as_text=True)
+        assert 'traces := [2,0,0,0,14];' in magma, magma
+        magma = self.tc.get('%s/80.5.h.a/magma' % base).get_data(as_text=True)
+        assert 'traces := [2,0,0,0,-50];' in magma, magma
+
+    def test_code_download_weight_one_has_no_magma(self):
+        """Magma cannot compute weight 1 newforms, so no Magma command snippet
+        should be offered anywhere for a weight 1 form."""
+        from lmfdb.classical_modular_forms.web_newform import WebNewform
+        from lmfdb import db
+
+        form = WebNewform(db.mf_newforms.lookup('23.1.b.a'))
+        code = form.code
+        # no section keeps a magma entry, and magma is gone from both language
+        # selectors ('prompt' is the one CodeSnippet uses, 'show' is the fallback)
+        assert 'magma' not in code['prompt'], code['prompt']
+        assert 'magma' not in code['show'], code['show']
+        assert not [key for key, val in code.items()
+                    if isinstance(val, dict) and 'magma' in val]
+        assert form.code_langs == ['pari', 'sage'], form.code_langs
+
+        base = '/ModularForm/GL2/Q/holomorphic/download_code_newform'
+        # a direct request is rejected rather than returning a header-only or
+        # failing script
+        assert self.tc.get('%s/23.1.b.a/magma' % base).status_code == 404
+        # ...while the other languages are unaffected
+        for lang in ['pari', 'gp', 'sage']:
+            res = self.tc.get('%s/23.1.b.a/%s' % (base, lang))
+            assert res.status_code == 200
+            assert 'mfinit(' in res.get_data(as_text=True) or 'cuspidal_submodule' in res.get_data(as_text=True)
+
+        page = self.tc.get('/ModularForm/GL2/Q/holomorphic/23/1/b/a/').get_data(as_text=True)
+        assert "show_code('magma'" not in page  # no Magma in the language selector
+        assert 'Magma commands' not in page     # nor in the downloads list
+        assert 'CuspForms' not in page          # no Magma code block is rendered
+        # the separate 'Modular form to Magma' data export is still offered
+        assert 'Modular form to Magma' in page
+        # and the other languages are still there
+        assert "show_code('pari'" in page and "show_code('sage'" in page
+        assert 'PariGP commands' in page and 'SageMath commands' in page
+
+        # weight at least 2 still offers Magma everywhere
+        page = self.tc.get('/ModularForm/GL2/Q/holomorphic/80/5/h/b/').get_data(as_text=True)
+        assert "show_code('magma'" in page
+        assert 'Magma commands' in page
+        assert 'NewformDecomposition(' in page
+
+    def test_trace_wrapping(self):
+        """Long trace lists are broken across lines; check that each language
+        gets a syntactically valid multiline assignment which reconstructs the
+        original list. 675.1.g has trace bound 46, so its trace list is long
+        enough to wrap."""
+        from lmfdb.classical_modular_forms.web_newform import wrap_traces
+        from lmfdb import db
+
+        traces = db.mf_newforms.lookup('675.1.g.a', ['traces'])['traces']
+        trace_bound = db.mf_newspaces.lucky({'label': '675.1.g'}, 'trace_bound')
+        traces_list = str(traces[0:trace_bound]).replace(" ", "")
+
+        sage = wrap_traces("traces = " + traces_list)
+        gp = wrap_traces("traces = " + traces_list, " \\")
+        magma = wrap_traces("traces := " + traces_list + ";")
+        for wrapped in [sage, gp, magma]:
+            assert '\n' in wrapped, wrapped  # the point of the test
+            assert all(len(line) <= 75 for line in wrapped.split('\n')), wrapped
+        # gp needs an explicit continuation on every line but the last
+        gp_lines = gp.split('\n')
+        assert all(line.endswith(', \\') for line in gp_lines[:-1]), gp
+        assert not gp_lines[-1].endswith('\\'), gp
+        # sage and magma continue a bracketed list, so they take no continuation
+        assert '\\' not in sage and '\\' not in magma
+        assert magma.endswith(';')
+        # every line break happens after a comma, so nothing is split mid-number
+        # and the assignment reconstructs the original list exactly
+        for wrapped, sep in [(sage, "traces = "), (gp, "traces = "), (magma, "traces := ")]:
+            rhs = wrapped.replace(" \\", "").replace("\n", "").replace(sep, "").rstrip(";")
+            assert literal_eval(rhs) == traces[0:trace_bound], wrapped
+
+        # a short trace list is left on a single line
+        assert wrap_traces("traces = [1,-2]", " \\") == "traces = [1,-2]"
+
+        # and the generated gp snippet really does contain the wrapped form
+        gp_code = self.tc.get(
+            '/ModularForm/GL2/Q/holomorphic/download_code_newform/675.1.g.a/gp').get_data(as_text=True)
+        assert gp in gp_code, gp_code

lmfdb/classical_modular_forms/web_newform.py

@@ -141,6 +141,29 @@ def parity_text(val):
     return 'odd' if val == -1 else 'even'
 
 
+def wrap_traces(traces_string, continuation=""):
+    """Break a trace assignment across lines so that it fits in the code box.
+
+    Line breaks are only inserted after a comma, so that each line is valid in
+    the target language once ``continuation`` (e.g. ``" \\"`` for gp, which
+    requires an explicit line continuation) is appended to it.  Some newspaces
+    have a large trace bound, and hence a long trace list: 3912.1.cp.a is an
+    extreme example.
+    """
+    line_length = 70
+    i = 0
+    wrapped = ""
+    while i < len(traces_string):
+        wrapped += traces_string[i:i + line_length]
+        i += line_length
+        while wrapped[-1] != "," and i < len(traces_string):
+            wrapped += traces_string[i]
+            i += 1
+        if i < len(traces_string):
+            wrapped += continuation + "\n"
+    return wrapped
+
+
 class WebNewform():
     def __init__(self, data, space=None, all_m=False, all_n=False, embedding_label=None):
         # TODO validate data
@@ -412,12 +435,14 @@ def downloads(self):
         else:
             label = '%s.%s' % (self.label, self.embedding_label)
             downloads.append(('Coefficient data to text', url_for('.download_embedded_newform', label=label)))
-        downloads.append(
-                ('Magma commands', url_for(".cmf_code_download", label=self.label, download_type='magma')))
-        downloads.append(
-                ('PariGP commands', url_for(".cmf_code_download", label=self.label, download_type='pari')))
-        downloads.append(
-                ('SageMath commands', url_for(".cmf_code_download", label=self.label, download_type='sage')))
+        # Only offer the languages that this newform actually has snippets for
+        # (weight 1 has no Magma, see WebNewform.code)
+        for lang, name in [('magma', 'Magma commands'),
+                           ('pari', 'PariGP commands'),
+                           ('sage', 'SageMath commands')]:
+            if lang in self.code_langs:
+                downloads.append(
+                    (name, url_for(".cmf_code_download", label=self.label, download_type=lang)))
 
         downloads.append(('Underlying data', url_for('.mf_data', label=label)))
         return downloads
@@ -1444,31 +1469,41 @@ def code(self):
         vals = conrey_chi.genvalues
         sage_genvalues = get_sage_genvalues(self.level, self.char_order, vals, sage_zeta_order)
         sage_trace_bound = self.ns_data.get('trace_bound')
-        traces_string = "traces = "+str(self.traces[0:sage_trace_bound]).replace(" ","")
-        #format string to look nice in the code box if it's long (check 3912/1/cp/a e.g.)
-        line_length = 70
-        i = 0
-        sage_traces_up_to_bound = ""
-        while i < len(traces_string):
-            sage_traces_up_to_bound += traces_string[i:i+line_length]
-            i += line_length
-            while sage_traces_up_to_bound[-1] != "," and i < len(traces_string):
-                sage_traces_up_to_bound += traces_string[i]
-                i += 1
-            if i < len(traces_string):
-                sage_traces_up_to_bound += "\n"
+        traces_list = str(self.traces[0:sage_trace_bound]).replace(" ","")
 
         data = { 'N': self.level,
                  'k': self.weight,
                  'conrey_index': self.conrey_index,
                  'sage_zeta_order': sage_zeta_order,
                  'sage_genvalues': sage_genvalues,
                  'sage_trace_bound': sage_trace_bound,
-                 'sage_traces': sage_traces_up_to_bound,
+                 'sage_traces': wrap_traces("traces = " + traces_list),
+                 'gp_traces': wrap_traces("traces = " + traces_list, " \\"),
+                 'magma_traces': wrap_traces("traces := " + traces_list + ";"),
                }
         for prop in code:
             if not isinstance(code[prop], dict):
                 continue
             for lang in code[prop]:
-                code[prop][lang] = code[prop][lang].format(**data)
+                if isinstance(code[prop][lang], str):
+                    code[prop][lang] = code[prop][lang].format(**data)
+        if self.weight == 1:
+            # Magma does not support weight 1 newforms (Newforms errors out and
+            # modular symbols require weight at least 2), so we cannot select
+            # the newform there.  Drop Magma from every section and from the
+            # language selectors, so that the page, the downloads list and
+            # download_code all agree that Magma is unavailable here.
+            for key, val in code.items():
+                if isinstance(val, dict):
+                    val.pop('magma', None)
+                elif isinstance(val, list):
+                    code[key] = [lang for lang in val if lang != 'magma']
         return code
+
+    @lazy_attribute
+    def code_langs(self):
+        """The languages that code snippets are available in for this newform.
+
+        Weight 1 newforms have no Magma snippet; see :meth:`code`.
+        """
+        return [lang for lang in ('magma', 'pari', 'sage') if lang in self.code['prompt']]

lmfdb/tests/generate_snippet_tests.py

@@ -32,7 +32,10 @@
              'sage_gap': 'sage --simple-prompt',   # Used for evaluating GAP-type group objects defined within Sage
              'magma': 'magma -b',
              'oscar': 'julia',
-             'gp': "sage -gp -D prompt='gp> ' -D breakloop=0 -D colors='no,no,no,no,no,no,no' -D readline=0 -q",
+             # parisizemax lets the PARI stack grow on demand; the 8MB default is
+             # not enough for some snippets (mfinit of a weight 1 newspace of
+             # moderate level already overflows it)
+             'gp': "sage -gp -D prompt='gp> ' -D breakloop=0 -D colors='no,no,no,no,no,no,no' -D readline=0 -D parisizemax=2G -q",
              'gap': """sage -gap -b -T -r -A -m 256m -o 512m -x 800 -c 'SetUserPreference("UseColorsInTerminal",false); SetUserPreference("UseColorPrompt",false); ColorPrompt(false);'""",
              }
 prompt_dict = {'sage': 'sage:', 'sage_gap': 'sage:', 'magma': 'magma> ', 'oscar': 'julia>', 'gp': 'gp> ', 'gap': 'gap> '}
@@ -111,6 +114,12 @@ def _start_snippet_procs(langs, chimp_spec=None):
             #     command never contains the full prompt.
             #
             #  SetColumns(0) disables line-wrapping, and SetAutoColumns(false) stops Magma from re-adjusting to the pty width, so that log files are stable across environments.
+            #
+            #  Note that Magma seeds its random number generator from the clock, and several of its algorithms are randomized
+            #  (e.g. the decomposition underlying NewformDecomposition, and the search underlying Generators of an elliptic
+            #  curve).  The answers are always correct, but a choice of generator can differ between runs, so a Magma log is
+            #  not reproducible: e.g. a q-expansion may come out in terms of a rather than -a.  Calling SetSeed here would fix
+            #  that, but the committed Magma logs were generated without it, so they would all have to be regenerated first.
 
             magma = pexpect.spawn(exec_dict['magma'],
                                   echo=False,
@@ -231,13 +240,17 @@ def create_snippet_tests(yaml_file_path=None, ignore_langs=[], test=False, only_
         contents = yaml.load(code_file.open(), Loader=yaml.FullLoader)
         if 'snippet_test' in contents:
             langs |= set(contents['prompt'].keys())
-    langs -= set(ignore_langs)
+    # The yaml files spell it 'pari', but the executable and the log files are
+    # named after 'gp'.  Normalize before applying --ignore/--only, so that
+    # either spelling works on the command line.
+    def normalize(names):
+        return {'gp' if name == 'pari' else name for name in names}
+
+    langs = normalize(langs)
+    langs -= normalize(ignore_langs)
     if only_langs is not None:
-        langs &= only_langs
+        langs &= normalize(only_langs)
 
-    if 'pari' in langs:
-        langs.remove('pari')
-        langs.add('gp')
     if len(langs) == 0:
         print("No valid languages selected")
         return 1