[ruby-dev:16843] Re: math.c 1.10

なかだです。

At Thu, 11 Apr 2002 15:23:40 +0900,
OZAWA -Crouton- Sakuro wrote:
> asinh(x) = log(x + √(x^2+1)) [-∞<x<∞]
> acosh(x) = log(x + √(x^2-1)) [x>=1]
>            NaN [otherwise] errno=EDOM
> atanh(x) = (1/2) * log((1+x)/(1-x)) [-1<x<1]
>            NaN [otherwise] errno=EDOM
> 
> だそうです。
> 
> x<0のときは、
> asinh(x) = -log(√(x^2+1)-x)
> として桁落ちを防ぎ、

これは、sinhは(0,0)に対して点対称なので、絶対値を取っておいて最
後に符号を戻せばいいですよね。

> x≒0のときは級数展開して
> asinh(x) = x - x^3/6 ...

Σ(-1^i*-x^(2i-1)/(2i-1)!) ?

> atanh(x) = x + x^3/3 ...

Σ(x^(2i-1)/(2i-1)) ?

> を使え、と奥村先生の本に載ってました。
> 数値計算上は、0と比較しちゃいけませんので極小のεを定義しなきゃなりませんが。

(1.0/(1<<DBL_MANT_DIG/2))くらいでいいですかね。

久しぶりだからこれでいいのかどうか…。


#include <errno.h>
#include <math.h>

#define BIG_CRITERIA (1.0*(1<<(DBL_MANT_DIG/2+2)))
#define SMALL_CRITERIA (1.0/(1<<DBL_MANT_DIG/2))

double
acosh(x)
    double x;
{
    if (x < 1)
	x = -1;			/* NaN */
    else if (x == 1)
	return 0;
    else if (x < BIG_CRITERIA)
	x += sqrt((x + 1) * (x - 1));
    else
	x += x;
    return log(x);
}

#include <errno.h>
#include <math.h>
#include <float.h>

#define BIG_CRITERIA (1.0*(1<<(DBL_MANT_DIG/2+2)))
#define SMALL_CRITERIA (1.0/(1<<DBL_MANT_DIG/2))

double
asinh(x)
    double x;
{
    int neg = x < 0;
    double z;

    if (x == 0) return x;
    z = fabs(x);
    if (z < SMALL_CRITERIA) {
	double f = 1, d = z, x2 = z * z, zn;
	int i = 1;

	while (f *= ++i, f *= ++i, (zn = z + (d = -d * x2 / f)) != z)
	    z = zn;
    }
    else {
	z = log(z + (z > BIG_CRITERIA ? z : sqrt(z * z + 1)));
    }
    if (neg) z = -z;
    return z;
}

#include <errno.h>
#include <math.h>
#include <float.h>

#define BIG_CRITERIA (1.0*(1<<(DBL_MANT_DIG/2+2)))
#define SMALL_CRITERIA (1.0/(1<<DBL_MANT_DIG/2))

double
atanh(x)
    double x;
{
    int neg = x < 0;
    double z;

    if (x == 0) return x;
    z = fabs(x);
    if (z < SMALL_CRITERIA) {
	double f = 1, d = z, x2 = z * z, zn;

	while ((zn = z + (d = d * x2 / (f += 2))) != z)
	    z = zn;
    }
    else {
	z = log(z > 1 ? 0 : (1 + z) / (1 - z)) / 2;
    }
    if (neg) z = -z;
    return z;
}


-- 
--- 僕の前にBugはない。
--- 僕の後ろにBugはできる。
    中田 伸悦