Compute decomposition of $S_2(\Gamma_0(N))$, $N \leq 10000$

def f(N,k):
    """
    INPUT:
        N -- positive integer
        k -- even integer >= 2
    OUTPUT:
        (N,k) -- level and weight
        (t0,t1,t2,t3,tall)  -- timings
        Modular symbols space -- with new cuspidal subspace decomposed.
    """
    tall = cputime(); t = cputime()
    M = ModularSymbols(N,k,sign=1)
    t0 = cputime(t); t = cputime()
    S = M.cuspidal_subspace()
    t1 = cputime(t); t = cputime()
    new = S.new_subspace()
    t2 = cputime(t); t = cputime()
    D = new.decomposition()
    t3 = cputime(t)
    return ( (N, k),  (t0,t1,t2,t3,cputime(tall)),  M)

We estimate that this calculation will take about 230 CPU days using exactly the code in Sage now with no further optimization.