erfinv

В Octave она реализована так:

function [y, iterations] = erfinv (x)

if (nargin != 1) print_usage (); endif

maxit = 100; if (isa (x, «single»)) tol = eps («single»); else tol = eps; endif

iterations = 0;

sz = size (x); nel = numel (x);

x = reshape (x, nel, 1); y = zeros (nel, 1);

## x < -1 or x > 1 ==> NaN y(abs (x) >= 1) = NaN; y(x == -1) = -Inf; y(x == +1) = +Inf;

i = find ((x > -1) & (x < 1)); if (any (i)) s = sqrt (pi) / 2; z = sqrt (-log (1 - abs (x(i)))) .* sign (x(i)); while (any (abs (erf (z) - x(i)) > tol * abs (x(i)))) z = z - (erf (z) - x(i)) .* exp (z.^2) * s; if (++iterations > maxit) warning («erfinv: iteration limit exceeded»); break; endif endwhile y(i) = z; endif

y = reshape (y, sz);

endfunction

Переписать на си или фортран не так уж и сложно. Где-то строчек шесть-семь получится. Тем более, что все вычисления внутри if идут. Также легко гуглятся и более оптимизированные варианты.

Если библиотеки смотреть, то многие коммерческие умеют erfinv. Из бесплатных не знаю, подумал на GSL, но там нет.

Для плюсцов erfinv есть, например, в Boost.

Пардон, форматирование слетело.

Вот октавовский код.

https://gist.github.com/1338250

Если надо, ЛОРовские аналитики всегда помогут с переводом на сишечку.

mkl же.

Столлман не одобряет.