anchor.diff:  double complex, intent(in), dimension(:,:):: mat
   double complex, intent(in), dimension(:,:):: with
   double complex, intent(out), dimension(:,:):: output
  -output = matmul(&
  -	with,&
  -	matmul(&
  -		mat , conjg(transpose(&
  -			with &
  -			))&
  -		))
  +output = matmul(with,matmul(mat , conjg(transpose(with))))
   return
   end subroutine sandwich_internal
   
  
positive.code: 		do t_counter = 1,size(time_dependence,2)
  			propagators(:,:,t_counter,p_counter) = 0
  			do op_counter = 1,size(time_dependence,1)
  				propagators(:,:,t_counter,p_counter) = &
  					propagators(:,:,t_counter,p_counter) &
  					+ p(op_counter,p_counter) &
  					* time_dependence(op_counter,t_counter) &
  					* operators(:,:,op_counter)
  			end do
  			call complexexp(propagators(:,:,t_counter,p_counter)) ! this passes a HERMITIAN argument, and gets back the exponent of the 2*pi*i* that argument
  			!call exact(propagators(:,:,t_counter,p_counter)) ! this passes a HERMITIAN argument, and gets back the exponent
  			! for now, i'm not actually propagating, just exponentiating
  			if (t_counter.gt.1) then
  				propagators(:,:,t_counter,p_counter) = &
  					matmul(&
  					propagators(:,:,t_counter,p_counter),&
  					propagators(:,:,t_counter-1,p_counter))
  			endif
  		enddo
  		!write(*,*) 'done propagating parameter ',p_counter,' out of ',size(p,2)
  
negative.code: 		do p_counter = 1,size(rho,4)
  			!write (*,*) 'multiply',t_counter,p_counter
  			if (size(C,3).gt.1) then
  				tempmat = matmul(rho(:,:,t_counter,p_counter),conjg(transpose(C(:,:,p_counter))))
  			else
  				tempmat = matmul(rho(:,:,t_counter,p_counter),conjg(transpose(C(:,:,1))))
  			endif
  			traceout(t_counter,p_counter) = 0
  			call trace_internal(tempmat,traceout(t_counter,p_counter))
  