One idea is to integrate once to get the sample points then compute the remaining integrals as sums:

     sample = Transpose@
        SortBy[First@Last@Reap[
           NIntegrate[x (c = Cos[10 x + Cos[x]]), {x, 0, 5}, 
                      EvaluationMonitor :> Sow[{x, c}]]], #[[1]] &];
     wt = ((#[[3]] - #[[1]])/2) & /@ Partition[Join[{0}, sample[[1]], {5}], 3, 1];
     wt.(sample[[2]] #) & /@ {sample[[1]], sample[[1]]^2, sample[[1]]^3}

> {0.0133333, 0.133275, 0.861541}

Note the accuracy is not terribly good, NIntegrate gives:

> {0.0125266, 0.131514, 0.855716}

Somethings a bit off in my quick&dirty trapezoid integration but i think this can be made to work.

## roll your own ##

For this example there really is little benefit to `NIntegrate`'s  adaptive sampling so we might as well just use a uniform sampling:

     np = 650;
     wt = 5/(np - 1) Join[{ 1/2},  ConstantArray[ 1 , np - 2], {1 /2}] // N;
     x = ( Range[0, np - 1] 5/(np - 1)) // N;
     fast = Cos[10 # + Cos[#]] & /@ x // N;
     (# fast) .wt & /@ {x, x^2, x^3}

> {0.0125239, 0.131539, 0.855964}