;;;; Implement spectralnorm algorithm from Computer Language Shootout

;;;; usage: csc -O3 spectalnorm.scm
;;;;        ./spectralnorm 500

;(require-extension srfi-4)

(declare (block) (disable-interrupts))

(define (A i j)
 (/ 1.0 (+ (* (+ i j)
              (/ (+ i j 1.0) 2.0))
           i 1.0)))

(define (A-times-u N u Au)
  (do ((i 0 (fx+ i 1)))
      ((fx>= i N))
    (vector-set! Au i
		    (let loop ((j 0) (sum 0))
		      (if (fx>= j N)
			  sum
			  (loop (fx+ j 1)
				(+ sum (* (A i j) (vector-ref u j)))))))))

(define (At-times-u N u Au)
  (do ((i 0 (fx+ i 1)))
      ((fx>= i N))
    (vector-set! Au i
		    (let loop ((j 0) (sum 0))
		      (if (fx>= j N)
			  sum
			  (loop (fx+ j 1)
				(+ sum (* (A j i) (vector-ref u j)))))))))

(define (AtA-times-u N u AtAu)
 (let ((v (make-vector N)))
   (A-times-u N u v)
   (At-times-u N v AtAu)))


#>!
double eval_A(int i, int j)
{
 return (1.0/((i+j)*(i+j+1)/2+i+1));
}

void eval_A_times_u(int N, const double u[], double Au[])
{
 int i,j;
 for(i=0;i<N;i++)
   {
     Au[i]=0;
     for(j=0;j<N;j++) Au[i]+=eval_A(i,j)*u[j];
   }
}

void eval_At_times_u(int N, const double u[], double Au[])
{
 int i,j;
 for(i=0;i<N;i++)
   {
     Au[i]=0;
     for(j=0;j<N;j++) Au[i]+=eval_A(j,i)*u[j];
   }
}

void eval_AtA_times_u(int N, const double u[], double AtAu[])
{ double v[N]; eval_A_times_u(N,u,v); eval_At_times_u(N,v,AtAu); }

double eval_dot(int N, const double a[], const double b[])
{
 int i;
 double res = 0.0;
 for (i = 0; i < N; ++i) {
   res += a[i] * b[i];
 }
 return (res);
}
<#

(define (dot N u v)
 (do ((i 0 (fx+ i 1))
      (sum 0 (+ sum (* (vector-ref u i)
		       (vector-ref v i)))))
     ((fx>= i N) sum) ))

(define (main-1)
 (let* ((N (with-input-from-string (car (command-line-arguments)) read))
        (u (make-vector N 1.0))
        (v (make-vector N)))
   (do ((i 0 (fx+ i 1)))
       ((fx>= i 10))
     (eval_AtA_times_u N u v)
     (eval_AtA_times_u N v u))
   (print (sqrt (/ (dot N v u) (dot N v v))))))

(define (main-2)
 (let* ((N (with-input-from-string (car (command-line-arguments)) read))
        (u (make-vector N 1.0))
        (v (make-vector N)))
   (do ((i 0 (fx+ i 1)))
       ((fx>= i 10))
     (AtA-times-u N u v)
     (AtA-times-u N v u))
   (print (sqrt (/ (dot N v u) (dot N v v))))))

;(time (main-1))
(time (main-2))