Your errors are dominated by overflow. These kind of functions are normally computed in a logorithmic maniford, not the real number manifold.
Consider: pow(x,254)
almost always in overflow for most of your x range.pow(M_E,-x)
at underflow for most of your x range
Those two can be combined: exp( log(x)*256-x), that reduces the over/under flow, and obviously doing the whole problem in the log space reduces it still further.
Your range from x to 100000 is going to do that unless you convert the integration to a log basis, if you calculated log_chi_cdf, you should be ok -- then just convert at the end