Re: Fixes for the new step methods in the range.c & numeric.c

まつもと ゆきひろです

以下のようなメールがruby-talkに来たのですが、私にはどこがどう違うのか
良く分かりませんでした。もしよければ知恵を貸していただけませんか？

In message "Fixes for the new step methods in the range.c & numeric.c"
    on 02/05/06, [e-mail] <[e-mail]> writes:
|
|This is a multi-part message in MIME format.
|
|------=_NextPart_000_0001_01C1F50B.6E551680
|Content-Type: text/plain;
|	charset="us-ascii"
|Content-Transfer-Encoding: 7bit
|
|Hi,
|
|an included patch for the step method in range.c is a simple
|bug-fix - for example
|
|---
|class Class
|  alias succ superclass
|end
|
|# currently throws an exception complaining that succ for nil
|# is undefined
|(Fixnum..Object).step(2) {|x| p x }
|---
|
|The patch for the step method in numeric.c adds a sequence 
|end point adjustment strategy - #step_fuzzyness measures the
|``fuzzyness'' of step-length additions around the end-point
|and adds a safety margin (there exist probably much better 
|step_fuzzyness measures) - in pure Ruby:  
|
|---
|class Float
|    def test_step(b,s)
|      n = ((b - self)/ s).round
|      last = n*s + self
|      
|      fuzzyness = step_fuzzyness(n,s)
|      if (s > 0 && last > b + fuzzyness) || \
|         (s < 0 && last < b - fuzzyness)
|         n-=1
|         last = n*s + self
|      end
|      
|      if (last - b).abs < fuzzyness
|        last = b
|      end
|        
|      (n-1).times {|i| yield i*s + self }
|      yield last
|	self
|    end
|    
|    EPSILON = 2.0**(-43)
|    def step_fuzzyness(n, s)
|      tmp = ((n-1)..(n+1)).collect do |i| 
|       ((i*s + self - ((i-1)*s + self)) - s).abs
|      end
|      tmp.push s.abs*EPSILON
|      2*tmp.max
|    end
|end
|
|
|# with this patch on my machine ...
|M = 52343
|T = 12418.979
|K = 20300012
|L = 15673
|
|p b = (K/T)                  # => 1634.595887
|p e = (-L*M + K)/T           # => -64423.31749
|p s = -M/T                   # => -4.214758717
|p err = M/T * (2.0**(-40))# => 3.677147842e-12
|
|# with this patch on my machine
|
|last = nil
|b.step(e,s) {|last| }
|p e - last       # => 0.0
|
|b.step(e+err,s) {|last| }
|p (e+err) - last # => -4.214758717
|
|b.step(e-err,s) {|last| }
|p (e-err) - last # => -7.275957614e-12
|
|
|# a bigger error is necessary to demonstrate
|# fuzzyness in pure Ruby 
|
|p err = M/T * (2.0**(-38.00))
|
|b.test_step(e,s) {|last| }
|p e - last       # => 0.0
|
|b.test_step(e+err,s) {|last| }
|p (e+err) - last # => -4.214758717
|
|b.test_step(e-err,s) {|last| }
|p (e-err) - last # => -1.455191523e-11
|---
|
|
|/Christoph
|
|------=_NextPart_000_0001_01C1F50B.6E551680
|Content-Type: application/octet-stream;
|	name="range.patch"
|Content-Transfer-Encoding: quoted-printable
|Content-Disposition: attachment;
|	filename="range.patch"
|
|--- range.c.1.34	Fri May  3 18:44:42 2002=0A=
|+++ range.c	Fri May  3 00:23:28 2002=0A=
|@@ -71,7 +71,7 @@=0A=
|     VALUE obj;=0A=
| {=0A=
|     VALUE beg, end, flags;=0A=
|-    =0A=
|+=0A=
|     rb_scan_args(argc, argv, "21", &beg, &end, &flags);=0A=
|     /* Ranges are immutable, so that they should be initialized only =
|once. */=0A=
|     if (rb_ivar_defined(obj, id_beg)) {=0A=
|@@ -304,26 +304,21 @@=0A=
|     else {			/* generic each */=0A=
| 	VALUE v =3D b;=0A=
| 	long lim =3D NUM2LONG(step);=0A=
|-	long i;=0A=
|+	long i =3D lim;=0A=
| =0A=
| 	if (lim <=3D 0) {=0A=
| 	    rb_raise(rb_eArgError, "step can't be <=3D 0");=0A=
| 	}=0A=
|-	if (EXCL(range)) {=0A=
|-	    while (r_lt(v, e)) {=0A=
|-		if (r_eq(v, e)) break;=0A=
|-		rb_yield(v);=0A=
|-		for (i=3D0; i<lim; i++) =0A=
|-		    v =3D rb_funcall(v, id_succ, 0, 0);=0A=
|-	    }=0A=
|-	}=0A=
|-	else {=0A=
|-	    while (r_le(v, e)) {=0A=
|+	for(;NUM2INT(rb_funcall(v, id_cmp, 1, e)) < 0; v =3D rb_funcall(v, =
|id_succ, 0, 0))=0A=
|+		if (i =3D=3D lim) {=0A=
|+			rb_yield(v);=0A=
|+			i =3D 1;=0A=
|+		} else {=0A=
|+			i++;=0A=
|+		}=0A=
|+	if (!(EXCL(range)) && (i =3D=3D lim) && \=0A=
|+	    (v =3D=3D e || NUM2INT(rb_funcall(v, id_cmp, 1, e)) =3D=3D 0)) {=0A=
| 		rb_yield(v);=0A=
|-		if (r_eq(v, e)) break;=0A=
|-		for (i=3D0; i<lim; i++) =0A=
|-		    v =3D rb_funcall(v, id_succ, 0, 0);=0A=
|-	    }=0A=
| 	}=0A=
|     }=0A=
|     return range;=0A=
|
|------=_NextPart_000_0001_01C1F50B.6E551680
|Content-Type: application/octet-stream;
|	name="numeric.patch"
|Content-Transfer-Encoding: quoted-printable
|Content-Disposition: attachment;
|	filename="numeric.patch"
|
|--- numeric.c.1.47	Fri May  3 19:02:33 2002=0A=
|+++ numeric.c	Sun May  5 12:17:11 2002=0A=
|@@ -224,7 +224,7 @@=0A=
| 	return rb_str_new2(value < 0 ? "-Infinity" : "Infinity");=0A=
|     else if(isnan(value))=0A=
| 	return rb_str_new2("NaN");=0A=
|-    =0A=
|+=0A=
|     avalue =3D fabs(value);=0A=
|     if (avalue =3D=3D 0.0) {=0A=
| 	fmt =3D "%.1f";=0A=
|@@ -234,13 +234,13 @@=0A=
| 	while (d1 < 1.0) d1 *=3D 10.0;=0A=
| 	d1 =3D modf(d1, &d2);=0A=
| 	if (d1 =3D=3D 0) fmt =3D "%.1e";=0A=
|-    }    =0A=
|+    }=0A=
|     else if (avalue >=3D 1.0e10) {=0A=
| 	d1 =3D avalue;=0A=
| 	while (d1 > 10.0) d1 /=3D 10.0;=0A=
| 	d1 =3D modf(d1, &d2);=0A=
| 	if (d1 =3D=3D 0) fmt =3D "%.1e";=0A=
|-    }    =0A=
|+    }=0A=
|     else if ((d1 =3D modf(value, &d2)) =3D=3D 0) {=0A=
| 	fmt =3D "%.1f";=0A=
|     }=0A=
|@@ -641,7 +641,7 @@=0A=
| =0A=
| static VALUE flo_is_nan_p(num)=0A=
|      VALUE num;=0A=
|-{     =0A=
|+{=0A=
| =0A=
|   double value =3D RFLOAT(num)->value;=0A=
| =0A=
|@@ -650,7 +650,7 @@=0A=
| =0A=
| static VALUE flo_is_infinite_p(num)=0A=
|      VALUE num;=0A=
|-{     =0A=
|+{=0A=
|   double value =3D RFLOAT(num)->value;=0A=
| =0A=
|   if (isinf(value)) {=0A=
|@@ -662,12 +662,12 @@=0A=
| =0A=
| static VALUE flo_is_finite_p(num)=0A=
|      VALUE num;=0A=
|-{     =0A=
|+{=0A=
|   double value =3D RFLOAT(num)->value;=0A=
| =0A=
|   if (isinf(value) || isnan(value))=0A=
|     return Qfalse;=0A=
|-  =0A=
|+=0A=
|   return Qtrue;=0A=
| }=0A=
| =0A=
|@@ -762,6 +762,25 @@=0A=
|     return flo_truncate(rb_Float(num));=0A=
| }=0A=
| =0A=
|+=0A=
|+static double=0A=
|+step_fuzzyness(n,a,b,c)=0A=
|+	long  n;=0A=
|+	double a,b,c;=0A=
|+{=0A=
|+	/* eps =3D 2^(-43)  */=0A=
|+	const double eps =3D 1.136868377e-13;=0A=
|+	double n0 =3D (n-2)*c + a;=0A=
|+	double n1 =3D (n-1)*c + a;=0A=
|+	double n2 =3D (n)*c + a;=0A=
|+	double n3 =3D ((n)+1)*c + a;=0A=
|+	/* max-macro is not universal */=0A=
|+	#define fuzzy_max(a,b)    (((a) > (b)) ? (a) : (b))=0A=
|+	return 2*fuzzy_max(fuzzy_max(fabs(n1-n0-c),fabs(n2-n1-c)),\=0A=
|+	                   fuzzy_max(fabs(n3-n2-c),fabs(c)*eps));=0A=
|+}=0A=
|+=0A=
|+=0A=
| static VALUE=0A=
| num_step(argc, argv, from)=0A=
|     int argc;=0A=
|@@ -798,17 +817,32 @@=0A=
| 	}=0A=
|     }=0A=
|     else if (TYPE(from) =3D=3D T_FLOAT || TYPE(to) =3D=3D T_FLOAT || =
|TYPE(step) =3D=3D T_FLOAT) {=0A=
|-	const double epsilon =3D 2.2204460492503131E-16;=0A=
| 	double beg =3D NUM2DBL(from);=0A=
| 	double end =3D NUM2DBL(to);=0A=
| 	double unit =3D NUM2DBL(step);=0A=
|-	double n =3D (end - beg)/unit;=0A=
|-	long i;=0A=
|+	double fuzzy,last;=0A=
|+	long   n,i;=0A=
|+	n =3D floor((end - beg)/unit + 0.5);=0A=
|+=0A=
|+	if (n < 0) {=0A=
|+		return from;=0A=
|+	}=0A=
|+	fuzzy =3D step_fuzzyness(n,beg,end,unit);=0A=
|+	last =3D n* unit + beg;=0A=
| =0A=
|-	n =3D floor(n + n*epsilon) + 1;=0A=
|-	for (i=3D0; i<n; i++) {=0A=
|-	    rb_yield(rb_float_new(i*unit+beg));=0A=
|+	if ((unit > 0.0 && last > end + fuzzy) || (unit < 0.0 && last < end - =
|fuzzy)) {=0A=
|+		n--;=0A=
|+		last =3D n*unit + beg;=0A=
| 	}=0A=
|+	if (fabs(last - end) < fuzzy) {=0A=
|+		last =3D end;=0A=
|+	}=0A=
|+=0A=
|+	for (i =3D 0; i<n; i++) {=0A=
|+		rb_yield(rb_float_new(i*unit + beg));=0A=
|+	}=0A=
|+	rb_yield(rb_float_new(last));=0A=
|+=0A=
|     }=0A=
|     else {=0A=
| 	VALUE i =3D from;=0A=
|
|------=_NextPart_000_0001_01C1F50B.6E551680--
|
|