Mercurial > octave-libgccjit

comparison liboctave/numeric/lo-specfun.cc @ 18928:161ebb78ac1b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use gnulib::log and gnulib::logf functions * Faddeeva.cc, lo-mappers.cc, lo-specfun.cc: Use gnulib::log and gnulib::logf instead of log and logf.
author John W. Eaton <jwe@octave.org>
date Wed, 16 Jul 2014 19:56:22 -0400
parents 6113e0c6920b
children
comparison
equal deleted inserted replaced
18927:09eb8a2ddb02 18928:161ebb78ac1b
629 s += (t *= u*u) / (i+1); 629 s += (t *= u*u) / (i+1);
630 630
631 retval = 2 * (s + 1) * u; 631 retval = 2 * (s + 1) * u;
632 } 632 }
633 else 633 else
634 retval = log (1 + x); 634 retval = gnulib::log (1 + x);
635 635
636 return retval; 636 return retval;
637 } 637 }
638 #endif 638 #endif
639 639
688 s += (t *= u*u) / (i+1); 688 s += (t *= u*u) / (i+1);
689 689
690 retval = 2 * (s + 1) * u; 690 retval = 2 * (s + 1) * u;
691 } 691 }
692 else 692 else
693 retval = log (1 + x); 693 retval = gnulib::logf (1 + x);
694 694
695 return retval; 695 return retval;
696 } 696 }
697 #endif 697 #endif
698 698
2974 2974
2975 2975
2976 Complex rc_log1p (double x) 2976 Complex rc_log1p (double x)
2977 { 2977 {
2978 const double pi = 3.14159265358979323846; 2978 const double pi = 3.14159265358979323846;
2979 return x < -1.0 ? Complex (log (-(1.0 + x)), pi) : Complex (log1p (x)); 2979 return (x < -1.0
2980 ? Complex (gnulib::log (-(1.0 + x)), pi)
2981 : Complex (log1p (x)));
2980 } 2982 }
2981 2983
2982 FloatComplex rc_log1p (float x) 2984 FloatComplex rc_log1p (float x)
2983 { 2985 {
2984 const float pi = 3.14159265358979323846f; 2986 const float pi = 3.14159265358979323846f;
2985 return x < -1.0f ? FloatComplex (logf (-(1.0f + x)), pi) 2987 return x < -1.0f ? FloatComplex (gnulib::logf (-(1.0f + x)), pi)
2986 : FloatComplex (log1pf (x)); 2988 : FloatComplex (log1pf (x));
2987 } 2989 }
2988 2990
2989 // This algorithm is due to P. J. Acklam. 2991 // This algorithm is due to P. J. Acklam.
2990 // See http://home.online.no/~pjacklam/notes/invnorm/ 2992 // See http://home.online.no/~pjacklam/notes/invnorm/
3034 y = yn / yd; 3036 y = yn / yd;
3035 } 3037 }
3036 else if (ax < 1.0) 3038 else if (ax < 1.0)
3037 { 3039 {
3038 // Tail region. 3040 // Tail region.
3039 const double q = sqrt (-2*log (0.5*(1-ax))); 3041 const double q = sqrt (-2*gnulib::log (0.5*(1-ax)));
3040 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; 3042 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3041 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; 3043 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3042 y = yn / yd * signum (-x); 3044 y = yn / yd * signum (-x);
3043 } 3045 }
3044 else if (ax == 1.0) 3046 else if (ax == 1.0)
3112 y = yn / yd; 3114 y = yn / yd;
3113 } 3115 }
3114 else if (x > 0.0 && x < 2.0) 3116 else if (x > 0.0 && x < 2.0)
3115 { 3117 {
3116 // Tail region. 3118 // Tail region.
3117 const double q = x < 1 ? sqrt (-2*log (0.5*x)) : sqrt (-2*log (0.5*(2-x))); 3119 const double q = x < 1 ? sqrt (-2*gnulib::log (0.5*x)) : sqrt (-2*gnulib::log (0.5*(2-x)));
3118 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]; 3120 const double yn = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5];
3119 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0; 3121 const double yd = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0;
3120 y = yn / yd; 3122 y = yn / yd;
3121 if (x < pbreak_lo) 3123 if (x < pbreak_lo)
3122 y = -y; 3124 y = -y;
3231 value = value + term; 3233 value = value + term;
3232 temp = fabs (term); 3234 temp = fabs (term);
3233 3235
3234 if (temp <= acu && temp <= acu * value) 3236 if (temp <= acu && temp <= acu * value)
3235 { 3237 {
3236 value = value * exp (pp * log (xx) 3238 value = value * exp (pp * gnulib::log (xx)
3237 + (qq - 1.0) * log (cx) - beta) / pp; 3239 + (qq - 1.0) * gnulib::log (cx) - beta) / pp;
3238 3240
3239 if (indx) 3241 if (indx)
3240 { 3242 {
3241 value = 1.0 - value; 3243 value = 1.0 - value;
3242 } 3244 }
3330 indx = false; 3332 indx = false;
3331 } 3333 }
3332 3334
3333 // Calculate the initial approximation. 3335 // Calculate the initial approximation.
3334 3336
3335 r = sqrt (- log (a * a)); 3337 r = sqrt (- gnulib::log (a * a));
3336 3338
3337 ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r); 3339 ycur = r - (2.30753 + 0.27061 * r) / (1.0 + (0.99229 + 0.04481 * r) * r);
3338 3340
3339 if (1.0 < pp && 1.0 < qq) 3341 if (1.0 < pp && 1.0 < qq)
3340 { 3342 {
3351 t = 1.0 / (9.0 * qq); 3353 t = 1.0 / (9.0 * qq);
3352 t = r * pow (1.0 - t + ycur * sqrt (t), 3); 3354 t = r * pow (1.0 - t + ycur * sqrt (t), 3);
3353 3355
3354 if (t <= 0.0) 3356 if (t <= 0.0)
3355 { 3357 {
3356 value = 1.0 - exp ((log ((1.0 - a) * qq) + beta) / qq); 3358 value = 1.0 - exp ((gnulib::log ((1.0 - a) * qq) + beta) / qq);
3357 } 3359 }
3358 else 3360 else
3359 { 3361 {
3360 t = (4.0 * pp + r - 2.0) / t; 3362 t = (4.0 * pp + r - 2.0) / t;
3361 3363
3362 if (t <= 1.0) 3364 if (t <= 1.0)
3363 { 3365 {
3364 value = exp ((log (a * pp) + beta) / pp); 3366 value = exp ((gnulib::log (a * pp) + beta) / pp);
3365 } 3367 }
3366 else 3368 else
3367 { 3369 {
3368 value = 1.0 - 2.0 / (t + 1.0); 3370 value = 1.0 - 2.0 / (t + 1.0);
3369 } 3371 }
3401 { 3403 {
3402 return value; 3404 return value;
3403 } 3405 }
3404 3406
3405 xin = value; 3407 xin = value;
3406 ycur = (ycur - a) * exp (beta + r * log (xin) + t * log (1.0 - xin)); 3408 ycur = (ycur - a) * exp (beta + r * gnulib::log (xin)
3409 + t * gnulib::log (1.0 - xin));
3407 3410
3408 if (ycur * yprev <= 0.0) 3411 if (ycur * yprev <= 0.0)
3409 { 3412 {
3410 prev = std::max (sq, fpu); 3413 prev = std::max (sq, fpu);
3411 } 3414 }