math.d

Revision 2293, 126.1 kB (checked in by andrei, 14 years ago)
unconditionalized unittest for round
Property svn:eol-style set to `native`

Line
1	// Written in the D programming language.
2
3	/**
4	* Elementary mathematical functions
5	*
6	* Contains the elementary mathematical functions (powers, roots,
7	* and trignometric functions), and low-level floating-point operations.
8	* Mathematical special functions are available in std.mathspecial.
9	*
10	* The functionality closely follows the IEEE754-2008 standard for
11	* floating-point arithmetic, including the use of camelCase names rather
12	* than C99-style lower case names. All of these functions behave correctly
13	* when presented with an infinity or NaN.
14	*
15	* Unlike C, there is no global 'errno' variable. Consequently, almost all of
16	* these functions are pure nothrow.
17	*
18	* Status:
19	* The gamma and error functions have been superceded by improved versions in
20	* std.mathspecial. They will be officially deprecated in std.math in DMD2.055.
21	* The semantics and names of feqrel and approxEqual will be revised.
22	*
23	* Source: $(PHOBOSSRC std/_math.d)
24	* Macros:
25	* WIKI = Phobos/StdMath
26	*
27	* TABLE_SV = <table border=1 cellpadding=4 cellspacing=0>
28	* <caption>Special Values</caption>
29	* $0</table>
30	* SVH = $(TR $(TH $1) $(TH $2))
31	* SV = $(TR $(TD $1) $(TD $2))
32	*
33	* NAN = $(RED NAN)
34	* SUP = <span style="vertical-align:super;font-size:smaller">$0</span>
35	* GAMMA = Γ
36	* THETA = θ
37	* INTEGRAL = ∫
38	* INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>)
39	* POWER = $1<sup>$2</sup>
40	* SUB = $1<sub>$2</sub>
41	* BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>)
42	* CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG ))
43	* PLUSMN = ±
44	* INFIN = ∞
45	* PLUSMNINF = ±∞
46	* PI = π
47	* LT = <
48	* GT = >
49	* SQRT = √
50	* HALF = ½
51	*
52	* Copyright: Copyright Digital Mars 2000 - 2009.
53	* License: <a href="http://www.boost.org/LICENSE_1_0.txt">Boost License 1.0</a>.
54	* Authors: $(WEB digitalmars.com, Walter Bright),
55	* Don Clugston
56	*/
57	module std.math;
58
59	import core.stdc.math;
60	import std.range, std.traits;
61
62	version(unittest) {
63	import std.typetuple;
64	}
65
66	version(LDC) {
67	import ldc.intrinsics;
68	}
69
70	version(DigitalMars){
71	version = INLINE_YL2X; // x87 has opcodes for these
72	}
73
74	version (X86){
75	version = X86_Any;
76	}
77
78	version (X86_64){
79	version = X86_Any;
80	}
81
82	version(D_InlineAsm_X86){
83	version = InlineAsm_X86_Any;
84	}
85	else version(D_InlineAsm_X86_64){
86	version = InlineAsm_X86_Any;
87	}
88
89
90
91
92	private:
93	/*
94	* The following IEEE 'real' formats are currently supported:
95	* 64 bit Big-endian 'double' (eg PowerPC)
96	* 128 bit Big-endian 'quadruple' (eg SPARC)
97	* 64 bit Little-endian 'double' (eg x86-SSE2)
98	* 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium).
99	* 128 bit Little-endian 'quadruple' (not implemented on any known processor!)
100	*
101	* Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support
102	*/
103	version(LittleEndian) {
104	static assert(real.mant_dig == 53 \|\| real.mant_dig==64
105	\|\| real.mant_dig == 113,
106	"Only 64-bit, 80-bit, and 128-bit reals"
107	" are supported for LittleEndian CPUs");
108	} else {
109	static assert(real.mant_dig == 53 \|\| real.mant_dig==106
110	\|\| real.mant_dig == 113,
111	"Only 64-bit and 128-bit reals are supported for BigEndian CPUs."
112	" double-double reals have partial support");
113	}
114
115	// Constants used for extracting the components of the representation.
116	// They supplement the built-in floating point properties.
117	template floatTraits(T) {
118	// EXPMASK is a ushort mask to select the exponent portion (without sign)
119	// EXPPOS_SHORT is the index of the exponent when represented as a ushort array.
120	// SIGNPOS_BYTE is the index of the sign when represented as a ubyte array.
121	// RECIP_EPSILON is the value such that (smallest_denormal) * RECIP_EPSILON == T.min_normal
122	enum T RECIP_EPSILON = (1/T.epsilon);
123	static if (T.mant_dig == 24)
124	{ // float
125	enum ushort EXPMASK = 0x7F80;
126	enum ushort EXPBIAS = 0x3F00;
127	enum uint EXPMASK_INT = 0x7F80_0000;
128	enum uint MANTISSAMASK_INT = 0x007F_FFFF;
129	version(LittleEndian) {
130	enum EXPPOS_SHORT = 1;
131	} else {
132	enum EXPPOS_SHORT = 0;
133	}
134	}
135	else static if (T.mant_dig == 53) // double, or real==double
136	{
137	enum ushort EXPMASK = 0x7FF0;
138	enum ushort EXPBIAS = 0x3FE0;
139	enum uint EXPMASK_INT = 0x7FF0_0000;
140	enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only
141	version(LittleEndian) {
142	enum EXPPOS_SHORT = 3;
143	enum SIGNPOS_BYTE = 7;
144	} else {
145	enum EXPPOS_SHORT = 0;
146	enum SIGNPOS_BYTE = 0;
147	}
148	}
149	else static if (T.mant_dig == 64) // real80
150	{
151	enum ushort EXPMASK = 0x7FFF;
152	enum ushort EXPBIAS = 0x3FFE;
153	version(LittleEndian)
154	{
155	enum EXPPOS_SHORT = 4;
156	enum SIGNPOS_BYTE = 9;
157	}
158	else
159	{
160	enum EXPPOS_SHORT = 0;
161	enum SIGNPOS_BYTE = 0;
162	}
163	} else static if (T.mant_dig == 113){ // quadruple
164	enum ushort EXPMASK = 0x7FFF;
165	version(LittleEndian)
166	{
167	enum EXPPOS_SHORT = 7;
168	enum SIGNPOS_BYTE = 15;
169	}
170	else
171	{
172	enum EXPPOS_SHORT = 0;
173	enum SIGNPOS_BYTE = 0;
174	}
175	} else static if (T.mant_dig == 106) { // doubledouble
176	enum ushort EXPMASK = 0x7FF0;
177	// the exponent byte is not unique
178	version(LittleEndian)
179	{
180	enum EXPPOS_SHORT = 7; // [3] is also an exp short
181	enum SIGNPOS_BYTE = 15;
182	}
183	else
184	{
185	enum EXPPOS_SHORT = 0; // [4] is also an exp short
186	enum SIGNPOS_BYTE = 0;
187	}
188	}
189	}
190
191	// These apply to all floating-point types
192	version(LittleEndian)
193	{
194	enum MANTISSA_LSB = 0;
195	enum MANTISSA_MSB = 1;
196	}
197	else
198	{
199	enum MANTISSA_LSB = 1;
200	enum MANTISSA_MSB = 0;
201	}
202
203	public:
204
205	enum real E = 2.7182818284590452354L; /** e */ // 0x1.5BF0A8B1_45769535_5FF5p+1L
206	enum real LOG2T = 0x1.a934f0979a3715fcp+1; /** $(SUB log, 2)10 */ // 3.32193 fldl2t
207	enum real LOG2E = 0x1.71547652b82fe178p+0; /** $(SUB log, 2)e */ // 1.4427 fldl2e
208	enum real LOG2 = 0x1.34413509f79fef32p-2; /** $(SUB log, 10)2 */ // 0.30103 fldlg2
209	enum real LOG10E = 0.43429448190325182765; /** $(SUB log, 10)e */
210	enum real LN2 = 0x1.62e42fefa39ef358p-1; /** ln 2 */ // 0.693147 fldln2
211	enum real LN10 = 2.30258509299404568402; /** ln 10 */
212	enum real PI = 0x1.921fb54442d1846ap+1; /** $(_PI) */ // 3.14159 fldpi
213	enum real PI_2 = 1.57079632679489661923; /** $(PI) / 2 */
214	enum real PI_4 = 0.78539816339744830962; /** $(PI) / 4 */
215	enum real M_1_PI = 0.31830988618379067154; /** 1 / $(PI) */
216	enum real M_2_PI = 0.63661977236758134308; /** 2 / $(PI) */
217	enum real M_2_SQRTPI = 1.12837916709551257390; /** 2 / $(SQRT)$(PI) */
218	enum real SQRT2 = 1.41421356237309504880; /** $(SQRT)2 */
219	enum real SQRT1_2 = 0.70710678118654752440; /** $(SQRT)$(HALF) */
220
221	/*
222	Octal versions:
223	PI/64800 0.00001 45530 36176 77347 02143 15351 61441 26767
224	PI/180 0.01073 72152 11224 72344 25603 54276 63351 22056
225	PI/8 0.31103 75524 21026 43021 51423 06305 05600 67016
226	SQRT(1/PI) 0.44067 27240 41233 33210 65616 51051 77327 77303
227	2/PI 0.50574 60333 44710 40522 47741 16537 21752 32335
228	PI/4 0.62207 73250 42055 06043 23046 14612 13401 56034
229	SQRT(2/PI) 0.63041 05147 52066 24106 41762 63612 00272 56161
230
231	PI 3.11037 55242 10264 30215 14230 63050 56006 70163
232	LOG2 0.23210 11520 47674 77674 61076 11263 26013 37111
233	*/
234
235	/***********************************
236	* Calculates the absolute value
237	*
238	* For complex numbers, abs(z) = sqrt( $(POWER z.re, 2) + $(POWER z.im, 2) )
239	* = hypot(z.re, z.im).
240	*/
241	Num abs(Num)(Num x) @safe pure nothrow
242	if (is(typeof(Num.init >= 0)) && is(typeof(-Num.init)) &&
243	!(is(Num* : const(ifloat)) \|\| is(Num : const(idouble*))
244	\|\| is(Num* : const(ireal*))))
245	{
246	static if (isFloatingPoint!(Num))
247	return fabs(x);
248	else
249	return x>=0 ? x : -x;
250	}
251
252	auto abs(Num)(Num z) @safe pure nothrow
253	if (is(Num* : const(cfloat)) \|\| is(Num : const(cdouble*))
254	\|\| is(Num* : const(creal*)))
255	{
256	return hypot(z.re, z.im);
257	}
258
259	/** ditto */
260	real abs(Num)(Num y) @safe pure nothrow
261	if (is(Num* : const(ifloat)) \|\| is(Num : const(idouble*))
262	\|\| is(Num* : const(ireal*)))
263	{
264	return fabs(y.im);
265	}
266
267
268	unittest
269	{
270	assert(isIdentical(abs(-0.0L), 0.0L));
271	assert(isNaN(abs(real.nan)));
272	assert(abs(-real.infinity) == real.infinity);
273	assert(abs(-3.2Li) == 3.2L);
274	assert(abs(71.6Li) == 71.6L);
275	assert(abs(-56) == 56);
276	assert(abs(2321312L) == 2321312L);
277	assert(abs(-1+1i) == sqrt(2.0));
278	}
279
280	/***********************************
281	* Complex conjugate
282	*
283	* conj(x + iy) = x - iy
284	*
285	* Note that z * conj(z) = $(POWER z.re, 2) - $(POWER z.im, 2)
286	* is always a real number
287	*/
288	creal conj(creal z) @safe pure nothrow
289	{
290	return z.re - z.im*1i;
291	}
292
293	/** ditto */
294	ireal conj(ireal y) @safe pure nothrow
295	{
296	return -y;
297	}
298
299
300	unittest
301	{
302	assert(conj(7 + 3i) == 7-3i);
303	ireal z = -3.2Li;
304	assert(conj(z) == -z);
305	}
306
307	/***********************************
308	* Returns cosine of x. x is in radians.
309	*
310	* $(TABLE_SV
311	* $(TR $(TH x) $(TH cos(x)) $(TH invalid?))
312	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) )
313	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes) )
314	* )
315	* Bugs:
316	* Results are undefined if \|x\| >= $(POWER 2,64).
317	*/
318
319	real cos(real x) @safe pure nothrow; /* intrinsic */
320
321	/***********************************
322	* Returns sine of x. x is in radians.
323	*
324	* $(TABLE_SV
325	* $(TR $(TH x) $(TH sin(x)) $(TH invalid?))
326	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
327	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
328	* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes))
329	* )
330	* Bugs:
331	* Results are undefined if \|x\| >= $(POWER 2,64).
332	*/
333
334	real sin(real x) @safe pure nothrow; /* intrinsic */
335
336
337	/***********************************
338	* sine, complex and imaginary
339	*
340	* sin(z) = sin(z.re)cosh(z.im) + cos(z.re)sinh(z.im)i
341	*
342	* If both sin($(THETA)) and cos($(THETA)) are required,
343	* it is most efficient to use expi($(THETA)).
344	*/
345	creal sin(creal z) @safe pure nothrow
346	{
347	creal cs = expi(z.re);
348	creal csh = coshisinh(z.im);
349	return cs.im * csh.re + cs.re * csh.im * 1i;
350	}
351
352	/** ditto */
353	ireal sin(ireal y) @safe pure nothrow
354	{
355	return cosh(y.im)*1i;
356	}
357
358	unittest
359	{
360	assert(sin(0.0+0.0i) == 0.0);
361	assert(sin(2.0+0.0i) == sin(2.0L) );
362	}
363
364	/***********************************
365	* cosine, complex and imaginary
366	*
367	* cos(z) = cos(z.re)cosh(z.im) - sin(z.re)sinh(z.im)i
368	*/
369	creal cos(creal z) @safe pure nothrow
370	{
371	creal cs = expi(z.re);
372	creal csh = coshisinh(z.im);
373	return cs.re * csh.re - cs.im * csh.im * 1i;
374	}
375
376	/** ditto */
377	real cos(ireal y) @safe pure nothrow
378	{
379	return cosh(y.im);
380	}
381
382	unittest{
383	assert(cos(0.0+0.0i)==1.0);
384	assert(cos(1.3L+0.0i)==cos(1.3L));
385	assert(cos(5.2Li)== cosh(5.2L));
386	}
387
388	/****************************************************************************
389	* Returns tangent of x. x is in radians.
390	*
391	* $(TABLE_SV
392	* $(TR $(TH x) $(TH tan(x)) $(TH invalid?))
393	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
394	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
395	* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD yes))
396	* )
397	*/
398
399	real tan(real x) @trusted pure nothrow
400	{
401	version(D_InlineAsm_X86) {
402	asm
403	{
404	fld x[EBP] ; // load theta
405	}
406	} else version(D_InlineAsm_X86_64) {
407	asm {
408	fld x[RBP] ; // load theta
409	}
410	}
411	version(InlineAsm_X86_Any) {
412	asm
413	{
414	fxam ; // test for oddball values
415	fstsw AX ;
416	sahf ;
417	jc trigerr ; // x is NAN, infinity, or empty
418	// 387's can handle denormals
419	SC18: fptan ;
420	fstp ST(0) ; // dump X, which is always 1
421	fstsw AX ;
422	sahf ;
423	jnp Lret ; // C2 = 1 (x is out of range)
424
425	// Do argument reduction to bring x into range
426	fldpi ;
427	fxch ;
428	SC17: fprem1 ;
429	fstsw AX ;
430	sahf ;
431	jp SC17 ;
432	fstp ST(1) ; // remove pi from stack
433	jmp SC18 ;
434
435	trigerr:
436	jnp Lret ; // if theta is NAN, return theta
437	fstp ST(0) ; // dump theta
438	}
439	return real.nan;
440
441	Lret:
442	;
443	} else {
444	return core.stdc.math.tanl(x);
445	}
446	}
447
448	unittest
449	{
450	static real vals[][2] = // angle,tan
451	[
452	[ 0, 0],
453	[ .5, .5463024898],
454	[ 1, 1.557407725],
455	[ 1.5, 14.10141995],
456	[ 2, -2.185039863],
457	[ 2.5,-.7470222972],
458	[ 3, -.1425465431],
459	[ 3.5, .3745856402],
460	[ 4, 1.157821282],
461	[ 4.5, 4.637332055],
462	[ 5, -3.380515006],
463	[ 5.5,-.9955840522],
464	[ 6, -.2910061914],
465	[ 6.5, .2202772003],
466	[ 10, .6483608275],
467
468	// special angles
469	[ PI_4, 1],
470	//[ PI_2, real.infinity], // PI_2 is not _exactly_ pi/2.
471	[ 3*PI_4, -1],
472	[ PI, 0],
473	[ 5*PI_4, 1],
474	//[ 3*PI_2, -real.infinity],
475	[ 7*PI_4, -1],
476	[ 2*PI, 0],
477	];
478	int i;
479
480	for (i = 0; i < vals.length; i++)
481	{
482	real x = vals[i][0];
483	real r = vals[i][1];
484	real t = tan(x);
485
486	//printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
487	if (!isIdentical(r, t)) assert(fabs(r-t) <= .0000001);
488
489	x = -x;
490	r = -r;
491	t = tan(x);
492	//printf("tan(%Lg) = %Lg, should be %Lg\n", x, t, r);
493	if (!isIdentical(r, t) && !(r!<>=0 && t!<>=0)) assert(fabs(r-t) <= .0000001);
494	}
495	// overflow
496	assert(isNaN(tan(real.infinity)));
497	assert(isNaN(tan(-real.infinity)));
498	// NaN propagation
499	assert(isIdentical( tan(NaN(0x0123L)), NaN(0x0123L) ));
500	}
501
502	/***************
503	* Calculates the arc cosine of x,
504	* returning a value ranging from 0 to $(PI).
505	*
506	* $(TABLE_SV
507	* $(TR $(TH x) $(TH acos(x)) $(TH invalid?))
508	* $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes))
509	* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes))
510	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
511	* )
512	*/
513	real acos(real x) @safe pure nothrow
514	{
515	return atan2(sqrt(1-x*x), x);
516	}
517
518	/// ditto
519	double acos(double x) @safe pure nothrow { return acos(cast(real)x); }
520	/// ditto
521	float acos(float x) @safe pure nothrow { return acos(cast(real)x); }
522
523	/***************
524	* Calculates the arc sine of x,
525	* returning a value ranging from -$(PI)/2 to $(PI)/2.
526	*
527	* $(TABLE_SV
528	* $(TR $(TH x) $(TH asin(x)) $(TH invalid?))
529	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
530	* $(TR $(TD $(GT)1.0) $(TD $(NAN)) $(TD yes))
531	* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD yes))
532	* )
533	*/
534	real asin(real x) @safe pure nothrow
535	{
536	return atan2(x, sqrt(1-x*x));
537	}
538	/// ditto
539	double asin(double x) @safe pure nothrow { return asin(cast(real)x); }
540	/// ditto
541	float asin(float x) @safe pure nothrow { return asin(cast(real)x); }
542
543	/***************
544	* Calculates the arc tangent of x,
545	* returning a value ranging from -$(PI)/2 to $(PI)/2.
546	*
547	* $(TABLE_SV
548	* $(TR $(TH x) $(TH atan(x)) $(TH invalid?))
549	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
550	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN)) $(TD yes))
551	* )
552	*/
553	real atan(real x) @safe pure nothrow { return atan2(x, 1.0L); }
554	/// ditto
555	double atan(double x) @safe pure nothrow { return atan(cast(real)x); }
556	/// ditto
557	float atan(float x) @safe pure nothrow { return atan(cast(real)x); }
558
559	/***************
560	* Calculates the arc tangent of y / x,
561	* returning a value ranging from -$(PI) to $(PI).
562	*
563	* $(TABLE_SV
564	* $(TR $(TH y) $(TH x) $(TH atan(y, x)))
565	* $(TR $(TD $(NAN)) $(TD anything) $(TD $(NAN)) )
566	* $(TR $(TD anything) $(TD $(NAN)) $(TD $(NAN)) )
567	* $(TR $(TD $(PLUSMN)0.0) $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) )
568	* $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) $(TD $(PLUSMN)0.0) )
569	* $(TR $(TD $(PLUSMN)0.0) $(TD $(LT)0.0) $(TD $(PLUSMN)$(PI)))
570	* $(TR $(TD $(PLUSMN)0.0) $(TD -0.0) $(TD $(PLUSMN)$(PI)))
571	* $(TR $(TD $(GT)0.0) $(TD $(PLUSMN)0.0) $(TD $(PI)/2) )
572	* $(TR $(TD $(LT)0.0) $(TD $(PLUSMN)0.0) $(TD -$(PI)/2) )
573	* $(TR $(TD $(GT)0.0) $(TD $(INFIN)) $(TD $(PLUSMN)0.0) )
574	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD anything) $(TD $(PLUSMN)$(PI)/2))
575	* $(TR $(TD $(GT)0.0) $(TD -$(INFIN)) $(TD $(PLUSMN)$(PI)) )
576	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(INFIN)) $(TD $(PLUSMN)$(PI)/4))
577	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD -$(INFIN)) $(TD $(PLUSMN)3$(PI)/4))
578	* )
579	*/
580	real atan2(real y, real x) @trusted pure nothrow
581	{
582	version(InlineAsm_X86_Any)
583	{
584	asm {
585	fld y;
586	fld x;
587	fpatan;
588	}
589	}
590	else
591	{
592	return core.stdc.math.atan2l(y,x);
593	}
594	}
595
596	/// ditto
597	double atan2(double y, double x) @safe pure nothrow
598	{
599	return atan2(cast(real)y, cast(real)x);
600	}
601
602	/// ditto
603	float atan2(float y, float x) @safe pure nothrow
604	{
605	return atan2(cast(real)y, cast(real)x);
606	}
607
608	/***********************************
609	* Calculates the hyperbolic cosine of x.
610	*
611	* $(TABLE_SV
612	* $(TR $(TH x) $(TH cosh(x)) $(TH invalid?))
613	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)0.0) $(TD no) )
614	* )
615	*/
616	real cosh(real x) @safe pure nothrow
617	{
618	// cosh = (exp(x)+exp(-x))/2.
619	// The naive implementation works correctly.
620	real y = exp(x);
621	return (y + 1.0/y) * 0.5;
622	}
623	/// ditto
624	double cosh(double x) @safe pure nothrow { return cosh(cast(real)x); }
625	/// ditto
626	float cosh(float x) @safe pure nothrow { return cosh(cast(real)x); }
627
628
629	/***********************************
630	* Calculates the hyperbolic sine of x.
631	*
632	* $(TABLE_SV
633	* $(TR $(TH x) $(TH sinh(x)) $(TH invalid?))
634	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no))
635	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no))
636	* )
637	*/
638	real sinh(real x) @safe pure nothrow
639	{
640	// sinh(x) = (exp(x)-exp(-x))/2;
641	// Very large arguments could cause an overflow, but
642	// the maximum value of x for which exp(x) + exp(-x)) != exp(x)
643	// is x = 0.5 * (real.mant_dig) * LN2. // = 22.1807 for real80.
644	if (fabs(x) > real.mant_dig * LN2) {
645	return copysign(0.5 * exp(fabs(x)), x);
646	}
647	real y = expm1(x);
648	return 0.5 * y / (y+1) * (y+2);
649	}
650	/// ditto
651	double sinh(double x) @safe pure nothrow { return sinh(cast(real)x); }
652	/// ditto
653	float sinh(float x) @safe pure nothrow { return sinh(cast(real)x); }
654
655
656	/***********************************
657	* Calculates the hyperbolic tangent of x.
658	*
659	* $(TABLE_SV
660	* $(TR $(TH x) $(TH tanh(x)) $(TH invalid?))
661	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) )
662	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)1.0) $(TD no))
663	* )
664	*/
665	real tanh(real x) @safe pure nothrow
666	{
667	// tanh(x) = (exp(x) - exp(-x))/(exp(x)+exp(-x))
668	if (fabs(x) > real.mant_dig * LN2) {
669	return copysign(1, x);
670	}
671	real y = expm1(2*x);
672	return y / (y + 2);
673	}
674	/// ditto
675	double tanh(double x) @safe pure nothrow { return tanh(cast(real)x); }
676	/// ditto
677	float tanh(float x) @safe pure nothrow { return tanh(cast(real)x); }
678
679	private:
680	/* Returns cosh(x) + I * sinh(x)
681	* Only one call to exp() is performed.
682	*/
683	creal coshisinh(real x) @safe pure nothrow
684	{
685	// See comments for cosh, sinh.
686	if (fabs(x) > real.mant_dig * LN2) {
687	real y = exp(fabs(x));
688	return y * 0.5 + 0.5i * copysign(y, x);
689	} else {
690	real y = expm1(x);
691	return (y + 1.0 + 1.0/(y + 1.0)) * 0.5 + 0.5i * y / (y+1) * (y+2);
692	}
693	}
694
695	unittest {
696	creal c = coshisinh(3.0L);
697	assert(c.re == cosh(3.0L));
698	assert(c.im == sinh(3.0L));
699	}
700
701	public:
702
703	/***********************************
704	* Calculates the inverse hyperbolic cosine of x.
705	*
706	* Mathematically, acosh(x) = log(x + sqrt( x*x - 1))
707	*
708	* $(TABLE_DOMRG
709	* $(DOMAIN 1..$(INFIN))
710	* $(RANGE 1..log(real.max), $(INFIN)) )
711	* $(TABLE_SV
712	* $(SVH x, acosh(x) )
713	* $(SV $(NAN), $(NAN) )
714	* $(SV $(LT)1, $(NAN) )
715	* $(SV 1, 0 )
716	* $(SV +$(INFIN),+$(INFIN))
717	* )
718	*/
719	real acosh(real x) @safe pure nothrow
720	{
721	if (x > 1/real.epsilon)
722	return LN2 + log(x);
723	else
724	return log(x + sqrt(x*x - 1));
725	}
726	/// ditto
727	double acosh(double x) @safe pure nothrow { return acosh(cast(real)x); }
728	/// ditto
729	float acosh(float x) @safe pure nothrow { return acosh(cast(real)x); }
730
731
732	unittest
733	{
734	assert(isNaN(acosh(0.9)));
735	assert(isNaN(acosh(real.nan)));
736	assert(acosh(1.0)==0.0);
737	assert(acosh(real.infinity) == real.infinity);
738	}
739
740	/***********************************
741	* Calculates the inverse hyperbolic sine of x.
742	*
743	* Mathematically,
744	* ---------------
745	* asinh(x) = log( x + sqrt( x*x + 1 )) // if x >= +0
746	* asinh(x) = -log(-x + sqrt( x*x + 1 )) // if x <= -0
747	* -------------
748	*
749	* $(TABLE_SV
750	* $(SVH x, asinh(x) )
751	* $(SV $(NAN), $(NAN) )
752	* $(SV $(PLUSMN)0, $(PLUSMN)0 )
753	* $(SV $(PLUSMN)$(INFIN),$(PLUSMN)$(INFIN))
754	* )
755	*/
756	real asinh(real x) @safe pure nothrow
757	{
758	return (fabs(x) > 1 / real.epsilon)
759	// beyond this point, xx + 1 == xx
760	? copysign(LN2 + log(fabs(x)), x)
761	// sqrt(xx + 1) == 1 + x x / ( 1 + sqrt(x*x + 1) )
762	: copysign(log1p(fabs(x) + xx / (1 + sqrt(xx + 1)) ), x);
763	}
764	/// ditto
765	double asinh(double x) @safe pure nothrow { return asinh(cast(real)x); }
766	/// ditto
767	float asinh(float x) @safe pure nothrow { return asinh(cast(real)x); }
768
769	unittest
770	{
771	assert(isIdentical(asinh(0.0), 0.0));
772	assert(isIdentical(asinh(-0.0), -0.0));
773	assert(asinh(real.infinity) == real.infinity);
774	assert(asinh(-real.infinity) == -real.infinity);
775	assert(isNaN(asinh(real.nan)));
776	}
777
778	/***********************************
779	* Calculates the inverse hyperbolic tangent of x,
780	* returning a value from ranging from -1 to 1.
781	*
782	* Mathematically, atanh(x) = log( (1+x)/(1-x) ) / 2
783	*
784	*
785	* $(TABLE_DOMRG
786	* $(DOMAIN -$(INFIN)..$(INFIN))
787	* $(RANGE -1..1) )
788	* $(TABLE_SV
789	* $(SVH x, acosh(x) )
790	* $(SV $(NAN), $(NAN) )
791	* $(SV $(PLUSMN)0, $(PLUSMN)0)
792	* $(SV -$(INFIN), -0)
793	* )
794	*/
795	real atanh(real x) @safe pure nothrow
796	{
797	// log( (1+x)/(1-x) ) == log ( 1 + (2*x)/(1-x) )
798	return 0.5 * log1p( 2 * x / (1 - x) );
799	}
800	/// ditto
801	double atanh(double x) @safe pure nothrow { return atanh(cast(real)x); }
802	/// ditto
803	float atanh(float x) @safe pure nothrow { return atanh(cast(real)x); }
804
805
806	unittest
807	{
808	assert(isIdentical(atanh(0.0), 0.0));
809	assert(isIdentical(atanh(-0.0),-0.0));
810	assert(isNaN(atanh(real.nan)));
811	assert(isNaN(atanh(-real.infinity)));
812	}
813
814	/*****************************************
815	* Returns x rounded to a long value using the current rounding mode.
816	* If the integer value of x is
817	* greater than long.max, the result is
818	* indeterminate.
819	*/
820	long rndtol(real x) @safe pure nothrow; /* intrinsic */
821
822
823	/*****************************************
824	* Returns x rounded to a long value using the FE_TONEAREST rounding mode.
825	* If the integer value of x is
826	* greater than long.max, the result is
827	* indeterminate.
828	*/
829	extern (C) real rndtonl(real x);
830
831	/***************************************
832	* Compute square root of x.
833	*
834	* $(TABLE_SV
835	* $(TR $(TH x) $(TH sqrt(x)) $(TH invalid?))
836	* $(TR $(TD -0.0) $(TD -0.0) $(TD no))
837	* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD yes))
838	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
839	* )
840	*/
841
842	@safe pure nothrow
843	{
844	float sqrt(float x); /* intrinsic */
845	double sqrt(double x); /* intrinsic */ /// ditto
846	real sqrt(real x); /* intrinsic */ /// ditto
847	}
848
849	@trusted pure nothrow { // Should be @safe. See bugs 4628, 4630.
850	// Create explicit overloads for integer sqrts. No ddoc for these because
851	// hopefully a more elegant solution will eventually be found, so we don't
852	// want people relying too heavily on the minutiae of this, for example,
853	// by taking the address of sqrt(int) or something.
854	real sqrt(byte x) { return sqrt(cast(real) x); }
855	real sqrt(ubyte x) { return sqrt(cast(real) x); }
856	real sqrt(short x) { return sqrt(cast(real) x); }
857	real sqrt(ushort x) { return sqrt(cast(real) x); }
858	real sqrt(int x) { return sqrt(cast(real) x); }
859	real sqrt(uint x) { return sqrt(cast(real) x); }
860	real sqrt(long x) { return sqrt(cast(real) x); }
861	real sqrt(ulong x) { return sqrt(cast(real) x); }
862	}
863
864	unittest {
865	alias TypeTuple!(byte, ubyte, short, ushort,
866	int, uint, long, ulong, float, double, real) Numerics;
867	foreach(T; Numerics) {
868	immutable T two = 2;
869	assert(approxEqual(sqrt(two), SQRT2),
870	"sqrt unittest failed on type " ~ T.stringof);
871	}
872	}
873
874	creal sqrt(creal z) @safe pure nothrow
875	{
876	creal c;
877	real x,y,w,r;
878
879	if (z == 0)
880	{
881	c = 0 + 0i;
882	}
883	else
884	{
885	real z_re = z.re;
886	real z_im = z.im;
887
888	x = fabs(z_re);
889	y = fabs(z_im);
890	if (x >= y)
891	{
892	r = y / x;
893	w = sqrt(x) * sqrt(0.5 * (1 + sqrt(1 + r * r)));
894	}
895	else
896	{
897	r = x / y;
898	w = sqrt(y) * sqrt(0.5 * (r + sqrt(1 + r * r)));
899	}
900
901	if (z_re >= 0)
902	{
903	c = w + (z_im / (w + w)) * 1.0i;
904	}
905	else
906	{
907	if (z_im < 0)
908	w = -w;
909	c = z_im / (w + w) + w * 1.0i;
910	}
911	}
912	return c;
913	}
914
915	/**
916	* Calculates e$(SUP x).
917	*
918	* $(TABLE_SV
919	* $(TR $(TH x) $(TH e$(SUP x)) )
920	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
921	* $(TR $(TD -$(INFIN)) $(TD +0.0) )
922	* $(TR $(TD $(NAN)) $(TD $(NAN)) )
923	* )
924	*/
925	real exp(real x) @safe pure nothrow
926	{
927	version(D_InlineAsm_X86)
928	{
929	// e^^x = 2^^(LOG2E*x)
930	// (This is valid because the overflow & underflow limits for exp
931	// and exp2 are so similar).
932	return exp2(LOG2E*x);
933	}
934	else version(D_InlineAsm_X86_64)
935	{
936	// e^^x = 2^^(LOG2E*x)
937	// (This is valid because the overflow & underflow limits for exp
938	// and exp2 are so similar).
939	return exp2(LOG2E*x);
940	} else {
941	return core.stdc.math.exp(x);
942	}
943	}
944	/// ditto
945	double exp(double x) @safe pure nothrow { return exp(cast(real)x); }
946	/// ditto
947	float exp(float x) @safe pure nothrow { return exp(cast(real)x); }
948
949
950	/**
951	* Calculates the value of the natural logarithm base (e)
952	* raised to the power of x, minus 1.
953	*
954	* For very small x, expm1(x) is more accurate
955	* than exp(x)-1.
956	*
957	* $(TABLE_SV
958	* $(TR $(TH x) $(TH e$(SUP x)-1) )
959	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
960	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
961	* $(TR $(TD -$(INFIN)) $(TD -1.0) )
962	* $(TR $(TD $(NAN)) $(TD $(NAN)) )
963	* )
964	*/
965	real expm1(real x) @trusted pure nothrow
966	{
967	version(D_InlineAsm_X86) {
968	enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
969	asm {
970	/* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
971	* Author: Don Clugston.
972	*
973	* expm1(x) = 2^^(rndint(y))* 2^^(y-rndint(y)) - 1 where y = LN2*x.
974	* = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^^(rndint(y))
975	* and 2ym1 = (2^^(y-rndint(y))-1).
976	* If 2rndy < 0.5*real.epsilon, result is -1.
977	* Implementation is otherwise the same as for exp2()
978	*/
979	naked;
980	fld real ptr [ESP+4] ; // x
981	mov AX, [ESP+4+8]; // AX = exponent and sign
982	sub ESP, 12+8; // Create scratch space on the stack
983	// [ESP,ESP+2] = scratchint
984	// [ESP+4..+6, +8..+10, +10] = scratchreal
985	// set scratchreal mantissa = 1.0
986	mov dword ptr [ESP+8], 0;
987	mov dword ptr [ESP+8+4], 0x80000000;
988	and AX, 0x7FFF; // drop sign bit
989	cmp AX, 0x401D; // avoid InvalidException in fist
990	jae L_extreme;
991	fldl2e;
992	fmulp ST(1), ST; // y = x*log2(e)
993	fist dword ptr [ESP]; // scratchint = rndint(y)
994	fisub dword ptr [ESP]; // y - rndint(y)
995	// and now set scratchreal exponent
996	mov EAX, [ESP];
997	add EAX, 0x3fff;
998	jle short L_largenegative;
999	cmp EAX,0x8000;
1000	jge short L_largepositive;
1001	mov [ESP+8+8],AX;
1002	f2xm1; // 2ym1 = 2^^(y-rndint(y)) -1
1003	fld real ptr [ESP+8] ; // 2rndy = 2^^rndint(y)
1004	fmul ST(1), ST; // ST=2rndy, ST(1)=2rndy*2ym1
1005	fld1;
1006	fsubp ST(1), ST; // ST = 2rndy-1, ST(1) = 2rndy * 2ym1 - 1
1007	faddp ST(1), ST; // ST = 2rndy * 2ym1 + 2rndy - 1
1008	add ESP,12+8;
1009	ret PARAMSIZE;
1010
1011	L_extreme: // Extreme exponent. X is very large positive, very
1012	// large negative, infinity, or NaN.
1013	fxam;
1014	fstsw AX;
1015	test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1016	jz L_was_nan; // if x is NaN, returns x
1017	test AX, 0x0200;
1018	jnz L_largenegative;
1019	L_largepositive:
1020	// Set scratchreal = real.max.
1021	// squaring it will create infinity, and set overflow flag.
1022	mov word ptr [ESP+8+8], 0x7FFE;
1023	fstp ST(0), ST;
1024	fld real ptr [ESP+8]; // load scratchreal
1025	fmul ST(0), ST; // square it, to create havoc!
1026	L_was_nan:
1027	add ESP,12+8;
1028	ret PARAMSIZE;
1029	L_largenegative:
1030	fstp ST(0), ST;
1031	fld1;
1032	fchs; // return -1. Underflow flag is not set.
1033	add ESP,12+8;
1034	ret PARAMSIZE;
1035	}
1036	} else version(D_InlineAsm_X86_64) {
1037	asm
1038	{
1039	/* expm1() for x87 80-bit reals, IEEE754-2008 conformant.
1040	* Author: Don Clugston.
1041	*
1042	* expm1(x) = 2^(rndint(y))* 2^(y-rndint(y)) - 1 where y = LN2*x.
1043	* = 2rndy * 2ym1 + 2rndy - 1, where 2rndy = 2^(rndint(y))
1044	* and 2ym1 = (2^(y-rndint(y))-1).
1045	* If 2rndy < 0.5*real.epsilon, result is -1.
1046	* Implementation is otherwise the same as for exp2()
1047	*/
1048	naked;
1049	fld real ptr [RSP+8] ; // x
1050	mov AX, [RSP+8+8]; // AX = exponent and sign
1051	sub RSP, 24; // Create scratch space on the stack
1052	// [RSP,RSP+2] = scratchint
1053	// [RSP+4..+6, +8..+10, +10] = scratchreal
1054	// set scratchreal mantissa = 1.0
1055	mov dword ptr [RSP+8], 0;
1056	mov dword ptr [RSP+8+4], 0x80000000;
1057	and AX, 0x7FFF; // drop sign bit
1058	cmp AX, 0x401D; // avoid InvalidException in fist
1059	jae L_extreme;
1060	fldl2e;
1061	fmul ; // y = x*log2(e)
1062	fist dword ptr [RSP]; // scratchint = rndint(y)
1063	fisub dword ptr [RSP]; // y - rndint(y)
1064	// and now set scratchreal exponent
1065	mov EAX, [RSP];
1066	add EAX, 0x3fff;
1067	jle short L_largenegative;
1068	cmp EAX,0x8000;
1069	jge short L_largepositive;
1070	mov [RSP+8+8],AX;
1071	f2xm1; // 2^(y-rndint(y)) -1
1072	fld real ptr [RSP+8] ; // 2^rndint(y)
1073	fmul ST(1), ST;
1074	fld1;
1075	fsubp ST(1), ST;
1076	fadd;
1077	add RSP,24;
1078	ret;
1079
1080	L_extreme: // Extreme exponent. X is very large positive, very
1081	// large negative, infinity, or NaN.
1082	fxam;
1083	fstsw AX;
1084	test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1085	jz L_was_nan; // if x is NaN, returns x
1086	test AX, 0x0200;
1087	jnz L_largenegative;
1088	L_largepositive:
1089	// Set scratchreal = real.max.
1090	// squaring it will create infinity, and set overflow flag.
1091	mov word ptr [RSP+8+8], 0x7FFE;
1092	fstp ST(0), ST;
1093	fld real ptr [RSP+8]; // load scratchreal
1094	fmul ST(0), ST; // square it, to create havoc!
1095	L_was_nan:
1096	add RSP,24;
1097	ret;
1098
1099	L_largenegative:
1100	fstp ST(0), ST;
1101	fld1;
1102	fchs; // return -1. Underflow flag is not set.
1103	add RSP,24;
1104	ret;
1105	}
1106	} else {
1107	return core.stdc.math.expm1(x);
1108	}
1109	}
1110
1111
1112
1113	/**
1114	* Calculates 2$(SUP x).
1115	*
1116	* $(TABLE_SV
1117	* $(TR $(TH x) $(TH exp2(x)) )
1118	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) )
1119	* $(TR $(TD -$(INFIN)) $(TD +0.0) )
1120	* $(TR $(TD $(NAN)) $(TD $(NAN)) )
1121	* )
1122	*/
1123	real exp2(real x) @trusted pure nothrow
1124	{
1125	version(D_InlineAsm_X86) {
1126	enum { PARAMSIZE = (real.sizeof+3)&(0xFFFF_FFFC) } // always a multiple of 4
1127	asm {
1128	/* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1129	* Author: Don Clugston.
1130	*
1131	* exp2(x) = 2^^(rndint(x))* 2^^(y-rndint(x))
1132	* The trick for high performance is to avoid the fscale(28cycles on core2),
1133	* frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1134	*
1135	* We can do frndint by using fist. BUT we can't use it for huge numbers,
1136	* because it will set the Invalid Operation flag if overflow or NaN occurs.
1137	* Fortunately, whenever this happens the result would be zero or infinity.
1138	*
1139	* We can perform fscale by directly poking into the exponent. BUT this doesn't
1140	* work for the (very rare) cases where the result is subnormal. So we fall back
1141	* to the slow method in that case.
1142	*/
1143	naked;
1144	fld real ptr [ESP+4] ; // x
1145	mov AX, [ESP+4+8]; // AX = exponent and sign
1146	sub ESP, 12+8; // Create scratch space on the stack
1147	// [ESP,ESP+2] = scratchint
1148	// [ESP+4..+6, +8..+10, +10] = scratchreal
1149	// set scratchreal mantissa = 1.0
1150	mov dword ptr [ESP+8], 0;
1151	mov dword ptr [ESP+8+4], 0x80000000;
1152	and AX, 0x7FFF; // drop sign bit
1153	cmp AX, 0x401D; // avoid InvalidException in fist
1154	jae L_extreme;
1155	fist dword ptr [ESP]; // scratchint = rndint(x)
1156	fisub dword ptr [ESP]; // x - rndint(x)
1157	// and now set scratchreal exponent
1158	mov EAX, [ESP];
1159	add EAX, 0x3fff;
1160	jle short L_subnormal;
1161	cmp EAX,0x8000;
1162	jge short L_overflow;
1163	mov [ESP+8+8],AX;
1164	L_normal:
1165	f2xm1;
1166	fld1;
1167	faddp ST(1), ST; // 2^^(x-rndint(x))
1168	fld real ptr [ESP+8] ; // 2^^rndint(x)
1169	add ESP,12+8;
1170	fmulp ST(1), ST;
1171	ret PARAMSIZE;
1172
1173	L_subnormal:
1174	// Result will be subnormal.
1175	// In this rare case, the simple poking method doesn't work.
1176	// The speed doesn't matter, so use the slow fscale method.
1177	fild dword ptr [ESP]; // scratchint
1178	fld1;
1179	fscale;
1180	fstp real ptr [ESP+8]; // scratchreal = 2^^scratchint
1181	fstp ST(0),ST; // drop scratchint
1182	jmp L_normal;
1183
1184	L_extreme: // Extreme exponent. X is very large positive, very
1185	// large negative, infinity, or NaN.
1186	fxam;
1187	fstsw AX;
1188	test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1189	jz L_was_nan; // if x is NaN, returns x
1190	// set scratchreal = real.min_normal
1191	// squaring it will return 0, setting underflow flag
1192	mov word ptr [ESP+8+8], 1;
1193	test AX, 0x0200;
1194	jnz L_waslargenegative;
1195	L_overflow:
1196	// Set scratchreal = real.max.
1197	// squaring it will create infinity, and set overflow flag.
1198	mov word ptr [ESP+8+8], 0x7FFE;
1199	L_waslargenegative:
1200	fstp ST(0), ST;
1201	fld real ptr [ESP+8]; // load scratchreal
1202	fmul ST(0), ST; // square it, to create havoc!
1203	L_was_nan:
1204	add ESP,12+8;
1205	ret PARAMSIZE;
1206	}
1207	} else version(D_InlineAsm_X86_64) {
1208	asm {
1209	/* exp2() for x87 80-bit reals, IEEE754-2008 conformant.
1210	* Author: Don Clugston.
1211	*
1212	* exp2(x) = 2^(rndint(x))* 2^(y-rndint(x))
1213	* The trick for high performance is to avoid the fscale(28cycles on core2),
1214	* frndint(19 cycles), leaving f2xm1(19 cycles) as the only slow instruction.
1215	*
1216	* We can do frndint by using fist. BUT we can't use it for huge numbers,
1217	* because it will set the Invalid Operation flag is overflow or NaN occurs.
1218	* Fortunately, whenever this happens the result would be zero or infinity.
1219	*
1220	* We can perform fscale by directly poking into the exponent. BUT this doesn't
1221	* work for the (very rare) cases where the result is subnormal. So we fall back
1222	* to the slow method in that case.
1223	*/
1224	naked;
1225	fld real ptr [RSP+8] ; // x
1226	mov AX, [RSP+8+8]; // AX = exponent and sign
1227	sub RSP, 24; // Create scratch space on the stack
1228	// [RSP,RSP+2] = scratchint
1229	// [RSP+4..+6, +8..+10, +10] = scratchreal
1230	// set scratchreal mantissa = 1.0
1231	mov dword ptr [RSP+8], 0;
1232	mov dword ptr [RSP+8+4], 0x80000000;
1233	and AX, 0x7FFF; // drop sign bit
1234	cmp AX, 0x401D; // avoid InvalidException in fist
1235	jae L_extreme;
1236	fist dword ptr [RSP]; // scratchint = rndint(x)
1237	fisub dword ptr [RSP]; // x - rndint(x)
1238	// and now set scratchreal exponent
1239	mov EAX, [RSP];
1240	add EAX, 0x3fff;
1241	jle short L_subnormal;
1242	cmp EAX,0x8000;
1243	jge short L_overflow;
1244	mov [RSP+8+8],AX;
1245	L_normal:
1246	f2xm1;
1247	fld1;
1248	fadd; // 2^(x-rndint(x))
1249	fld real ptr [RSP+8] ; // 2^rndint(x)
1250	add RSP,24;
1251	fmulp ST(1), ST;
1252	ret;
1253
1254	L_subnormal:
1255	// Result will be subnormal.
1256	// In this rare case, the simple poking method doesn't work.
1257	// The speed doesn't matter, so use the slow fscale method.
1258	fild dword ptr [RSP]; // scratchint
1259	fld1;
1260	fscale;
1261	fstp real ptr [RSP+8]; // scratchreal = 2^scratchint
1262	fstp ST(0),ST; // drop scratchint
1263	jmp L_normal;
1264
1265	L_extreme: // Extreme exponent. X is very large positive, very
1266	// large negative, infinity, or NaN.
1267	fxam;
1268	fstsw AX;
1269	test AX, 0x0400; // NaN_or_zero, but we already know x!=0
1270	jz L_was_nan; // if x is NaN, returns x
1271	// set scratchreal = real.min
1272	// squaring it will return 0, setting underflow flag
1273	mov word ptr [RSP+8+8], 1;
1274	test AX, 0x0200;
1275	jnz L_waslargenegative;
1276	L_overflow:
1277	// Set scratchreal = real.max.
1278	// squaring it will create infinity, and set overflow flag.
1279	mov word ptr [RSP+8+8], 0x7FFE;
1280	L_waslargenegative:
1281	fstp ST(0), ST;
1282	fld real ptr [RSP+8]; // load scratchreal
1283	fmul ST(0), ST; // square it, to create havoc!
1284	L_was_nan:
1285	add RSP,24;
1286	ret;
1287	}
1288	} else {
1289	return core.stdc.math.exp2(x);
1290	}
1291	}
1292
1293	unittest{
1294	assert(exp2(0.5L)== SQRT2);
1295	assert(exp2(8.0L) == 256.0);
1296	assert(exp2(-9.0L)== 1.0L/512.0);
1297	assert(exp(3.0L) == EEE);
1298	}
1299
1300	unittest
1301	{
1302	FloatingPointControl ctrl;
1303	ctrl.disableExceptions(FloatingPointControl.allExceptions);
1304	ctrl.rounding = FloatingPointControl.roundToNearest;
1305
1306	// @@BUG@@: Non-immutable array literals are ridiculous.
1307	// Note that these are only valid for 80-bit reals: overflow will be different for 64-bit reals.
1308	static const real [2][] exptestpoints =
1309	[ // x, exp(x)
1310	[1.0L, E ],
1311	[0.5L, 0x1.A612_98E1_E069_BC97p+0L ],
1312	[3.0L, EEE ],
1313	[0x1.1p13L, 0x1.29aeffefc8ec645p+12557L ], // near overflow
1314	[-0x1.18p13L, 0x1.5e4bf54b4806db9p-12927L ], // near underflow
1315	[-0x1.625p13L, 0x1.a6bd68a39d11f35cp-16358L],
1316	[-0x1p30L, 0 ], // underflow - subnormal
1317	[-0x1.62DAFp13L, 0x1.96c53d30277021dp-16383L ],
1318	[-0x1.643p13L, 0x1p-16444L ],
1319	[-0x1.645p13L, 0 ], // underflow to zero
1320	[0x1p80L, real.infinity ], // far overflow
1321	[real.infinity, real.infinity ],
1322	[0x1.7p13L, real.infinity ] // close overflow
1323	];
1324	real x;
1325	IeeeFlags f;
1326	for (int i=0; i<exptestpoints.length;++i) {
1327	resetIeeeFlags();
1328	x = exp(exptestpoints[i][0]);
1329	f = ieeeFlags();
1330	assert(x == exptestpoints[i][1]);
1331	// Check the overflow bit
1332	assert(f.overflow() == (fabs(x) == real.infinity));
1333	// Check the underflow bit
1334	assert(f.underflow() == (fabs(x) < real.min_normal));
1335	// Invalid and div by zero shouldn't be affected.
1336	assert(!f.invalid);
1337	assert(!f.divByZero);
1338	}
1339	// Ideally, exp(0) would not set the inexact flag.
1340	// Unfortunately, fldl2e sets it!
1341	// So it's not realistic to avoid setting it.
1342	assert(exp(0.0L) == 1.0);
1343
1344	// NaN propagation. Doesn't set flags, bcos was already NaN.
1345	resetIeeeFlags();
1346	x = exp(real.nan);
1347	f = ieeeFlags();
1348	assert(isIdentical(x,real.nan));
1349	assert(f.flags == 0);
1350
1351	resetIeeeFlags();
1352	x = exp(-real.nan);
1353	f = ieeeFlags;
1354	assert(isIdentical(x, -real.nan));
1355	assert(f.flags == 0);
1356
1357	x = exp(NaN(0x123));
1358	assert(isIdentical(x, NaN(0x123)));
1359
1360	// High resolution test
1361	assert(exp(0.5L) == 0x1.A612_98E1_E069_BC97_2DFE_FAB6D_33Fp+0L);
1362
1363	}
1364
1365
1366	/**
1367	* Calculate cos(y) + i sin(y).
1368	*
1369	* On many CPUs (such as x86), this is a very efficient operation;
1370	* almost twice as fast as calculating sin(y) and cos(y) separately,
1371	* and is the preferred method when both are required.
1372	*/
1373	creal expi(real y) @trusted pure nothrow
1374	{
1375	version(InlineAsm_X86_Any)
1376	{
1377	asm
1378	{
1379	fld y;
1380	fsincos;
1381	fxch ST(1), ST(0);
1382	}
1383	}
1384	else
1385	{
1386	return cos(y) + sin(y)*1i;
1387	}
1388	}
1389
1390	unittest
1391	{
1392	assert(expi(1.3e5L) == cos(1.3e5L) + sin(1.3e5L) * 1i);
1393	assert(expi(0.0L) == 1L + 0.0Li);
1394	}
1395
1396	/*********************************************************************
1397	* Separate floating point value into significand and exponent.
1398	*
1399	* Returns:
1400	* Calculate and return $(I x) and $(I exp) such that
1401	* value =$(I x)*2$(SUP exp) and
1402	* .5 $(LT)= \|$(I x)\| $(LT) 1.0
1403	*
1404	* $(I x) has same sign as value.
1405	*
1406	* $(TABLE_SV
1407	* $(TR $(TH value) $(TH returns) $(TH exp))
1408	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD 0))
1409	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD int.max))
1410	* $(TR $(TD -$(INFIN)) $(TD -$(INFIN)) $(TD int.min))
1411	* $(TR $(TD $(PLUSMN)$(NAN)) $(TD $(PLUSMN)$(NAN)) $(TD int.min))
1412	* )
1413	*/
1414	real frexp(real value, out int exp) @trusted pure nothrow
1415	{
1416	ushort* vu = cast(ushort*)&value;
1417	long* vl = cast(long*)&value;
1418	uint ex;
1419	alias floatTraits!(real) F;
1420
1421	ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1422	static if (real.mant_dig == 64) { // real80
1423	if (ex) { // If exponent is non-zero
1424	if (ex == F.EXPMASK) { // infinity or NaN
1425	if (*vl & 0x7FFF_FFFF_FFFF_FFFF) { // NaN
1426	*vl \|= 0xC000_0000_0000_0000; // convert NaNS to NaNQ
1427	exp = int.min;
1428	} else if (vu[F.EXPPOS_SHORT] & 0x8000) { // negative infinity
1429	exp = int.min;
1430	} else { // positive infinity
1431	exp = int.max;
1432	}
1433	} else {
1434	exp = ex - F.EXPBIAS;
1435	vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FFE;
1436	}
1437	} else if (!*vl) {
1438	// value is +-0.0
1439	exp = 0;
1440	} else {
1441	// denormal
1442	value *= F.RECIP_EPSILON;
1443	ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1444	exp = ex - F.EXPBIAS - real.mant_dig + 1;
1445	vu[F.EXPPOS_SHORT] = (0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FFE;
1446	}
1447	} else static if (real.mant_dig == 113) { // quadruple
1448	if (ex) { // If exponent is non-zero
1449	if (ex == F.EXPMASK) { // infinity or NaN
1450	if (vl[MANTISSA_LSB] \|
1451	( vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) { // NaN
1452	// convert NaNS to NaNQ
1453	vl[MANTISSA_MSB] \|= 0x0000_8000_0000_0000;
1454	exp = int.min;
1455	} else if (vu[F.EXPPOS_SHORT] & 0x8000) { // negative infinity
1456	exp = int.min;
1457	} else { // positive infinity
1458	exp = int.max;
1459	}
1460	} else {
1461	exp = ex - F.EXPBIAS;
1462	vu[F.EXPPOS_SHORT] =
1463	cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FFE);
1464	}
1465	} else if ((vl[MANTISSA_LSB]
1466	\|(vl[MANTISSA_MSB] & 0x0000_FFFF_FFFF_FFFF)) == 0) {
1467	// value is +-0.0
1468	exp = 0;
1469	} else {
1470	// denormal
1471	value *= F.RECIP_EPSILON;
1472	ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1473	exp = ex - F.EXPBIAS - real.mant_dig + 1;
1474	vu[F.EXPPOS_SHORT] =
1475	cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FFE);
1476	}
1477	} else static if (real.mant_dig==53) { // real is double
1478	if (ex) { // If exponent is non-zero
1479	if (ex == F.EXPMASK) { // infinity or NaN
1480	if (*vl == 0x7FF0_0000_0000_0000) { // positive infinity
1481	exp = int.max;
1482	} else if (*vl == 0xFFF0_0000_0000_0000) { // negative infinity
1483	exp = int.min;
1484	} else { // NaN
1485	*vl \|= 0x0008_0000_0000_0000; // convert NaNS to NaNQ
1486	exp = int.min;
1487	}
1488	} else {
1489	exp = (ex - F.EXPBIAS) >> 4;
1490	vu[F.EXPPOS_SHORT] = cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FE0);
1491	}
1492	} else if (!(*vl & 0x7FFF_FFFF_FFFF_FFFF)) {
1493	// value is +-0.0
1494	exp = 0;
1495	} else {
1496	// denormal
1497	value *= F.RECIP_EPSILON;
1498	ex = vu[F.EXPPOS_SHORT] & F.EXPMASK;
1499	exp = ((ex - F.EXPBIAS)>> 4) - real.mant_dig + 1;
1500	vu[F.EXPPOS_SHORT] =
1501	cast(ushort)((0x8000 & vu[F.EXPPOS_SHORT]) \| 0x3FE0);
1502	}
1503	} else { //static if(real.mant_dig==106) // doubledouble
1504	assert (0, "frexp not implemented");
1505	}
1506	return value;
1507	}
1508
1509
1510	unittest
1511	{
1512	static real vals[][3] = // x,frexp,exp
1513	[
1514	[0.0, 0.0, 0],
1515	[-0.0, -0.0, 0],
1516	[1.0, .5, 1],
1517	[-1.0, -.5, 1],
1518	[2.0, .5, 2],
1519	[double.min_normal/2.0, .5, -1022],
1520	[real.infinity,real.infinity,int.max],
1521	[-real.infinity,-real.infinity,int.min],
1522	[real.nan,real.nan,int.min],
1523	[-real.nan,-real.nan,int.min],
1524	];
1525
1526	int i;
1527
1528	for (i = 0; i < vals.length; i++) {
1529	real x = vals[i][0];
1530	real e = vals[i][1];
1531	int exp = cast(int)vals[i][2];
1532	int eptr;
1533	real v = frexp(x, eptr);
1534	assert(isIdentical(e, v));
1535	assert(exp == eptr);
1536
1537	}
1538	static if (real.mant_dig == 64) {
1539	static real extendedvals[][3] = [ // x,frexp,exp
1540	[0x1.a5f1c2eb3fe4efp+73L, 0x1.A5F1C2EB3FE4EFp-1L, 74], // normal
1541	[0x1.fa01712e8f0471ap-1064L, 0x1.fa01712e8f0471ap-1L, -1063],
1542	[real.min_normal, .5, -16381],
1543	[real.min_normal/2.0L, .5, -16382] // denormal
1544	];
1545
1546	for (i = 0; i < extendedvals.length; i++) {
1547	real x = extendedvals[i][0];
1548	real e = extendedvals[i][1];
1549	int exp = cast(int)extendedvals[i][2];
1550	int eptr;
1551	real v = frexp(x, eptr);
1552	assert(isIdentical(e, v));
1553	assert(exp == eptr);
1554
1555	}
1556	}
1557	}
1558
1559	/******************************************
1560	* Extracts the exponent of x as a signed integral value.
1561	*
1562	* If x is not a special value, the result is the same as
1563	* $(D cast(int)logb(x)).
1564	*
1565	* $(TABLE_SV
1566	* $(TR $(TH x) $(TH ilogb(x)) $(TH Range error?))
1567	* $(TR $(TD 0) $(TD FP_ILOGB0) $(TD yes))
1568	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD int.max) $(TD no))
1569	* $(TR $(TD $(NAN)) $(TD FP_ILOGBNAN) $(TD no))
1570	* )
1571	*/
1572	int ilogb(real x) @trusted nothrow { return core.stdc.math.ilogbl(x); }
1573
1574	alias core.stdc.math.FP_ILOGB0 FP_ILOGB0;
1575	alias core.stdc.math.FP_ILOGBNAN FP_ILOGBNAN;
1576
1577
1578	/*******************************************
1579	* Compute n * 2$(SUP exp)
1580	* References: frexp
1581	*/
1582
1583	real ldexp(real n, int exp) @safe pure nothrow; /* intrinsic */
1584
1585	unittest {
1586	assert(ldexp(1, -16384) == 0x1p-16384L);
1587	assert(ldexp(1, -16382) == 0x1p-16382L);
1588	int x;
1589	real n = frexp(0x1p-16384L, x);
1590	assert(n==0.5L);
1591	assert(x==-16383);
1592	assert(ldexp(n, x)==0x1p-16384L);
1593
1594	}
1595
1596	/**************************************
1597	* Calculate the natural logarithm of x.
1598	*
1599	* $(TABLE_SV
1600	* $(TR $(TH x) $(TH log(x)) $(TH divide by 0?) $(TH invalid?))
1601	* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
1602	* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
1603	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
1604	* )
1605	*/
1606
1607	real log(real x) @safe pure nothrow
1608	{
1609	version (INLINE_YL2X)
1610	return yl2x(x, LN2);
1611	else
1612	return core.stdc.math.logl(x);
1613	}
1614
1615	unittest
1616	{
1617	assert(log(E) == 1);
1618	}
1619
1620	/**************************************
1621	* Calculate the base-10 logarithm of x.
1622	*
1623	* $(TABLE_SV
1624	* $(TR $(TH x) $(TH log10(x)) $(TH divide by 0?) $(TH invalid?))
1625	* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
1626	* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes))
1627	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no))
1628	* )
1629	*/
1630
1631	real log10(real x) @safe pure nothrow
1632	{
1633	version (INLINE_YL2X)
1634	return yl2x(x, LOG2);
1635	else
1636	return core.stdc.math.log10l(x);
1637	}
1638
1639	unittest
1640	{
1641	//printf("%Lg\n", log10(1000) - 3);
1642	assert(fabs(log10(1000) - 3) < .000001);
1643	}
1644
1645	/******************************************
1646	* Calculates the natural logarithm of 1 + x.
1647	*
1648	* For very small x, log1p(x) will be more accurate than
1649	* log(1 + x).
1650	*
1651	* $(TABLE_SV
1652	* $(TR $(TH x) $(TH log1p(x)) $(TH divide by 0?) $(TH invalid?))
1653	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) $(TD no))
1654	* $(TR $(TD -1.0) $(TD -$(INFIN)) $(TD yes) $(TD no))
1655	* $(TR $(TD $(LT)-1.0) $(TD $(NAN)) $(TD no) $(TD yes))
1656	* $(TR $(TD +$(INFIN)) $(TD -$(INFIN)) $(TD no) $(TD no))
1657	* )
1658	*/
1659
1660	real log1p(real x) @safe pure nothrow
1661	{
1662	version(INLINE_YL2X)
1663	{
1664	// On x87, yl2xp1 is valid if and only if -0.5 <= lg(x) <= 0.5,
1665	// ie if -0.29<=x<=0.414
1666	return (fabs(x) <= 0.25) ? yl2xp1(x, LN2) : yl2x(x+1, LN2);
1667	}
1668	else
1669	{
1670	return core.stdc.math.log1pl(x);
1671	}
1672	}
1673
1674	/***************************************
1675	* Calculates the base-2 logarithm of x:
1676	* $(SUB log, 2)x
1677	*
1678	* $(TABLE_SV
1679	* $(TR $(TH x) $(TH log2(x)) $(TH divide by 0?) $(TH invalid?))
1680	* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) $(TD no) )
1681	* $(TR $(TD $(LT)0.0) $(TD $(NAN)) $(TD no) $(TD yes) )
1682	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no) $(TD no) )
1683	* )
1684	*/
1685	real log2(real x) @safe pure nothrow
1686	{
1687	version (INLINE_YL2X)
1688	return yl2x(x, 1);
1689	else
1690	return core.stdc.math.log2l(x);
1691	}
1692
1693	/*****************************************
1694	* Extracts the exponent of x as a signed integral value.
1695	*
1696	* If x is subnormal, it is treated as if it were normalized.
1697	* For a positive, finite x:
1698	*
1699	* 1 $(LT)= $(I x) * FLT_RADIX$(SUP -logb(x)) $(LT) FLT_RADIX
1700	*
1701	* $(TABLE_SV
1702	* $(TR $(TH x) $(TH logb(x)) $(TH divide by 0?) )
1703	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
1704	* $(TR $(TD $(PLUSMN)0.0) $(TD -$(INFIN)) $(TD yes) )
1705	* )
1706	*/
1707	real logb(real x) @trusted nothrow { return core.stdc.math.logbl(x); }
1708
1709	/************************************
1710	* Calculates the remainder from the calculation x/y.
1711	* Returns:
1712	* The value of x - i * y, where i is the number of times that y can
1713	* be completely subtracted from x. The result has the same sign as x.
1714	*
1715	* $(TABLE_SV
1716	* $(TR $(TH x) $(TH y) $(TH modf(x, y)) $(TH invalid?))
1717	* $(TR $(TD $(PLUSMN)0.0) $(TD not 0.0) $(TD $(PLUSMN)0.0) $(TD no))
1718	* $(TR $(TD $(PLUSMNINF)) $(TD anything) $(TD $(NAN)) $(TD yes))
1719	* $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD $(NAN)) $(TD yes))
1720	* $(TR $(TD !=$(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x) $(TD no))
1721	* )
1722	*/
1723	real modf(real x, ref real y) @trusted nothrow { return core.stdc.math.modfl(x,&y); }
1724
1725	/*************************************
1726	* Efficiently calculates x * 2$(SUP n).
1727	*
1728	* scalbn handles underflow and overflow in
1729	* the same fashion as the basic arithmetic operators.
1730	*
1731	* $(TABLE_SV
1732	* $(TR $(TH x) $(TH scalb(x)))
1733	* $(TR $(TD $(PLUSMNINF)) $(TD $(PLUSMNINF)) )
1734	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) )
1735	* )
1736	*/
1737	real scalbn(real x, int n) @trusted nothrow
1738	{
1739	version(InlineAsm_X86_Any) {
1740	// scalbnl is not supported on DMD-Windows, so use asm.
1741	asm {
1742	fild n;
1743	fld x;
1744	fscale;
1745	fstp ST(1), ST;
1746	}
1747	} else {
1748	return core.stdc.math.scalbnl(x, n);
1749	}
1750	}
1751
1752	unittest {
1753	assert(scalbn(-real.infinity, 5) == -real.infinity);
1754	}
1755
1756	/***************
1757	* Calculates the cube root of x.
1758	*
1759	* $(TABLE_SV
1760	* $(TR $(TH $(I x)) $(TH cbrt(x)) $(TH invalid?))
1761	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMN)0.0) $(TD no) )
1762	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes) )
1763	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD $(PLUSMN)$(INFIN)) $(TD no) )
1764	* )
1765	*/
1766	real cbrt(real x) @trusted nothrow { return core.stdc.math.cbrtl(x); }
1767
1768
1769	/*******************************
1770	* Returns \|x\|
1771	*
1772	* $(TABLE_SV
1773	* $(TR $(TH x) $(TH fabs(x)))
1774	* $(TR $(TD $(PLUSMN)0.0) $(TD +0.0) )
1775	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) )
1776	* )
1777	*/
1778	real fabs(real x) @safe pure nothrow; /* intrinsic */
1779
1780
1781	/***********************************************************************
1782	* Calculates the length of the
1783	* hypotenuse of a right-angled triangle with sides of length x and y.
1784	* The hypotenuse is the value of the square root of
1785	* the sums of the squares of x and y:
1786	*
1787	* sqrt($(POW x, 2) + $(POW y, 2))
1788	*
1789	* Note that hypot(x, y), hypot(y, x) and
1790	* hypot(x, -y) are equivalent.
1791	*
1792	* $(TABLE_SV
1793	* $(TR $(TH x) $(TH y) $(TH hypot(x, y)) $(TH invalid?))
1794	* $(TR $(TD x) $(TD $(PLUSMN)0.0) $(TD \|x\|) $(TD no))
1795	* $(TR $(TD $(PLUSMNINF)) $(TD y) $(TD +$(INFIN)) $(TD no))
1796	* $(TR $(TD $(PLUSMNINF)) $(TD $(NAN)) $(TD +$(INFIN)) $(TD no))
1797	* )
1798	*/
1799
1800	real hypot(real x, real y) @safe pure nothrow
1801	{
1802	// Scale x and y to avoid underflow and overflow.
1803	// If one is huge and the other tiny, return the larger.
1804	// If both are huge, avoid overflow by scaling by 1/sqrt(real.max/2).
1805	// If both are tiny, avoid underflow by scaling by sqrt(real.min_normal*real.epsilon).
1806
1807	enum real SQRTMIN = 0.5*sqrt(real.min_normal); // This is a power of 2.
1808	enum real SQRTMAX = 1.0L/SQRTMIN; // 2^^((max_exp)/2) = nextUp(sqrt(real.max))
1809
1810	static assert(2(SQRTMAX/2)(SQRTMAX/2) <= real.max);
1811	static assert(real.min_normalreal.max>2 && real.min_normalreal.max<=4); // Proves that sqrt(real.max) ~~ 0.5/sqrt(real.min_normal)
1812
1813	real u = fabs(x);
1814	real v = fabs(y);
1815	if (!(u >= v)) // check for NaN as well.
1816	{
1817	v = u;
1818	u = fabs(y);
1819	if (u == real.infinity) return u; // hypot(inf, nan) == inf
1820	if (v == real.infinity) return v; // hypot(nan, inf) == inf
1821	}
1822	// Now u >= v, or else one is NaN.
1823	if (v >= SQRTMAX*0.5)
1824	{
1825	// hypot(huge, huge) -- avoid overflow
1826	u = SQRTMIN0.5;
1827	v = SQRTMIN0.5;
1828	return sqrt(uu + vv) * SQRTMAX * 2.0;
1829	}
1830	if (u <= SQRTMIN)
1831	{
1832	// hypot (tiny, tiny) -- avoid underflow
1833	// This is only necessary to avoid setting the underflow
1834	// flag.
1835	u *= SQRTMAX / real.epsilon;
1836	v *= SQRTMAX / real.epsilon;
1837	return sqrt(uu + vv) * SQRTMIN * real.epsilon;
1838	}
1839	if (u * real.epsilon > v)
1840	{
1841	// hypot (huge, tiny) = huge
1842	return u;
1843	}
1844
1845	// both are in the normal range
1846	return sqrt(uu + vv);
1847	}
1848
1849	unittest
1850	{
1851	static real vals[][3] = // x,y,hypot
1852	[
1853	[ 0.0, 0.0, 0.0],
1854	[ 0.0, -0.0, 0.0],
1855	[ -0.0, -0.0, 0.0],
1856	[ 3.0, 4.0, 5.0],
1857	[ -300, -400, 500],
1858	[0.0, 7.0, 7.0],
1859	[9.0, 9*real.epsilon, 9.0],
1860	[88/(64sqrt(real.min_normal)), 105/(64sqrt(real.min_normal)), 137/(64*sqrt(real.min_normal))],
1861	[88/(128sqrt(real.min_normal)), 105/(128sqrt(real.min_normal)), 137/(128*sqrt(real.min_normal))],
1862	[3real.min_normalreal.epsilon, 4real.min_normalreal.epsilon, 5real.min_normalreal.epsilon],
1863	[ real.min_normal, real.min_normal, sqrt(2.0L)*real.min_normal],
1864	[ real.max/sqrt(2.0L), real.max/sqrt(2.0L), real.max],
1865	[ real.infinity, real.nan, real.infinity],
1866	[ real.nan, real.infinity, real.infinity],
1867	[ real.nan, real.nan, real.nan],
1868	[ real.nan, real.max, real.nan],
1869	[ real.max, real.nan, real.nan],
1870	];
1871	for (int i = 0; i < vals.length; i++)
1872	{
1873	real x = vals[i][0];
1874	real y = vals[i][1];
1875	real z = vals[i][2];
1876	real h = hypot(x, y);
1877	assert(isIdentical(z, h));
1878	}
1879	}
1880
1881	/**********************************
1882	* Returns the error function of x.
1883	*
1884	* <img src="erf.gif" alt="error function">
1885	*/
1886	real erf(real x) @trusted nothrow { return core.stdc.math.erfl(x); }
1887
1888	/**********************************
1889	* Returns the complementary error function of x, which is 1 - erf(x).
1890	*
1891	* <img src="erfc.gif" alt="complementary error function">
1892	*/
1893	real erfc(real x) @trusted nothrow { return core.stdc.math.erfcl(x); }
1894
1895	/***********************************
1896	* Natural logarithm of gamma function.
1897	*
1898	* Returns the base e (2.718...) logarithm of the absolute
1899	* value of the gamma function of the argument.
1900	*
1901	* For reals, lgamma is equivalent to log(fabs(gamma(x))).
1902	*
1903	* $(TABLE_SV
1904	* $(TR $(TH x) $(TH lgamma(x)) $(TH invalid?))
1905	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
1906	* $(TR $(TD integer $(LT)= 0) $(TD +$(INFIN)) $(TD yes))
1907	* $(TR $(TD $(PLUSMN)$(INFIN)) $(TD +$(INFIN)) $(TD no))
1908	* )
1909	*/
1910	real lgamma(real x) @trusted nothrow
1911	{
1912	return core.stdc.math.lgammal(x);
1913
1914	// Use etc.gamma.lgamma for those C systems that are missing it
1915	}
1916
1917	/***********************************
1918	* The Gamma function, $(GAMMA)(x)
1919	*
1920	* $(GAMMA)(x) is a generalisation of the factorial function
1921	* to real and complex numbers.
1922	* Like x!, $(GAMMA)(x+1) = x*$(GAMMA)(x).
1923	*
1924	* Mathematically, if z.re > 0 then
1925	* $(GAMMA)(z) = $(INTEGRATE 0, $(INFIN)) $(POWER t, z-1)$(POWER e, -t) dt
1926	*
1927	* $(TABLE_SV
1928	* $(TR $(TH x) $(TH $(GAMMA)(x)) $(TH invalid?))
1929	* $(TR $(TD $(NAN)) $(TD $(NAN)) $(TD yes))
1930	* $(TR $(TD $(PLUSMN)0.0) $(TD $(PLUSMNINF)) $(TD yes))
1931	* $(TR $(TD integer $(GT)0) $(TD (x-1)!) $(TD no))
1932	* $(TR $(TD integer $(LT)0) $(TD $(NAN)) $(TD yes))
1933	* $(TR $(TD +$(INFIN)) $(TD +$(INFIN)) $(TD no))
1934	* $(TR $(TD -$(INFIN)) $(TD $(NAN)) $(TD yes))
1935	* )
1936	*
1937	* References:
1938	* $(LINK http://en.wikipedia.org/wiki/Gamma_function),
1939	* $(LINK http://www.netlib.org/cephes/ldoubdoc.html#gamma)
1940	*/
1941	real tgamma(real x) @trusted nothrow
1942	{
1943	return core.stdc.math.tgammal(x);
1944
1945	// Use etc.gamma.tgamma for those C systems that are missing it
1946	}
1947
1948	/**************************************
1949	* Returns the value of x rounded upward to the next integer
1950	* (toward positive infinity).
1951	*/
1952	real ceil(real x) @trusted nothrow { return core.stdc.math.ceill(x); }
1953
1954	/**************************************
1955	* Returns the value of x rounded downward to the next integer
1956	* (toward negative infinity).
1957	*/
1958	real floor(real x) @trusted nothrow { return core.stdc.math.floorl(x); }
1959
1960	/******************************************
1961	* Rounds x to the nearest integer value, using the current rounding
1962	* mode.
1963	*
1964	* Unlike the rint functions, nearbyint does not raise the
1965	* FE_INEXACT exception.
1966	*/
1967	real nearbyint(real x) @trusted nothrow { return core.stdc.math.nearbyintl(x); }
1968
1969	/**********************************
1970	* Rounds x to the nearest integer value, using the current rounding
1971	* mode.
1972	* If the return value is not equal to x, the FE_INEXACT
1973	* exception is raised.
1974	* $(B nearbyint) performs
1975	* the same operation, but does not set the FE_INEXACT exception.
1976	*/
1977	real rint(real x) @safe pure nothrow; /* intrinsic */
1978
1979	/***************************************
1980	* Rounds x to the nearest integer value, using the current rounding
1981	* mode.
1982	*
1983	* This is generally the fastest method to convert a floating-point number
1984	* to an integer. Note that the results from this function
1985	* depend on the rounding mode, if the fractional part of x is exactly 0.5.
1986	* If using the default rounding mode (ties round to even integers)
1987	* lrint(4.5) == 4, lrint(5.5)==6.
1988	*/
1989	long lrint(real x) @trusted pure nothrow
1990	{
1991	version(InlineAsm_X86_Any)
1992	{
1993	long n;
1994	asm
1995	{
1996	fld x;
1997	fistp n;
1998	}
1999	return n;
2000	} else {
2001	return core.stdc.math.llrintl(x);
2002	}
2003	}
2004
2005	/*******************************************
2006	* Return the value of x rounded to the nearest integer.
2007	* If the fractional part of x is exactly 0.5, the return value is rounded to
2008	* the even integer.
2009	*/
2010	real round(real x) @trusted nothrow { return core.stdc.math.roundl(x); }
2011
2012	/**********************************************
2013	* Return the value of x rounded to the nearest integer.
2014	*
2015	* If the fractional part of x is exactly 0.5, the return value is rounded
2016	* away from zero.
2017	*/
2018	long lround(real x) @trusted nothrow
2019	{
2020	version (Posix)
2021	return core.stdc.math.llroundl(x);
2022	else
2023	assert (0, "lround not implemented");
2024	}
2025
2026	version(Posix)
2027	{
2028	unittest
2029	{
2030	assert(lround(0.49) == 0);
2031	assert(lround(0.5) == 1);
2032	assert(lround(1.5) == 2);
2033	}
2034	}
2035
2036	/****************************************************
2037	* Returns the integer portion of x, dropping the fractional portion.
2038	*
2039	* This is also known as "chop" rounding.
2040	*/
2041	real trunc(real x) @trusted nothrow { return core.stdc.math.truncl(x); }
2042
2043	/****************************************************
2044	* Calculate the remainder x REM y, following IEC 60559.
2045	*
2046	* REM is the value of x - y * n, where n is the integer nearest the exact
2047	* value of x / y.
2048	* If \|n - x / y\| == 0.5, n is even.
2049	* If the result is zero, it has the same sign as x.
2050	* Otherwise, the sign of the result is the sign of x / y.
2051	* Precision mode has no effect on the remainder functions.
2052	*
2053	* remquo returns n in the parameter n.
2054	*
2055	* $(TABLE_SV
2056	* $(TR $(TH x) $(TH y) $(TH remainder(x, y)) $(TH n) $(TH invalid?))
2057	* $(TR $(TD $(PLUSMN)0.0) $(TD not 0.0) $(TD $(PLUSMN)0.0) $(TD 0.0) $(TD no))
2058	* $(TR $(TD $(PLUSMNINF)) $(TD anything) $(TD $(NAN)) $(TD ?) $(TD yes))
2059	* $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD $(NAN)) $(TD ?) $(TD yes))
2060	* $(TR $(TD != $(PLUSMNINF)) $(TD $(PLUSMNINF)) $(TD x) $(TD ?) $(TD no))
2061	* )
2062	*
2063	* Note: remquo not supported on windows
2064	*/
2065	real remainder(real x, real y) @trusted nothrow { return core.stdc.math.remainderl(x, y); }
2066
2067	real remquo(real x, real y, out int n) @trusted nothrow /// ditto
2068	{
2069	version (Posix)
2070	return core.stdc.math.remquol(x, y, &n);
2071	else
2072	assert (0, "remquo not implemented");
2073	}
2074
2075	/** IEEE exception status flags ('sticky bits')
2076
2077	These flags indicate that an exceptional floating-point condition has occurred.
2078	They indicate that a NaN or an infinity has been generated, that a result
2079	is inexact, or that a signalling NaN has been encountered. If floating-point
2080	exceptions are enabled (unmasked), a hardware exception will be generated
2081	instead of setting these flags.
2082
2083	Example:
2084	----
2085	real a=3.5;
2086	// Set all the flags to zero
2087	resetIeeeFlags();
2088	assert(!ieeeFlags.divByZero);
2089	// Perform a division by zero.
2090	a/=0.0L;
2091	assert(a==real.infinity);
2092	assert(ieeeFlags.divByZero);
2093	// Create a NaN
2094	a*=0.0L;
2095	assert(ieeeFlags.invalid);
2096	assert(isNaN(a));
2097
2098	// Check that calling func() has no effect on the
2099	// status flags.
2100	IeeeFlags f = ieeeFlags;
2101	func();
2102	assert(ieeeFlags == f);
2103
2104	----
2105	*/
2106	struct IeeeFlags
2107	{
2108	private:
2109	// The x87 FPU status register is 16 bits.
2110	// The Pentium SSE2 status register is 32 bits.
2111	uint flags;
2112	version (X86_Any) {
2113	// Applies to both x87 status word (16 bits) and SSE2 status word(32 bits).
2114	enum : int {
2115	INEXACT_MASK = 0x20,
2116	UNDERFLOW_MASK = 0x10,
2117	OVERFLOW_MASK = 0x08,
2118	DIVBYZERO_MASK = 0x04,
2119	INVALID_MASK = 0x01
2120	}
2121	// Don't bother about denormals, they are not supported on most CPUs.
2122	// DENORMAL_MASK = 0x02;
2123	} else version (PPC) {
2124	// PowerPC FPSCR is a 32-bit register.
2125	enum : int {
2126	INEXACT_MASK = 0x600,
2127	UNDERFLOW_MASK = 0x010,
2128	OVERFLOW_MASK = 0x008,
2129	DIVBYZERO_MASK = 0x020,
2130	INVALID_MASK = 0xF80 // PowerPC has five types of invalid exceptions.
2131	}
2132	} else version(SPARC) { // SPARC FSR is a 32bit register
2133	//(64 bits for Sparc 7 & 8, but high 32 bits are uninteresting).
2134	enum : int {
2135	INEXACT_MASK = 0x020,
2136	UNDERFLOW_MASK = 0x080,
2137	OVERFLOW_MASK = 0x100,
2138	DIVBYZERO_MASK = 0x040,
2139	INVALID_MASK = 0x200
2140	}
2141	} else
2142	static assert(0, "Not implemented");
2143	private:
2144	static uint getIeeeFlags()
2145	{
2146	version(D_InlineAsm_X86) {
2147	asm {
2148	fstsw AX;
2149	// NOTE: If compiler supports SSE2, need to OR the result with
2150	// the SSE2 status register.
2151	// Clear all irrelevant bits
2152	and EAX, 0x03D;
2153	}
2154	} else version(D_InlineAsm_X86_64) {
2155	asm {
2156	fstsw AX;
2157	// NOTE: If compiler supports SSE2, need to OR the result with
2158	// the SSE2 status register.
2159	// Clear all irrelevant bits
2160	and RAX, 0x03D;
2161	}
2162	} else version (SPARC) {
2163	/*
2164	int retval;
2165	asm { st %fsr, retval; }
2166	return retval;
2167	*/
2168	assert(0, "Not yet supported");
2169	} else
2170	assert(0, "Not yet supported");
2171	}
2172	static void resetIeeeFlags()
2173	{
2174	version(InlineAsm_X86_Any) {
2175	asm {
2176	fnclex;
2177	}
2178	} else {
2179	/* SPARC:
2180	int tmpval;
2181	asm { st %fsr, tmpval; }
2182	tmpval &=0xFFFF_FC00;
2183	asm { ld tmpval, %fsr; }
2184	*/
2185	assert(0, "Not yet supported");
2186	}
2187	}
2188	public:
2189	/// The result cannot be represented exactly, so rounding occured.
2190	/// (example: x = sin(0.1); )
2191	bool inexact() { return (flags & INEXACT_MASK) != 0; }
2192	/// A zero was generated by underflow (example: x = real.min*real.epsilon/2;)
2193	bool underflow() { return (flags & UNDERFLOW_MASK) != 0; }
2194	/// An infinity was generated by overflow (example: x = real.max*2;)
2195	bool overflow() { return (flags & OVERFLOW_MASK) != 0; }
2196	/// An infinity was generated by division by zero (example: x = 3/0.0; )
2197	bool divByZero() { return (flags & DIVBYZERO_MASK) != 0; }
2198	/// A machine NaN was generated. (example: x = real.infinity * 0.0; )
2199	bool invalid() { return (flags & INVALID_MASK) != 0; }
2200	}
2201
2202
2203	/// Set all of the floating-point status flags to false.
2204	void resetIeeeFlags() { IeeeFlags.resetIeeeFlags; }
2205
2206	/// Return a snapshot of the current state of the floating-point status flags.
2207	IeeeFlags ieeeFlags()
2208	{
2209	return IeeeFlags(IeeeFlags.getIeeeFlags());
2210	}
2211
2212	/** Control the Floating point hardware
2213
2214	Change the IEEE754 floating-point rounding mode and the floating-point
2215	hardware exceptions.
2216
2217	By default, the rounding mode is roundToNearest and all hardware exceptions
2218	are disabled. For most applications, debugging is easier if the $(I division
2219	by zero), $(I overflow), and $(I invalid operation) exceptions are enabled.
2220	These three are combined into a $(I severeExceptions) value for convenience.
2221	Note in particular that if $(I invalidException) is enabled, a hardware trap
2222	will be generated whenever an uninitialized floating-point variable is used.
2223
2224	All changes are temporary. The previous state is restored at the
2225	end of the scope.
2226
2227
2228	Example:
2229	----
2230	{
2231	// Enable hardware exceptions for division by zero, overflow to infinity,
2232	// invalid operations, and uninitialized floating-point variables.
2233
2234	FloatingPointControl fpctrl;
2235	fpctrl.enableExceptions(FloatingPointControl.severeExceptions);
2236
2237	double y = x*3.0; // will generate a hardware exception, if x is uninitialized.
2238	//
2239	fpctrl.rounding = FloatingPointControl.roundUp;
2240
2241	// The hardware exceptions will be disabled when leaving this scope.
2242	// The original rounding mode will also be restored.
2243	}
2244
2245	----
2246
2247	*/
2248	struct FloatingPointControl
2249	{
2250	alias uint RoundingMode;
2251
2252	/** IEEE rounding modes.
2253	* The default mode is roundToNearest.
2254	*/
2255	enum : RoundingMode
2256	{
2257	roundToNearest = 0x0000,
2258	roundDown = 0x0400,
2259	roundUp = 0x0800,
2260	roundToZero = 0x0C00
2261	};
2262	/** IEEE hardware exceptions.
2263	* By default, all exceptions are masked (disabled).
2264	*/
2265	enum : uint
2266	{
2267	inexactException = 0x20,
2268	underflowException = 0x10,
2269	overflowException = 0x08,
2270	divByZeroException = 0x04,
2271	invalidException = 0x01,
2272	/// Severe = The overflow, division by zero, and invalid exceptions.
2273	severeExceptions = overflowException \| divByZeroException
2274	\| invalidException,
2275	allExceptions = severeExceptions \| underflowException
2276	\| inexactException,
2277	};
2278	private:
2279	enum ushort EXCEPTION_MASK = 0x3F;
2280	enum ushort ROUNDING_MASK = 0xC00;
2281	public:
2282	/// Enable (unmask) specific hardware exceptions. Multiple exceptions may be ORed together.
2283	void enableExceptions(uint exceptions)
2284	{
2285	initialize();
2286	setControlState(getControlState() & ~(exceptions & EXCEPTION_MASK));
2287	}
2288	/// Disable (mask) specific hardware exceptions. Multiple exceptions may be ORed together.
2289	void disableExceptions(uint exceptions)
2290	{
2291	initialize();
2292	setControlState(getControlState() \| (exceptions & EXCEPTION_MASK));
2293	}
2294	//// Change the floating-point hardware rounding mode
2295	void rounding(RoundingMode newMode)
2296	{
2297	ushort old = getControlState();
2298	setControlState((old & ~ROUNDING_MASK) \| (newMode & ROUNDING_MASK));
2299	}
2300	/// Return the exceptions which are currently enabled (unmasked)
2301	static uint enabledExceptions()
2302	{
2303	return (getControlState() & EXCEPTION_MASK) ^ EXCEPTION_MASK;
2304	}
2305	/// Return the currently active rounding mode
2306	static RoundingMode rounding()
2307	{
2308	return cast(RoundingMode)(getControlState() & ROUNDING_MASK);
2309	}
2310	/// Clear all pending exceptions, then restore the original exception state and rounding mode.
2311	~this()
2312	{
2313	clearExceptions();
2314	setControlState(savedState);
2315	}
2316	private:
2317	ushort savedState;
2318
2319	bool initialized=false;
2320	void initialize()
2321	{
2322	// BUG: This works around the absence of this() constructors.
2323	if (initialized) return;
2324	clearExceptions();
2325	savedState = getControlState();
2326	initialized=true;
2327	}
2328	// Clear all pending exceptions
2329	static void clearExceptions()
2330	{
2331	version (InlineAsm_X86_Any)
2332	{
2333	asm
2334	{
2335	fclex;
2336	}
2337	}
2338	else
2339	assert(0, "Not yet supported");
2340	}
2341	// Read from the control register
2342	static ushort getControlState()
2343	{
2344	version (D_InlineAsm_X86)
2345	{
2346	short cont;
2347	asm
2348	{
2349	xor EAX, EAX;
2350	fstcw cont;
2351	}
2352	return cont;
2353	}
2354	else
2355	version (D_InlineAsm_X86_64)
2356	{
2357	short cont;
2358	asm
2359	{
2360	xor RAX, RAX;
2361	fstcw cont;
2362	}
2363	return cont;
2364	}
2365	else
2366	assert(0, "Not yet supported");
2367	}
2368	// Set the control register
2369	static void setControlState(ushort newState)
2370	{
2371	version (InlineAsm_X86_Any)
2372	{
2373	asm
2374	{
2375	fclex;
2376	fldcw newState;
2377	}
2378	}
2379	else
2380	assert(0, "Not yet supported");
2381	}
2382	}
2383
2384	unittest
2385	{
2386	{
2387	FloatingPointControl ctrl;
2388	ctrl.enableExceptions(FloatingPointControl.divByZeroException
2389	\| FloatingPointControl.overflowException);
2390	assert(ctrl.enabledExceptions() ==
2391	(FloatingPointControl.divByZeroException
2392	\| FloatingPointControl.overflowException));
2393
2394	ctrl.rounding = FloatingPointControl.roundUp;
2395	assert(FloatingPointControl.rounding == FloatingPointControl.roundUp);
2396	}
2397	assert(FloatingPointControl.rounding
2398	== FloatingPointControl.roundToNearest);
2399	assert(FloatingPointControl.enabledExceptions() ==0);
2400	}
2401
2402
2403	/*********************************
2404	* Returns !=0 if e is a NaN.
2405	*/
2406
2407	bool isNaN(real x) @trusted pure nothrow
2408	{
2409	alias floatTraits!(real) F;
2410	static if (real.mant_dig==53) { // double
2411	ulong* p = cast(ulong *)&x;
2412	return ((*p & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000)
2413	&& *p & 0x000F_FFFF_FFFF_FFFF;
2414	} else static if (real.mant_dig==64) { // real80
2415	ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2416	ulong* ps = cast(ulong *)&x;
2417	return e == F.EXPMASK &&
2418	*ps & 0x7FFF_FFFF_FFFF_FFFF; // not infinity
2419	} else static if (real.mant_dig==113) { // quadruple
2420	ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2421	ulong* ps = cast(ulong *)&x;
2422	return e == F.EXPMASK &&
2423	(ps[MANTISSA_LSB] \| (ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))!=0;
2424	} else {
2425	return x!=x;
2426	}
2427	}
2428
2429
2430	unittest
2431	{
2432	assert(isNaN(float.nan));
2433	assert(isNaN(-double.nan));
2434	assert(isNaN(real.nan));
2435
2436	assert(!isNaN(53.6));
2437	assert(!isNaN(float.infinity));
2438	}
2439
2440	/*********************************
2441	* Returns !=0 if e is finite (not infinite or $(NAN)).
2442	*/
2443
2444	int isFinite(real e) @trusted pure nothrow
2445	{
2446	alias floatTraits!(real) F;
2447	ushort* pe = cast(ushort *)&e;
2448	return (pe[F.EXPPOS_SHORT] & F.EXPMASK) != F.EXPMASK;
2449	}
2450
2451	unittest
2452	{
2453	assert(isFinite(1.23));
2454	assert(!isFinite(double.infinity));
2455	assert(!isFinite(float.nan));
2456	}
2457
2458
2459	/*********************************
2460	* Returns !=0 if x is normalized (not zero, subnormal, infinite, or $(NAN)).
2461	*/
2462
2463	/* Need one for each format because subnormal floats might
2464	* be converted to normal reals.
2465	*/
2466
2467	int isNormal(X)(X x) @trusted pure nothrow
2468	{
2469	alias floatTraits!(X) F;
2470
2471	static if(real.mant_dig==106) { // doubledouble
2472	// doubledouble is normal if the least significant part is normal.
2473	return isNormal((cast(double*)&x)[MANTISSA_LSB]);
2474	} else {
2475	ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2476	return (e != F.EXPMASK && e!=0);
2477	}
2478	}
2479
2480
2481	unittest
2482	{
2483	float f = 3;
2484	double d = 500;
2485	real e = 10e+48;
2486
2487	assert(isNormal(f));
2488	assert(isNormal(d));
2489	assert(isNormal(e));
2490	f = d = e = 0;
2491	assert(!isNormal(f));
2492	assert(!isNormal(d));
2493	assert(!isNormal(e));
2494	assert(!isNormal(real.infinity));
2495	assert(isNormal(-real.max));
2496	assert(!isNormal(real.min_normal/4));
2497
2498	}
2499
2500	/*********************************
2501	* Is number subnormal? (Also called "denormal".)
2502	* Subnormals have a 0 exponent and a 0 most significant mantissa bit.
2503	*/
2504
2505	/* Need one for each format because subnormal floats might
2506	* be converted to normal reals.
2507	*/
2508
2509	int isSubnormal(float f) @trusted pure nothrow
2510	{
2511	uint p = cast(uint )&f;
2512	return (p & 0x7F80_0000) == 0 && p & 0x007F_FFFF;
2513	}
2514
2515	unittest
2516	{
2517	float f = 3.0;
2518
2519	for (f = 1.0; !isSubnormal(f); f /= 2)
2520	assert(f != 0);
2521	}
2522
2523	/// ditto
2524
2525	int isSubnormal(double d) @trusted pure nothrow
2526	{
2527	uint p = cast(uint )&d;
2528	return (p[MANTISSA_MSB] & 0x7FF0_0000) == 0
2529	&& (p[MANTISSA_LSB] \|\| p[MANTISSA_MSB] & 0x000F_FFFF);
2530	}
2531
2532	unittest
2533	{
2534	double f;
2535
2536	for (f = 1; !isSubnormal(f); f /= 2)
2537	assert(f != 0);
2538	}
2539
2540	/// ditto
2541
2542	int isSubnormal(real x) @trusted pure nothrow
2543	{
2544	alias floatTraits!(real) F;
2545	static if (real.mant_dig == 53) { // double
2546	return isSubnormal(cast(double)x);
2547	} else static if (real.mant_dig == 113) { // quadruple
2548	ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2549	long* ps = cast(long *)&x;
2550	return (e == 0 &&
2551	(((ps[MANTISSA_LSB]\|(ps[MANTISSA_MSB]& 0x0000_FFFF_FFFF_FFFF))) !=0));
2552	} else static if (real.mant_dig==64) { // real80
2553	ushort* pe = cast(ushort *)&x;
2554	long* ps = cast(long *)&x;
2555
2556	return (pe[F.EXPPOS_SHORT] & F.EXPMASK) == 0 && *ps > 0;
2557	} else { // double double
2558	return isSubnormal((cast(double*)&x)[MANTISSA_MSB]);
2559	}
2560	}
2561
2562	unittest
2563	{
2564	real f;
2565
2566	for (f = 1; !isSubnormal(f); f /= 2)
2567	assert(f != 0);
2568	}
2569
2570	/*********************************
2571	* Return !=0 if e is $(PLUSMN)$(INFIN).
2572	*/
2573
2574	bool isInfinity(real x) @trusted pure nothrow
2575	{
2576	alias floatTraits!(real) F;
2577	static if (real.mant_dig == 53) { // double
2578	return ((cast(ulong )&x) & 0x7FFF_FFFF_FFFF_FFFF)
2579	== 0x7FF8_0000_0000_0000;
2580	} else static if(real.mant_dig == 106) { //doubledouble
2581	return (((cast(ulong *)&x)[MANTISSA_MSB]) & 0x7FFF_FFFF_FFFF_FFFF)
2582	== 0x7FF8_0000_0000_0000;
2583	} else static if (real.mant_dig == 113) { // quadruple
2584	long* ps = cast(long *)&x;
2585	return (ps[MANTISSA_LSB] == 0)
2586	&& (ps[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_0000_0000_0000;
2587	} else { // real80
2588	ushort e = cast(ushort)(F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT]);
2589	ulong* ps = cast(ulong *)&x;
2590	// On Motorola 68K, infinity can have hidden bit=1 or 0. On x86, it is always 1.
2591	return e == F.EXPMASK && (*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0;
2592	}
2593	}
2594
2595	unittest
2596	{
2597	assert(isInfinity(float.infinity));
2598	assert(!isInfinity(float.nan));
2599	assert(isInfinity(double.infinity));
2600	assert(isInfinity(-real.infinity));
2601
2602	assert(isInfinity(-1.0 / 0.0));
2603	}
2604
2605	/*********************************
2606	* Is the binary representation of x identical to y?
2607	*
2608	* Same as ==, except that positive and negative zero are not identical,
2609	* and two $(NAN)s are identical if they have the same 'payload'.
2610	*/
2611
2612	bool isIdentical(real x, real y) @trusted pure nothrow
2613	{
2614	// We're doing a bitwise comparison so the endianness is irrelevant.
2615	long* pxs = cast(long *)&x;
2616	long* pys = cast(long *)&y;
2617	static if (real.mant_dig == 53)
2618	{ //double
2619	return pxs[0] == pys[0];
2620	}
2621	else static if (real.mant_dig == 113 \|\| real.mant_dig==106)
2622	{
2623	// quadruple or doubledouble
2624	return pxs[0] == pys[0] && pxs[1] == pys[1];
2625	}
2626	else
2627	{ // real80
2628	ushort* pxe = cast(ushort *)&x;
2629	ushort* pye = cast(ushort *)&y;
2630	return pxe[4] == pye[4] && pxs[0] == pys[0];
2631	}
2632	}
2633
2634	/*********************************
2635	* Return 1 if sign bit of e is set, 0 if not.
2636	*/
2637
2638	int signbit(real x) @trusted pure nothrow
2639	{
2640	return ((cast(ubyte *)&x)[floatTraits!(real).SIGNPOS_BYTE] & 0x80) != 0;
2641	}
2642
2643	unittest
2644	{
2645	debug (math) printf("math.signbit.unittest\n");
2646	assert(!signbit(float.nan));
2647	assert(signbit(-float.nan));
2648	assert(!signbit(168.1234));
2649	assert(signbit(-168.1234));
2650	assert(!signbit(0.0));
2651	assert(signbit(-0.0));
2652	assert(signbit(-double.max));
2653	assert(!signbit(double.max));
2654	}
2655
2656	/*********************************
2657	* Return a value composed of to with from's sign bit.
2658	*/
2659
2660	real copysign(real to, real from) @trusted pure nothrow
2661	{
2662	ubyte* pto = cast(ubyte *)&to;
2663	const ubyte* pfrom = cast(ubyte *)&from;
2664
2665	alias floatTraits!(real) F;
2666	pto[F.SIGNPOS_BYTE] &= 0x7F;
2667	pto[F.SIGNPOS_BYTE] \|= pfrom[F.SIGNPOS_BYTE] & 0x80;
2668	return to;
2669	}
2670
2671	unittest
2672	{
2673	real e;
2674
2675	e = copysign(21, 23.8);
2676	assert(e == 21);
2677
2678	e = copysign(-21, 23.8);
2679	assert(e == 21);
2680
2681	e = copysign(21, -23.8);
2682	assert(e == -21);
2683
2684	e = copysign(-21, -23.8);
2685	assert(e == -21);
2686
2687	e = copysign(real.nan, -23.8);
2688	assert(isNaN(e) && signbit(e));
2689	}
2690
2691	/*********************************
2692	Returns $(D -1) if $(D x < 0), $(D x) if $(D x == 0), $(D 1) if
2693	$(D x > 0), and $(NAN) if x==$(NAN).
2694	*/
2695	F sgn(F)(F x) @safe pure nothrow
2696	{
2697	// @@@TODO@@@: make this faster
2698	return x > 0 ? 1 : x < 0 ? -1 : x;
2699	}
2700
2701	unittest
2702	{
2703	debug (math) printf("math.sgn.unittest\n");
2704	assert(sgn(168.1234) == 1);
2705	assert(sgn(-168.1234) == -1);
2706	assert(sgn(0.0) == 0);
2707	assert(sgn(-0.0) == 0);
2708	}
2709
2710	// Functions for NaN payloads
2711	/*
2712	* A 'payload' can be stored in the significand of a $(NAN). One bit is required
2713	* to distinguish between a quiet and a signalling $(NAN). This leaves 22 bits
2714	* of payload for a float; 51 bits for a double; 62 bits for an 80-bit real;
2715	* and 111 bits for a 128-bit quad.
2716	*/
2717	/**
2718	* Create a quiet $(NAN), storing an integer inside the payload.
2719	*
2720	* For floats, the largest possible payload is 0x3F_FFFF.
2721	* For doubles, it is 0x3_FFFF_FFFF_FFFF.
2722	* For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
2723	*/
2724	real NaN(ulong payload) @trusted pure nothrow
2725	{
2726	static if (real.mant_dig == 64) { //real80
2727	ulong v = 3; // implied bit = 1, quiet bit = 1
2728	} else {
2729	ulong v = 2; // no implied bit. quiet bit = 1
2730	}
2731
2732	ulong a = payload;
2733
2734	// 22 Float bits
2735	ulong w = a & 0x3F_FFFF;
2736	a -= w;
2737
2738	v <<=22;
2739	v \|= w;
2740	a >>=22;
2741
2742	// 29 Double bits
2743	v <<=29;
2744	w = a & 0xFFF_FFFF;
2745	v \|= w;
2746	a -= w;
2747	a >>=29;
2748
2749	static if (real.mant_dig == 53) { // double
2750	v \|=0x7FF0_0000_0000_0000;
2751	real x;
2752	* cast(ulong *)(&x) = v;
2753	return x;
2754	} else {
2755	v <<=11;
2756	a &= 0x7FF;
2757	v \|= a;
2758	real x = real.nan;
2759	// Extended real bits
2760	static if (real.mant_dig==113) { //quadruple
2761	v<<=1; // there's no implicit bit
2762	version(LittleEndian) {
2763	cast(ulong)(6+cast(ubyte*)(&x)) = v;
2764	} else {
2765	cast(ulong)(2+cast(ubyte*)(&x)) = v;
2766	}
2767	} else { // real80
2768	* cast(ulong *)(&x) = v;
2769	}
2770	return x;
2771	}
2772	}
2773
2774	/**
2775	* Extract an integral payload from a $(NAN).
2776	*
2777	* Returns:
2778	* the integer payload as a ulong.
2779	*
2780	* For floats, the largest possible payload is 0x3F_FFFF.
2781	* For doubles, it is 0x3_FFFF_FFFF_FFFF.
2782	* For 80-bit or 128-bit reals, it is 0x3FFF_FFFF_FFFF_FFFF.
2783	*/
2784	ulong getNaNPayload(real x) @trusted pure nothrow
2785	{
2786	// assert(isNaN(x));
2787	static if (real.mant_dig == 53) {
2788	ulong m = cast(ulong )(&x);
2789	// Make it look like an 80-bit significand.
2790	// Skip exponent, and quiet bit
2791	m &= 0x0007_FFFF_FFFF_FFFF;
2792	m <<= 10;
2793	} else static if (real.mant_dig==113) { // quadruple
2794	version(LittleEndian) {
2795	ulong m = cast(ulong)(6+cast(ubyte*)(&x));
2796	} else {
2797	ulong m = cast(ulong)(2+cast(ubyte*)(&x));
2798	}
2799	m>>=1; // there's no implicit bit
2800	} else {
2801	ulong m = cast(ulong )(&x);
2802	}
2803	// ignore implicit bit and quiet bit
2804	ulong f = m & 0x3FFF_FF00_0000_0000L;
2805	ulong w = f >>> 40;
2806	w \|= (m & 0x00FF_FFFF_F800L) << (22 - 11);
2807	w \|= (m & 0x7FF) << 51;
2808	return w;
2809	}
2810
2811	debug(UnitTest) {
2812	unittest {
2813	real nan4 = NaN(0x789_ABCD_EF12_3456);
2814	static if (real.mant_dig == 64 \|\| real.mant_dig==113) {
2815	assert (getNaNPayload(nan4) == 0x789_ABCD_EF12_3456);
2816	} else {
2817	assert (getNaNPayload(nan4) == 0x1_ABCD_EF12_3456);
2818	}
2819	double nan5 = nan4;
2820	assert (getNaNPayload(nan5) == 0x1_ABCD_EF12_3456);
2821	float nan6 = nan4;
2822	assert (getNaNPayload(nan6) == 0x12_3456);
2823	nan4 = NaN(0xFABCD);
2824	assert (getNaNPayload(nan4) == 0xFABCD);
2825	nan6 = nan4;
2826	assert (getNaNPayload(nan6) == 0xFABCD);
2827	nan5 = NaN(0x100_0000_0000_3456);
2828	assert(getNaNPayload(nan5) == 0x0000_0000_3456);
2829	}
2830	}
2831
2832	/**
2833	* Calculate the next largest floating point value after x.
2834	*
2835	* Return the least number greater than x that is representable as a real;
2836	* thus, it gives the next point on the IEEE number line.
2837	*
2838	* $(TABLE_SV
2839	* $(SVH x, nextUp(x) )
2840	* $(SV -$(INFIN), -real.max )
2841	* $(SV $(PLUSMN)0.0, real.min_normal*real.epsilon )
2842	* $(SV real.max, $(INFIN) )
2843	* $(SV $(INFIN), $(INFIN) )
2844	* $(SV $(NAN), $(NAN) )
2845	* )
2846	*/
2847	real nextUp(real x) @trusted pure nothrow
2848	{
2849	alias floatTraits!(real) F;
2850	static if (real.mant_dig == 53) { // double
2851	return nextUp(cast(double)x);
2852	} else static if(real.mant_dig==113) { // quadruple
2853	ushort e = F.EXPMASK & (cast(ushort *)&x)[F.EXPPOS_SHORT];
2854	if (e == F.EXPMASK) { // NaN or Infinity
2855	if (x == -real.infinity) return -real.max;
2856	return x; // +Inf and NaN are unchanged.
2857	}
2858	ulong* ps = cast(ulong *)&e;
2859	if (ps[MANTISSA_LSB] & 0x8000_0000_0000_0000) { // Negative number
2860	if (ps[MANTISSA_LSB] == 0
2861	&& ps[MANTISSA_MSB] == 0x8000_0000_0000_0000) {
2862	// it was negative zero, change to smallest subnormal
2863	ps[MANTISSA_LSB] = 0x0000_0000_0000_0001;
2864	ps[MANTISSA_MSB] = 0;
2865	return x;
2866	}
2867	--*ps;
2868	if (ps[MANTISSA_LSB]==0) --ps[MANTISSA_MSB];
2869	} else { // Positive number
2870	++ps[MANTISSA_LSB];
2871	if (ps[MANTISSA_LSB]==0) ++ps[MANTISSA_MSB];
2872	}
2873	return x;
2874
2875	} else static if(real.mant_dig==64){ // real80
2876	// For 80-bit reals, the "implied bit" is a nuisance...
2877	ushort pe = cast(ushort )&x;
2878	ulong ps = cast(ulong )&x;
2879
2880	if ((pe[F.EXPPOS_SHORT] & F.EXPMASK) == F.EXPMASK) {
2881	// First, deal with NANs and infinity
2882	if (x == -real.infinity) return -real.max;
2883	return x; // +Inf and NaN are unchanged.
2884	}
2885	if (pe[F.EXPPOS_SHORT] & 0x8000) {
2886	// Negative number -- need to decrease the significand
2887	--*ps;
2888	// Need to mask with 0x7FFF... so subnormals are treated correctly.
2889	if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0x7FFF_FFFF_FFFF_FFFF) {
2890	if (pe[F.EXPPOS_SHORT] == 0x8000) { // it was negative zero
2891	*ps = 1;
2892	pe[F.EXPPOS_SHORT] = 0; // smallest subnormal.
2893	return x;
2894	}
2895	--pe[F.EXPPOS_SHORT];
2896	if (pe[F.EXPPOS_SHORT] == 0x8000) {
2897	return x; // it's become a subnormal, implied bit stays low.
2898	}
2899	*ps = 0xFFFF_FFFF_FFFF_FFFF; // set the implied bit
2900	return x;
2901	}
2902	return x;
2903	} else {
2904	// Positive number -- need to increase the significand.
2905	// Works automatically for positive zero.
2906	++*ps;
2907	if ((*ps & 0x7FFF_FFFF_FFFF_FFFF) == 0) {
2908	// change in exponent
2909	++pe[F.EXPPOS_SHORT];
2910	*ps = 0x8000_0000_0000_0000; // set the high bit
2911	}
2912	}
2913	return x;
2914	} // doubledouble is not supported
2915	}
2916
2917	/** ditto */
2918	double nextUp(double x) @trusted pure nothrow
2919	{
2920	ulong ps = cast(ulong )&x;
2921
2922	if ((*ps & 0x7FF0_0000_0000_0000) == 0x7FF0_0000_0000_0000) {
2923	// First, deal with NANs and infinity
2924	if (x == -x.infinity) return -x.max;
2925	return x; // +INF and NAN are unchanged.
2926	}
2927	if (*ps & 0x8000_0000_0000_0000) { // Negative number
2928	if (*ps == 0x8000_0000_0000_0000) { // it was negative zero
2929	*ps = 0x0000_0000_0000_0001; // change to smallest subnormal
2930	return x;
2931	}
2932	--*ps;
2933	} else { // Positive number
2934	++*ps;
2935	}
2936	return x;
2937	}
2938
2939	/** ditto */
2940	float nextUp(float x) @trusted pure nothrow
2941	{
2942	uint ps = cast(uint )&x;
2943
2944	if ((*ps & 0x7F80_0000) == 0x7F80_0000) {
2945	// First, deal with NANs and infinity
2946	if (x == -x.infinity) return -x.max;
2947	return x; // +INF and NAN are unchanged.
2948	}
2949	if (*ps & 0x8000_0000) { // Negative number
2950	if (*ps == 0x8000_0000) { // it was negative zero
2951	*ps = 0x0000_0001; // change to smallest subnormal
2952	return x;
2953	}
2954	--*ps;
2955	} else { // Positive number
2956	++*ps;
2957	}
2958	return x;
2959	}
2960
2961	/**
2962	* Calculate the next smallest floating point value before x.
2963	*
2964	* Return the greatest number less than x that is representable as a real;
2965	* thus, it gives the previous point on the IEEE number line.
2966	*
2967	* $(TABLE_SV
2968	* $(SVH x, nextDown(x) )
2969	* $(SV $(INFIN), real.max )
2970	* $(SV $(PLUSMN)0.0, -real.min_normal*real.epsilon )
2971	* $(SV -real.max, -$(INFIN) )
2972	* $(SV -$(INFIN), -$(INFIN) )
2973	* $(SV $(NAN), $(NAN) )
2974	* )
2975	*/
2976	real nextDown(real x) @safe pure nothrow
2977	{
2978	return -nextUp(-x);
2979	}
2980
2981	/** ditto */
2982	double nextDown(double x) @safe pure nothrow
2983	{
2984	return -nextUp(-x);
2985	}
2986
2987	/** ditto */
2988	float nextDown(float x) @safe pure nothrow
2989	{
2990	return -nextUp(-x);
2991	}
2992
2993	unittest {
2994	assert( nextDown(1.0 + real.epsilon) == 1.0);
2995	}
2996
2997	unittest {
2998	static if (real.mant_dig == 64) {
2999
3000	// Tests for 80-bit reals
3001	assert(isIdentical(nextUp(NaN(0xABC)), NaN(0xABC)));
3002	// negative numbers
3003	assert( nextUp(-real.infinity) == -real.max );
3004	assert( nextUp(-1.0L-real.epsilon) == -1.0 );
3005	assert( nextUp(-2.0L) == -2.0 + real.epsilon);
3006	// denormals and zero
3007	assert( nextUp(-real.min_normal) == -real.min_normal*(1-real.epsilon) );
3008	assert( nextUp(-real.min_normal(1-real.epsilon)) == -real.min_normal(1-2*real.epsilon) );
3009	assert( isIdentical(-0.0L, nextUp(-real.min_normal*real.epsilon)) );
3010	assert( nextUp(-0.0L) == real.min_normal*real.epsilon );
3011	assert( nextUp(0.0L) == real.min_normal*real.epsilon );
3012	assert( nextUp(real.min_normal*(1-real.epsilon)) == real.min_normal );
3013	assert( nextUp(real.min_normal) == real.min_normal*(1+real.epsilon) );
3014	// positive numbers
3015	assert( nextUp(1.0L) == 1.0 + real.epsilon );
3016	assert( nextUp(2.0L-real.epsilon) == 2.0 );
3017	assert( nextUp(real.max) == real.infinity );
3018	assert( nextUp(real.infinity)==real.infinity );
3019	}
3020
3021	double n = NaN(0xABC);
3022	assert(isIdentical(nextUp(n), n));
3023	// negative numbers
3024	assert( nextUp(-double.infinity) == -double.max );
3025	assert( nextUp(-1-double.epsilon) == -1.0 );
3026	assert( nextUp(-2.0) == -2.0 + double.epsilon);
3027	// denormals and zero
3028
3029	assert( nextUp(-double.min_normal) == -double.min_normal*(1-double.epsilon) );
3030	assert( nextUp(-double.min_normal(1-double.epsilon)) == -double.min_normal(1-2*double.epsilon) );
3031	assert( isIdentical(-0.0, nextUp(-double.min_normal*double.epsilon)) );
3032	assert( nextUp(0.0) == double.min_normal*double.epsilon );
3033	assert( nextUp(-0.0) == double.min_normal*double.epsilon );
3034	assert( nextUp(double.min_normal*(1-double.epsilon)) == double.min_normal );
3035	assert( nextUp(double.min_normal) == double.min_normal*(1+double.epsilon) );
3036	// positive numbers
3037	assert( nextUp(1.0) == 1.0 + double.epsilon );
3038	assert( nextUp(2.0-double.epsilon) == 2.0 );
3039	assert( nextUp(double.max) == double.infinity );
3040
3041	float fn = NaN(0xABC);
3042	assert(isIdentical(nextUp(fn), fn));
3043	float f = -float.min_normal*(1-float.epsilon);
3044	float f1 = -float.min_normal;
3045	assert( nextUp(f1) == f);
3046	f = 1.0f+float.epsilon;
3047	f1 = 1.0f;
3048	assert( nextUp(f1) == f );
3049	f1 = -0.0f;
3050	assert( nextUp(f1) == float.min_normal*float.epsilon);
3051	assert( nextUp(float.infinity)==float.infinity );
3052
3053	assert(nextDown(1.0L+real.epsilon)==1.0);
3054	assert(nextDown(1.0+double.epsilon)==1.0);
3055	f = 1.0f+float.epsilon;
3056	assert(nextDown(f)==1.0);
3057	assert(nextafter(1.0+real.epsilon, -real.infinity)==1.0);
3058	}
3059
3060
3061
3062	/******************************************
3063	* Calculates the next representable value after x in the direction of y.
3064	*
3065	* If y > x, the result will be the next largest floating-point value;
3066	* if y < x, the result will be the next smallest value.
3067	* If x == y, the result is y.
3068	*
3069	* Remarks:
3070	* This function is not generally very useful; it's almost always better to use
3071	* the faster functions nextUp() or nextDown() instead.
3072	*
3073	* The FE_INEXACT and FE_OVERFLOW exceptions will be raised if x is finite and
3074	* the function result is infinite. The FE_INEXACT and FE_UNDERFLOW
3075	* exceptions will be raised if the function value is subnormal, and x is
3076	* not equal to y.
3077	*/
3078	T nextafter(T)(T x, T y) @safe pure nothrow
3079	{
3080	if (x==y) return y;
3081	return ((y>x) ? nextUp(x) : nextDown(x));
3082	}
3083
3084	unittest
3085	{
3086	float a = 1;
3087	assert(is(typeof(nextafter(a, a)) == float));
3088	assert(nextafter(a, a.infinity) > a);
3089
3090	double b = 2;
3091	assert(is(typeof(nextafter(b, b)) == double));
3092	assert(nextafter(b, b.infinity) > b);
3093
3094	real c = 3;
3095	assert(is(typeof(nextafter(c, c)) == real));
3096	assert(nextafter(c, c.infinity) > c);
3097	}
3098
3099	//real nexttoward(real x, real y) { return core.stdc.math.nexttowardl(x, y); }
3100
3101	/*******************************************
3102	* Returns the positive difference between x and y.
3103	* Returns:
3104	* $(TABLE_SV
3105	* $(TR $(TH x, y) $(TH fdim(x, y)))
3106	* $(TR $(TD x $(GT) y) $(TD x - y))
3107	* $(TR $(TD x $(LT)= y) $(TD +0.0))
3108	* )
3109	*/
3110	real fdim(real x, real y) @safe pure nothrow { return (x > y) ? x - y : +0.0; }
3111
3112	/****************************************
3113	* Returns the larger of x and y.
3114	*/
3115	real fmax(real x, real y) @safe pure nothrow { return x > y ? x : y; }
3116
3117	/****************************************
3118	* Returns the smaller of x and y.
3119	*/
3120	real fmin(real x, real y) @safe pure nothrow { return x < y ? x : y; }
3121
3122	/**************************************
3123	* Returns (x * y) + z, rounding only once according to the
3124	* current rounding mode.
3125	*
3126	* BUGS: Not currently implemented - rounds twice.
3127	*/
3128	real fma(real x, real y, real z) @safe pure nothrow { return (x * y) + z; }
3129
3130	/*******************************************************************
3131	* Compute the value of x $(SUP n), where n is an integer
3132	*/
3133	Unqual!F pow(F, G)(F x, G n) @trusted pure nothrow
3134	if (isFloatingPoint!(F) && isIntegral!(G))
3135	{
3136	real p = 1.0, v = void;
3137	Unsigned!(Unqual!G) m = n;
3138	if (n < 0)
3139	{
3140	switch (n)
3141	{
3142	case -1:
3143	return 1 / x;
3144	case -2:
3145	return 1 / (x * x);
3146	default:
3147	}
3148
3149	m = -n;
3150	v = p / x;
3151	}
3152	else
3153	{
3154	switch (n)
3155	{
3156	case 0:
3157	return 1.0;
3158	case 1:
3159	return x;
3160	case 2:
3161	return x * x;
3162	default:
3163	}
3164
3165	v = x;
3166	}
3167
3168	while (1)
3169	{
3170	if (m & 1)
3171	p *= v;
3172	m >>= 1;
3173	if (!m)
3174	break;
3175	v *= v;
3176	}
3177	return p;
3178	}
3179
3180	unittest
3181	{
3182	// Make sure it instantiates and works properly on immutable values and
3183	// with various integer and float types.
3184	immutable real x = 46;
3185	immutable float xf = x;
3186	immutable double xd = x;
3187	immutable uint one = 1;
3188	immutable ushort two = 2;
3189	immutable ubyte three = 3;
3190	immutable ulong eight = 8;
3191
3192	immutable int neg1 = -1;
3193	immutable short neg2 = -2;
3194	immutable byte neg3 = -3;
3195	immutable long neg8 = -8;
3196
3197
3198	assert(pow(x,0) == 1.0);
3199	assert(pow(xd,one) == x);
3200	assert(pow(xf,two) == x * x);
3201	assert(pow(x,three) == x * x * x);
3202	assert(pow(x,eight) == (x * x) * (x * x) * (x * x) * (x * x));
3203
3204	assert(pow(x, neg1) == 1 / x);
3205	assert(pow(xd, neg2) == 1 / (x * x));
3206	assert(pow(x, neg3) == 1 / (x * x * x));
3207	assert(pow(xf, neg8) == 1 / ((x * x) * (x * x) * (x * x) * (x * x)));
3208	}
3209
3210	/** Compute the value of an integer x, raised to the power of a positive
3211	* integer n.
3212	*
3213	* If both x and n are 0, the result is 1.
3214	* If n is negative, an integer divide error will occur at runtime,
3215	* regardless of the value of x.
3216	*/
3217
3218	typeof(Unqual!(F).init * Unqual!(G).init) pow(F, G)(F x, G n) @trusted pure nothrow
3219	if (isIntegral!(F) && isIntegral!(G))
3220	{
3221	if (n<0) return x/0; // Only support positive powers
3222	typeof(return) p, v = void;
3223	Unqual!G m = n;
3224
3225	switch (m)
3226	{
3227	case 0:
3228	p = 1;
3229	break;
3230
3231	case 1:
3232	p = x;
3233	break;
3234
3235	case 2:
3236	p = x * x;
3237	break;
3238
3239	default:
3240	v = x;
3241	p = 1;
3242	while (1){
3243	if (m & 1)
3244	p *= v;
3245	m >>= 1;
3246	if (!m)
3247	break;
3248	v *= v;
3249	}
3250	break;
3251	}
3252	return p;
3253	}
3254
3255	unittest
3256	{
3257	immutable int one = 1;
3258	immutable byte two = 2;
3259	immutable ubyte three = 3;
3260	immutable short four = 4;
3261	immutable long ten = 10;
3262
3263	assert(pow(two, three) == 8);
3264	assert(pow(two, ten) == 1024);
3265	assert(pow(one, ten) == 1);
3266	assert(pow(ten, four) == 10_000);
3267	assert(pow(four, 10) == 1_048_576);
3268	assert(pow(three, four) == 81);
3269
3270	}
3271
3272	/*Computes integer to floating point powers./
3273	real pow(I, F)(I x, F y) @trusted pure nothrow
3274	if(isIntegral!I && isFloatingPoint!F)
3275	{
3276	return pow(cast(real) x, cast(Unqual!F) y);
3277	}
3278
3279	/*********************************************
3280	* Calculates x$(SUP y).
3281	*
3282	* $(TABLE_SV
3283	* $(TR $(TH x) $(TH y) $(TH pow(x, y))
3284	* $(TH div 0) $(TH invalid?))
3285	* $(TR $(TD anything) $(TD $(PLUSMN)0.0) $(TD 1.0)
3286	* $(TD no) $(TD no) )
3287	* $(TR $(TD \|x\| $(GT) 1) $(TD +$(INFIN)) $(TD +$(INFIN))
3288	* $(TD no) $(TD no) )
3289	* $(TR $(TD \|x\| $(LT) 1) $(TD +$(INFIN)) $(TD +0.0)
3290	* $(TD no) $(TD no) )
3291	* $(TR $(TD \|x\| $(GT) 1) $(TD -$(INFIN)) $(TD +0.0)
3292	* $(TD no) $(TD no) )
3293	* $(TR $(TD \|x\| $(LT) 1) $(TD -$(INFIN)) $(TD +$(INFIN))
3294	* $(TD no) $(TD no) )
3295	* $(TR $(TD +$(INFIN)) $(TD $(GT) 0.0) $(TD +$(INFIN))
3296	* $(TD no) $(TD no) )
3297	* $(TR $(TD +$(INFIN)) $(TD $(LT) 0.0) $(TD +0.0)
3298	* $(TD no) $(TD no) )
3299	* $(TR $(TD -$(INFIN)) $(TD odd integer $(GT) 0.0) $(TD -$(INFIN))
3300	* $(TD no) $(TD no) )
3301	* $(TR $(TD -$(INFIN)) $(TD $(GT) 0.0, not odd integer) $(TD +$(INFIN))
3302	* $(TD no) $(TD no))
3303	* $(TR $(TD -$(INFIN)) $(TD odd integer $(LT) 0.0) $(TD -0.0)
3304	* $(TD no) $(TD no) )
3305	* $(TR $(TD -$(INFIN)) $(TD $(LT) 0.0, not odd integer) $(TD +0.0)
3306	* $(TD no) $(TD no) )
3307	* $(TR $(TD $(PLUSMN)1.0) $(TD $(PLUSMN)$(INFIN)) $(TD $(NAN))
3308	* $(TD no) $(TD yes) )
3309	* $(TR $(TD $(LT) 0.0) $(TD finite, nonintegral) $(TD $(NAN))
3310	* $(TD no) $(TD yes))
3311	* $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(LT) 0.0) $(TD $(PLUSMNINF))
3312	* $(TD yes) $(TD no) )
3313	* $(TR $(TD $(PLUSMN)0.0) $(TD $(LT) 0.0, not odd integer) $(TD +$(INFIN))
3314	* $(TD yes) $(TD no))
3315	* $(TR $(TD $(PLUSMN)0.0) $(TD odd integer $(GT) 0.0) $(TD $(PLUSMN)0.0)
3316	* $(TD no) $(TD no) )
3317	* $(TR $(TD $(PLUSMN)0.0) $(TD $(GT) 0.0, not odd integer) $(TD +0.0)
3318	* $(TD no) $(TD no) )
3319	* )
3320	*/
3321
3322	Unqual!(Largest!(F, G)) pow(F, G)(F x, G y) @trusted pure nothrow
3323	if (isFloatingPoint!(F) && isFloatingPoint!(G))
3324	{
3325	alias typeof(return) Float;
3326
3327	static real impl(real x, real y) pure nothrow
3328	{
3329	if (isNaN(y))
3330	return y;
3331
3332	if (y == 0)
3333	return 1; // even if x is $(NAN)
3334	if (isNaN(x) && y != 0)
3335	return x;
3336	if (isInfinity(y))
3337	{
3338	if (fabs(x) > 1)
3339	{
3340	if (signbit(y))
3341	return +0.0;
3342	else
3343	return F.infinity;
3344	}
3345	else if (fabs(x) == 1)
3346	{
3347	return y * 0; // generate NaN.
3348	}
3349	else // < 1
3350	{
3351	if (signbit(y))
3352	return F.infinity;
3353	else
3354	return +0.0;
3355	}
3356	}
3357	if (isInfinity(x))
3358	{
3359	if (signbit(x))
3360	{ long i;
3361
3362	i = cast(long)y;
3363	if (y > 0)
3364	{
3365	if (i == y && i & 1)
3366	return -F.infinity;
3367	else
3368	return F.infinity;
3369	}
3370	else if (y < 0)
3371	{
3372	if (i == y && i & 1)
3373	return -0.0;
3374	else
3375	return +0.0;
3376	}
3377	}
3378	else
3379	{
3380	if (y > 0)
3381	return F.infinity;
3382	else if (y < 0)
3383	return +0.0;
3384	}
3385	}
3386
3387	if (x == 0.0)
3388	{
3389	if (signbit(x))
3390	{ long i;
3391
3392	i = cast(long)y;
3393	if (y > 0)
3394	{
3395	if (i == y && i & 1)
3396	return -0.0;
3397	else
3398	return +0.0;
3399	}
3400	else if (y < 0)
3401	{
3402	if (i == y && i & 1)
3403	return -F.infinity;
3404	else
3405	return F.infinity;
3406	}
3407	}
3408	else
3409	{
3410	if (y > 0)
3411	return +0.0;
3412	else if (y < 0)
3413	return F.infinity;
3414	}
3415	}
3416	double sign = 1.0;
3417	if (x < 0) {
3418	// Result is real only if y is an integer
3419	// Check for a non-zero fractional part
3420	if (y > -1.0 / real.epsilon && y < 1.0 / real.epsilon)
3421	{
3422	long w = cast(long)y;
3423	if (w != y)
3424	return sqrt(x); // Complex result -- create a NaN
3425	if (w & 1) sign = -1.0;
3426	}
3427	x = -x;
3428	}
3429	version(INLINE_YL2X) {
3430	// If x > 0, x ^^ y == 2 ^^ ( y * log2(x) )
3431	// TODO: This is not accurate in practice. A fast and accurate
3432	// (though complicated) method is described in:
3433	// "An efficient rounding boundary test for pow(x, y)
3434	// in double precision", C.Q. Lauter and V. Lefèvre, INRIA (2007).
3435	return sign * exp2( yl2x(x, y) );
3436	} else {
3437	return sign * core.stdc.math.powl(x, y);
3438	}
3439	}
3440	return impl(x, y);
3441	}
3442
3443	unittest
3444	{
3445	// Test all the special values. These unittests can be run on Windows
3446	// by temporarily changing the version(linux) to version(all).
3447	immutable float zero = 0;
3448	immutable real one = 1;
3449	immutable double two = 2;
3450	immutable float three = 3;
3451	immutable float fnan = float.nan;
3452	immutable double dnan = double.nan;
3453	immutable real rnan = real.nan;
3454	immutable dinf = double.infinity;
3455	immutable rninf = -real.infinity;
3456
3457	assert(pow(fnan, zero) == 1);
3458	assert(pow(dnan, zero) == 1);
3459	assert(pow(rnan, zero) == 1);
3460
3461	assert(pow(two, dinf) == double.infinity);
3462	assert(isIdentical(pow(0.2f, dinf), +0.0));
3463	assert(pow(0.99999999L, rninf) == real.infinity);
3464	assert(isIdentical(pow(1.000000001, rninf), +0.0));
3465	assert(pow(dinf, 0.001) == dinf);
3466	assert(isIdentical(pow(dinf, -0.001), +0.0));
3467	assert(pow(rninf, 3.0L) == rninf);
3468	assert(pow(rninf, 2.0L) == real.infinity);
3469	assert(isIdentical(pow(rninf, -3.0), -0.0));
3470	assert(isIdentical(pow(rninf, -2.0), +0.0));
3471
3472	// @@@BUG@@@ somewhere
3473	version(OSX) {} else assert(isNaN(pow(one, dinf)));
3474	version(OSX) {} else assert(isNaN(pow(-one, dinf)));
3475	assert(isNaN(pow(-0.2, PI)));
3476	// boundary cases. Note that epsilon == 2^^-n for some n,
3477	// so 1/epsilon == 2^^n is always even.
3478	assert(pow(-1.0L, 1/real.epsilon - 1.0L) == -1.0L);
3479	assert(pow(-1.0L, 1/real.epsilon) == 1.0L);
3480	assert(isNaN(pow(-1.0L, 1/real.epsilon-0.5L)));
3481	assert(isNaN(pow(-1.0L, -1/real.epsilon+0.5L)));
3482
3483	assert(pow(0.0, -3.0) == double.infinity);
3484	assert(pow(-0.0, -3.0) == -double.infinity);
3485	assert(pow(0.0, -PI) == double.infinity);
3486	assert(pow(-0.0, -PI) == double.infinity);
3487	assert(isIdentical(pow(0.0, 5.0), 0.0));
3488	assert(isIdentical(pow(-0.0, 5.0), -0.0));
3489	assert(isIdentical(pow(0.0, 6.0), 0.0));
3490	assert(isIdentical(pow(-0.0, 6.0), 0.0));
3491
3492	// Now, actual numbers.
3493	assert(approxEqual(pow(two, three), 8.0));
3494	assert(approxEqual(pow(two, -2.5), 0.1767767));
3495
3496	// Test integer to float power.
3497	immutable uint twoI = 2;
3498	assert(approxEqual(pow(twoI, three), 8.0));
3499	}
3500
3501	/**************************************
3502	* To what precision is x equal to y?
3503	*
3504	* Returns: the number of mantissa bits which are equal in x and y.
3505	* eg, 0x1.F8p+60 and 0x1.F1p+60 are equal to 5 bits of precision.
3506	*
3507	* $(TABLE_SV
3508	* $(TR $(TH x) $(TH y) $(TH feqrel(x, y)))
3509	* $(TR $(TD x) $(TD x) $(TD real.mant_dig))
3510	* $(TR $(TD x) $(TD $(GT)= 2*x) $(TD 0))
3511	* $(TR $(TD x) $(TD $(LT)= x/2) $(TD 0))
3512	* $(TR $(TD $(NAN)) $(TD any) $(TD 0))
3513	* $(TR $(TD any) $(TD $(NAN)) $(TD 0))
3514	* )
3515	*/
3516	int feqrel(X)(X x, X y) @trusted pure nothrow
3517	if (isFloatingPoint!(X))
3518	{
3519	/* Public Domain. Author: Don Clugston, 18 Aug 2005.
3520	*/
3521	static if (X.mant_dig == 106) { // doubledouble.
3522	if (cast(double)(&x)[MANTISSA_MSB] == cast(double)(&y)[MANTISSA_MSB]) {
3523	return double.mant_dig
3524	+ feqrel(cast(double*)(&x)[MANTISSA_LSB],
3525	cast(double*)(&y)[MANTISSA_LSB]);
3526	} else {
3527	return feqrel(cast(double*)(&x)[MANTISSA_MSB],
3528	cast(double*)(&y)[MANTISSA_MSB]);
3529	}
3530	} else static if (X.mant_dig==64 \|\| X.mant_dig==113 \|\| X.mant_dig==53) {
3531
3532	if (x == y) return X.mant_dig; // ensure diff!=0, cope with INF.
3533
3534	X diff = fabs(x - y);
3535
3536	ushort pa = cast(ushort )(&x);
3537	ushort pb = cast(ushort )(&y);
3538	ushort pd = cast(ushort )(&diff);
3539
3540	alias floatTraits!(X) F;
3541
3542	// The difference in abs(exponent) between x or y and abs(x-y)
3543	// is equal to the number of significand bits of x which are
3544	// equal to y. If negative, x and y have different exponents.
3545	// If positive, x and y are equal to 'bitsdiff' bits.
3546	// AND with 0x7FFF to form the absolute value.
3547	// To avoid out-by-1 errors, we subtract 1 so it rounds down
3548	// if the exponents were different. This means 'bitsdiff' is
3549	// always 1 lower than we want, except that if bitsdiff==0,
3550	// they could have 0 or 1 bits in common.
3551
3552	static if (X.mant_dig==64 \|\| X.mant_dig==113) { // real80 or quadruple
3553	int bitsdiff = ( ((pa[F.EXPPOS_SHORT] & F.EXPMASK)
3554	+ (pb[F.EXPPOS_SHORT] & F.EXPMASK) - 1) >> 1)
3555	- pd[F.EXPPOS_SHORT];
3556	} else static if (X.mant_dig==53) { // double
3557	int bitsdiff = (( ((pa[F.EXPPOS_SHORT]&0x7FF0)
3558	+ (pb[F.EXPPOS_SHORT]&0x7FF0)-0x10)>>1)
3559	- (pd[F.EXPPOS_SHORT]&0x7FF0))>>4;
3560	}
3561	if (pd[F.EXPPOS_SHORT] == 0)
3562	{ // Difference is denormal
3563	// For denormals, we need to add the number of zeros that
3564	// lie at the start of diff's significand.
3565	// We do this by multiplying by 2^^real.mant_dig
3566	diff *= F.RECIP_EPSILON;
3567	return bitsdiff + X.mant_dig - pd[F.EXPPOS_SHORT];
3568	}
3569
3570	if (bitsdiff > 0)
3571	return bitsdiff + 1; // add the 1 we subtracted before
3572
3573	// Avoid out-by-1 errors when factor is almost 2.
3574	static if (X.mant_dig==64 \|\| X.mant_dig==113) { // real80 or quadruple
3575	return (bitsdiff == 0) ? (pa[F.EXPPOS_SHORT] == pb[F.EXPPOS_SHORT]) : 0;
3576	} else static if (X.mant_dig==53) { // double
3577	if (bitsdiff == 0
3578	&& !((pa[F.EXPPOS_SHORT] ^ pb[F.EXPPOS_SHORT])& F.EXPMASK)) {
3579	return 1;
3580	} else return 0;
3581	}
3582	}
3583	}
3584
3585	unittest
3586	{
3587	// Exact equality
3588	assert(feqrel(real.max,real.max)==real.mant_dig);
3589	assert(feqrel(0.0L,0.0L)==real.mant_dig);
3590	assert(feqrel(7.1824L,7.1824L)==real.mant_dig);
3591	assert(feqrel(real.infinity,real.infinity)==real.mant_dig);
3592
3593	// a few bits away from exact equality
3594	real w=1;
3595	for (int i=1; i<real.mant_dig-1; ++i) {
3596	assert(feqrel(1+w*real.epsilon,1.0L)==real.mant_dig-i);
3597	assert(feqrel(1-w*real.epsilon,1.0L)==real.mant_dig-i);
3598	assert(feqrel(1.0L,1+(w-1)*real.epsilon)==real.mant_dig-i+1);
3599	w*=2;
3600	}
3601	assert(feqrel(1.5+real.epsilon,1.5L)==real.mant_dig-1);
3602	assert(feqrel(1.5-real.epsilon,1.5L)==real.mant_dig-1);
3603	assert(feqrel(1.5-real.epsilon,1.5+real.epsilon)==real.mant_dig-2);
3604
3605	assert(feqrel(real.min_normal/8,real.min_normal/17)==3);;
3606
3607	// Numbers that are close
3608	assert(feqrel(0x1.Bp+84, 0x1.B8p+84)==5);
3609	assert(feqrel(0x1.8p+10, 0x1.Cp+10)==2);
3610	assert(feqrel(1.5*(1-real.epsilon), 1.0L)==2);
3611	assert(feqrel(1.5, 1.0)==1);
3612	assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
3613
3614	// Factors of 2
3615	assert(feqrel(real.max,real.infinity)==0);
3616	assert(feqrel(2*(1-real.epsilon), 1.0L)==1);
3617	assert(feqrel(1.0, 2.0)==0);
3618	assert(feqrel(4.0, 1.0)==0);
3619
3620	// Extreme inequality
3621	assert(feqrel(real.nan,real.nan)==0);
3622	assert(feqrel(0.0L,-real.nan)==0);
3623	assert(feqrel(real.nan,real.infinity)==0);
3624	assert(feqrel(real.infinity,-real.infinity)==0);
3625	assert(feqrel(-real.max,real.infinity)==0);
3626	assert(feqrel(real.max,-real.max)==0);
3627	}
3628
3629	package: // Not public yet
3630	/* Return the value that lies halfway between x and y on the IEEE number line.
3631	*
3632	* Formally, the result is the arithmetic mean of the binary significands of x
3633	* and y, multiplied by the geometric mean of the binary exponents of x and y.
3634	* x and y must have the same sign, and must not be NaN.
3635	* Note: this function is useful for ensuring O(log n) behaviour in algorithms
3636	* involving a 'binary chop'.
3637	*
3638	* Special cases:
3639	* If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value
3640	* is the arithmetic mean (x + y) / 2.
3641	* If x and y are even powers of 2, the return value is the geometric mean,
3642	* ieeeMean(x, y) = sqrt(x * y).
3643	*
3644	*/
3645	T ieeeMean(T)(T x, T y) @trusted pure nothrow
3646	in {
3647	// both x and y must have the same sign, and must not be NaN.
3648	assert(signbit(x) == signbit(y));
3649	assert(x<>=0 && y<>=0);
3650	}
3651	body {
3652	// Runtime behaviour for contract violation:
3653	// If signs are opposite, or one is a NaN, return 0.
3654	if (!((x>=0 && y>=0) \|\| (x<=0 && y<=0))) return 0.0;
3655
3656	// The implementation is simple: cast x and y to integers,
3657	// average them (avoiding overflow), and cast the result back to a floating-point number.
3658
3659	alias floatTraits!(real) F;
3660	T u;
3661	static if (T.mant_dig==64) { // real80
3662	// There's slight additional complexity because they are actually
3663	// 79-bit reals...
3664	ushort ue = cast(ushort )&u;
3665	ulong ul = cast(ulong )&u;
3666	ushort xe = cast(ushort )&x;
3667	ulong xl = cast(ulong )&x;
3668	ushort ye = cast(ushort )&y;
3669	ulong yl = cast(ulong )&y;
3670	// Ignore the useless implicit bit. (Bonus: this prevents overflows)
3671	ulong m = ((xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((yl) & 0x7FFF_FFFF_FFFF_FFFFL);
3672
3673	// @@@ BUG? @@@
3674	// Cast shouldn't be here
3675	ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK)
3676	+ (ye[F.EXPPOS_SHORT] & F.EXPMASK));
3677	if (m & 0x8000_0000_0000_0000L) {
3678	++e;
3679	m &= 0x7FFF_FFFF_FFFF_FFFFL;
3680	}
3681	// Now do a multi-byte right shift
3682	uint c = e & 1; // carry
3683	e >>= 1;
3684	m >>>= 1;
3685	if (c) m \|= 0x4000_0000_0000_0000L; // shift carry into significand
3686	if (e) *ul = m \| 0x8000_0000_0000_0000L; // set implicit bit...
3687	else *ul = m; // ... unless exponent is 0 (denormal or zero).
3688	ue[4]= e \| (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit
3689	} else static if(T.mant_dig == 113) { //quadruple
3690	// This would be trivial if 'ucent' were implemented...
3691	ulong ul = cast(ulong )&u;
3692	ulong xl = cast(ulong )&x;
3693	ulong yl = cast(ulong )&y;
3694	// Multi-byte add, then multi-byte right shift.
3695	ulong mh = ((xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL)
3696	+ (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL));
3697	// Discard the lowest bit (to avoid overflow)
3698	ulong ml = (xl[MANTISSA_LSB]>>>1) + (yl[MANTISSA_LSB]>>>1);
3699	// add the lowest bit back in, if necessary.
3700	if (xl[MANTISSA_LSB] & yl[MANTISSA_LSB] & 1) {
3701	++ml;
3702	if (ml==0) ++mh;
3703	}
3704	mh >>>=1;
3705	ul[MANTISSA_MSB] = mh \| (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000);
3706	ul[MANTISSA_LSB] = ml;
3707	} else static if (T.mant_dig == double.mant_dig) {
3708	ulong ul = cast(ulong )&u;
3709	ulong xl = cast(ulong )&x;
3710	ulong yl = cast(ulong )&y;
3711	ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL)
3712	+ ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1;
3713	m \|= ((*xl) & 0x8000_0000_0000_0000L);
3714	*ul = m;
3715	} else static if (T.mant_dig == float.mant_dig) {
3716	uint ul = cast(uint )&u;
3717	uint xl = cast(uint )&x;
3718	uint yl = cast(uint )&y;
3719	uint m = (((xl) & 0x7FFF_FFFF) + ((yl) & 0x7FFF_FFFF)) >>> 1;
3720	m \|= ((*xl) & 0x8000_0000);
3721	*ul = m;
3722	} else {
3723	assert(0, "Not implemented");
3724	}
3725	return u;
3726	}
3727
3728	unittest {
3729	assert(ieeeMean(-0.0,-1e-20)<0);
3730	assert(ieeeMean(0.0,1e-20)>0);
3731
3732	assert(ieeeMean(1.0L,4.0L)==2L);
3733	assert(ieeeMean(2.01.013,8.01.013)==4*1.013);
3734	assert(ieeeMean(-1.0L,-4.0L)==-2L);
3735	assert(ieeeMean(-1.0,-4.0)==-2);
3736	assert(ieeeMean(-1.0f,-4.0f)==-2f);
3737	assert(ieeeMean(-1.0,-2.0)==-1.5);
3738	assert(ieeeMean(-1(1+8real.epsilon),-2(1+8real.epsilon))
3739	==-1.5(1+5real.epsilon));
3740	assert(ieeeMean(0x1p60,0x1p-10)==0x1p25);
3741	static if (real.mant_dig==64) { // x87, 80-bit reals
3742	assert(ieeeMean(1.0L,real.infinity)==0x1p8192L);
3743	assert(ieeeMean(0.0L,real.infinity)==1.5);
3744	}
3745	assert(ieeeMean(0.5real.min_normal(1-4real.epsilon),0.5real.min_normal)
3746	== 0.5real.min_normal(1-2*real.epsilon));
3747	}
3748
3749	public:
3750
3751
3752	/***********************************
3753	* Evaluate polynomial A(x) = $(SUB a, 0) + $(SUB a, 1)x + $(SUB a, 2)$(POWER x,2)
3754	* + $(SUB a,3)$(POWER x,3); ...
3755	*
3756	* Uses Horner's rule A(x) = $(SUB a, 0) + x($(SUB a, 1) + x($(SUB a, 2)
3757	* + x($(SUB a, 3) + ...)))
3758	* Params:
3759	* A = array of coefficients $(SUB a, 0), $(SUB a, 1), etc.
3760	*/
3761	real poly(real x, const real[] A) @trusted pure nothrow
3762	in
3763	{
3764	assert(A.length > 0);
3765	}
3766	body
3767	{
3768	version (D_InlineAsm_X86)
3769	{
3770	version (Windows)
3771	{
3772	// BUG: This code assumes a frame pointer in EBP.
3773	asm // assembler by W. Bright
3774	{
3775	// EDX = (A.length - 1) * real.sizeof
3776	mov ECX,A[EBP] ; // ECX = A.length
3777	dec ECX ;
3778	lea EDX,[ECX][ECX*8] ;
3779	add EDX,ECX ;
3780	add EDX,A+4[EBP] ;
3781	fld real ptr [EDX] ; // ST0 = coeff[ECX]
3782	jecxz return_ST ;
3783	fld x[EBP] ; // ST0 = x
3784	fxch ST(1) ; // ST1 = x, ST0 = r
3785	align 4 ;
3786	L2: fmul ST,ST(1) ; // r *= x
3787	fld real ptr -10[EDX] ;
3788	sub EDX,10 ; // deg--
3789	faddp ST(1),ST ;
3790	dec ECX ;
3791	jne L2 ;
3792	fxch ST(1) ; // ST1 = r, ST0 = x
3793	fstp ST(0) ; // dump x
3794	align 4 ;
3795	return_ST: ;
3796	;
3797	}
3798	}
3799	else version (linux)
3800	{
3801	asm // assembler by W. Bright
3802	{
3803	// EDX = (A.length - 1) * real.sizeof
3804	mov ECX,A[EBP] ; // ECX = A.length
3805	dec ECX ;
3806	lea EDX,[ECX*8] ;
3807	lea EDX,[EDX][ECX*4] ;
3808	add EDX,A+4[EBP] ;
3809	fld real ptr [EDX] ; // ST0 = coeff[ECX]
3810	jecxz return_ST ;
3811	fld x[EBP] ; // ST0 = x
3812	fxch ST(1) ; // ST1 = x, ST0 = r
3813	align 4 ;
3814	L2: fmul ST,ST(1) ; // r *= x
3815	fld real ptr -12[EDX] ;
3816	sub EDX,12 ; // deg--
3817	faddp ST(1),ST ;
3818	dec ECX ;
3819	jne L2 ;
3820	fxch ST(1) ; // ST1 = r, ST0 = x
3821	fstp ST(0) ; // dump x
3822	align 4 ;
3823	return_ST: ;
3824	;
3825	}
3826	}
3827	else version (OSX)
3828	{
3829	asm // assembler by W. Bright
3830	{
3831	// EDX = (A.length - 1) * real.sizeof
3832	mov ECX,A[EBP] ; // ECX = A.length
3833	dec ECX ;
3834	lea EDX,[ECX*8] ;
3835	add EDX,EDX ;
3836	add EDX,A+4[EBP] ;
3837	fld real ptr [EDX] ; // ST0 = coeff[ECX]
3838	jecxz return_ST ;
3839	fld x[EBP] ; // ST0 = x
3840	fxch ST(1) ; // ST1 = x, ST0 = r
3841	align 4 ;
3842	L2: fmul ST,ST(1) ; // r *= x
3843	fld real ptr -16[EDX] ;
3844	sub EDX,16 ; // deg--
3845	faddp ST(1),ST ;
3846	dec ECX ;
3847	jne L2 ;
3848	fxch ST(1) ; // ST1 = r, ST0 = x
3849	fstp ST(0) ; // dump x
3850	align 4 ;
3851	return_ST: ;
3852	;
3853	}
3854	}
3855	else version (FreeBSD)
3856	{
3857	asm // assembler by W. Bright
3858	{
3859	// EDX = (A.length - 1) * real.sizeof
3860	mov ECX,A[EBP] ; // ECX = A.length
3861	dec ECX ;
3862	lea EDX,[ECX*8] ;
3863	lea EDX,[EDX][ECX*4] ;
3864	add EDX,A+4[EBP] ;
3865	fld real ptr [EDX] ; // ST0 = coeff[ECX]
3866	jecxz return_ST ;
3867	fld x[EBP] ; // ST0 = x
3868	fxch ST(1) ; // ST1 = x, ST0 = r
3869	align 4 ;
3870	L2: fmul ST,ST(1) ; // r *= x
3871	fld real ptr -12[EDX] ;
3872	sub EDX,12 ; // deg--
3873	faddp ST(1),ST ;
3874	dec ECX ;
3875	jne L2 ;
3876	fxch ST(1) ; // ST1 = r, ST0 = x
3877	fstp ST(0) ; // dump x
3878	align 4 ;
3879	return_ST: ;
3880	;
3881	}
3882	}
3883	else
3884	{
3885	static assert(0);
3886	}
3887	}
3888	else
3889	{
3890	sizediff_t i = A.length - 1;
3891	real r = A[i];
3892	while (--i >= 0)
3893	{
3894	r *= x;
3895	r += A[i];
3896	}
3897	return r;
3898	}
3899	}
3900
3901	unittest
3902	{
3903	debug (math) printf("math.poly.unittest\n");
3904	real x = 3.1;
3905	static real pp[] = [56.1, 32.7, 6];
3906
3907	assert( poly(x, pp) == (56.1L + (32.7L + 6L * x) * x) );
3908	}
3909
3910	/**
3911	Computes whether $(D lhs) is approximately equal to $(D rhs)
3912	admitting a maximum relative difference $(D maxRelDiff) and a
3913	maximum absolute difference $(D maxAbsDiff).
3914
3915	If the two inputs are ranges, $(D approxEqual) returns true if and
3916	only if the ranges have the same number of elements and if $(D
3917	approxEqual) evaluates to $(D true) for each pair of elements.
3918	*/
3919	bool approxEqual(T, U, V)(T lhs, U rhs, V maxRelDiff, V maxAbsDiff = 1e-5)
3920	{
3921	static if (isInputRange!T)
3922	{
3923	static if (isInputRange!U)
3924	{
3925	// Two ranges
3926	for (;; lhs.popFront, rhs.popFront)
3927	{
3928	if (lhs.empty) return rhs.empty;
3929	if (rhs.empty) return lhs.empty;
3930	if (!approxEqual(lhs.front, rhs.front, maxRelDiff, maxAbsDiff))
3931	return false;
3932	}
3933	}
3934	else
3935	{
3936	// lhs is range, rhs is number
3937	for (; !lhs.empty; lhs.popFront)
3938	{
3939	if (!approxEqual(lhs.front, rhs, maxRelDiff, maxAbsDiff))
3940	return false;
3941	}
3942	return true;
3943	}
3944	}
3945	else
3946	{
3947	static if (isInputRange!U)
3948	{
3949	// lhs is number, rhs is array
3950	return approxEqual(rhs, lhs, maxRelDiff, maxAbsDiff);
3951	}
3952	else
3953	{
3954	// two numbers
3955	//static assert(is(T : real) && is(U : real));
3956	if (rhs == 0)
3957	{
3958	return fabs(lhs) <= maxAbsDiff;
3959	}
3960	static if (is(typeof(lhs.infinity)) && is(typeof(rhs.infinity)))
3961	{
3962	if (lhs == lhs.infinity && rhs == rhs.infinity \|\|
3963	lhs == -lhs.infinity && rhs == -rhs.infinity) return true;
3964	}
3965	return fabs((lhs - rhs) / rhs) <= maxRelDiff
3966	\|\| maxAbsDiff != 0 && fabs(lhs - rhs) <= maxAbsDiff;
3967	}
3968	}
3969	}
3970
3971	/**
3972	Returns $(D approxEqual(lhs, rhs, 1e-2, 1e-5)).
3973	*/
3974	bool approxEqual(T, U)(T lhs, U rhs)
3975	{
3976	return approxEqual(lhs, rhs, 1e-2, 1e-5);
3977	}
3978
3979	unittest
3980	{
3981	assert(approxEqual(1.0, 1.0099));
3982	assert(!approxEqual(1.0, 1.011));
3983	float[] arr1 = [ 1.0, 2.0, 3.0 ];
3984	double[] arr2 = [ 1.001, 1.999, 3 ];
3985	assert(approxEqual(arr1, arr2));
3986
3987	real num = real.infinity;
3988	assert(num == real.infinity); // Passes.
3989	assert(approxEqual(num, real.infinity)); // Fails.
3990	num = -real.infinity;
3991	assert(num == -real.infinity); // Passes.
3992	assert(approxEqual(num, -real.infinity)); // Fails.
3993	}
3994
3995	// Included for backwards compatibility with Phobos1
3996	alias isNaN isnan;
3997	alias isFinite isfinite;
3998	alias isNormal isnormal;
3999	alias isSubnormal issubnormal;
4000	alias isInfinity isinf;
4001
4002	/* **********************************
4003	* Building block functions, they
4004	* translate to a single x87 instruction.
4005	*/
4006
4007	real yl2x(real x, real y) @safe pure nothrow; // y * log2(x)
4008	real yl2xp1(real x, real y) @safe pure nothrow; // y * log2(x + 1)
4009
4010	unittest
4011	{
4012	version (INLINE_YL2X)
4013	{
4014	assert(yl2x(1024, 1) == 10);
4015	assert(yl2xp1(1023, 1) == 10);
4016	}
4017	}
4018
4019	unittest
4020	{
4021	real num = real.infinity;
4022	assert(num == real.infinity); // Passes.
4023	assert(approxEqual(num, real.infinity)); // Fails.
4024	}

Note: See TracBrowser for help on using the browser.

phobos

root/trunk/phobos/std/math.d

Download in other formats: