erf, erfc (誤差関数, 相補誤差関数) を求める PHP の関数

PHP には erf(), erfc() がなかったので、Ruby に Math.erf(), Math.erfc() として実装されている C のソース を PHP に移植しました。プログラムは自由に使用ください。

PHP function

<?php

/*
 reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten
 (New Algorithm handbook in C language) (Gijyutsu hyouronsha, Tokyo, 1991) p.227 [in Japanese]

 C source (for Ruby Interpreter)
 http://www.google.com/codesearch?hl=ja&q=+show:i7gZ8RAuObo:ErgYZrNE2TI:1koF_ZMUv9Q&  sa=N&cd=2&ct=rc&cs_p=ftp://ftp.sunfreeware.com/pub/freeware/SOURCES/ruby-1.8.6.tar.gz&  cs_f=ruby-1.8.6/missing/erf.c
*/

/* Incomplete gamma function 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt  */
function p_gamma($a,$x,$loggamma_a){
 if($x>=1+$a){ return(1-q_gamma($a,$x,$loggamma_a)); }
 if($x==0){ return(0); }
 $result=$term=exp($a*log($x)-$x-$loggamma_a)/$a;
  for($k=1;$k<1000;$k++){
   $term*=$x/($a+$k);
   $previous=$result;
   $result+=$term;
   if($result==$previous){ return($result); }
 }
 trigger_error("d:p_gamma() could not converge.");
 return($result);
}

/* Incomplete gamma function 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt  */
function q_gamma($a,$x,$loggamma_a){
 $la=1;
 $lb=1+$x-$a; /* Laguerre polynomial */

 if($x<1+$a){ return(1-p_gamma($a,$x,$loggamma_a)); }
 $w=exp($a*log($x)-$x-$loggamma_a);
 $result=$w/$lb;
 for ($k=2;$k<1000;$k++) {
  $temp=(($k-1-$a)*($lb-$la)+($k+$x)*$lb)/$k;
  $la=$lb;
  $lb=$temp;
  $w*=($k-1-$a)/$k;
  $temp=$w/($la*$lb);
  $previous=$result;
  $result+=$temp;
  if($result==$previous){ return($result); }
 }
 trigger_error("d:q_gamma() could not converge.");
 return($result);
}

define(LOG_PI_OVER_2,0.572364942924700087071713675675); /* log_e(PI)/2 */

function erf($x){
if(!is_finite($x)){
  if(is_nan($x)){ return($x); } /* erf(NaN) = NaN */
  return($x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */
}
if($x>=0){ return(p_gamma(0.5,$x*$x, LOG_PI_OVER_2)); }
else{ return(-1*p_gamma(0.5,$x*$x, LOG_PI_OVER_2)); }
}

function erfc($x){
 if(!is_finite($x)){
  if(is_nan($x)){return($x); } /* erfc(NaN) = NaN */
  return($x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */
 }
 if($x>=0){ return(q_gamma(0.5,$x*$x,LOG_PI_OVER_2)); }
 else{ return(1+p_gamma(0.5,$x*$x,LOG_PI_OVER_2)); }
}

/* 実行結果 */
echo("erf(-2.0)=".erf(-2.0)."<br />"); // erf(-2.0)=-0.99532226501895
echo("erf(-1.6)=".erf(-1.6)."<br />"); // erf(-1.6)=-0.97634838334464
echo("erf(-1.2)=".erf(-1.2)."<br />"); // erf(-1.2)=-0.91031397822964
echo("erf(-0.8)=".erf(-0.8)."<br />"); // erf(-0.8)=-0.74210096470766
echo("erf(-0.4)=".erf(-0.4)."<br />"); // erf(-0.4)=-0.42839235504667
echo("erf(0.0)=".erf(0.0)."<br />"); // erf(0.0)=0
echo("erf(0.4)=".erf(0.4)."<br />"); // erf(0.4)=0.42839235504667
echo("erf(0.8)=".erf(0.8)."<br />"); // erf(0.8)=0.74210096470766
echo("erf(1.2)=".erf(1.2)."<br />"); // erf(1.2)=0.91031397822964
echo("erf(1.6)=".erf(1.6)."<br />"); // erf(1.6)=0.97634838334464
echo("erf(2.0)=".erf(2.0)."<br />"); // erf(2.0)=0.99532226501895

echo("erfc(-2.0)=".erfc(-2.0)."<br />"); // erfc(-2.0)=1.995322265019
echo("erfc(-1.6)=".erfc(-1.6)."<br />"); // erfc(-1.6)=1.9763483833446
echo("erfc(-1.2)=".erfc(-1.2)."<br />"); // erfc(-1.2)=1.9103139782296
echo("erfc(-0.8)=".erfc(-0.8)."<br />"); // erfc(-0.8)=1.7421009647077
echo("erfc(-0.4)=".erfc(-0.4)."<br />"); // erfc(-0.4)=1.4283923550467
echo("erfc(0.0)=".erfc(0.0)."<br />"); // erfc(0.0)=1
echo("erfc(0.4)=".erfc(0.4)."<br />"); // erfc(0.4)=0.57160764495333
echo("erfc(0.8)=".erfc(0.8)."<br />"); // erfc(0.8)=0.25789903529234
echo("erfc(1.2)=".erfc(1.2)."<br />"); // erfc(1.2)=0.089686021770365
echo("erfc(1.6)=".erfc(1.6)."<br />"); // erfc(1.6)=0.023651616655356
echo("erfc(2.0)=".erfc(2.0)."<br />"); // erfc(2.0)=0.0046777349810473

?>

Link

• Takekatsu HIRAMURA's page