Class DD

1	/*
2	* Copyright (c) 2016 Martin Davis.
3	*
4	* All rights reserved. This program and the accompanying materials
5	* are made available under the terms of the Eclipse Public License 2.0
6	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
7	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
8	* and the Eclipse Distribution License is available at
9	*
10	* http://www.eclipse.org/org/documents/edl-v10.php.
11	*/
12	package org.locationtech.jts.math;
13
14	import java.io.Serializable;
15
16	/**
17	* Implements extended-precision floating-point numbers
18	* which maintain 106 bits (approximately 30 decimal digits) of precision.
19	* <p>
20	* A DoubleDouble uses a representation containing two double-precision values.
21	* A number x is represented as a pair of doubles, x.hi and x.lo,
22	* such that the number represented by x is x.hi + x.lo, where
23	* <pre>
24	* \|x.lo\| <= 0.5*ulp(x.hi)
25	* </pre>
26	* and ulp(y) means "unit in the last place of y".
27	* The basic arithmetic operations are implemented using
28	* convenient properties of IEEE-754 floating-point arithmetic.
29	* <p>
30	* The range of values which can be represented is the same as in IEEE-754.
31	* The precision of the representable numbers
32	* is twice as great as IEEE-754 double precision.
33	* <p>
34	* The correctness of the arithmetic algorithms relies on operations
35	* being performed with standard IEEE-754 double precision and rounding.
36	* This is the Java standard arithmetic model, but for performance reasons
37	* Java implementations are not
38	* constrained to using this standard by default.
39	* Some processors (notably the Intel Pentium architecture) perform
40	* floating point operations in (non-IEEE-754-standard) extended-precision.
41	* A JVM implementation may choose to use the non-standard extended-precision
42	* as its default arithmetic mode.
43	* To prevent this from happening, this code uses the
44	* Java <tt>strictfp</tt> modifier,
45	* which forces all operations to take place in the standard IEEE-754 rounding model.
46	* <p>
47	* The API provides both a set of value-oriented operations
48	* and a set of mutating operations.
49	* Value-oriented operations treat DoubleDouble values as
50	* immutable; operations on them return new objects carrying the result
51	* of the operation. This provides a simple and safe semantics for
52	* writing DoubleDouble expressions. However, there is a performance
53	* penalty for the object allocations required.
54	* The mutable interface updates object values in-place.
55	* It provides optimum memory performance, but requires
56	* care to ensure that aliasing errors are not created
57	* and constant values are not changed.
58	* <p>
59	* For example, the following code example constructs three DD instances:
60	* two to hold the input values and one to hold the result of the addition.
61	* <pre>
62	* DD a = new DD(2.0);
63	* DD b = new DD(3.0);
64	* DD c = a.add(b);
65	* </pre>
66	* In contrast, the following approach uses only one object:
67	* <pre>
68	* DD a = new DD(2.0);
69	* a.selfAdd(3.0);
70	* </pre>
71	* <p>
72	* This implementation uses algorithms originally designed variously by
73	* Knuth, Kahan, Dekker, and Linnainmaa.
74	* Douglas Priest developed the first C implementation of these techniques.
75	* Other more recent C++ implementation are due to Keith M. Briggs and David Bailey et al.
76	*
77	* <h3>References</h3>
78	* <ul>
79	* <li>Priest, D., <i>Algorithms for Arbitrary Precision Floating Point Arithmetic</i>,
80	* in P. Kornerup and D. Matula, Eds., Proc. 10th Symposium on Computer Arithmetic,
81	* IEEE Computer Society Press, Los Alamitos, Calif., 1991.
82	* <li>Yozo Hida, Xiaoye S. Li and David H. Bailey,
83	* <i>Quad-Double Arithmetic: Algorithms, Implementation, and Application</i>,
84	* manuscript, Oct 2000; Lawrence Berkeley National Laboratory Report BNL-46996.
85	* <li>David Bailey, <i>High Precision Software Directory</i>;
86	* <tt>http://crd.lbl.gov/~dhbailey/mpdist/index.html</tt>
87	* </ul>
88	*
89	*
90	* @author Martin Davis
91	*
92	*/
93	public strictfp final class DD
94	implements Serializable, Comparable, Cloneable
95	{
96	/**
97	* The value nearest to the constant Pi.
98	*/
99	public static final DD PI = new DD(
100	3.141592653589793116e+00,
101	1.224646799147353207e-16);
102
103	/**
104	* The value nearest to the constant 2 * Pi.
105	*/
106	public static final DD TWO_PI = new DD(
107	6.283185307179586232e+00,
108	2.449293598294706414e-16);
109
110	/**
111	* The value nearest to the constant Pi / 2.
112	*/
113	public static final DD PI_2 = new DD(
114	1.570796326794896558e+00,
115	6.123233995736766036e-17);
116
117	/**
118	* The value nearest to the constant e (the natural logarithm base).
119	*/
120	public static final DD E = new DD(
121	2.718281828459045091e+00,
122	1.445646891729250158e-16);
123
124	/**
125	* A value representing the result of an operation which does not return a valid number.
126	*/
127	public static final DD NaN = new DD(Double.NaN, Double.NaN);
128
129	/**
130	* The smallest representable relative difference between two {link @ DoubleDouble} values
131	*/
132	public static final double EPS = 1.23259516440783e-32; /* = 2^-106 */
133
134	private static DD createNaN()
135	{
136	return new DD(Double.NaN, Double.NaN);
137	}
138
139	/**
140	* Converts the string argument to a DoubleDouble number.
141	*
142	* @param str a string containing a representation of a numeric value
143	* @return the extended precision version of the value
144	* @throws NumberFormatException if <tt>s</tt> is not a valid representation of a number
145	*/
146	public static DD valueOf(String str)
147	throws NumberFormatException
148	{
149	return parse(str);
150	}
151
152	/**
153	* Converts the <tt>double</tt> argument to a DoubleDouble number.
154	*
155	* @param x a numeric value
156	* @return the extended precision version of the value
157	*/
158	public static DD valueOf(double x) { return new DD(x); }
159
160	/**
161	* The value to split a double-precision value on during multiplication
162	*/
163	private static final double SPLIT = 134217729.0D; // 2^27+1, for IEEE double
164
165	/**
166	* The high-order component of the double-double precision value.
167	*/
168	private double hi = 0.0;
169
170	/**
171	* The low-order component of the double-double precision value.
172	*/
173	private double lo = 0.0;
174
175	/**
176	* Creates a new DoubleDouble with value 0.0.
177	*/
178	public DD()
179	{
180	init(0.0);
181	}
182
183	/**
184	* Creates a new DoubleDouble with value x.
185	*
186	* @param x the value to initialize
187	*/
188	public DD(double x)
189	{
190	init(x);
191	}
192
193	/**
194	* Creates a new DoubleDouble with value (hi, lo).
195	*
196	* @param hi the high-order component
197	* @param lo the high-order component
198	*/
199	public DD(double hi, double lo)
200	{
201	init(hi, lo);
202	}
203
204	/**
205	* Creates a new DoubleDouble with value equal to the argument.
206	*
207	* @param dd the value to initialize
208	*/
209	public DD(DD dd)
210	{
211	init(dd);
212	}
213
214	/**
215	* Creates a new DoubleDouble with value equal to the argument.
216	*
217	* @param str the value to initialize by
218	* @throws NumberFormatException if <tt>str</tt> is not a valid representation of a number
219	*/
220	public DD(String str)
221	throws NumberFormatException
222	{
223	this(parse(str));
224	}
225
226	/**
227	* Creates a new DoubleDouble with the value of the argument.
228	*
229	* @param dd the DoubleDouble value to copy
230	* @return a copy of the input value
231	*/
232	public static DD copy(DD dd)
233	{
234	return new DD(dd);
235	}
236
237	/**
238	* Creates and returns a copy of this value.
239	*
240	* @return a copy of this value
241	*/
242	public Object clone()
243	{
244	try {
245	return super.clone();
246	}
247	catch (CloneNotSupportedException ex) {
248	// should never reach here
249	return null;
250	}
251	}
252
253	private final void init(double x)
254	{
255	this.hi = x;
256	this.lo = 0.0;
257	}
258
259	private final void init(double hi, double lo)
260	{
261	this.hi = hi;
262	this.lo = lo;
263	}
264
265	private final void init(DD dd)
266	{
267	hi = dd.hi;
268	lo = dd.lo;
269	}
270
271	/*
272	double getHighComponent() { return hi; }
273
274	double getLowComponent() { return lo; }
275	*/
276
277	// Testing only - should not be public
278	/*
279	public void RENORM()
280	{
281	double s = hi + lo;
282	double err = lo - (s - hi);
283	hi = s;
284	lo = err;
285	}
286	*/
287
288	/**
289	* Set the value for the DD object. This method supports the mutating
290	* operations concept described in the class documentation (see above).
291	* @param value a DD instance supplying an extended-precision value.
292	* @return a self-reference to the DD instance.
293	*/
294	public DD setValue(DD value) {
295	init(value);
296	return this;
297	}
298
299	/**
300	* Set the value for the DD object. This method supports the mutating
301	* operations concept described in the class documentation (see above).
302	* @param value a floating point value to be stored in the instance.
303	* @return a self-reference to the DD instance.
304	*/
305	public DD setValue(double value) {
306	init(value);
307	return this;
308	}
309
310
311	/**
312	* Returns a new DoubleDouble whose value is <tt>(this + y)</tt>.
313	*
314	* @param y the addend
315	* @return <tt>(this + y)</tt>
316	*/
317	public final DD add(DD y)
318	{
319	return copy(this).selfAdd(y);
320	}
321
322	/**
323	* Returns a new DoubleDouble whose value is <tt>(this + y)</tt>.
324	*
325	* @param y the addend
326	* @return <tt>(this + y)</tt>
327	*/
328	public final DD add(double y)
329	{
330	return copy(this).selfAdd(y);
331	}
332
333	/**
334	* Adds the argument to the value of <tt>this</tt>.
335	* To prevent altering constants,
336	* this method <b>must only</b> be used on values known to
337	* be newly created.
338	*
339	* @param y the addend
340	* @return this object, increased by y
341	*/
342	public final DD selfAdd(DD y)
343	{
344	return selfAdd(y.hi, y.lo);
345	}
346
347	/**
348	* Adds the argument to the value of <tt>this</tt>.
349	* To prevent altering constants,
350	* this method <b>must only</b> be used on values known to
351	* be newly created.
352	*
353	* @param y the addend
354	* @return this object, increased by y
355	*/
356	public final DD selfAdd(double y)
357	{
358	double H, h, S, s, e, f;
359	S = hi + y;
360	e = S - hi;
361	s = S - e;
362	s = (y - e) + (hi - s);
363	f = s + lo;
364	H = S + f;
365	h = f + (S - H);
366	hi = H + h;
367	lo = h + (H - hi);
368	return this;
369	// return selfAdd(y, 0.0);
370	}
371
372	private final DD selfAdd(double yhi, double ylo)
373	{
374	double H, h, T, t, S, s, e, f;
375	S = hi + yhi;
376	T = lo + ylo;
377	e = S - hi;
378	f = T - lo;
379	s = S-e;
380	t = T-f;
381	s = (yhi-e)+(hi-s);
382	t = (ylo-f)+(lo-t);
383	e = s+T; H = S+e; h = e+(S-H); e = t+h;
384
385	double zhi = H + e;
386	double zlo = e + (H - zhi);
387	hi = zhi;
388	lo = zlo;
389	return this;
390	}
391
392	/**
393	* Computes a new DoubleDouble object whose value is <tt>(this - y)</tt>.
394	*
395	* @param y the subtrahend
396	* @return <tt>(this - y)</tt>
397	*/
398	public final DD subtract(DD y)
399	{
400	return add(y.negate());
401	}
402
403	/**
404	* Computes a new DoubleDouble object whose value is <tt>(this - y)</tt>.
405	*
406	* @param y the subtrahend
407	* @return <tt>(this - y)</tt>
408	*/
409	public final DD subtract(double y)
410	{
411	return add(-y);
412	}
413
414
415	/**
416	* Subtracts the argument from the value of <tt>this</tt>.
417	* To prevent altering constants,
418	* this method <b>must only</b> be used on values known to
419	* be newly created.
420	*
421	* @param y the addend
422	* @return this object, decreased by y
423	*/
424	public final DD selfSubtract(DD y)
425	{
426	if (isNaN()) return this;
427	return selfAdd(-y.hi, -y.lo);
428	}
429
430	/**
431	* Subtracts the argument from the value of <tt>this</tt>.
432	* To prevent altering constants,
433	* this method <b>must only</b> be used on values known to
434	* be newly created.
435	*
436	* @param y the addend
437	* @return this object, decreased by y
438	*/
439	public final DD selfSubtract(double y)
440	{
441	if (isNaN()) return this;
442	return selfAdd(-y, 0.0);
443	}
444
445	/**
446	* Returns a new DoubleDouble whose value is <tt>-this</tt>.
447	*
448	* @return <tt>-this</tt>
449	*/
450	public final DD negate()
451	{
452	if (isNaN()) return this;
453	return new DD(-hi, -lo);
454	}
455
456	/**
457	* Returns a new DoubleDouble whose value is <tt>(this * y)</tt>.
458	*
459	* @param y the multiplicand
460	* @return <tt>(this * y)</tt>
461	*/
462	public final DD multiply(DD y)
463	{
464	if (y.isNaN()) return createNaN();
465	return copy(this).selfMultiply(y);
466	}
467
468	/**
469	* Returns a new DoubleDouble whose value is <tt>(this * y)</tt>.
470	*
471	* @param y the multiplicand
472	* @return <tt>(this * y)</tt>
473	*/
474	public final DD multiply(double y)
475	{
476	if (Double.isNaN(y)) return createNaN();
477	return copy(this).selfMultiply(y, 0.0);
478	}
479
480	/**
481	* Multiplies this object by the argument, returning <tt>this</tt>.
482	* To prevent altering constants,
483	* this method <b>must only</b> be used on values known to
484	* be newly created.
485	*
486	* @param y the value to multiply by
487	* @return this object, multiplied by y
488	*/
489	public final DD selfMultiply(DD y)
490	{
491	return selfMultiply(y.hi, y.lo);
492	}
493
494	/**
495	* Multiplies this object by the argument, returning <tt>this</tt>.
496	* To prevent altering constants,
497	* this method <b>must only</b> be used on values known to
498	* be newly created.
499	*
500	* @param y the value to multiply by
501	* @return this object, multiplied by y
502	*/
503	public final DD selfMultiply(double y)
504	{
505	return selfMultiply(y, 0.0);
506	}
507
508	private final DD selfMultiply(double yhi, double ylo)
509	{
510	double hx, tx, hy, ty, C, c;
511	C = SPLIT * hi; hx = C-hi; c = SPLIT * yhi;
512	hx = C-hx; tx = hi-hx; hy = c-yhi;
513	C = hi*yhi; hy = c-hy; ty = yhi-hy;
514	c = ((((hxhy-C)+hxty)+txhy)+txty)+(hiylo+loyhi);
515	double zhi = C+c; hx = C-zhi;
516	double zlo = c+hx;
517	hi = zhi;
518	lo = zlo;
519	return this;
520	}
521
522	/**
523	* Computes a new DoubleDouble whose value is <tt>(this / y)</tt>.
524	*
525	* @param y the divisor
526	* @return a new object with the value <tt>(this / y)</tt>
527	*/
528	public final DD divide(DD y)
529	{
530	double hc, tc, hy, ty, C, c, U, u;
531	C = hi/y.hi; c = SPLITC; hc =c-C; u = SPLITy.hi; hc = c-hc;
532	tc = C-hc; hy = u-y.hi; U = C * y.hi; hy = u-hy; ty = y.hi-hy;
533	u = (((hchy-U)+hcty)+tchy)+tcty;
534	c = ((((hi-U)-u)+lo)-C*y.lo)/y.hi;
535	u = C+c;
536
537	double zhi = u;
538	double zlo = (C-u)+c;
539	return new DD(zhi, zlo);
540	}
541
542	/**
543	* Computes a new DoubleDouble whose value is <tt>(this / y)</tt>.
544	*
545	* @param y the divisor
546	* @return a new object with the value <tt>(this / y)</tt>
547	*/
548	public final DD divide(double y)
549	{
550	if (Double.isNaN(y)) return createNaN();
551	return copy(this).selfDivide(y, 0.0);
552	}
553
554	/**
555	* Divides this object by the argument, returning <tt>this</tt>.
556	* To prevent altering constants,
557	* this method <b>must only</b> be used on values known to
558	* be newly created.
559	*
560	* @param y the value to divide by
561	* @return this object, divided by y
562	*/
563	public final DD selfDivide(DD y)
564	{
565	return selfDivide(y.hi, y.lo);
566	}
567
568	/**
569	* Divides this object by the argument, returning <tt>this</tt>.
570	* To prevent altering constants,
571	* this method <b>must only</b> be used on values known to
572	* be newly created.
573	*
574	* @param y the value to divide by
575	* @return this object, divided by y
576	*/
577	public final DD selfDivide(double y)
578	{
579	return selfDivide(y, 0.0);
580	}
581
582	private final DD selfDivide(double yhi, double ylo)
583	{
584	double hc, tc, hy, ty, C, c, U, u;
585	C = hi/yhi; c = SPLITC; hc =c-C; u = SPLITyhi; hc = c-hc;
586	tc = C-hc; hy = u-yhi; U = C * yhi; hy = u-hy; ty = yhi-hy;
587	u = (((hchy-U)+hcty)+tchy)+tcty;
588	c = ((((hi-U)-u)+lo)-C*ylo)/yhi;
589	u = C+c;
590
591	hi = u;
592	lo = (C-u)+c;
593	return this;
594	}
595
596	/**
597	* Returns a DoubleDouble whose value is <tt>1 / this</tt>.
598	*
599	* @return the reciprocal of this value
600	*/
601	public final DD reciprocal()
602	{
603	double hc, tc, hy, ty, C, c, U, u;
604	C = 1.0/hi;
605	c = SPLIT*C;
606	hc =c-C;
607	u = SPLIT*hi;
608	hc = c-hc; tc = C-hc; hy = u-hi; U = C*hi; hy = u-hy; ty = hi-hy;
609	u = (((hchy-U)+hcty)+tchy)+tcty;
610	c = ((((1.0-U)-u))-C*lo)/hi;
611
612	double zhi = C+c;
613	double zlo = (C-zhi)+c;
614	return new DD(zhi, zlo);
615	}
616
617	/**
618	* Returns the largest (closest to positive infinity)
619	* value that is not greater than the argument
620	* and is equal to a mathematical integer.
621	* Special cases:
622	* <ul>
623	* <li>If this value is NaN, returns NaN.
624	* </ul>
625	*
626	* @return the largest (closest to positive infinity)
627	* value that is not greater than the argument
628	* and is equal to a mathematical integer.
629	*/
630	public DD floor()
631	{
632	if (isNaN()) return NaN;
633	double fhi=Math.floor(hi);
634	double flo = 0.0;
635	// Hi is already integral. Floor the low word
636	if (fhi == hi) {
637	flo = Math.floor(lo);
638	}
639	// do we need to renormalize here?
640	return new DD(fhi, flo);
641	}
642
643	/**
644	* Returns the smallest (closest to negative infinity) value
645	* that is not less than the argument and is equal to a mathematical integer.
646	* Special cases:
647	* <ul>
648	* <li>If this value is NaN, returns NaN.
649	* </ul>
650	*
651	* @return the smallest (closest to negative infinity) value
652	* that is not less than the argument and is equal to a mathematical integer.
653	*/
654	public DD ceil()
655	{
656	if (isNaN()) return NaN;
657	double fhi=Math.ceil(hi);
658	double flo = 0.0;
659	// Hi is already integral. Ceil the low word
660	if (fhi == hi) {
661	flo = Math.ceil(lo);
662	// do we need to renormalize here?
663	}
664	return new DD(fhi, flo);
665	}
666
667	/**
668	* Returns an integer indicating the sign of this value.
669	* <ul>
670	* <li>if this value is > 0, returns 1
671	* <li>if this value is < 0, returns -1
672	* <li>if this value is = 0, returns 0
673	* <li>if this value is NaN, returns 0
674	* </ul>
675	*
676	* @return an integer indicating the sign of this value
677	*/
678	public int signum()
679	{
680	if (hi > 0) return 1;
681	if (hi < 0) return -1;
682	if (lo > 0) return 1;
683	if (lo < 0) return -1;
684	return 0;
685	}
686
687	/**
688	* Rounds this value to the nearest integer.
689	* The value is rounded to an integer by adding 1/2 and taking the floor of the result.
690	* Special cases:
691	* <ul>
692	* <li>If this value is NaN, returns NaN.
693	* </ul>
694	*
695	* @return this value rounded to the nearest integer
696	*/
697	public DD rint()
698	{
699	if (isNaN()) return this;
700	// may not be 100% correct
701	DD plus5 = this.add(0.5);
702	return plus5.floor();
703	}
704
705	/**
706	* Returns the integer which is largest in absolute value and not further
707	* from zero than this value.
708	* Special cases:
709	* <ul>
710	* <li>If this value is NaN, returns NaN.
711	* </ul>
712	*
713	* @return the integer which is largest in absolute value and not further from zero than this value
714	*/
715	public DD trunc()
716	{
717	if (isNaN()) return NaN;
718	if (isPositive())
719	return floor();
720	else
721	return ceil();
722	}
723
724	/**
725	* Returns the absolute value of this value.
726	* Special cases:
727	* <ul>
728	* <li>If this value is NaN, it is returned.
729	* </ul>
730	*
731	* @return the absolute value of this value
732	*/
733	public DD abs()
734	{
735	if (isNaN()) return NaN;
736	if (isNegative())
737	return negate();
738	return new DD(this);
739	}
740
741	/**
742	* Computes the square of this value.
743	*
744	* @return the square of this value.
745	*/
746	public DD sqr()
747	{
748	return this.multiply(this);
749	}
750
751	/**
752	* Squares this object.
753	* To prevent altering constants,
754	* this method <b>must only</b> be used on values known to
755	* be newly created.
756	*
757	* @return the square of this value.
758	*/
759	public DD selfSqr()
760	{
761	return this.selfMultiply(this);
762	}
763
764	/**
765	* Computes the square of this value.
766	*
767	* @return the square of this value.
768	*/
769	public static DD sqr(double x)
770	{
771	return valueOf(x).selfMultiply(x);
772	}
773
774	/**
775	* Computes the positive square root of this value.
776	* If the number is NaN or negative, NaN is returned.
777	*
778	* @return the positive square root of this number.
779	* If the argument is NaN or less than zero, the result is NaN.
780	*/
781	public DD sqrt()
782	{
783	/* Strategy: Use Karp's trick: if x is an approximation
784	to sqrt(a), then
785
786	sqrt(a) = ax + [a - (ax)^2] * x / 2 (approx)
787
788	The approximation is accurate to twice the accuracy of x.
789	Also, the multiplication (ax) and [-]x can be done with
790	only half the precision.
791	*/
792
793	if (isZero())
794	return valueOf(0.0);
795
796	if (isNegative()) {
797	return NaN;
798	}
799
800	double x = 1.0 / Math.sqrt(hi);
801	double ax = hi * x;
802
803	DD axdd = valueOf(ax);
804	DD diffSq = this.subtract(axdd.sqr());
805	double d2 = diffSq.hi * (x * 0.5);
806
807	return axdd.add(d2);
808	}
809
810	public static DD sqrt(double x)
811	{
812	return valueOf(x).sqrt();
813	}
814
815	/**
816	* Computes the value of this number raised to an integral power.
817	* Follows semantics of Java Math.pow as closely as possible.
818	*
819	* @param exp the integer exponent
820	* @return x raised to the integral power exp
821	*/
822	public DD pow(int exp)
823	{
824	if (exp == 0.0)
825	return valueOf(1.0);
826
827	DD r = new DD(this);
828	DD s = valueOf(1.0);
829	int n = Math.abs(exp);
830
831	if (n > 1) {
832	/* Use binary exponentiation */
833	while (n > 0) {
834	if (n % 2 == 1) {
835	s.selfMultiply(r);
836	}
837	n /= 2;
838	if (n > 0)
839	r = r.sqr();
840	}
841	} else {
842	s = r;
843	}
844
845	/* Compute the reciprocal if n is negative. */
846	if (exp < 0)
847	return s.reciprocal();
848	return s;
849	}
850
851	/**
852	* Computes the determinant of the 2x2 matrix with the given entries.
853	*
854	* @param x1 a double value
855	* @param y1 a double value
856	* @param x2 a double value
857	* @param y2 a double value
858	* @return the determinant of the values
859	*/
860	public static DD determinant(double x1, double y1, double x2, double y2)
861	{
862	return determinant(valueOf(x1), valueOf(y1), valueOf(x2), valueOf(y2) );
863	}
864
865	/**
866	* Computes the determinant of the 2x2 matrix with the given entries.
867	*
868	* @param x1 a matrix entry
869	* @param y1 a matrix entry
870	* @param x2 a matrix entry
871	* @param y2 a matrix entry
872	* @return the determinant of the matrix of values
873	*/
874	public static DD determinant(DD x1, DD y1, DD x2, DD y2)
875	{
876	DD det = x1.multiply(y2).selfSubtract(y1.multiply(x2));
877	return det;
878	}
879
880	/*------------------------------------------------------------
881	* Ordering Functions
882	*------------------------------------------------------------
883	*/
884
885	/**
886	* Computes the minimum of this and another DD number.
887	*
888	* @param x a DD number
889	* @return the minimum of the two numbers
890	*/
891	public DD min(DD x) {
892	if (this.le(x)) {
893	return this;
894	}
895	else {
896	return x;
897	}
898	}
899
900	/**
901	* Computes the maximum of this and another DD number.
902	*
903	* @param x a DD number
904	* @return the maximum of the two numbers
905	*/
906	public DD max(DD x) {
907	if (this.ge(x)) {
908	return this;
909	}
910	else {
911	return x;
912	}
913	}
914
915	/*------------------------------------------------------------
916	* Conversion Functions
917	*------------------------------------------------------------
918	*/
919
920	/**
921	* Converts this value to the nearest double-precision number.
922	*
923	* @return the nearest double-precision number to this value
924	*/
925	public double doubleValue()
926	{
927	return hi + lo;
928	}
929
930	/**
931	* Converts this value to the nearest integer.
932	*
933	* @return the nearest integer to this value
934	*/
935	public int intValue()
936	{
937	return (int) hi;
938	}
939
940	/*------------------------------------------------------------
941	* Predicates
942	*------------------------------------------------------------
943	*/
944
945	/**
946	* Tests whether this value is equal to 0.
947	*
948	* @return true if this value is equal to 0
949	*/
950	public boolean isZero()
951	{
952	return hi == 0.0 && lo == 0.0;
953	}
954
955	/**
956	* Tests whether this value is less than 0.
957	*
958	* @return true if this value is less than 0
959	*/
960	public boolean isNegative()
961	{
962	return hi < 0.0 \|\| (hi == 0.0 && lo < 0.0);
963	}
964
965	/**
966	* Tests whether this value is greater than 0.
967	*
968	* @return true if this value is greater than 0
969	*/
970	public boolean isPositive()
971	{
972	return hi > 0.0 \|\| (hi == 0.0 && lo > 0.0);
973	}
974
975	/**
976	* Tests whether this value is NaN.
977	*
978	* @return true if this value is NaN
979	*/
980	public boolean isNaN() { return Double.isNaN(hi); }
981
982	/**
983	* Tests whether this value is equal to another <tt>DoubleDouble</tt> value.
984	*
985	* @param y a DoubleDouble value
986	* @return true if this value = y
987	*/
988	public boolean equals(DD y)
989	{
990	return hi == y.hi && lo == y.lo;
991	}
992
993	/**
994	* Tests whether this value is greater than another <tt>DoubleDouble</tt> value.
995	* @param y a DoubleDouble value
996	* @return true if this value > y
997	*/
998	public boolean gt(DD y)
999	{
1000	return (hi > y.hi) \|\| (hi == y.hi && lo > y.lo);
1001	}
1002	/**
1003	* Tests whether this value is greater than or equals to another <tt>DoubleDouble</tt> value.
1004	* @param y a DoubleDouble value
1005	* @return true if this value >= y
1006	*/
1007	public boolean ge(DD y)
1008	{
1009	return (hi > y.hi) \|\| (hi == y.hi && lo >= y.lo);
1010	}
1011	/**
1012	* Tests whether this value is less than another <tt>DoubleDouble</tt> value.
1013	* @param y a DoubleDouble value
1014	* @return true if this value < y
1015	*/
1016	public boolean lt(DD y)
1017	{
1018	return (hi < y.hi) \|\| (hi == y.hi && lo < y.lo);
1019	}
1020	/**
1021	* Tests whether this value is less than or equal to another <tt>DoubleDouble</tt> value.
1022	* @param y a DoubleDouble value
1023	* @return true if this value <= y
1024	*/
1025	public boolean le(DD y)
1026	{
1027	return (hi < y.hi) \|\| (hi == y.hi && lo <= y.lo);
1028	}
1029
1030	/**
1031	* Compares two DoubleDouble objects numerically.
1032	*
1033	* @return -1,0 or 1 depending on whether this value is less than, equal to
1034	* or greater than the value of <tt>o</tt>
1035	*/
1036	public int compareTo(Object o)
1037	{
1038	DD other = (DD) o;
1039
1040	if (hi < other.hi) return -1;
1041	if (hi > other.hi) return 1;
1042	if (lo < other.lo) return -1;
1043	if (lo > other.lo) return 1;
1044	return 0;
1045	}
1046
1047
1048	/*------------------------------------------------------------
1049	* Output
1050	*------------------------------------------------------------
1051	*/
1052
1053	private static final int MAX_PRINT_DIGITS = 32;
1054	private static final DD TEN = DD.valueOf(10.0);
1055	private static final DD ONE = DD.valueOf(1.0);
1056	private static final String SCI_NOT_EXPONENT_CHAR = "E";
1057	private static final String SCI_NOT_ZERO = "0.0E0";
1058
1059	/**
1060	* Dumps the components of this number to a string.
1061	*
1062	* @return a string showing the components of the number
1063	*/
1064	public String dump()
1065	{
1066	return "DD<" + hi + ", " + lo + ">";
1067	}
1068
1069	/**
1070	* Returns a string representation of this number, in either standard or scientific notation.
1071	* If the magnitude of the number is in the range [ 10<sup>-3</sup>, 10<sup>8</sup> ]
1072	* standard notation will be used. Otherwise, scientific notation will be used.
1073	*
1074	* @return a string representation of this number
1075	*/
1076	public String toString()
1077	{
1078	int mag = magnitude(hi);
1079	if (mag >= -3 && mag <= 20)
1080	return toStandardNotation();
1081	return toSciNotation();
1082	}
1083
1084	/**
1085	* Returns the string representation of this value in standard notation.
1086	*
1087	* @return the string representation in standard notation
1088	*/
1089	public String toStandardNotation()
1090	{
1091	String specialStr = getSpecialNumberString();
1092	if (specialStr != null)
1093	return specialStr;
1094
1095	int[] magnitude = new int[1];
1096	String sigDigits = extractSignificantDigits(true, magnitude);
1097	int decimalPointPos = magnitude[0] + 1;
1098
1099	String num = sigDigits;
1100	// add a leading 0 if the decimal point is the first char
1101	if (sigDigits.charAt(0) == '.') {
1102	num = "0" + sigDigits;
1103	}
1104	else if (decimalPointPos < 0) {
1105	num = "0." + stringOfChar('0', -decimalPointPos) + sigDigits;
1106	}
1107	else if (sigDigits.indexOf('.') == -1) {
1108	// no point inserted - sig digits must be smaller than magnitude of number
1109	// add zeroes to end to make number the correct size
1110	int numZeroes = decimalPointPos - sigDigits.length();
1111	String zeroes = stringOfChar('0', numZeroes);
1112	num = sigDigits + zeroes + ".0";
1113	}
1114
1115	if (this.isNegative())
1116	return "-" + num;
1117	return num;
1118	}
1119
1120	/**
1121	* Returns the string representation of this value in scientific notation.
1122	*
1123	* @return the string representation in scientific notation
1124	*/
1125	public String toSciNotation()
1126	{
1127	// special case zero, to allow as
1128	if (isZero())
1129	return SCI_NOT_ZERO;
1130
1131	String specialStr = getSpecialNumberString();
1132	if (specialStr != null)
1133	return specialStr;
1134
1135	int[] magnitude = new int[1];
1136	String digits = extractSignificantDigits(false, magnitude);
1137	String expStr = SCI_NOT_EXPONENT_CHAR + magnitude[0];
1138
1139	// should never have leading zeroes
1140	// MD - is this correct? Or should we simply strip them if they are present?
1141	if (digits.charAt(0) == '0') {
1142	throw new IllegalStateException("Found leading zero: " + digits);
1143	}
1144
1145	// add decimal point
1146	String trailingDigits = "";
1147	if (digits.length() > 1)
1148	trailingDigits = digits.substring(1);
1149	String digitsWithDecimal = digits.charAt(0) + "." + trailingDigits;
1150
1151	if (this.isNegative())
1152	return "-" + digitsWithDecimal + expStr;
1153	return digitsWithDecimal + expStr;
1154	}
1155
1156
1157	/**
1158	* Extracts the significant digits in the decimal representation of the argument.
1159	* A decimal point may be optionally inserted in the string of digits
1160	* (as long as its position lies within the extracted digits
1161	* - if not, the caller must prepend or append the appropriate zeroes and decimal point).
1162	*
1163	* @param y the number to extract ( >= 0)
1164	* @param decimalPointPos the position in which to insert a decimal point
1165	* @return the string containing the significant digits and possibly a decimal point
1166	*/
1167	private String extractSignificantDigits(boolean insertDecimalPoint, int[] magnitude)
1168	{
1169	DD y = this.abs();
1170	// compute correct magnitude of y
1171	int mag = magnitude(y.hi);
1172	DD scale = TEN.pow(mag);
1173	y = y.divide(scale);
1174
1175	// fix magnitude if off by one
1176	if (y.gt(TEN)) {
1177	y = y.divide(TEN);
1178	mag += 1;
1179	}
1180	else if (y.lt(ONE)) {
1181	y = y.multiply(TEN);
1182	mag -= 1;
1183	}
1184
1185	int decimalPointPos = mag + 1;
1186	StringBuffer buf = new StringBuffer();
1187	int numDigits = MAX_PRINT_DIGITS - 1;
1188	for (int i = 0; i <= numDigits; i++) {
1189	if (insertDecimalPoint && i == decimalPointPos) {
1190	buf.append('.');
1191	}
1192	int digit = (int) y.hi;
1193	// System.out.println("printDump: [" + i + "] digit: " + digit + " y: " + y.dump() + " buf: " + buf);
1194
1195	/**
1196	* This should never happen, due to heuristic checks on remainder below
1197	*/
1198	if (digit < 0 \|\| digit > 9) {
1199	// System.out.println("digit > 10 : " + digit);
1200	// throw new IllegalStateException("Internal errror: found digit = " + digit);
1201	}
1202	/**
1203	* If a negative remainder is encountered, simply terminate the extraction.
1204	* This is robust, but maybe slightly inaccurate.
1205	* My current hypothesis is that negative remainders only occur for very small lo components,
1206	* so the inaccuracy is tolerable
1207	*/
1208	if (digit < 0) {
1209	break;
1210	// throw new IllegalStateException("Internal errror: found digit = " + digit);
1211	}
1212	boolean rebiasBy10 = false;
1213	char digitChar = 0;
1214	if (digit > 9) {
1215	// set flag to re-bias after next 10-shift
1216	rebiasBy10 = true;
1217	// output digit will end up being '9'
1218	digitChar = '9';
1219	}
1220	else {
1221	digitChar = (char) ('0' + digit);
1222	}
1223	buf.append(digitChar);
1224	y = (y.subtract(DD.valueOf(digit))
1225	.multiply(TEN));
1226	if (rebiasBy10)
1227	y.selfAdd(TEN);
1228
1229	boolean continueExtractingDigits = true;
1230	/**
1231	* Heuristic check: if the remaining portion of
1232	* y is non-positive, assume that output is complete
1233	*/
1234	// if (y.hi <= 0.0)
1235	// if (y.hi < 0.0)
1236	// continueExtractingDigits = false;
1237	/**
1238	* Check if remaining digits will be 0, and if so don't output them.
1239	* Do this by comparing the magnitude of the remainder with the expected precision.
1240	*/
1241	int remMag = magnitude(y.hi);
1242	if (remMag < 0 && Math.abs(remMag) >= (numDigits - i))
1243	continueExtractingDigits = false;
1244	if (! continueExtractingDigits)
1245	break;
1246	}
1247	magnitude[0] = mag;
1248	return buf.toString();
1249	}
1250
1251
1252	/**
1253	* Creates a string of a given length containing the given character
1254	*
1255	* @param ch the character to be repeated
1256	* @param len the len of the desired string
1257	* @return the string
1258	*/
1259	private static String stringOfChar(char ch, int len)
1260	{
1261	StringBuffer buf = new StringBuffer();
1262	for (int i = 0; i < len; i++) {
1263	buf.append(ch);
1264	}
1265	return buf.toString();
1266	}
1267
1268	/**
1269	* Returns the string for this value if it has a known representation.
1270	* (E.g. NaN or 0.0)
1271	*
1272	* @return the string for this special number
1273	* or null if the number is not a special number
1274	*/
1275	private String getSpecialNumberString()
1276	{
1277	if (isZero()) return "0.0";
1278	if (isNaN()) return "NaN ";
1279	return null;
1280	}
1281
1282
1283
1284	/**
1285	* Determines the decimal magnitude of a number.
1286	* The magnitude is the exponent of the greatest power of 10 which is less than
1287	* or equal to the number.
1288	*
1289	* @param x the number to find the magnitude of
1290	* @return the decimal magnitude of x
1291	*/
1292	private static int magnitude(double x)
1293	{
1294	double xAbs = Math.abs(x);
1295	double xLog10 = Math.log(xAbs) / Math.log(10);
1296	int xMag = (int) Math.floor(xLog10);
1297	/**
1298	* Since log computation is inexact, there may be an off-by-one error
1299	* in the computed magnitude.
1300	* Following tests that magnitude is correct, and adjusts it if not
1301	*/
1302	double xApprox = Math.pow(10, xMag);
1303	if (xApprox * 10 <= xAbs)
1304	xMag += 1;
1305
1306	return xMag;
1307	}
1308
1309
1310	/*------------------------------------------------------------
1311	* Input
1312	*------------------------------------------------------------
1313	*/
1314
1315	/**
1316	* Converts a string representation of a real number into a DoubleDouble value.
1317	* The format accepted is similar to the standard Java real number syntax.
1318	* It is defined by the following regular expression:
1319	* <pre>
1320	* [<tt>+</tt>\|<tt>-</tt>] {<i>digit</i>} [ <tt>.</tt> {<i>digit</i>} ] [ ( <tt>e</tt> \| <tt>E</tt> ) [<tt>+</tt>\|<tt>-</tt>] {<i>digit</i>}+
1321	* </pre>
1322	*
1323	* @param str the string to parse
1324	* @return the value of the parsed number
1325	* @throws NumberFormatException if <tt>str</tt> is not a valid representation of a number
1326	*/
1327	public static DD parse(String str)
1328	throws NumberFormatException
1329	{
1330	int i = 0;
1331	int strlen = str.length();
1332
1333	// skip leading whitespace
1334	while (Character.isWhitespace(str.charAt(i)))
1335	i++;
1336
1337	// check for sign
1338	boolean isNegative = false;
1339	if (i < strlen) {
1340	char signCh = str.charAt(i);
1341	if (signCh == '-' \|\| signCh == '+') {
1342	i++;
1343	if (signCh == '-') isNegative = true;
1344	}
1345	}
1346
1347	// scan all digits and accumulate into an integral value
1348	// Keep track of the location of the decimal point (if any) to allow scaling later
1349	DD val = new DD();
1350
1351	int numDigits = 0;
1352	int numBeforeDec = 0;
1353	int exp = 0;
1354	boolean hasDecimalChar = false;
1355	while (true) {
1356	if (i >= strlen)
1357	break;
1358	char ch = str.charAt(i);
1359	i++;
1360	if (Character.isDigit(ch)) {
1361	double d = ch - '0';
1362	val.selfMultiply(TEN);
1363	// MD: need to optimize this
1364	val.selfAdd(d);
1365	numDigits++;
1366	continue;
1367	}
1368	if (ch == '.') {
1369	numBeforeDec = numDigits;
1370	hasDecimalChar = true;
1371	continue;
1372	}
1373	if (ch == 'e' \|\| ch == 'E') {
1374	String expStr = str.substring(i);
1375	// this should catch any format problems with the exponent
1376	try {
1377	exp = Integer.parseInt(expStr);
1378	}
1379	catch (NumberFormatException ex) {
1380	throw new NumberFormatException("Invalid exponent " + expStr + " in string " + str);
1381	}
1382	break;
1383	}
1384	throw new NumberFormatException("Unexpected character '" + ch
1385	+ "' at position " + i
1386	+ " in string " + str);
1387	}
1388	DD val2 = val;
1389
1390	// correct number of digits before decimal sign if we don't have a decimal sign in the string
1391	if (!hasDecimalChar) numBeforeDec = numDigits;
1392
1393	// scale the number correctly
1394	int numDecPlaces = numDigits - numBeforeDec - exp;
1395	if (numDecPlaces == 0) {
1396	val2 = val;
1397	}
1398	else if (numDecPlaces > 0) {
1399	DD scale = TEN.pow(numDecPlaces);
1400	val2 = val.divide(scale);
1401	}
1402	else if (numDecPlaces < 0) {
1403	DD scale = TEN.pow(-numDecPlaces);
1404	val2 = val.multiply(scale);
1405	}
1406	// apply leading sign, if any
1407	if (isNegative) {
1408	return val2.negate();
1409	}
1410	return val2;
1411
1412	}
1413	}
1414