I. Python 笑传之吃吃饼¶

分数：100 分

背景描述¶

在日常逛 arxiv 的时候发现了这样一篇论文，称其实现了一种名为 HIApy（Human Implicit Automated Python）的 python 解释器，能在兼容 cphyton 的同时将纯 python 代码的运行速度提高 10 倍。刚把标题发到课题组大群，正当要仔细研读其理论基数和具体实现时，浏览器突然崩溃了，重新查找那篇论文和关键词时却怎么也找不到，就好像从来就没有那篇论文一样。

由于阅读时间过短，只记得他采用了人工智能（Labor Intelligence）技术分析 Python 代码，好巧不巧导师私信说有个和大佬交谈得到的课题，希望能用这篇文献中的方法实现一下，并甩了一段 gemini 味浓到溢出的代码过来。显然，这段代码使用 autograd 的自动微分功能对函数 L 求偏导，换成 jax 库的 jit 技术能显著提高运行速度，但导师立刻拒绝了这个提议，并且说要比 jax 做的更好。

百般无奈，你只能实践 HIApy，来完成这个任务。

任务描述¶

优化下面这段 python 代码，并尽可能提高 cal 函数的运行速度。

名词解释：

n：自由度，正整数
q：广义变量，长度为 n 的 64 位浮点 numpy 数组
qdot：广义变量随时间的导数，长度为 n 的 64 位浮点 numpy 数组
qddot：广义变量随时间的二次导数，长度为 n 的 64 位浮点 numpy 数组
L：一个简单的纯函数，执行一系列固定的计算，这些计算包含加，减，乘，取负，以及 sin，cos 这两种三角函数

代码：

import numpy
import autograd
def gc(L,n):
    ff_=autograd.grad(L,0)
    def ff(a,b):
        return ff_(a,b,autograd.numpy)
    mf_=autograd.hessian(L, 1)
    def mf(a,b):
        return mf_(a,b,autograd.numpy)
    cf_=autograd.jacobian(autograd.grad(L,1),0)
    def cf(a,b):
        return cf_(a,b,autograd.numpy)
    def cal(q,qdot,qddot):
        F = ff(q,qdot)
        M = mf(q,qdot)
        C = cf(q,qdot)
        qddot[:]=numpy.linalg.solve(M, F - C.dot(qdot))
    return cal

你可以使用类似下面的代码测试你的 gc 函数。这里有一个简单的 L 函数和相应的时间演化计算实现，实际评测会使用不同的 L 并可能根据 L 的特性选择不一样的时间演化算法。

import matplotlib.pyplot
import time
def check(a,b):
    def l(q,qdot,usermath):
        t=0.0
        t=t+qdot[0]*qdot[0]
        t=t+qdot[1]*qdot[1]+qdot[0]*qdot[0]+qdot[0]*qdot[1]*2.0*usermath.cos(q[0]-q[1])
        v=-2*usermath.cos(q[0])-usermath.cos(q[1])
        return t-v
    n=2
    cal=gc(l,n)
    q=numpy.zeros(n,dtype=numpy.float64)
    qdot=numpy.zeros(n,dtype=numpy.float64)
    k1=numpy.zeros(n,dtype=numpy.float64)
    q[0]=a
    qdot[0]=b
    dt=0.02
    recn=100
    r=numpy.zeros((recn,n),dtype=numpy.float64)
    for i in range(recn):
        for j in range(10):
            cal(q,qdot,k1)
            qdot_=qdot+0.5*dt*k1
            q_=q+0.5*dt*qdot_
            cal(q_,qdot_,k1)
            qdot_=qdot+0.5*dt*k1
            q_=q+0.5*dt*qdot_
            cal(q_,qdot_,k1)
            qdot_=qdot+dt*k1
            q=q+0.5*dt*(qdot_+qdot)
            qdot=qdot_
        r[i]=q
    return r
yl=check(0,0.7)
matplotlib.pyplot.plot(yl)

你会看到一张这样的图片：

测试¶

只需要提交一份 Python 代码，内容为优化后的 gc 函数，文件大小限制为 100 KB
这个任务确实无法有效利用多线程并行计算，因此本题将在 Intel Xeon Platinum 8358 CPU 的单个核心上运行。
仅有一个测试样例，自由度 n = 20
计算总时限 300 s
目标时间 1.0s
分数正比于运行时间的倒数，score = 100/time

注意¶

时间演化的固有的运行时间也会计算在内，这部分大致在 0.25s 左右。
Python 本身已经安装的包不多，但可以使用系统自带的工具和库。
不需要考虑稀疏性的事。
不要试图拿到 L 函数并在提交中硬编码相关内容。
gc 函数的执行时间和后续计算共享时限，没有额外的限制，不参与评分。
不建议在 gc函数中访问本机外的内容，比如pip install