(define d 4)
(define e 4)
(define (iter n)(lambda(p)(if (positive? n)
  (begin (p (- n 1))((iter (- n 1)) p)))))
(define itr4 (iter d))
(define (g1 p) (let ((nv (make-vector d)))
     (itr4 (lambda (j) (vector-set! nv j (p j))))
     (lambda (x) (vector-ref nv x))))
(define (gn g) (lambda(p)(let ((a (g1 (lambda(w)
   (g (lambda q (apply p (cons w q))))))))
   (lambda (x . y) (apply (a x) y)))))
(define g2 (gn g1))
(define g3 (gn g2))
(define g4 (gn g3))

(define (eta i j) (if (= i j) (if (zero? i) -1 1) 0))
(define ran (let ((seed 4429)(md (expt 2 20)))(lambda ()
  (set! seed (modulo (* seed 17) md)) (modulo seed 11))))
(define (pgen d n gen z) (if (zero? d) z
  (let ((v (make-vector n)))
    ((iter n)(lambda (j)(vector-set! v j
      (pgen (- d 1) (- n j) gen (if (zero? j) z (gen)))))) v)))
(define (paddto n k) (lambda (p1 p2)
    ((iter k)
      (if (= n 1)
        (lambda(j) (vector-set! p2 j
             (+ (vector-ref p1 j)( vector-ref p2 j))))
        (lambda(j) ((paddto (- n 1) (- k j))
          (vector-ref p1 j) (vector-ref p2 j)))))))
(define (padd n k) (lambda (p1 p2) (let ((v (make-vector k)))
   ((iter k) (lambda(j) (vector-set! v j
      ((if (= n 1) + (padd (- n 1) (- k j)))
       (vector-ref p1 j) (vector-ref p2 j))))) v)))
; The multiply operation required by the current job requires specifying
; the degree of the product and it is assured by the outer logic that
; enough terms of the factors will be present.
; For that matter the code "(padd 2 k)" will yield a polynomial
; of degree k-1 even from inputs of greater degrees.
(define (pmul n k) (lambda (p1 p2) (let ((v (make-vector k)))
   ((iter k) (lambda(j) (if (= n 1)

(let ((s 0))((iter (+ j 1))(lambda(i) (set! s (+ s (* (vector-ref p1 i)
    (vector-ref p2 (- j i))))))) (vector-set! v j s)) 

(let ((s (pgen (- n 1) (- k j) (lambda()0) 0))) ((iter (+ j 1)) (lambda (i)
   ((paddto (- n 1) (- k j)) ((pmul (- n 1)(- k j)) (vector-ref p1 i)
        (vector-ref p2 (- j i))) s))) (vector-set! v j s))

))) v)))

(define (mape ms x)(if (vector? x) (let* ((k (vector-length x))
   (nv (make-vector k)))
   ((iter k) (lambda(k)(vector-set! nv k (mape ms (vector-ref x k))))) nv)
   (ms x)))

(define (ip e a b)(let ((sum (pgen d e (lambda()0) 0))
    (ad (padd d e))
    (ml (pmul d e)))
  (itr4(lambda (k) (set! sum (ad sum (ml (a k)(b k))))))
 sum))

(define mm (lambda (x y) (g2 (lambda (i k)
      (ip e (lambda( j)(x i j))(lambda (j)(y j k)))))))
(define (inv y)(let* ((zp (pgen d e (lambda()0) 0))
  (op (pgen d e (lambda()0) 1))
  (mop (pgen d e (lambda()0) -1))
  (et (g2 (lambda (i j)
      (case (eta i j) ((0) zp) ((1) op) ((-1) mop)
        (else (pgen d e (lambda()0) (eta i j)))))))
  (ma (lambda (x y) (g2 (lambda (i j) ((padd d e) (x i j) (y i j))))))
  (s0 (g2 (lambda (i j) zp)))
  ; s0 is a mutable reference to an immutable array of mvps.
  ; Successive values are sums of ever more terms.
  (x (g2 (lambda (i j) ((padd d e) (et i j) (mape - (y i j))))))
  ; y = et - x => x = et - y
  (st et) ; will take on term values
  (etx (g2 (lambda (i k)(ip e (lambda( j)(et i j))(lambda (j)(x j k)))))))
 ((iter e) (lambda (n)(set! s0 (ma s0 st))
    (if (not (zero? n))(set! st (mm etx st)))))s0))
(define (ftl g) (let z ((n d)) (if (zero? n) '()
    (cons (g (- d n)) (z (- n 1))))))
(define (curry g)(lambda(x) (lambda y (apply g (cons x y)))))
(define ml (lambda (k) (lambda (g) (if (zero? k) (g)
   (map (ml (- k 1)) (ftl (curry g)))))))

; To take the partial derivative of a mvp with respect to one of its variables.
; Below (deriv p e j) calls for the derivative of mvp p wrt variable j.
; The input is of degree e (or greater) and the yield will be of degree e-1.
; This is the first function that knows which is the first variable of a mvp!
(define (deriv p e j) (let ((nv (make-vector (- e 1))))
     ((iter ( - e 1))(lambda (k) (vector-set! nv k
         (if (zero? j) (if (zero? k) (vector-ref p 1)
            (mape (lambda (z) (* z (+ k 1))) (vector-ref p (+ k 1))))
           (deriv (vector-ref p k) (- e k) (- j 1)))))) nv))

(define (zt d e n v) (let ((zz (pgen d e (lambda()0) 0)))
  (let z ((vl ((ml n) v))) (or 
    (equal? vl zz)
    (and (pair? vl) (z (car vl)) (z (cdr vl)))
    (null? vl)))))

;test:
(define (test j k)
  (let ((acc (pgen j k (lambda()0) 0))) ((paddto j k)(pgen j k
  (let ((x 0))(lambda()(set! x (+ x 1))x)) 0) acc) acc))
(let ((a (pgen 4 5 ran (ran)))(b (pgen 4 5 ran (ran)))) 
  (let ((x ((padd 4 5) a b))(ac (pgen 4 5 (lambda()0) 0)))
    ((paddto 4 5) a ac)((paddto 4 5) b ac)(equal? ac x)))
(let ((x (pgen 4 5 (lambda () 3) 5))(y (pgen 4 5 (lambda () 3) 5)))
   (equal? x y))
(define (gx n) (pgen 4 5 (let ((x 0))(lambda()(set! x (+ x 1))
   (if (= x n) 1 0))) 0))
(let* ((n (lambda ()(pgen 4 5 ran (ran))))(a (n))(b (n))(c (n))
         (t (pmul 4 5))(p (padd 4 5)))
  (list (equal? (t a (p b c))(p (t a b)(t a c)))
     (equal? (t a (t b c)) (t (t a b) c))))

  
(define Rs
(let* ((gg (g2 (lambda (i j) (pgen d e ran (eta i j)))))
  (g (lambda (i j)(if (< i j) (gg i j) (gg j i))))
  (gi (inv g))
  (gder (g3 (lambda (i j k) (deriv (g j k) 4 i))))
  (2chrst1 (let ((da (padd d (- e 1))))
    (g3 (lambda (i j k) (da (mape - (gder i j k))
          (da (gder j k i) (gder k i j)))))))
  (2chrst (g3 (lambda (i j k)(ip (- e 1)
    (lambda(n) (gi n i)) (lambda (n) (2chrst1 n j k))))))
  (2chder (g4 (lambda (i j k l) (deriv (2chrst j k l) (- e 1) i))))
  (4Riemann (let ((2t (lambda (a)(mape (lambda(x)(* 2 x)) a)))
      (sub (lambda (x y)((padd d (- e 2)) x (mape - y)))))
    (g4 (lambda (m n a b)
     ((padd d (- e 2))(2t (sub (2chder a m n b) (2chder b m n a)))
    (sub (ip (- e 2) (lambda(s)(2chrst m s a))(lambda(s)(2chrst s n b)))
       (ip (- e 2) (lambda(s)(2chrst m s b))(lambda(s)(2chrst s n a)))))))))
  (4R (g4 (lambda (i j k l)
    (ip (- e 2) (lambda(x)(g x i))(lambda(x)(4Riemann x j k l)))))))
(write (let* ((e (- e 1))(pa (padd d e)))
   (zt d e 3 (g3 (lambda (i j k)
      (pa (2chrst i j k) (mape - (2chrst i k j))))))))
(list 2chrst 4R 4Riemann)))

;tests
; Riemann symetries:
(let* ((R (cadr Rs))(e (- e 2))(pa (padd d e))) (list
  (zt d e 4 (g4 (lambda (i j k l) (pa (R i j k l) (R i j l k)))))
  (zt d e 4 (g4 (lambda (i j k l) (pa (R i j k l) (R j i k l)))))
  (zt d e 4 (g4 (lambda (i j k l) (pa (R i j k l) (mape - (R k l i j))))))
  (zt d e 4 (g4 (lambda (i j k l) (pa (pa
     (R i j k l) (R i k l j)) (R i l j k)))))))
; Bianchi identities
(write "Bianchi:")
(let* ((R (caddr Rs))(pa (padd 4 1))
   (Rd ((gn g4) (lambda(i j k l m)(deriv (R i j k l) 2 m)))))
   (zt d 1 5 ((gn g4) (lambda(i j k l m)
     (pa(pa(Rd i j k l m) (Rd i j m k l)) (Rd i j l m k))))))
(write "better Bianchi:")
(let* ((R (caddr Rs))(pa (padd 4 1))(g5 (gn g4))(Gam (car Rs))
   (Rd (g5 (lambda(i j k l m)
     (mape (lambda(a)(* a -2)) (deriv (R i j k l) 2 m))))))
   ((iter 4)(lambda(q)
     (set! Rd (g5 (lambda (i j k l m) (pa (Rd i j k l m)
        (ip 1 (lambda (x)
        (Gam x m (cond ((= q 0) i)((= q 1) j)((= q 2) k)((= q 3) l))))
          (lambda (x)
      (Rd (if (= q 0) x i)(if (= q 1) x j)(if (= q 2) x k)(if (= q 3) x l)  
          m)))))))))
   (zt d 1 5 (g5 (lambda(i j k l m)
     (pa(pa(Rd i j k l m) (Rd i j m k l)) (Rd i j l m k))))))