{OperatorNorm(x,tol) = 
    n=length(x); b=vector(2*n,l, 1+(l^2-l)/2);
    lambda=0; lambda0=1; 
    while(abs(lambda-lambda0)>tol, 
        xhat=x/sqrt(sum(k=-n,-1, x[-k]^2));
        y=vector(n,k, sum(l=-n,-1, xhat[-l]/(b[k-l]+l)));
        x=vector(n,j, sum(k=-n,-1,  y[-k]/(b[j-k]-j)));
        lambda0=lambda; lambda=sum(k=-n,-1,xhat[-k]*x[-k]);
    ); 
    [sqrt(lambda),x]
}