Skip to content

Commit 9cfe52b

Browse files
samir-nasiblioleksandr-pavlyk
authored andcommitted
FIX: mkl_random.chisquare when df==1
1 parent 09deb55 commit 9cfe52b

File tree

1 file changed

+4
-15
lines changed

1 file changed

+4
-15
lines changed

mkl_random/src/mkl_distributions.cpp

Lines changed: 4 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -659,25 +659,14 @@ irk_noncentral_chisquare_vec(irk_state *state, npy_intp len, double *res, const
659659

660660
mkl_free(pvec);
661661
} else {
662-
float *fuvec = NULL;
663-
664-
/* noncentral_chisquare(1, nonc) ~ sqrt(nonc)*(-1)^[U<0.5] + Z */
665-
err = vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, state->stream, len, res, d_zero, d_one);
662+
/* noncentral_chisquare(1, nonc) ~ (Z + sqrt(nonc))**2 for df == 1 */
666663
loc = sqrt(nonc);
667-
668-
fuvec = (float *) mkl_malloc(len*sizeof(float), 64);
669-
assert(fuvec != NULL);
670-
671-
err = vsRngUniform(VSL_RNG_METHOD_UNIFORM_STD, state->stream, len, fuvec, (const float) d_zero, (const float) d_one);
664+
err = vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, state->stream, len, res, loc, d_one);
672665
assert(err == VSL_STATUS_OK);
673-
674-
#pragma ivdep
675-
for(i=0; i<len; i++) res[i] += (fuvec[i] < 0.5) ? -loc : loc;
676-
677-
mkl_free(fuvec);
666+
/* squaring could result in an overflow */
667+
vmdSqr(len, res, res, VML_HA);
678668
}
679669
}
680-
681670
}
682671

683672
void

0 commit comments

Comments
 (0)