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

なかだです。

At Fri, 12 Apr 2002 01:54:23 +0900,
OZAWA -Crouton- Sakuro wrote:
> > > acosh がなくて、asinhなどがあるというのは考えにくいんで、こ
> > > れら全部をひとつのファイル(missing/acos.c)に入れて、
> 
> ファイル名、acos.cでいいんですか?
> archyp.c くらいにはしといたほうがいいんじゃないかと。

たしかにそっちのほうが適切っぽいんですが、AC_REPLACE_FUNCSだと
関数名をそのままファイル名にしちゃうんですよね。こういう感じか
なぁ。

それはそうと、アルゴリズム事典を発掘して確認してみたら、思いっ
きり間違ってました。これでだいたい最低精度は12桁くらいです。


Index: Makefile.in
===================================================================
RCS file: /cvs/ruby/src/ruby/Makefile.in,v
retrieving revision 1.32
diff -u -2 -p -r1.32 Makefile.in
--- Makefile.in	2002/04/02 06:48:41	1.32
+++ Makefile.in	2002/04/11 18:13:33
@@ -169,4 +169,7 @@ alloca.@OBJEXT@: $(srcdir)/missing/alloc
 	$(CC) -I. $(CFLAGS) $(CPPFLAGS) -c $(srcdir)/missing/alloca.c
 
+archyp.@OBJEXT@: $(srcdir)/missing/archyp.c
+	$(CC) -I. $(CFLAGS) $(CPPFLAGS) -c $(srcdir)/missing/archyp.c
+
 crypt.@OBJEXT@: $(srcdir)/missing/crypt.c
 	$(CC) -I. $(CFLAGS) $(CPPFLAGS) -c $(srcdir)/missing/crypt.c
@@ -174,4 +177,7 @@ crypt.@OBJEXT@: $(srcdir)/missing/crypt.
 dup2.@OBJEXT@: $(srcdir)/missing/dup2.c
 	$(CC) -I. $(CFLAGS) $(CPPFLAGS) -c $(srcdir)/missing/dup2.c
+
+fileblocks.@OBJEXT@: $(srcdir)/missing/fileblocks.c
+	$(CC) -I. $(CFLAGS) $(CPPFLAGS) -c $(srcdir)/missing/fileblocks.c
 
 finite.@OBJEXT@: $(srcdir)/missing/finite.c
Index: configure.in
===================================================================
RCS file: /cvs/ruby/src/ruby/configure.in,v
retrieving revision 1.116
diff -u -2 -p -r1.116 configure.in
--- configure.in	2002/04/11 10:02:59	1.116
+++ configure.in	2002/04/11 17:58:03
@@ -309,5 +309,8 @@ AC_CHECK_FUNCS(fmod killpg drand48 rando
 	      getpgrp setpgrp getpgid setpgid getgroups getpriority getrlimit\
 	      dlopen sigprocmask sigaction _setjmp setsid telldir seekdir fchmod\
-	      mktime timegm cosh sinh tanh)
+	      mktime timegm cosh sinh tanh acosh asinh atanh)
+case $ac_cv_func_acosh$ac_cv_func_asinh$ac_cv_func_atanh in
+*no*) AC_LIBOBJ(archyp);;
+esac
 AC_STRUCT_TIMEZONE
 AC_CACHE_CHECK(for struct tm.tm_gmtoff, rb_cv_member_struct_tm_tm_gmtoff,
Index: missing/archyp.c
===================================================================
RCS file: archyp.c
diff -N archyp.c
--- /dev/null	Thu Jan  1 09:00:00 1970
+++ archyp.c	Fri Apr 12 03:16:10 2002
@@ -0,0 +1,88 @@
+/**********************************************************************
+
+  archyp.c -
+
+  $Author$
+  $Date$
+  created at: Fri Apr 12 03:05:29 JST 2002
+
+ public domain rewrite of acosh(3), asinh(3) and atanh(3)
+
+**********************************************************************/
+
+#include <errno.h>
+#include <float.h>
+#include <math.h>
+
+/* DBL_MANT_DIG must be less than 4 times of bits of int */
+#ifdef DBL_MANT_DIG
+#define DBL_MANT_DIG 53		/* in this case, at least 12 digit precision */
+#endif
+#define BIG_CRITERIA_BIT (1<<DBL_MANT_DIG/2)
+#if BIG_CRITERIA_BIT > 0
+#define BIG_CRITERIA (1.0*BIG_CRITERIA_BIT)
+#else
+#define BIG_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/2+1-DBL_MANT_DIG/4)))
+#endif
+#define SMALL_CRITERIA_BIT (1<<(DBL_MANT_DIG/3))
+#if SMALL_CRITERIA_BIT > 0
+#define SMALL_CRITERIA (1.0/SMALL_CRITERIA_BIT)
+#else
+#define SMALL_CRITERIA (1.0*(1<<DBL_MANT_DIG/4)*(1<<(DBL_MANT_DIG/3+1-DBL_MANT_DIG/4)))
+#endif
+
+#ifndef HAVE_ACOSH
+double
+acosh(x)
+    double x;
+{
+    if (x < 1)
+	x = -1;			/* NaN */
+    else if (x == 1)
+	return 0;
+    else if (x > BIG_CRITERIA)
+	x += x;
+    else
+	x += sqrt((x + 1) * (x - 1));
+    return log(x);
+}
+#endif
+
+#ifndef HAVE_ASINH
+double
+asinh(x)
+    double x;
+{
+    int neg = x < 0;
+    double z = fabs(x);
+
+    if (z < SMALL_CRITERIA) return x;
+    if (z < (1.0/(1<<DBL_MANT_DIG/5))) {
+	double x2 = z * z;
+	z *= 1 + x2 * (-1.0/6.0 + x2 * 3.0/40.0);
+    }
+    else if (z > BIG_CRITERIA) {
+	z = log(z + z);
+    }
+    else {
+	z = log(z + sqrt(z * z + 1));
+    }
+    if (neg) z = -z;
+    return z;
+}
+#endif
+
+#ifndef HAVE_ATANH
+double
+atanh(x)
+    double x;
+{
+    int neg = x < 0;
+    double z = fabs(x);
+
+    if (z < SMALL_CRITERIA) return x;
+    z = log(z > 1 ? -1 : (1 + z) / (1 - z)) / 2;
+    if (neg) z = -z;
+    return z;
+}
+#endif


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