Changeset 1889

Timestamp:

03/13/07 03:23:59 (2 years ago)

Author:

Don Clugston

Message:

Minor refactoring: added a trivial function rationalPoly() which reduces code size in many of the math modules, and creates an opportunity for optimisation.

Files:

trunk/tango/math/Bessel.d (modified) (6 diffs)
trunk/tango/math/Core.d (modified) (1 diff)
trunk/tango/math/ErrorFunction.d (modified) (5 diffs)
trunk/tango/math/GammaFunction.d (modified) (1 diff)

Legend:

: Unmodified
: Added
: Removed
: Modified
: Copied
: Moved

trunk/tango/math/Bessel.d

r1435	r1889
132	132	const real [] JZ = [5.783185962946784521176L, 30.47126234366208639908L, 7.488700679069518344489e1L];
133	133	y = (xx - JZ[0]) * (xx - JZ[1]) * (xx - JZ[2]);
134		y *= ~~poly( xx, j0n) / poly( xx~~, j0d);
	134	y *= rationalPoly( xx, j0n, j0d);
135	135	return y;
136	136	}
…	…
138	138	y = fabs(x);
139	139	xx = 1.0/xx;
140		phase = ~~poly( xx, j0phasen) / poly( xx~~, j0phased);
	140	phase = rationalPoly( xx, j0phasen, j0phased);
141	141
142	142	z = 1.0/y;
143		modulus = ~~poly( z, j0modulusn) / poly( z~~, j0modulusd);
	143	modulus = rationalPoly( z, j0modulusn, j0modulusd);
144	144
145	145	y = modulus * cos( y - PI_4 + z*phase) / sqrt(y);
…	…
246	246	if ( xx < 20.25L ) {
247	247	y = M_2_PI * log(x) * cylBessel_j0(x);
248		y += ~~poly( xx, y0n) / poly( xx~~, y0d);
	248	y += rationalPoly( xx, y0n, y0d);
249	249	} else {
250	250	const real [] Y0Z = [3.957678419314857868376e0L, 7.086051060301772697624e0L,
251	251	1.022234504349641701900e1L, 1.336109747387276347827e1L];
252	252	y = (x - Y0Z[0])(x - Y0Z[1])(x - Y0Z[2])*(x - Y0Z[3]);
253		y *= ~~poly( x, y059n) / poly( x~~, y059d);
	253	y *= rationalPoly( x, y059n, y059d);
254	254	}
255	255	return y;
…	…
258	258	y = fabs(x);
259	259	xx = 1.0/xx;
260		phase = ~~poly( xx, j0phasen) / poly( xx~~, j0phased);
	260	phase = rationalPoly( xx, j0phasen, j0phased);
261	261
262	262	z = 1.0/y;
263		modulus = ~~poly( z, j0modulusn) / poly( z~~, j0modulusd);
	263	modulus = rationalPoly( z, j0modulusn, j0modulusd);
264	264
265	265	y = modulus * sin( y - PI_4 + z*phase) / sqrt(y);
…	…
458	458	y = fabs(x);
459	459	xx = 1.0/xx;
460		phase = ~~poly( xx, j1phasen) / poly( xx~~, j1phased);
	460	phase = rationalPoly( xx, j1phasen, j1phased);
461	461	z = 1.0/y;
462		modulus = ~~poly( z, j1modulusn) / poly( z~~, j1modulusd);
	462	modulus = rationalPoly( z, j1modulusn, j1modulusd);
463	463
464	464	const real M_3PI_4 = 3 * PI_4;
…	…
584	584	8.59600586833116892643e0L, 1.17491548308398812434e1L];
585	585	y = (x - Y1Z[0])(x - Y1Z[1])(x - Y1Z[2])*(x - Y1Z[3]);
586		y *= ~~poly( x, y159n) / poly( x~~, y159d);
	586	y *= rationalPoly( x, y159n, y159d);
587	587	}
588	588	return y;
589	589	}
590	590	xx = 1.0/xx;
591		phase = ~~poly( xx, j1phasen) / poly( xx~~, j1phased);
592		modulus = ~~poly( z, j1modulusn) / poly( z~~, j1modulusd);
	591	phase = rationalPoly( xx, j1phasen, j1phased);
	592	modulus = rationalPoly( z, j1modulusn, j1modulusd);
593	593
594	594	const real M_3PI_4 = 3 * PI_4;

trunk/tango/math/Core.d

r1587	r1889
1677	1677	}
1678	1678
	1679	package {
	1680	real rationalPoly(real x, real [] numerator, real [] denominator)
	1681	{
	1682	return poly(x, numerator)/poly(x, denominator);
	1683	}
	1684	}
	1685
1679	1686	private enum : int { MANTDIG_2 = real.mant_dig/2 } // Compiler workaround
1680	1687

trunk/tango/math/ErrorFunction.d

r1528	r1889
167	167	real p, q;
168	168
169		if( x < 8.0 ) {
170		p = poly( y, P);
171		q = poly( y, Q);
172		} else {
173		q = y * y;
174		p = y * poly( q, R);
175		q = poly( q, S);
176		}
177		y = (z * p)/q;
	169	if( x < 8.0 ) y = z * rationalPoly(y, P, Q);
	170	else y = z * y * rationalPoly(y * y, R, S);
178	171
179	172	if (a < 0.0L)
…	…
203	196
204	197	if (x < 8.0) {
205		p = poly( y, P);
206		q = poly( y, Q);
	198	return rationalPoly( y, P, Q);
207	199	} else {
208		q = y * y;
209		p = y * poly( q, R);
210		q = poly( q, S);
211		}
212		return p/q;
	200	return y * rationalPoly(y*y, R, S);
	201	}
213	202	}
214	203
…	…
246	235
247	236	real z = x * x;
248		return x * ~~poly(z, T) / poly(z~~, U);
	237	return x * rationalPoly(z, T, U);
249	238	}
250	239
…	…
523	512	y = y - 0.5L;
524	513	y2 = y * y;
525		x = y + y * (y2 * ~~poly( y2, P0)/poly( y2~~, Q0));
	514	x = y + y * (y2 * rationalPoly( y2, P0, Q0));
526	515	return x * SQRT2PI;
527	516	}
…	…
531	520	z = 1.0L/x;
532	521	if ( x < 8.0L ) {
533		x1 = z * ~~poly( z, P1)/poly( z~~, Q1);
	522	x1 = z * rationalPoly( z, P1, Q1);
534	523	} else if( x < 32.0L ) {
535		x1 = z * ~~poly( z, P2)/poly( z~~, Q2);
	524	x1 = z * rationalPoly( z, P2, Q2);
536	525	} else {
537		x1 = z * ~~poly( z, P3)/poly( z~~, Q3);
	526	x1 = z * rationalPoly( z, P3, Q3);
538	527	}
539	528	x = x0 - x1;

trunk/tango/math/GammaFunction.d

r1528	r1889
429	429	return log(z);
430	430	x = (nx - 2.0L) + f;
431		real p = x * ~~poly( x, logGammaNumerator ) / poly( x~~, logGammaDenominator);
	431	real p = x * rationalPoly( x, logGammaNumerator, logGammaDenominator);
432	432	return log(z) + p;
433	433	}

Navigation