ux=c(0,0,1)
Sx=matrix(c(1,0,-1,0,2,3,-1,3,7),byrow=T,ncol=3)
Sx

b=c(-1,2,1)
A=matrix(c(1,-1,2,3,-3,-1,-1,0,1),byrow=T,ncol=3)

A%*%ux+b   #mean of AX+b
A%*%Sx%*%t(A) #covariance of AX+b

M=matrix(c(1,-1,2,-1,-3,-1,2,-1,-1),byrow=T,ncol=3)
Sx%*%M # SxM
sum(diag(Sx%*%M)) # trace of SxM

#mean of X'MX+17
t(ux)%*%M%*%ux+sum(diag(Sx%*%M))+17 #requres M to be symmetric

c=c(1,0,0)
B=matrix(c(1,2,3,0,1,0,-1,-2,-3),byrow=T,ncol=3)
CovXY=matrix(c(1,0,-1,2,1,-1,-1,3,4),byrow=T,ncol=3)
uy=c(1,2,3)

CovXY[1,2] #covariance of X1 with Y2 is 1,2 element of CovXY

A%*%ux+B%*%uy+b-c  #mean of AX+BY+b-c
A%*%CovXY%*%t(B) #COV(AX+b,BY+c)