anchor.diff:          (fuse*
            (reshape u selection)
            (α #'+
  -            (α #'* 0.25d0   (reshape u interior))
  -            (α #'* 0.125d0  (reshape u (τ (i j) ((1+ i) j)) interior))
  -            (α #'* 0.125d0  (reshape u (τ (i j) ((1- i) j)) interior))
  -            (α #'* 0.125d0  (reshape u (τ (i j) (i (1+ j))) interior))
  -            (α #'* 0.125d0  (reshape u (τ (i j) (i (1- j))) interior))
  -            (α #'* 0.0625d0 (reshape u (τ (i j) ((1+ i) (1+ j))) interior))
  -            (α #'* 0.0625d0 (reshape u (τ (i j) ((1- i) (1+ j))) interior))
  -            (α #'* 0.0625d0 (reshape u (τ (i j) ((1+ i) (1- j))) interior))
  -            (α #'* 0.0625d0 (reshape u (τ (i j) ((1- i) (1- j))) interior))))
  +            (α #'* 1/4
  +               (reshape u interior))
  +            (α #'* 1/8
  +               (α #'+
  +                  (reshape u (τ (i j) ((1+ i) j)) interior)
  +                  (reshape u (τ (i j) ((1- i) j)) interior)
  +                  (reshape u (τ (i j) (i (1+ j))) interior)
  +                  (reshape u (τ (i j) (i (1- j))) interior)))
  +            (α #'* 1/16
  +               (reshape u (τ (i j) ((1+ i) (1+ j))) interior)
  +               (reshape u (τ (i j) ((1- i) (1+ j))) interior)
  +               (reshape u (τ (i j) ((1+ i) (1- j))) interior)
  +               (reshape u (τ (i j) ((1- i) (1- j))) interior))))
           (τ (i j) ((+ start-1 (/ (- i start-1) 2))
                     (+ start-2 (/ (- j start-2) 2)))))))
       ((~ start-1 1 end-1 ~ start-2 1 end-2 ~ start-3 1 end-3)
  
positive.code:         (fuse*
           (reshape u selection)
           (α #'+
              (α #'* (float 1/8 0d0) (reshape u interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) ((1+ i) j k)) interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) (i (1+ j) k)) interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) (i j (1+ k))) interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) ((1- i) j k)) interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) (i (1- j) k)) interior))
              (α #'* (float 1/16 0d0) (reshape u (τ (i j k) (i j (1- k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) (i (1+ j) (1+ k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1+ i) j (1+ k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1+ i) (1+ j) k)) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) (i (1- j) (1+ k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1- i) j (1+ k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1- i) (1+ j) k)) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) (i (1+ j) (1- k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1+ i) j (1- k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1+ i) (1- j) k)) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) (i (1- j) (1- k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1- i) j (1- k))) interior))
              (α #'* (float 1/32 0d0) (reshape u (τ (i j k) ((1- i) (1- j) k)) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1+ i) (1+ j) (1+ k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1+ i) (1+ j) (1- k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1+ i) (1- j) (1+ k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1+ i) (1- j) (1- k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1- i) (1+ j) (1+ k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1- i) (1+ j) (1- k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1- i) (1- j) (1+ k))) interior))
              (α #'* (float 1/64 0d0) (reshape u (τ (i j k) ((1- i) (1- j) (1- k))) interior))))
          (τ (i j k) ((+ start-1 (/ (- i start-1) 2))
                      (+ start-2 (/ (- j start-2) 2))
                      (+ start-3 (/ (- k start-3) 2)))))))))
  
negative.code:       (rbgs u f h 3) ; solve "exactly"
        (let* ((x (rbgs u f h v1))
               (r (restrict (residual x f h)))
               (c (v-cycle (reshape 0 (shape r)) r (* 2 h) v1 v2)))
          (rbgs (α #'+ x (prolongate c)) f h v2))))
  
  
  