What are good ways to numerically calculate integrals? Background (skip, skim, or read depending on your patience):

I was trying to estimate the value of a difficult arithmetic sum recently and I'm not quite sure how to go about it.

The summands are ~ 1/xlog^2 x and so clearly converge, but a finite sum as far as I was able to manage (a few million terms, though I could reach a billion with a bit more time) aren't enough to get a value I can trust to even two decimal places, and it's surely meaningless beyond four.

The obvious next step was to take the sum as far as I cared to go and use an integral for the rest. Fortunately the terms are just well-behaved enough that I have viable upper and lower bounds. With some effort I managed to generate three integrals that give upper and lower bounds (the lower bound is two piecewise integrals).

So I typed all this up in my programming language of choice and wrote a wrapper function that let me set precision and how far to set the sum/integral cutoff. It was only then that I realized that the integral was apparently rather numerically unstable -- even though the bounding functions were smooth in the domain of consideration. Every time I raised the precision I got a new answer, to the limit of my patience.

My best though so far was to subtract the main term 1/xlog^2 x out of the integral, since it has closed-form integral -1/log x. But this hasn't helped much so far.

I'm just looking for general strategies, either better ways to calculate the integrals or a different approach to the sum. Because the sum includes primes it's a pain to work with numerically... but Pierre Dusart's bounds are tight enough at large numbers that I think I can make the integral approach work. But first I have to know that what I think is the integral is, in fact, the integral!