michael@0: // The Great Computer Language Shootout
michael@0: // http://shootout.alioth.debian.org/
michael@0: //
michael@0: // contributed by Ian Osgood
michael@0: 
michael@0: function A(i,j) {
michael@0:   return 1/((i+j)*(i+j+1)/2+i+1);
michael@0: }
michael@0: 
michael@0: function Au(u,v) {
michael@0:   for (var i=0; i<u.length; ++i) {
michael@0:     var t = 0;
michael@0:     for (var j=0; j<u.length; ++j)
michael@0:       t += A(i,j) * u[j];
michael@0:     v[i] = t;
michael@0:   }
michael@0: }
michael@0: 
michael@0: function Atu(u,v) {
michael@0:   for (var i=0; i<u.length; ++i) {
michael@0:     var t = 0;
michael@0:     for (var j=0; j<u.length; ++j)
michael@0:       t += A(j,i) * u[j];
michael@0:     v[i] = t;
michael@0:   }
michael@0: }
michael@0: 
michael@0: function AtAu(u,v,w) {
michael@0:   Au(u,w);
michael@0:   Atu(w,v);
michael@0: }
michael@0: 
michael@0: function spectralnorm(n) {
michael@0:   var i, u=[], v=[], w=[], vv=0, vBv=0;
michael@0:   for (i=0; i<n; ++i) {
michael@0:     u[i] = 1; v[i] = w[i] = 0;
michael@0:   }
michael@0:   for (i=0; i<10; ++i) {
michael@0:     AtAu(u,v,w);
michael@0:     AtAu(v,u,w);
michael@0:   }
michael@0:   for (i=0; i<n; ++i) {
michael@0:     vBv += u[i]*v[i];
michael@0:     vv  += v[i]*v[i];
michael@0:   }
michael@0:   return Math.sqrt(vBv/vv);
michael@0: }
michael@0: 
michael@0: var actual = '';
michael@0: for (var i = 6; i <= 48; i *= 2) {
michael@0:   actual += spectralnorm(i) + ',';
michael@0: }
michael@0: assertEq(actual, "1.2657786149754053,1.2727355112619148,1.273989979775574,1.274190125290389,");