Come impedisco più efficacemente che la mia variabile casuale normalmente distribuita sia zero?

Question

Sto scrivendo un algoritmo Monte Carlo, in cui ad un certo punto ho bisogno di dividere per una variabile casuale. Più precisamente: la variabile casuale viene utilizzata come larghezza del passaggio per un quoziente di differenza, quindi in realtà in realtà moltiplico qualcosa per la variabile e poi di nuovo dividi da una funzione localmente lineare di questa espressione. Piace

double f(double);

std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<> >
  r( std::tr1::mt19937(time(NULL)),
     std::tr1::normal_distribution<>(0) );

double h = r();
double a = ( f(x+h) - f(x) ) / h;

Funziona bene per la maggior parte del tempo, ma fallisce quando h=0. Matematicamente, questo non è una preoccupazione perché in nessuna selezione finita (o, in effetti, numerabile) di variabili casuali normalmente distribuite, tutte saranno diverse con probabilità 1. Ma nell'implementazione digitale incontrerò un h==0 Ogni funzione di ≈2³² chiama (indipendentemente dal twister di Mersenne che ha un periodo più lungo dell'universo, produce ancora ordinario longS!).

È abbastanza semplice evitare questo problema, al momento sto facendo

double h = r();
while (h==0) h=r();

Ma non lo considero particolarmente elegante. C'è un modo migliore?

---

La funzione che sto valutando non è in realtà solo un semplice ℝ-> ℝ come f è, ma un ℝᵐxℝⁿ -> ℝ in cui calcolo il gradiente nelle variabili ℝᵐ mentre integrano numericamente sulle variabili ℝⁿ. L'intera funzione è sovrapposta con rumore imprevedibile (ma "coerente"), a volte con frequenze in sospeso specifiche (ma sconosciute), questo è ciò che mi mette nei guai quando lo provo con valori fissi per h.