実際の実装というのはわかりませんが,
C++のdistribution系のソースを見るとかできれば何かわかるのかも??
オフトピック
自分が思いつく方法としては,
ある特定の分布を対象にした実装をするのであれば,積分区間を解くには
「例えばある程度の荒さのヒストグラム的なものを用意しておいて,そこから初期値を決定し,
後はNewton法みたいなのを使ってたかだか数回iteration計算すれば十分」とかなんかそんな感じのを試すかも?
(この時点でもう既に遅いという話もありそうですが…)
分布が2次元でも,1次元の問題×2に分割するんじゃないでしょうか.
(2次元の分布をスライスした断面図の面積を確率とする1次元の分布に関してやる感じで)