超椭圆曲线上的Frobenius矩阵

Sage对Harvey-Kedlaya算法进行了高度优化,用于计算有限域上与曲线相关的Frobenius矩阵。这是davidharvey的一个实现,它是GPL的,只依赖于NTL和zn_poly(Sage中用于快速算法的C库) \((\ZZ/n\ZZ)[x]\)

我们导入hypellfrob函数,并在多项式上调用它 \(\ZZ\) .

sage: from sage.schemes.hyperelliptic_curves.hypellfrob import hypellfrob
sage: R.<x> = PolynomialRing(ZZ)
sage: f = x^5 + 2*x^2 + x + 1; p = 101
sage: M = hypellfrob(p, 1, f); M
[     O(101)      O(101) 93 + O(101) 62 + O(101)]
[     O(101)      O(101) 55 + O(101) 19 + O(101)]
[     O(101)      O(101) 65 + O(101) 42 + O(101)]
[     O(101)      O(101) 89 + O(101) 29 + O(101)]

我们做同样的计算但是 \(\ZZ/101^4\ZZ\) 这就给了足够的精度来识别精确的特征多项式 \(\ZZ[x]\) 作为自同态环的一个元素。这个计算仍然很快,只需要一秒的时间。

sage: M = hypellfrob(p, 4, f)   # about 0.25 seconds
sage: M[0,0]
91844754 + O(101^4)

Frobenius的特征多项式是 \(x^4 + 7x^3 + 167x^2 + 707x + 10201\) ,它决定了 \(\zeta\) 曲线函数 \(y^2= f(x)\) .

sage: M.charpoly()
(1 + O(101^4))*x^4 + (7 + O(101^4))*x^3 + (167 + O(101^4))*x^2 + (707 + O(101^4))*x + 10201 + O(101^4)

注解

练习:显式地写下zeta函数,在一些有限的域上计算点数,看看是否匹配。

模块化符号

模符号在用模形式计算的算法中起着关键作用 \(L\) -函数,椭圆曲线和模交换变种。Sage拥有最普遍的模块化符号实现,这要归功于我,jordiquer(巴塞罗那)和Craig Citro(Hida的学生)的工作。此外,用模符号计算是迄今为止我最喜欢的计算数学部分。对于Sage中的模块化符号,仍有许多调整和优化工作要做,以便使其成为世界上最快的实现,因为我的Magma实现在某些重要情况下仍然更好。

注解

我将在周五的讲座中更多地讨论模符号,这将是关于模的形式和相关的对象。

我们创造空间 \(M\) 重量的 \(4\) 一类同余子群的模符号 \(\Gamma_H(13)\) 水平 \(13\) . 然后我们计算这个空间的基,用Manin符号表示。最后,我们计算了Hecke算子 \(T_2\) 作用于 \(M\) 求其特征多项式,并将其因子化。我们还计算了尖点子空间的维数。

sage: M = ModularSymbols(GammaH(13,[3]), weight=4)
sage: M
Modular Symbols space of dimension 14 for Congruence Subgroup
Gamma_H(13) with H generated by [3] of weight 4 with sign 0
and over Rational Field
sage: M.basis()
([X^2,(0,4)],
[X^2,(0,7)],
[X^2,(4,10)],
[X^2,(4,11)],
[X^2,(4,12)],
[X^2,(7,3)],
[X^2,(7,5)],
[X^2,(7,6)],
[X^2,(7,7)],
[X^2,(7,8)],
[X^2,(7,9)],
[X^2,(7,10)],
[X^2,(7,11)],
[X^2,(7,12)])
sage: factor(charpoly(M.T(2)))
(x - 7) * (x + 7) * (x - 9)^2 * (x + 5)^2
        * (x^2 - x - 4)^2 * (x^2 + 9)^2
sage: dimension(M.cuspidal_subspace())
10

{Cremona's Modular Symbols Library} Sage includes John Cremona's specialized and insanely fast implementation of modular symbols for weight 2 and trivial character. We illustrate below computing the space of modular symbols of level 20014, which has dimension \(5005\), along with a Hecke operator on this space. The whole computation below takes only a few seconds; a similar computation takes a few minutes using Sage's generic modular symbols code. Moreover, Cremona has done computations at levels over 200,000 using his library, so the code is known to scale well to large problems. The new code in Sage for modular symbols is much more general, but doesn't scale nearly so well (yet).

sage: M = CremonaModularSymbols(20014)       # few seconds
sage: M
Cremona Modular Symbols space of dimension 5005 for
Gamma_0(20014) of weight 2 with sign 0
sage: t = M.hecke_matrix(3)             # few seconds

枚举全实数字段

作为他枚举Shimura曲线项目的一部分,johnvoight为Sage提供了枚举完全实数字段的代码。该算法并不复杂,但它涉及一些“内部循环”,必须对其进行编码才能运行得非常快。使用Cython,Voight能够精确地实现牛顿迭代的变体,这是他解决问题所需要的。

函数 enumerate_totallyreal_fields_prim(n, B, ...) 不使用数据库(!)枚举本原(没有真子域)完全实的度域 \(n>1\) 有区别的 \(d \leq B\) .

我们计算了判别式的全实二次域 \(\leq 50\) . 下面的计算几乎是即时的,是实时完成的,不是表查找。

sage: enumerate_totallyreal_fields_prim(2,50)
[[5, x^2 - x - 1], [8, x^2 - 2], [12, x^2 - 3], [13, x^2 - x - 3],
 [17, x^2 - x - 4], [21, x^2 - x - 5], [24, x^2 - 6], [28, x^2 - 7],
 [29, x^2 - x - 7], [33, x^2 - x - 8], [37, x^2 - x - 9],
 [40, x^2 - 10], [41, x^2 - x - 10], [44, x^2 - 11]]

我们计算所有的全实五次判别域 \(\leq 10^5\) . 同样,这是实时完成的-这不是一个表查找!

sage: enumerate_totallyreal_fields_prim(5,10^5)
[[14641, x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1],
 [24217, x^5 - 5*x^3 - x^2 + 3*x + 1],
 [36497, x^5 - 2*x^4 - 3*x^3 + 5*x^2 + x - 1],
 [38569, x^5 - 5*x^3 + 4*x - 1],
 [65657, x^5 - x^4 - 5*x^3 + 2*x^2 + 5*x + 1],
 [70601, x^5 - x^4 - 5*x^3 + 2*x^2 + 3*x - 1],
 [81509, x^5 - x^4 - 5*x^3 + 3*x^2 + 5*x - 2],
 [81589, x^5 - 6*x^3 + 8*x - 1],
 [89417, x^5 - 6*x^3 - x^2 + 8*x + 3]]

伯努利数

Mathematica和Pari

来自Mathematica网站:

“今天我们打破了伯努利记录:从分析引擎到Mathematica 2008年4月29日,Oleksandr Pavlyk,核心技术,一周前,我拿到了Mathematica的最新开发版本,我输入了 BernoulliB[10^7] . 然后我等着。昨天——5天23小时51分37秒之后——我得到了结果!”

tomboothby在Sage中做了同样的计算,Sage使用Pari的bernfrac命令,该命令使用 \(\zeta\) 高精度阶乘,耗时2天12小时。

大卫哈维的伯恩姆

然后大卫哈维提出了一个全新的算法,可以很好地并行化。他给出了计算的时间安排 \(B_{{10^7}}\) 在他的机器上(在我的16核1.8ghz Opteron box上需要59分57秒):

PARI: 75 h, Mathematica: 142 h

bernmm (1 core) = 11.1 h, bernmm (10 cores) = 1.3 h

“在10个核心上运行5.5天,我 [大卫哈维] 计算 [伯努利数] \(B_k\) 对于 \(k = 10^8\) ,我相信这是一个新纪录。本质上,这是我之前在这个线程中建议的多模式算法,但是我想出了一些技巧来优化 \(B_k \text{{mod}} p\) ."

伯努利是世界上最快的数字。以下时间安排在24核2.66Ghz Xeon机箱上。

sage: w1 = bernoulli(100000, num_threads=16)     # 0.9 seconds wall time
sage: w2 = bernoulli(100000, algorithm='pari')   # long time (6s on sage.math, 2011)
sage: w1 == w2  # long time
True

多项式算法

一元多项式算法

Sage在中使用了billhart和davidharvey的GPL的Flint C库进行算术运算 \(\ZZ[x]\) . 它的主要名声是它是世界上多项式乘法最快的,例如,在下面的基准中,它在某些系统上比NTL和Magma快(尽管这些基准当然随着软件的改进而改变)。在幕后,Flint包含一些精心调整的离散傅立叶变换代码。

sage: Rflint = PolynomialRing(ZZ, 'x')
sage: f = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
sage: g = Rflint([ZZ.random_element(2^64) for _ in [1..32]])
sage: timeit('f*g')               # random output
625 loops, best of 3: 105 microseconds per loop
sage: Rntl = PolynomialRing(ZZ, 'x', implementation='NTL')
sage: f = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
sage: g = Rntl([ZZ.random_element(2^64) for _ in [1..32]])
sage: timeit('f*g')               # random output
625 loops, best of 3: 310 microseconds per loop
sage: ff = magma(f); gg = magma(g)  #optional - magma
sage: s = 'time v := [%s * %s : i in [1..10^5]];'%(ff.name(), gg.name()) #optional - magma
sage: magma.eval(s)     #optional - magma
'Time: ...'

奇异:多元多项式算法

在许多情况下,多元多项式算法在库模式中使用奇异(由于Martin Albrecht),这是相当快的。例如,下面我们对32003阶有限域进行法特曼基准测试,并将时间与岩浆进行比较。

sage: P.<x,y,z> = GF(32003)[]
sage: p = (x+y+z+1)^10
sage: q = p+1
sage: timeit('p*q')   # random output
125 loops, best of 3: 1.53 ms per loop

sage: p = (x+y+z+1)^20
sage: q = p+1
sage: timeit('p*q')   # not tested - timeout if SAGE_DEBUG=yes
5 loops, best of 3: 384 ms per loop

sage: pp = magma(p); qq = magma(q) #optional - magma
sage: s = 'time w := %s*%s;'%(pp.name(),qq.name()) #optional - magma
sage: magma.eval(s) #optional - magma
'Time: ...'