http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55645
--- Comment #4 from vincenzo Innocente <vincenzo.innocente at cern dot ch> 2012-12-11 10:43:21 UTC --- sure use c[i]=std::sqrt(a[i]/b[i]); Recent literature is plenty of examples mostly related to GPU code see for instance Random Variable Recycling - Numerical Algorithms Group www.nag.co.uk/.../gpus/random_variable_recycling_note.pdf In particular the erfinv implementations by Mike Giles Approximating the erfinv function - Mathematical Institute people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf (very simple, copy paste from last page or look below!) the Acklam's one http://home.online.no/~pjacklam/notes/invnorm/ (code at the bottom of http://home.online.no/~pjacklam/notes/invnorm/impl/lea/lea.c) or Shaw et al (http://arxiv.org/abs/0901.0638) I'm currently experimenting with Giles' code #include "approx_exp.h" #include "approx_log.h" #include <cmath> #define likely(x) (__builtin_expect(x, true)) inline float erfinv_like(float w) { w = w - 2.500000f; float p = 2.81022636e-08f; p = 3.43273939e-07f + p*w; p = -3.5233877e-06f + p*w; p = -4.39150654e-06f + p*w; p = 0.00021858087f + p*w; p = -0.00125372503f + p*w; p = -0.00417768164f + p*w; p = 0.246640727f + p*w; p = 1.50140941f + p*w; return p; } inline float erfinv_unlike(float w) { w = std::sqrt(w) - 3.000000f; float p = -0.000200214257f; p = 0.000100950558f + p*w; p = 0.00134934322f + p*w; p = -0.00367342844f + p*w; p = 0.00573950773f + p*w; p = -0.0076224613f + p*w; p = 0.00943887047f + p*w; p = 1.00167406f + p*w; p = 2.83297682f + p*w; return p; } inline float erfinv(float x) { float w, p; w = - unsafe_logf<8>((1.0f-x)*(1.0f+x)); if likely( w < 5.000000f ) p = erfinv_like(w); else erfinv_unlike(w); } return p*x; } float a[1024]; float b[1024]; void compute() { for (int i=0;i!=1024;++i) b[i]=erfinv(a[i]); }
