001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.dfp;
019    
020    import java.util.Arrays;
021    
022    import org.apache.commons.math.FieldElement;
023    
024    /**
025     *  Decimal floating point library for Java
026     *
027     *  <p>Another floating point class.  This one is built using radix 10000
028     *  which is 10<sup>4</sup>, so its almost decimal.</p>
029     *
030     *  <p>The design goals here are:
031     *  <ol>
032     *    <li>Decimal math, or close to it</li>
033     *    <li>Settable precision (but no mix between numbers using different settings)</li>
034     *    <li>Portability.  Code should be keep as portable as possible.</li>
035     *    <li>Performance</li>
036     *    <li>Accuracy  - Results should always be +/- 1 ULP for basic
037     *         algebraic operation</li>
038     *    <li>Comply with IEEE 854-1987 as much as possible.
039     *         (See IEEE 854-1987 notes below)</li>
040     *  </ol></p>
041     *
042     *  <p>Trade offs:
043     *  <ol>
044     *    <li>Memory foot print.  I'm using more memory than necessary to
045     *         represent numbers to get better performance.</li>
046     *    <li>Digits are bigger, so rounding is a greater loss.  So, if you
047     *         really need 12 decimal digits, better use 4 base 10000 digits
048     *         there can be one partially filled.</li>
049     *  </ol></p>
050     *
051     *  <p>Numbers are represented  in the following form:
052     *  <pre>
053     *  n  =  sign &times; mant &times; (radix)<sup>exp</sup>;</p>
054     *  </pre>
055     *  where sign is &plusmn;1, mantissa represents a fractional number between
056     *  zero and one.  mant[0] is the least significant digit.
057     *  exp is in the range of -32767 to 32768</p>
058     *
059     *  <p>IEEE 854-1987  Notes and differences</p>
060     *
061     *  <p>IEEE 854 requires the radix to be either 2 or 10.  The radix here is
062     *  10000, so that requirement is not met, but  it is possible that a
063     *  subclassed can be made to make it behave as a radix 10
064     *  number.  It is my opinion that if it looks and behaves as a radix
065     *  10 number then it is one and that requirement would be met.</p>
066     *
067     *  <p>The radix of 10000 was chosen because it should be faster to operate
068     *  on 4 decimal digits at once instead of one at a time.  Radix 10 behavior
069     *  can be realized by add an additional rounding step to ensure that
070     *  the number of decimal digits represented is constant.</p>
071     *
072     *  <p>The IEEE standard specifically leaves out internal data encoding,
073     *  so it is reasonable to conclude that such a subclass of this radix
074     *  10000 system is merely an encoding of a radix 10 system.</p>
075     *
076     *  <p>IEEE 854 also specifies the existence of "sub-normal" numbers.  This
077     *  class does not contain any such entities.  The most significant radix
078     *  10000 digit is always non-zero.  Instead, we support "gradual underflow"
079     *  by raising the underflow flag for numbers less with exponent less than
080     *  expMin, but don't flush to zero until the exponent reaches MIN_EXP-digits.
081     *  Thus the smallest number we can represent would be:
082     *  1E(-(MIN_EXP-digits-1)*4),  eg, for digits=5, MIN_EXP=-32767, that would
083     *  be 1e-131092.</p>
084     *
085     *  <p>IEEE 854 defines that the implied radix point lies just to the right
086     *  of the most significant digit and to the left of the remaining digits.
087     *  This implementation puts the implied radix point to the left of all
088     *  digits including the most significant one.  The most significant digit
089     *  here is the one just to the right of the radix point.  This is a fine
090     *  detail and is really only a matter of definition.  Any side effects of
091     *  this can be rendered invisible by a subclass.</p>
092     * @see DfpField
093     * @version $Revision: 1003889 $ $Date: 2010-10-02 23:11:55 +0200 (sam. 02 oct. 2010) $
094     * @since 2.2
095     */
096    public class Dfp implements FieldElement<Dfp> {
097    
098        /** The radix, or base of this system.  Set to 10000 */
099        public static final int RADIX = 10000;
100    
101        /** The minimum exponent before underflow is signaled.  Flush to zero
102         *  occurs at minExp-DIGITS */
103        public static final int MIN_EXP = -32767;
104    
105        /** The maximum exponent before overflow is signaled and results flushed
106         *  to infinity */
107        public static final int MAX_EXP =  32768;
108    
109        /** The amount under/overflows are scaled by before going to trap handler */
110        public static final int ERR_SCALE = 32760;
111    
112        /** Indicator value for normal finite numbers. */
113        public static final byte FINITE = 0;
114    
115        /** Indicator value for Infinity. */
116        public static final byte INFINITE = 1;
117    
118        /** Indicator value for signaling NaN. */
119        public static final byte SNAN = 2;
120    
121        /** Indicator value for quiet NaN. */
122        public static final byte QNAN = 3;
123    
124        /** String for NaN representation. */
125        private static final String NAN_STRING = "NaN";
126    
127        /** String for positive infinity representation. */
128        private static final String POS_INFINITY_STRING = "Infinity";
129    
130        /** String for negative infinity representation. */
131        private static final String NEG_INFINITY_STRING = "-Infinity";
132    
133        /** Name for traps triggered by addition. */
134        private static final String ADD_TRAP = "add";
135    
136        /** Name for traps triggered by multiplication. */
137        private static final String MULTIPLY_TRAP = "multiply";
138    
139        /** Name for traps triggered by division. */
140        private static final String DIVIDE_TRAP = "divide";
141    
142        /** Name for traps triggered by square root. */
143        private static final String SQRT_TRAP = "sqrt";
144    
145        /** Name for traps triggered by alignment. */
146        private static final String ALIGN_TRAP = "align";
147    
148        /** Name for traps triggered by truncation. */
149        private static final String TRUNC_TRAP = "trunc";
150    
151        /** Name for traps triggered by nextAfter. */
152        private static final String NEXT_AFTER_TRAP = "nextAfter";
153    
154        /** Name for traps triggered by lessThan. */
155        private static final String LESS_THAN_TRAP = "lessThan";
156    
157        /** Name for traps triggered by greaterThan. */
158        private static final String GREATER_THAN_TRAP = "greaterThan";
159    
160        /** Name for traps triggered by newInstance. */
161        private static final String NEW_INSTANCE_TRAP = "newInstance";
162    
163        /** Mantissa. */
164        protected int[] mant;
165    
166        /** Sign bit: & for positive, -1 for negative. */
167        protected byte sign;
168    
169        /** Exponent. */
170        protected int exp;
171    
172        /** Indicator for non-finite / non-number values. */
173        protected byte nans;
174    
175        /** Factory building similar Dfp's. */
176        private final DfpField field;
177    
178        /** Makes an instance with a value of zero.
179         * @param field field to which this instance belongs
180         */
181        protected Dfp(final DfpField field) {
182            mant = new int[field.getRadixDigits()];
183            sign = 1;
184            exp = 0;
185            nans = FINITE;
186            this.field = field;
187        }
188    
189        /** Create an instance from a byte value.
190         * @param field field to which this instance belongs
191         * @param x value to convert to an instance
192         */
193        protected Dfp(final DfpField field, byte x) {
194            this(field, (long) x);
195        }
196    
197        /** Create an instance from an int value.
198         * @param field field to which this instance belongs
199         * @param x value to convert to an instance
200         */
201        protected Dfp(final DfpField field, int x) {
202            this(field, (long) x);
203        }
204    
205        /** Create an instance from a long value.
206         * @param field field to which this instance belongs
207         * @param x value to convert to an instance
208         */
209        protected Dfp(final DfpField field, long x) {
210    
211            // initialize as if 0
212            mant = new int[field.getRadixDigits()];
213            nans = FINITE;
214            this.field = field;
215    
216            boolean isLongMin = false;
217            if (x == Long.MIN_VALUE) {
218                // special case for Long.MIN_VALUE (-9223372036854775808)
219                // we must shift it before taking its absolute value
220                isLongMin = true;
221                ++x;
222            }
223    
224            // set the sign
225            if (x < 0) {
226                sign = -1;
227                x = -x;
228            } else {
229                sign = 1;
230            }
231    
232            exp = 0;
233            while (x != 0) {
234                System.arraycopy(mant, mant.length - exp, mant, mant.length - 1 - exp, exp);
235                mant[mant.length - 1] = (int) (x % RADIX);
236                x /= RADIX;
237                exp++;
238            }
239    
240            if (isLongMin) {
241                // remove the shift added for Long.MIN_VALUE
242                // we know in this case that fixing the last digit is sufficient
243                for (int i = 0; i < mant.length - 1; i++) {
244                    if (mant[i] != 0) {
245                        mant[i]++;
246                        break;
247                    }
248                }
249            }
250        }
251    
252        /** Create an instance from a double value.
253         * @param field field to which this instance belongs
254         * @param x value to convert to an instance
255         */
256        protected Dfp(final DfpField field, double x) {
257    
258            // initialize as if 0
259            mant = new int[field.getRadixDigits()];
260            sign = 1;
261            exp = 0;
262            nans = FINITE;
263            this.field = field;
264    
265            long bits = Double.doubleToLongBits(x);
266            long mantissa = bits & 0x000fffffffffffffL;
267            int exponent = (int) ((bits & 0x7ff0000000000000L) >> 52) - 1023;
268    
269            if (exponent == -1023) {
270                // Zero or sub-normal
271                if (x == 0) {
272                    return;
273                }
274    
275                exponent++;
276    
277                // Normalize the subnormal number
278                while ( (mantissa & 0x0010000000000000L) == 0) {
279                    exponent--;
280                    mantissa <<= 1;
281                }
282                mantissa &= 0x000fffffffffffffL;
283            }
284    
285            if (exponent == 1024) {
286                // infinity or NAN
287                if (x != x) {
288                    sign = (byte) 1;
289                    nans = QNAN;
290                } else if (x < 0) {
291                    sign = (byte) -1;
292                    nans = INFINITE;
293                } else {
294                    sign = (byte) 1;
295                    nans = INFINITE;
296                }
297                return;
298            }
299    
300            Dfp xdfp = new Dfp(field, mantissa);
301            xdfp = xdfp.divide(new Dfp(field, 4503599627370496l)).add(field.getOne());  // Divide by 2^52, then add one
302            xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
303    
304            if ((bits & 0x8000000000000000L) != 0) {
305                xdfp = xdfp.negate();
306            }
307    
308            System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
309            sign = xdfp.sign;
310            exp  = xdfp.exp;
311            nans = xdfp.nans;
312    
313        }
314    
315        /** Copy constructor.
316         * @param d instance to copy
317         */
318        public Dfp(final Dfp d) {
319            mant  = d.mant.clone();
320            sign  = d.sign;
321            exp   = d.exp;
322            nans  = d.nans;
323            field = d.field;
324        }
325    
326        /** Create an instance from a String representation.
327         * @param field field to which this instance belongs
328         * @param s string representation of the instance
329         */
330        protected Dfp(final DfpField field, final String s) {
331    
332            // initialize as if 0
333            mant = new int[field.getRadixDigits()];
334            sign = 1;
335            exp = 0;
336            nans = FINITE;
337            this.field = field;
338    
339            boolean decimalFound = false;
340            final int rsize = 4;   // size of radix in decimal digits
341            final int offset = 4;  // Starting offset into Striped
342            final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
343    
344            // Check some special cases
345            if (s.equals(POS_INFINITY_STRING)) {
346                sign = (byte) 1;
347                nans = INFINITE;
348                return;
349            }
350    
351            if (s.equals(NEG_INFINITY_STRING)) {
352                sign = (byte) -1;
353                nans = INFINITE;
354                return;
355            }
356    
357            if (s.equals(NAN_STRING)) {
358                sign = (byte) 1;
359                nans = QNAN;
360                return;
361            }
362    
363            // Check for scientific notation
364            int p = s.indexOf("e");
365            if (p == -1) { // try upper case?
366                p = s.indexOf("E");
367            }
368    
369            final String fpdecimal;
370            int sciexp = 0;
371            if (p != -1) {
372                // scientific notation
373                fpdecimal = s.substring(0, p);
374                String fpexp = s.substring(p+1);
375                boolean negative = false;
376    
377                for (int i=0; i<fpexp.length(); i++)
378                {
379                    if (fpexp.charAt(i) == '-')
380                    {
381                        negative = true;
382                        continue;
383                    }
384                    if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9')
385                        sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
386                }
387    
388                if (negative) {
389                    sciexp = -sciexp;
390                }
391            } else {
392                // normal case
393                fpdecimal = s;
394            }
395    
396            // If there is a minus sign in the number then it is negative
397            if (fpdecimal.indexOf("-") !=  -1) {
398                sign = -1;
399            }
400    
401            // First off, find all of the leading zeros, trailing zeros, and significant digits
402            p = 0;
403    
404            // Move p to first significant digit
405            int decimalPos = 0;
406            for (;;) {
407                if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
408                    break;
409                }
410    
411                if (decimalFound && fpdecimal.charAt(p) == '0') {
412                    decimalPos--;
413                }
414    
415                if (fpdecimal.charAt(p) == '.') {
416                    decimalFound = true;
417                }
418    
419                p++;
420    
421                if (p == fpdecimal.length()) {
422                    break;
423                }
424            }
425    
426            // Copy the string onto Stripped
427            int q = offset;
428            striped[0] = '0';
429            striped[1] = '0';
430            striped[2] = '0';
431            striped[3] = '0';
432            int significantDigits=0;
433            for(;;) {
434                if (p == (fpdecimal.length())) {
435                    break;
436                }
437    
438                // Don't want to run pass the end of the array
439                if (q == mant.length*rsize+offset+1) {
440                    break;
441                }
442    
443                if (fpdecimal.charAt(p) == '.') {
444                    decimalFound = true;
445                    decimalPos = significantDigits;
446                    p++;
447                    continue;
448                }
449    
450                if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
451                    p++;
452                    continue;
453                }
454    
455                striped[q] = fpdecimal.charAt(p);
456                q++;
457                p++;
458                significantDigits++;
459            }
460    
461    
462            // If the decimal point has been found then get rid of trailing zeros.
463            if (decimalFound && q != offset) {
464                for (;;) {
465                    q--;
466                    if (q == offset) {
467                        break;
468                    }
469                    if (striped[q] == '0') {
470                        significantDigits--;
471                    } else {
472                        break;
473                    }
474                }
475            }
476    
477            // special case of numbers like "0.00000"
478            if (decimalFound && significantDigits == 0) {
479                decimalPos = 0;
480            }
481    
482            // Implicit decimal point at end of number if not present
483            if (!decimalFound) {
484                decimalPos = q-offset;
485            }
486    
487            // Find the number of significant trailing zeros
488            q = offset;  // set q to point to first sig digit
489            p = significantDigits-1+offset;
490    
491            int trailingZeros = 0;
492            while (p > q) {
493                if (striped[p] != '0') {
494                    break;
495                }
496                trailingZeros++;
497                p--;
498            }
499    
500            // Make sure the decimal is on a mod 10000 boundary
501            int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
502            q -= i;
503            decimalPos += i;
504    
505            // Make the mantissa length right by adding zeros at the end if necessary
506            while ((p - q) < (mant.length * rsize)) {
507                for (i = 0; i < rsize; i++) {
508                    striped[++p] = '0';
509                }
510            }
511    
512            // Ok, now we know how many trailing zeros there are,
513            // and where the least significant digit is
514            for (i = mant.length - 1; i >= 0; i--) {
515                mant[i] = (striped[q]   - '0') * 1000 +
516                          (striped[q+1] - '0') * 100  +
517                          (striped[q+2] - '0') * 10   +
518                          (striped[q+3] - '0');
519                q += 4;
520            }
521    
522    
523            exp = (decimalPos+sciexp) / rsize;
524    
525            if (q < striped.length) {
526                // Is there possible another digit?
527                round((striped[q] - '0')*1000);
528            }
529    
530        }
531    
532        /** Creates an instance with a non-finite value.
533         * @param field field to which this instance belongs
534         * @param sign sign of the Dfp to create
535         * @param nans code of the value, must be one of {@link #INFINITE},
536         * {@link #SNAN},  {@link #QNAN}
537         */
538        protected Dfp(final DfpField field, final byte sign, final byte nans) {
539            this.field = field;
540            this.mant    = new int[field.getRadixDigits()];
541            this.sign    = sign;
542            this.exp     = 0;
543            this.nans    = nans;
544        }
545    
546        /** Create an instance with a value of 0.
547         * Use this internally in preference to constructors to facilitate subclasses
548         * @return a new instance with a value of 0
549         */
550        public Dfp newInstance() {
551            return new Dfp(getField());
552        }
553    
554        /** Create an instance from a byte value.
555         * @param x value to convert to an instance
556         * @return a new instance with value x
557         */
558        public Dfp newInstance(final byte x) {
559            return new Dfp(getField(), x);
560        }
561    
562        /** Create an instance from an int value.
563         * @param x value to convert to an instance
564         * @return a new instance with value x
565         */
566        public Dfp newInstance(final int x) {
567            return new Dfp(getField(), x);
568        }
569    
570        /** Create an instance from a long value.
571         * @param x value to convert to an instance
572         * @return a new instance with value x
573         */
574        public Dfp newInstance(final long x) {
575            return new Dfp(getField(), x);
576        }
577    
578        /** Create an instance from a double value.
579         * @param x value to convert to an instance
580         * @return a new instance with value x
581         */
582        public Dfp newInstance(final double x) {
583            return new Dfp(getField(), x);
584        }
585    
586        /** Create an instance by copying an existing one.
587         * Use this internally in preference to constructors to facilitate subclasses.
588         * @param d instance to copy
589         * @return a new instance with the same value as d
590         */
591        public Dfp newInstance(final Dfp d) {
592    
593            // make sure we don't mix number with different precision
594            if (field.getRadixDigits() != d.field.getRadixDigits()) {
595                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
596                final Dfp result = newInstance(getZero());
597                result.nans = QNAN;
598                return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
599            }
600    
601            return new Dfp(d);
602    
603        }
604    
605        /** Create an instance from a String representation.
606         * Use this internally in preference to constructors to facilitate subclasses.
607         * @param s string representation of the instance
608         * @return a new instance parsed from specified string
609         */
610        public Dfp newInstance(final String s) {
611            return new Dfp(field, s);
612        }
613    
614        /** Creates an instance with a non-finite value.
615         * @param sig sign of the Dfp to create
616         * @param code code of the value, must be one of {@link #INFINITE},
617         * {@link #SNAN},  {@link #QNAN}
618         * @return a new instance with a non-finite value
619         */
620        public Dfp newInstance(final byte sig, final byte code) {
621            return field.newDfp(sig, code);
622        }
623    
624        /** Get the {@link org.apache.commons.math.Field Field} (really a {@link DfpField}) to which the instance belongs.
625         * <p>
626         * The field is linked to the number of digits and acts as a factory
627         * for {@link Dfp} instances.
628         * </p>
629         * @return {@link org.apache.commons.math.Field Field} (really a {@link DfpField}) to which the instance belongs
630         */
631        public DfpField getField() {
632            return field;
633        }
634    
635        /** Get the number of radix digits of the instance.
636         * @return number of radix digits
637         */
638        public int getRadixDigits() {
639            return field.getRadixDigits();
640        }
641    
642        /** Get the constant 0.
643         * @return a Dfp with value zero
644         */
645        public Dfp getZero() {
646            return field.getZero();
647        }
648    
649        /** Get the constant 1.
650         * @return a Dfp with value one
651         */
652        public Dfp getOne() {
653            return field.getOne();
654        }
655    
656        /** Get the constant 2.
657         * @return a Dfp with value two
658         */
659        public Dfp getTwo() {
660            return field.getTwo();
661        }
662    
663        /** Shift the mantissa left, and adjust the exponent to compensate.
664         */
665        protected void shiftLeft() {
666            for (int i = mant.length - 1; i > 0; i--) {
667                mant[i] = mant[i-1];
668            }
669            mant[0] = 0;
670            exp--;
671        }
672    
673        /* Note that shiftRight() does not call round() as that round() itself
674         uses shiftRight() */
675        /** Shift the mantissa right, and adjust the exponent to compensate.
676         */
677        protected void shiftRight() {
678            for (int i = 0; i < mant.length - 1; i++) {
679                mant[i] = mant[i+1];
680            }
681            mant[mant.length - 1] = 0;
682            exp++;
683        }
684    
685        /** Make our exp equal to the supplied one, this may cause rounding.
686         *  Also causes de-normalized numbers.  These numbers are generally
687         *  dangerous because most routines assume normalized numbers.
688         *  Align doesn't round, so it will return the last digit destroyed
689         *  by shifting right.
690         *  @param e desired exponent
691         *  @return last digit destroyed by shifting right
692         */
693        protected int align(int e) {
694            int lostdigit = 0;
695            boolean inexact = false;
696    
697            int diff = exp - e;
698    
699            int adiff = diff;
700            if (adiff < 0) {
701                adiff = -adiff;
702            }
703    
704            if (diff == 0) {
705                return 0;
706            }
707    
708            if (adiff > (mant.length + 1)) {
709                // Special case
710                Arrays.fill(mant, 0);
711                exp = e;
712    
713                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
714                dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
715    
716                return 0;
717            }
718    
719            for (int i = 0; i < adiff; i++) {
720                if (diff < 0) {
721                    /* Keep track of loss -- only signal inexact after losing 2 digits.
722                     * the first lost digit is returned to add() and may be incorporated
723                     * into the result.
724                     */
725                    if (lostdigit != 0) {
726                        inexact = true;
727                    }
728    
729                    lostdigit = mant[0];
730    
731                    shiftRight();
732                } else {
733                    shiftLeft();
734                }
735            }
736    
737            if (inexact) {
738                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
739                dotrap(DfpField.FLAG_INEXACT, ALIGN_TRAP, this, this);
740            }
741    
742            return lostdigit;
743    
744        }
745    
746        /** Check if instance is less than x.
747         * @param x number to check instance against
748         * @return true if instance is less than x and neither are NaN, false otherwise
749         */
750        public boolean lessThan(final Dfp x) {
751    
752            // make sure we don't mix number with different precision
753            if (field.getRadixDigits() != x.field.getRadixDigits()) {
754                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
755                final Dfp result = newInstance(getZero());
756                result.nans = QNAN;
757                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
758                return false;
759            }
760    
761            /* if a nan is involved, signal invalid and return false */
762            if (isNaN() || x.isNaN()) {
763                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
764                dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
765                return false;
766            }
767    
768            return compare(this, x) < 0;
769        }
770    
771        /** Check if instance is greater than x.
772         * @param x number to check instance against
773         * @return true if instance is greater than x and neither are NaN, false otherwise
774         */
775        public boolean greaterThan(final Dfp x) {
776    
777            // make sure we don't mix number with different precision
778            if (field.getRadixDigits() != x.field.getRadixDigits()) {
779                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
780                final Dfp result = newInstance(getZero());
781                result.nans = QNAN;
782                dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
783                return false;
784            }
785    
786            /* if a nan is involved, signal invalid and return false */
787            if (isNaN() || x.isNaN()) {
788                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
789                dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
790                return false;
791            }
792    
793            return compare(this, x) > 0;
794        }
795    
796        /** Check if instance is infinite.
797         * @return true if instance is infinite
798         */
799        public boolean isInfinite() {
800            return nans == INFINITE;
801        }
802    
803        /** Check if instance is not a number.
804         * @return true if instance is not a number
805         */
806        public boolean isNaN() {
807            return (nans == QNAN) || (nans == SNAN);
808        }
809    
810        /** Check if instance is equal to x.
811         * @param other object to check instance against
812         * @return true if instance is equal to x and neither are NaN, false otherwise
813         */
814        @Override
815        public boolean equals(final Object other) {
816    
817            if (other instanceof Dfp) {
818                final Dfp x = (Dfp) other;
819                if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
820                    return false;
821                }
822    
823                return compare(this, x) == 0;
824            }
825    
826            return false;
827    
828        }
829    
830        /**
831         * Gets a hashCode for the instance.
832         * @return a hash code value for this object
833         */
834        @Override
835        public int hashCode() {
836            return 17 + (sign << 8) + (nans << 16) + exp + Arrays.hashCode(mant);
837        }
838    
839        /** Check if instance is not equal to x.
840         * @param x number to check instance against
841         * @return true if instance is not equal to x and neither are NaN, false otherwise
842         */
843        public boolean unequal(final Dfp x) {
844            if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
845                return false;
846            }
847    
848            return greaterThan(x) || lessThan(x);
849        }
850    
851        /** Compare two instances.
852         * @param a first instance in comparison
853         * @param b second instance in comparison
854         * @return -1 if a<b, 1 if a>b and 0 if a==b
855         *  Note this method does not properly handle NaNs or numbers with different precision.
856         */
857        private static int compare(final Dfp a, final Dfp b) {
858            // Ignore the sign of zero
859            if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
860                a.nans == FINITE && b.nans == FINITE) {
861                return 0;
862            }
863    
864            if (a.sign != b.sign) {
865                if (a.sign == -1) {
866                    return -1;
867                } else {
868                    return 1;
869                }
870            }
871    
872            // deal with the infinities
873            if (a.nans == INFINITE && b.nans == FINITE) {
874                return a.sign;
875            }
876    
877            if (a.nans == FINITE && b.nans == INFINITE) {
878                return -b.sign;
879            }
880    
881            if (a.nans == INFINITE && b.nans == INFINITE) {
882                return 0;
883            }
884    
885            // Handle special case when a or b is zero, by ignoring the exponents
886            if (b.mant[b.mant.length-1] != 0 && a.mant[b.mant.length-1] != 0) {
887                if (a.exp < b.exp) {
888                    return -a.sign;
889                }
890    
891                if (a.exp > b.exp) {
892                    return a.sign;
893                }
894            }
895    
896            // compare the mantissas
897            for (int i = a.mant.length - 1; i >= 0; i--) {
898                if (a.mant[i] > b.mant[i]) {
899                    return a.sign;
900                }
901    
902                if (a.mant[i] < b.mant[i]) {
903                    return -a.sign;
904                }
905            }
906    
907            return 0;
908    
909        }
910    
911        /** Round to nearest integer using the round-half-even method.
912         *  That is round to nearest integer unless both are equidistant.
913         *  In which case round to the even one.
914         *  @return rounded value
915         */
916        public Dfp rint() {
917            return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
918        }
919    
920        /** Round to an integer using the round floor mode.
921         * That is, round toward -Infinity
922         *  @return rounded value
923         */
924        public Dfp floor() {
925            return trunc(DfpField.RoundingMode.ROUND_FLOOR);
926        }
927    
928        /** Round to an integer using the round ceil mode.
929         * That is, round toward +Infinity
930         *  @return rounded value
931         */
932        public Dfp ceil() {
933            return trunc(DfpField.RoundingMode.ROUND_CEIL);
934        }
935    
936        /** Returns the IEEE remainder.
937         * @param d divisor
938         * @return this less n &times; d, where n is the integer closest to this/d
939         */
940        public Dfp remainder(final Dfp d) {
941    
942            final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
943    
944            // IEEE 854-1987 says that if the result is zero, then it carries the sign of this
945            if (result.mant[mant.length-1] == 0) {
946                result.sign = sign;
947            }
948    
949            return result;
950    
951        }
952    
953        /** Does the integer conversions with the specified rounding.
954         * @param rmode rounding mode to use
955         * @return truncated value
956         */
957        protected Dfp trunc(final DfpField.RoundingMode rmode) {
958            boolean changed = false;
959    
960            if (isNaN()) {
961                return newInstance(this);
962            }
963    
964            if (nans == INFINITE) {
965                return newInstance(this);
966            }
967    
968            if (mant[mant.length-1] == 0) {
969                // a is zero
970                return newInstance(this);
971            }
972    
973            /* If the exponent is less than zero then we can certainly
974             * return zero */
975            if (exp < 0) {
976                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
977                Dfp result = newInstance(getZero());
978                result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
979                return result;
980            }
981    
982            /* If the exponent is greater than or equal to digits, then it
983             * must already be an integer since there is no precision left
984             * for any fractional part */
985    
986            if (exp >= mant.length) {
987                return newInstance(this);
988            }
989    
990            /* General case:  create another dfp, result, that contains the
991             * a with the fractional part lopped off.  */
992    
993            Dfp result = newInstance(this);
994            for (int i = 0; i < mant.length-result.exp; i++) {
995                changed |= result.mant[i] != 0;
996                result.mant[i] = 0;
997            }
998    
999            if (changed) {
1000                switch (rmode) {
1001                    case ROUND_FLOOR:
1002                        if (result.sign == -1) {
1003                            // then we must increment the mantissa by one
1004                            result = result.add(newInstance(-1));
1005                        }
1006                        break;
1007    
1008                    case ROUND_CEIL:
1009                        if (result.sign == 1) {
1010                            // then we must increment the mantissa by one
1011                            result = result.add(getOne());
1012                        }
1013                        break;
1014    
1015                    case ROUND_HALF_EVEN:
1016                    default:
1017                        final Dfp half = newInstance("0.5");
1018                        Dfp a = subtract(result);  // difference between this and result
1019                        a.sign = 1;            // force positive (take abs)
1020                        if (a.greaterThan(half)) {
1021                            a = newInstance(getOne());
1022                            a.sign = sign;
1023                            result = result.add(a);
1024                        }
1025    
1026                        /** If exactly equal to 1/2 and odd then increment */
1027                        if (a.equals(half) && result.exp > 0 && (result.mant[mant.length-result.exp]&1) != 0) {
1028                            a = newInstance(getOne());
1029                            a.sign = sign;
1030                            result = result.add(a);
1031                        }
1032                        break;
1033                }
1034    
1035                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);  // signal inexact
1036                result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
1037                return result;
1038            }
1039    
1040            return result;
1041        }
1042    
1043        /** Convert this to an integer.
1044         * If greater than 2147483647, it returns 2147483647. If less than -2147483648 it returns -2147483648.
1045         * @return converted number
1046         */
1047        public int intValue() {
1048            Dfp rounded;
1049            int result = 0;
1050    
1051            rounded = rint();
1052    
1053            if (rounded.greaterThan(newInstance(2147483647))) {
1054                return 2147483647;
1055            }
1056    
1057            if (rounded.lessThan(newInstance(-2147483648))) {
1058                return -2147483648;
1059            }
1060    
1061            for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
1062                result = result * RADIX + rounded.mant[i];
1063            }
1064    
1065            if (rounded.sign == -1) {
1066                result = -result;
1067            }
1068    
1069            return result;
1070        }
1071    
1072        /** Get the exponent of the greatest power of 10000 that is
1073         *  less than or equal to the absolute value of this.  I.E.  if
1074         *  this is 10<sup>6</sup> then log10K would return 1.
1075         *  @return integer base 10000 logarithm
1076         */
1077        public int log10K() {
1078            return exp - 1;
1079        }
1080    
1081        /** Get the specified  power of 10000.
1082         * @param e desired power
1083         * @return 10000<sup>e</sup>
1084         */
1085        public Dfp power10K(final int e) {
1086            Dfp d = newInstance(getOne());
1087            d.exp = e + 1;
1088            return d;
1089        }
1090    
1091        /** Get the exponent of the greatest power of 10 that is less than or equal to abs(this).
1092         *  @return integer base 10 logarithm
1093         */
1094        public int log10()  {
1095            if (mant[mant.length-1] > 1000) {
1096                return exp * 4 - 1;
1097            }
1098            if (mant[mant.length-1] > 100) {
1099                return exp * 4 - 2;
1100            }
1101            if (mant[mant.length-1] > 10) {
1102                return exp * 4 - 3;
1103            }
1104            return exp * 4 - 4;
1105        }
1106    
1107        /** Return the specified  power of 10.
1108         * @param e desired power
1109         * @return 10<sup>e</sup>
1110         */
1111        public Dfp power10(final int e) {
1112            Dfp d = newInstance(getOne());
1113    
1114            if (e >= 0) {
1115                d.exp = e / 4 + 1;
1116            } else {
1117                d.exp = (e + 1) / 4;
1118            }
1119    
1120            switch ((e % 4 + 4) % 4) {
1121                case 0:
1122                    break;
1123                case 1:
1124                    d = d.multiply(10);
1125                    break;
1126                case 2:
1127                    d = d.multiply(100);
1128                    break;
1129                default:
1130                    d = d.multiply(1000);
1131            }
1132    
1133            return d;
1134        }
1135    
1136        /** Negate the mantissa of this by computing the complement.
1137         *  Leaves the sign bit unchanged, used internally by add.
1138         *  Denormalized numbers are handled properly here.
1139         *  @param extra ???
1140         *  @return ???
1141         */
1142        protected int complement(int extra) {
1143    
1144            extra = RADIX-extra;
1145            for (int i = 0; i < mant.length; i++) {
1146                mant[i] = RADIX-mant[i]-1;
1147            }
1148    
1149            int rh = extra / RADIX;
1150            extra = extra - rh * RADIX;
1151            for (int i = 0; i < mant.length; i++) {
1152                final int r = mant[i] + rh;
1153                rh = r / RADIX;
1154                mant[i] = r - rh * RADIX;
1155            }
1156    
1157            return extra;
1158        }
1159    
1160        /** Add x to this.
1161         * @param x number to add
1162         * @return sum of this and x
1163         */
1164        public Dfp add(final Dfp x) {
1165    
1166            // make sure we don't mix number with different precision
1167            if (field.getRadixDigits() != x.field.getRadixDigits()) {
1168                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1169                final Dfp result = newInstance(getZero());
1170                result.nans = QNAN;
1171                return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1172            }
1173    
1174            /* handle special cases */
1175            if (nans != FINITE || x.nans != FINITE) {
1176                if (isNaN()) {
1177                    return this;
1178                }
1179    
1180                if (x.isNaN()) {
1181                    return x;
1182                }
1183    
1184                if (nans == INFINITE && x.nans == FINITE) {
1185                    return this;
1186                }
1187    
1188                if (x.nans == INFINITE && nans == FINITE) {
1189                    return x;
1190                }
1191    
1192                if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
1193                    return x;
1194                }
1195    
1196                if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
1197                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1198                    Dfp result = newInstance(getZero());
1199                    result.nans = QNAN;
1200                    result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1201                    return result;
1202                }
1203            }
1204    
1205            /* copy this and the arg */
1206            Dfp a = newInstance(this);
1207            Dfp b = newInstance(x);
1208    
1209            /* initialize the result object */
1210            Dfp result = newInstance(getZero());
1211    
1212            /* Make all numbers positive, but remember their sign */
1213            final byte asign = a.sign;
1214            final byte bsign = b.sign;
1215    
1216            a.sign = 1;
1217            b.sign = 1;
1218    
1219            /* The result will be signed like the arg with greatest magnitude */
1220            byte rsign = bsign;
1221            if (compare(a, b) > 0) {
1222                rsign = asign;
1223            }
1224    
1225            /* Handle special case when a or b is zero, by setting the exponent
1226           of the zero number equal to the other one.  This avoids an alignment
1227           which would cause catastropic loss of precision */
1228            if (b.mant[mant.length-1] == 0) {
1229                b.exp = a.exp;
1230            }
1231    
1232            if (a.mant[mant.length-1] == 0) {
1233                a.exp = b.exp;
1234            }
1235    
1236            /* align number with the smaller exponent */
1237            int aextradigit = 0;
1238            int bextradigit = 0;
1239            if (a.exp < b.exp) {
1240                aextradigit = a.align(b.exp);
1241            } else {
1242                bextradigit = b.align(a.exp);
1243            }
1244    
1245            /* complement the smaller of the two if the signs are different */
1246            if (asign != bsign) {
1247                if (asign == rsign) {
1248                    bextradigit = b.complement(bextradigit);
1249                } else {
1250                    aextradigit = a.complement(aextradigit);
1251                }
1252            }
1253    
1254            /* add the mantissas */
1255            int rh = 0; /* acts as a carry */
1256            for (int i = 0; i < mant.length; i++) {
1257                final int r = a.mant[i]+b.mant[i]+rh;
1258                rh = r / RADIX;
1259                result.mant[i] = r - rh * RADIX;
1260            }
1261            result.exp = a.exp;
1262            result.sign = rsign;
1263    
1264            /* handle overflow -- note, when asign!=bsign an overflow is
1265             * normal and should be ignored.  */
1266    
1267            if (rh != 0 && (asign == bsign)) {
1268                final int lostdigit = result.mant[0];
1269                result.shiftRight();
1270                result.mant[mant.length-1] = rh;
1271                final int excp = result.round(lostdigit);
1272                if (excp != 0) {
1273                    result = dotrap(excp, ADD_TRAP, x, result);
1274                }
1275            }
1276    
1277            /* normalize the result */
1278            for (int i = 0; i < mant.length; i++) {
1279                if (result.mant[mant.length-1] != 0) {
1280                    break;
1281                }
1282                result.shiftLeft();
1283                if (i == 0) {
1284                    result.mant[0] = aextradigit+bextradigit;
1285                    aextradigit = 0;
1286                    bextradigit = 0;
1287                }
1288            }
1289    
1290            /* result is zero if after normalization the most sig. digit is zero */
1291            if (result.mant[mant.length-1] == 0) {
1292                result.exp = 0;
1293    
1294                if (asign != bsign) {
1295                    // Unless adding 2 negative zeros, sign is positive
1296                    result.sign = 1;  // Per IEEE 854-1987 Section 6.3
1297                }
1298            }
1299    
1300            /* Call round to test for over/under flows */
1301            final int excp = result.round(aextradigit + bextradigit);
1302            if (excp != 0) {
1303                result = dotrap(excp, ADD_TRAP, x, result);
1304            }
1305    
1306            return result;
1307        }
1308    
1309        /** Returns a number that is this number with the sign bit reversed.
1310         * @return the opposite of this
1311         */
1312        public Dfp negate() {
1313            Dfp result = newInstance(this);
1314            result.sign = (byte) - result.sign;
1315            return result;
1316        }
1317    
1318        /** Subtract x from this.
1319         * @param x number to subtract
1320         * @return difference of this and a
1321         */
1322        public Dfp subtract(final Dfp x) {
1323            return add(x.negate());
1324        }
1325    
1326        /** Round this given the next digit n using the current rounding mode.
1327         * @param n ???
1328         * @return the IEEE flag if an exception occurred
1329         */
1330        protected int round(int n) {
1331            boolean inc = false;
1332            switch (field.getRoundingMode()) {
1333                case ROUND_DOWN:
1334                    inc = false;
1335                    break;
1336    
1337                case ROUND_UP:
1338                    inc = n != 0;       // round up if n!=0
1339                    break;
1340    
1341                case ROUND_HALF_UP:
1342                    inc = n >= 5000;  // round half up
1343                    break;
1344    
1345                case ROUND_HALF_DOWN:
1346                    inc = n > 5000;  // round half down
1347                    break;
1348    
1349                case ROUND_HALF_EVEN:
1350                    inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1);  // round half-even
1351                    break;
1352    
1353                case ROUND_HALF_ODD:
1354                    inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0);  // round half-odd
1355                    break;
1356    
1357                case ROUND_CEIL:
1358                    inc = sign == 1 && n != 0;  // round ceil
1359                    break;
1360    
1361                case ROUND_FLOOR:
1362                default:
1363                    inc = sign == -1 && n != 0;  // round floor
1364                    break;
1365            }
1366    
1367            if (inc) {
1368                // increment if necessary
1369                int rh = 1;
1370                for (int i = 0; i < mant.length; i++) {
1371                    final int r = mant[i] + rh;
1372                    rh = r / RADIX;
1373                    mant[i] = r - rh * RADIX;
1374                }
1375    
1376                if (rh != 0) {
1377                    shiftRight();
1378                    mant[mant.length-1] = rh;
1379                }
1380            }
1381    
1382            // check for exceptional cases and raise signals if necessary
1383            if (exp < MIN_EXP) {
1384                // Gradual Underflow
1385                field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
1386                return DfpField.FLAG_UNDERFLOW;
1387            }
1388    
1389            if (exp > MAX_EXP) {
1390                // Overflow
1391                field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
1392                return DfpField.FLAG_OVERFLOW;
1393            }
1394    
1395            if (n != 0) {
1396                // Inexact
1397                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1398                return DfpField.FLAG_INEXACT;
1399            }
1400    
1401            return 0;
1402    
1403        }
1404    
1405        /** Multiply this by x.
1406         * @param x multiplicand
1407         * @return product of this and x
1408         */
1409        public Dfp multiply(final Dfp x) {
1410    
1411            // make sure we don't mix number with different precision
1412            if (field.getRadixDigits() != x.field.getRadixDigits()) {
1413                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1414                final Dfp result = newInstance(getZero());
1415                result.nans = QNAN;
1416                return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1417            }
1418    
1419            Dfp result = newInstance(getZero());
1420    
1421            /* handle special cases */
1422            if (nans != FINITE || x.nans != FINITE) {
1423                if (isNaN()) {
1424                    return this;
1425                }
1426    
1427                if (x.isNaN()) {
1428                    return x;
1429                }
1430    
1431                if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] != 0) {
1432                    result = newInstance(this);
1433                    result.sign = (byte) (sign * x.sign);
1434                    return result;
1435                }
1436    
1437                if (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] != 0) {
1438                    result = newInstance(x);
1439                    result.sign = (byte) (sign * x.sign);
1440                    return result;
1441                }
1442    
1443                if (x.nans == INFINITE && nans == INFINITE) {
1444                    result = newInstance(this);
1445                    result.sign = (byte) (sign * x.sign);
1446                    return result;
1447                }
1448    
1449                if ( (x.nans == INFINITE && nans == FINITE && mant[mant.length-1] == 0) ||
1450                        (nans == INFINITE && x.nans == FINITE && x.mant[mant.length-1] == 0) ) {
1451                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1452                    result = newInstance(getZero());
1453                    result.nans = QNAN;
1454                    result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1455                    return result;
1456                }
1457            }
1458    
1459            int[] product = new int[mant.length*2];  // Big enough to hold even the largest result
1460    
1461            for (int i = 0; i < mant.length; i++) {
1462                int rh = 0;  // acts as a carry
1463                for (int j=0; j<mant.length; j++) {
1464                    int r = mant[i] * x.mant[j];    // multiply the 2 digits
1465                    r = r + product[i+j] + rh;  // add to the product digit with carry in
1466    
1467                    rh = r / RADIX;
1468                    product[i+j] = r - rh * RADIX;
1469                }
1470                product[i+mant.length] = rh;
1471            }
1472    
1473            // Find the most sig digit
1474            int md = mant.length * 2 - 1;  // default, in case result is zero
1475            for (int i = mant.length * 2 - 1; i >= 0; i--) {
1476                if (product[i] != 0) {
1477                    md = i;
1478                    break;
1479                }
1480            }
1481    
1482            // Copy the digits into the result
1483            for (int i = 0; i < mant.length; i++) {
1484                result.mant[mant.length - i - 1] = product[md - i];
1485            }
1486    
1487            // Fixup the exponent.
1488            result.exp = exp + x.exp + md - 2 * mant.length + 1;
1489            result.sign = (byte)((sign == x.sign)?1:-1);
1490    
1491            if (result.mant[mant.length-1] == 0) {
1492                // if result is zero, set exp to zero
1493                result.exp = 0;
1494            }
1495    
1496            final int excp;
1497            if (md > (mant.length-1)) {
1498                excp = result.round(product[md-mant.length]);
1499            } else {
1500                excp = result.round(0); // has no effect except to check status
1501            }
1502    
1503            if (excp != 0) {
1504                result = dotrap(excp, MULTIPLY_TRAP, x, result);
1505            }
1506    
1507            return result;
1508    
1509        }
1510    
1511        /** Multiply this by a single digit 0&lt;=x&lt;radix.
1512         * There are speed advantages in this special case
1513         * @param x multiplicand
1514         * @return product of this and x
1515         */
1516        public Dfp multiply(final int x) {
1517            Dfp result = newInstance(this);
1518    
1519            /* handle special cases */
1520            if (nans != FINITE) {
1521                if (isNaN()) {
1522                    return this;
1523                }
1524    
1525                if (nans == INFINITE && x != 0) {
1526                    result = newInstance(this);
1527                    return result;
1528                }
1529    
1530                if (nans == INFINITE && x == 0) {
1531                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1532                    result = newInstance(getZero());
1533                    result.nans = QNAN;
1534                    result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
1535                    return result;
1536                }
1537            }
1538    
1539            /* range check x */
1540            if (x < 0 || x >= RADIX) {
1541                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1542                result = newInstance(getZero());
1543                result.nans = QNAN;
1544                result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
1545                return result;
1546            }
1547    
1548            int rh = 0;
1549            for (int i = 0; i < mant.length; i++) {
1550                final int r = mant[i] * x + rh;
1551                rh = r / RADIX;
1552                result.mant[i] = r - rh * RADIX;
1553            }
1554    
1555            int lostdigit = 0;
1556            if (rh != 0) {
1557                lostdigit = result.mant[0];
1558                result.shiftRight();
1559                result.mant[mant.length-1] = rh;
1560            }
1561    
1562            if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
1563                result.exp = 0;
1564            }
1565    
1566            final int excp = result.round(lostdigit);
1567            if (excp != 0) {
1568                result = dotrap(excp, MULTIPLY_TRAP, result, result);
1569            }
1570    
1571            return result;
1572        }
1573    
1574        /** Divide this by divisor.
1575         * @param divisor divisor
1576         * @return quotient of this by divisor
1577         */
1578        public Dfp divide(Dfp divisor) {
1579            int dividend[]; // current status of the dividend
1580            int quotient[]; // quotient
1581            int remainder[];// remainder
1582            int qd;         // current quotient digit we're working with
1583            int nsqd;       // number of significant quotient digits we have
1584            int trial=0;    // trial quotient digit
1585            int minadj;     // minimum adjustment
1586            boolean trialgood; // Flag to indicate a good trail digit
1587            int md=0;       // most sig digit in result
1588            int excp;       // exceptions
1589    
1590            // make sure we don't mix number with different precision
1591            if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
1592                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1593                final Dfp result = newInstance(getZero());
1594                result.nans = QNAN;
1595                return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1596            }
1597    
1598            Dfp result = newInstance(getZero());
1599    
1600            /* handle special cases */
1601            if (nans != FINITE || divisor.nans != FINITE) {
1602                if (isNaN()) {
1603                    return this;
1604                }
1605    
1606                if (divisor.isNaN()) {
1607                    return divisor;
1608                }
1609    
1610                if (nans == INFINITE && divisor.nans == FINITE) {
1611                    result = newInstance(this);
1612                    result.sign = (byte) (sign * divisor.sign);
1613                    return result;
1614                }
1615    
1616                if (divisor.nans == INFINITE && nans == FINITE) {
1617                    result = newInstance(getZero());
1618                    result.sign = (byte) (sign * divisor.sign);
1619                    return result;
1620                }
1621    
1622                if (divisor.nans == INFINITE && nans == INFINITE) {
1623                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1624                    result = newInstance(getZero());
1625                    result.nans = QNAN;
1626                    result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1627                    return result;
1628                }
1629            }
1630    
1631            /* Test for divide by zero */
1632            if (divisor.mant[mant.length-1] == 0) {
1633                field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1634                result = newInstance(getZero());
1635                result.sign = (byte) (sign * divisor.sign);
1636                result.nans = INFINITE;
1637                result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
1638                return result;
1639            }
1640    
1641            dividend = new int[mant.length+1];  // one extra digit needed
1642            quotient = new int[mant.length+2];  // two extra digits needed 1 for overflow, 1 for rounding
1643            remainder = new int[mant.length+1]; // one extra digit needed
1644    
1645            /* Initialize our most significant digits to zero */
1646    
1647            dividend[mant.length] = 0;
1648            quotient[mant.length] = 0;
1649            quotient[mant.length+1] = 0;
1650            remainder[mant.length] = 0;
1651    
1652            /* copy our mantissa into the dividend, initialize the
1653           quotient while we are at it */
1654    
1655            for (int i = 0; i < mant.length; i++) {
1656                dividend[i] = mant[i];
1657                quotient[i] = 0;
1658                remainder[i] = 0;
1659            }
1660    
1661            /* outer loop.  Once per quotient digit */
1662            nsqd = 0;
1663            for (qd = mant.length+1; qd >= 0; qd--) {
1664                /* Determine outer limits of our quotient digit */
1665    
1666                // r =  most sig 2 digits of dividend
1667                final int divMsb = dividend[mant.length]*RADIX+dividend[mant.length-1];
1668                int min = divMsb       / (divisor.mant[mant.length-1]+1);
1669                int max = (divMsb + 1) / divisor.mant[mant.length-1];
1670    
1671                trialgood = false;
1672                while (!trialgood) {
1673                    // try the mean
1674                    trial = (min+max)/2;
1675    
1676                    /* Multiply by divisor and store as remainder */
1677                    int rh = 0;
1678                    for (int i = 0; i < mant.length + 1; i++) {
1679                        int dm = (i<mant.length)?divisor.mant[i]:0;
1680                        final int r = (dm * trial) + rh;
1681                        rh = r / RADIX;
1682                        remainder[i] = r - rh * RADIX;
1683                    }
1684    
1685                    /* subtract the remainder from the dividend */
1686                    rh = 1;  // carry in to aid the subtraction
1687                    for (int i = 0; i < mant.length + 1; i++) {
1688                        final int r = ((RADIX-1) - remainder[i]) + dividend[i] + rh;
1689                        rh = r / RADIX;
1690                        remainder[i] = r - rh * RADIX;
1691                    }
1692    
1693                    /* Lets analyze what we have here */
1694                    if (rh == 0) {
1695                        // trial is too big -- negative remainder
1696                        max = trial-1;
1697                        continue;
1698                    }
1699    
1700                    /* find out how far off the remainder is telling us we are */
1701                    minadj = (remainder[mant.length] * RADIX)+remainder[mant.length-1];
1702                    minadj = minadj / (divisor.mant[mant.length-1]+1);
1703    
1704                    if (minadj >= 2) {
1705                        min = trial+minadj;  // update the minimum
1706                        continue;
1707                    }
1708    
1709                    /* May have a good one here, check more thoroughly.  Basically
1710               its a good one if it is less than the divisor */
1711                    trialgood = false;  // assume false
1712                    for (int i = mant.length - 1; i >= 0; i--) {
1713                        if (divisor.mant[i] > remainder[i]) {
1714                            trialgood = true;
1715                        }
1716                        if (divisor.mant[i] < remainder[i]) {
1717                            break;
1718                        }
1719                    }
1720    
1721                    if (remainder[mant.length] != 0) {
1722                        trialgood = false;
1723                    }
1724    
1725                    if (trialgood == false) {
1726                        min = trial+1;
1727                    }
1728                }
1729    
1730                /* Great we have a digit! */
1731                quotient[qd] = trial;
1732                if (trial != 0 || nsqd != 0) {
1733                    nsqd++;
1734                }
1735    
1736                if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
1737                    // We have enough for this mode
1738                    break;
1739                }
1740    
1741                if (nsqd > mant.length) {
1742                    // We have enough digits
1743                    break;
1744                }
1745    
1746                /* move the remainder into the dividend while left shifting */
1747                dividend[0] = 0;
1748                for (int i = 0; i < mant.length; i++) {
1749                    dividend[i + 1] = remainder[i];
1750                }
1751            }
1752    
1753            /* Find the most sig digit */
1754            md = mant.length;  // default
1755            for (int i = mant.length + 1; i >= 0; i--) {
1756                if (quotient[i] != 0) {
1757                    md = i;
1758                    break;
1759                }
1760            }
1761    
1762            /* Copy the digits into the result */
1763            for (int i=0; i<mant.length; i++) {
1764                result.mant[mant.length-i-1] = quotient[md-i];
1765            }
1766    
1767            /* Fixup the exponent. */
1768            result.exp = exp - divisor.exp + md - mant.length;
1769            result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
1770    
1771            if (result.mant[mant.length-1] == 0) { // if result is zero, set exp to zero
1772                result.exp = 0;
1773            }
1774    
1775            if (md > (mant.length-1)) {
1776                excp = result.round(quotient[md-mant.length]);
1777            } else {
1778                excp = result.round(0);
1779            }
1780    
1781            if (excp != 0) {
1782                result = dotrap(excp, DIVIDE_TRAP, divisor, result);
1783            }
1784    
1785            return result;
1786        }
1787    
1788        /** Divide by a single digit less than radix.
1789         *  Special case, so there are speed advantages. 0 &lt;= divisor &lt; radix
1790         * @param divisor divisor
1791         * @return quotient of this by divisor
1792         */
1793        public Dfp divide(int divisor) {
1794    
1795            // Handle special cases
1796            if (nans != FINITE) {
1797                if (isNaN()) {
1798                    return this;
1799                }
1800    
1801                if (nans == INFINITE) {
1802                    return newInstance(this);
1803                }
1804            }
1805    
1806            // Test for divide by zero
1807            if (divisor == 0) {
1808                field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1809                Dfp result = newInstance(getZero());
1810                result.sign = sign;
1811                result.nans = INFINITE;
1812                result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
1813                return result;
1814            }
1815    
1816            // range check divisor
1817            if (divisor < 0 || divisor >= RADIX) {
1818                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1819                Dfp result = newInstance(getZero());
1820                result.nans = QNAN;
1821                result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
1822                return result;
1823            }
1824    
1825            Dfp result = newInstance(this);
1826    
1827            int rl = 0;
1828            for (int i = mant.length-1; i >= 0; i--) {
1829                final int r = rl*RADIX + result.mant[i];
1830                final int rh = r / divisor;
1831                rl = r - rh * divisor;
1832                result.mant[i] = rh;
1833            }
1834    
1835            if (result.mant[mant.length-1] == 0) {
1836                // normalize
1837                result.shiftLeft();
1838                final int r = rl * RADIX;        // compute the next digit and put it in
1839                final int rh = r / divisor;
1840                rl = r - rh * divisor;
1841                result.mant[0] = rh;
1842            }
1843    
1844            final int excp = result.round(rl * RADIX / divisor);  // do the rounding
1845            if (excp != 0) {
1846                result = dotrap(excp, DIVIDE_TRAP, result, result);
1847            }
1848    
1849            return result;
1850    
1851        }
1852    
1853        /** Compute the square root.
1854         * @return square root of the instance
1855         */
1856        public Dfp sqrt() {
1857    
1858            // check for unusual cases
1859            if (nans == FINITE && mant[mant.length-1] == 0) {
1860                // if zero
1861                return newInstance(this);
1862            }
1863    
1864            if (nans != FINITE) {
1865                if (nans == INFINITE && sign == 1) {
1866                    // if positive infinity
1867                    return newInstance(this);
1868                }
1869    
1870                if (nans == QNAN) {
1871                    return newInstance(this);
1872                }
1873    
1874                if (nans == SNAN) {
1875                    Dfp result;
1876    
1877                    field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1878                    result = newInstance(this);
1879                    result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
1880                    return result;
1881                }
1882            }
1883    
1884            if (sign == -1) {
1885                // if negative
1886                Dfp result;
1887    
1888                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1889                result = newInstance(this);
1890                result.nans = QNAN;
1891                result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
1892                return result;
1893            }
1894    
1895            Dfp x = newInstance(this);
1896    
1897            /* Lets make a reasonable guess as to the size of the square root */
1898            if (x.exp < -1 || x.exp > 1) {
1899                x.exp = this.exp / 2;
1900            }
1901    
1902            /* Coarsely estimate the mantissa */
1903            switch (x.mant[mant.length-1] / 2000) {
1904                case 0:
1905                    x.mant[mant.length-1] = x.mant[mant.length-1]/2+1;
1906                    break;
1907                case 2:
1908                    x.mant[mant.length-1] = 1500;
1909                    break;
1910                case 3:
1911                    x.mant[mant.length-1] = 2200;
1912                    break;
1913                default:
1914                    x.mant[mant.length-1] = 3000;
1915            }
1916    
1917            Dfp dx = newInstance(x);
1918    
1919            /* Now that we have the first pass estimate, compute the rest
1920           by the formula dx = (y - x*x) / (2x); */
1921    
1922            Dfp px  = getZero();
1923            Dfp ppx = getZero();
1924            while (x.unequal(px)) {
1925                dx = newInstance(x);
1926                dx.sign = -1;
1927                dx = dx.add(this.divide(x));
1928                dx = dx.divide(2);
1929                ppx = px;
1930                px = x;
1931                x = x.add(dx);
1932    
1933                if (x.equals(ppx)) {
1934                    // alternating between two values
1935                    break;
1936                }
1937    
1938                // if dx is zero, break.  Note testing the most sig digit
1939                // is a sufficient test since dx is normalized
1940                if (dx.mant[mant.length-1] == 0) {
1941                    break;
1942                }
1943            }
1944    
1945            return x;
1946    
1947        }
1948    
1949        /** Get a string representation of the instance.
1950         * @return string representation of the instance
1951         */
1952        @Override
1953        public String toString() {
1954            if (nans != FINITE) {
1955                // if non-finite exceptional cases
1956                if (nans == INFINITE) {
1957                    return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
1958                } else {
1959                    return NAN_STRING;
1960                }
1961            }
1962    
1963            if (exp > mant.length || exp < -1) {
1964                return dfp2sci();
1965            }
1966    
1967            return dfp2string();
1968    
1969        }
1970    
1971        /** Convert an instance to a string using scientific notation.
1972         * @return string representation of the instance in scientific notation
1973         */
1974        protected String dfp2sci() {
1975            char rawdigits[]    = new char[mant.length * 4];
1976            char outputbuffer[] = new char[mant.length * 4 + 20];
1977            int p;
1978            int q;
1979            int e;
1980            int ae;
1981            int shf;
1982    
1983            // Get all the digits
1984            p = 0;
1985            for (int i = mant.length - 1; i >= 0; i--) {
1986                rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
1987                rawdigits[p++] = (char) (((mant[i] / 100) %10) + '0');
1988                rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
1989                rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
1990            }
1991    
1992            // Find the first non-zero one
1993            for (p = 0; p < rawdigits.length; p++) {
1994                if (rawdigits[p] != '0') {
1995                    break;
1996                }
1997            }
1998            shf = p;
1999    
2000            // Now do the conversion
2001            q = 0;
2002            if (sign == -1) {
2003                outputbuffer[q++] = '-';
2004            }
2005    
2006            if (p != rawdigits.length) {
2007                // there are non zero digits...
2008                outputbuffer[q++] = rawdigits[p++];
2009                outputbuffer[q++] = '.';
2010    
2011                while (p<rawdigits.length) {
2012                    outputbuffer[q++] = rawdigits[p++];
2013                }
2014            } else {
2015                outputbuffer[q++] = '0';
2016                outputbuffer[q++] = '.';
2017                outputbuffer[q++] = '0';
2018                outputbuffer[q++] = 'e';
2019                outputbuffer[q++] = '0';
2020                return new String(outputbuffer, 0, 5);
2021            }
2022    
2023            outputbuffer[q++] = 'e';
2024    
2025            // Find the msd of the exponent
2026    
2027            e = exp * 4 - shf - 1;
2028            ae = e;
2029            if (e < 0) {
2030                ae = -e;
2031            }
2032    
2033            // Find the largest p such that p < e
2034            for (p = 1000000000; p > ae; p /= 10) {
2035                // nothing to do
2036            }
2037    
2038            if (e < 0) {
2039                outputbuffer[q++] = '-';
2040            }
2041    
2042            while (p > 0) {
2043                outputbuffer[q++] = (char)(ae / p + '0');
2044                ae = ae % p;
2045                p = p / 10;
2046            }
2047    
2048            return new String(outputbuffer, 0, q);
2049    
2050        }
2051    
2052        /** Convert an instance to a string using normal notation.
2053         * @return string representation of the instance in normal notation
2054         */
2055        protected String dfp2string() {
2056            char buffer[] = new char[mant.length*4 + 20];
2057            int p = 1;
2058            int q;
2059            int e = exp;
2060            boolean pointInserted = false;
2061    
2062            buffer[0] = ' ';
2063    
2064            if (e <= 0) {
2065                buffer[p++] = '0';
2066                buffer[p++] = '.';
2067                pointInserted = true;
2068            }
2069    
2070            while (e < 0) {
2071                buffer[p++] = '0';
2072                buffer[p++] = '0';
2073                buffer[p++] = '0';
2074                buffer[p++] = '0';
2075                e++;
2076            }
2077    
2078            for (int i = mant.length - 1; i >= 0; i--) {
2079                buffer[p++] = (char) ((mant[i] / 1000) + '0');
2080                buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
2081                buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
2082                buffer[p++] = (char) (((mant[i]) % 10) + '0');
2083                if (--e == 0) {
2084                    buffer[p++] = '.';
2085                    pointInserted = true;
2086                }
2087            }
2088    
2089            while (e > 0) {
2090                buffer[p++] = '0';
2091                buffer[p++] = '0';
2092                buffer[p++] = '0';
2093                buffer[p++] = '0';
2094                e--;
2095            }
2096    
2097            if (!pointInserted) {
2098                // Ensure we have a radix point!
2099                buffer[p++] = '.';
2100            }
2101    
2102            // Suppress leading zeros
2103            q = 1;
2104            while (buffer[q] == '0') {
2105                q++;
2106            }
2107            if (buffer[q] == '.') {
2108                q--;
2109            }
2110    
2111            // Suppress trailing zeros
2112            while (buffer[p-1] == '0') {
2113                p--;
2114            }
2115    
2116            // Insert sign
2117            if (sign < 0) {
2118                buffer[--q] = '-';
2119            }
2120    
2121            return new String(buffer, q, p - q);
2122    
2123        }
2124    
2125        /** Raises a trap.  This does not set the corresponding flag however.
2126         *  @param type the trap type
2127         *  @param what - name of routine trap occurred in
2128         *  @param oper - input operator to function
2129         *  @param result - the result computed prior to the trap
2130         *  @return The suggested return value from the trap handler
2131         */
2132        public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
2133            Dfp def = result;
2134    
2135            switch (type) {
2136                case DfpField.FLAG_INVALID:
2137                    def = newInstance(getZero());
2138                    def.sign = result.sign;
2139                    def.nans = QNAN;
2140                    break;
2141    
2142                case DfpField.FLAG_DIV_ZERO:
2143                    if (nans == FINITE && mant[mant.length-1] != 0) {
2144                        // normal case, we are finite, non-zero
2145                        def = newInstance(getZero());
2146                        def.sign = (byte)(sign*oper.sign);
2147                        def.nans = INFINITE;
2148                    }
2149    
2150                    if (nans == FINITE && mant[mant.length-1] == 0) {
2151                        //  0/0
2152                        def = newInstance(getZero());
2153                        def.nans = QNAN;
2154                    }
2155    
2156                    if (nans == INFINITE || nans == QNAN) {
2157                        def = newInstance(getZero());
2158                        def.nans = QNAN;
2159                    }
2160    
2161                    if (nans == INFINITE || nans == SNAN) {
2162                        def = newInstance(getZero());
2163                        def.nans = QNAN;
2164                    }
2165                    break;
2166    
2167                case DfpField.FLAG_UNDERFLOW:
2168                    if ( (result.exp+mant.length) < MIN_EXP) {
2169                        def = newInstance(getZero());
2170                        def.sign = result.sign;
2171                    } else {
2172                        def = newInstance(result);  // gradual underflow
2173                    }
2174                    result.exp = result.exp + ERR_SCALE;
2175                    break;
2176    
2177                case DfpField.FLAG_OVERFLOW:
2178                    result.exp = result.exp - ERR_SCALE;
2179                    def = newInstance(getZero());
2180                    def.sign = result.sign;
2181                    def.nans = INFINITE;
2182                    break;
2183    
2184                default: def = result; break;
2185            }
2186    
2187            return trap(type, what, oper, def, result);
2188    
2189        }
2190    
2191        /** Trap handler.  Subclasses may override this to provide trap
2192         *  functionality per IEEE 854-1987.
2193         *
2194         *  @param type  The exception type - e.g. FLAG_OVERFLOW
2195         *  @param what  The name of the routine we were in e.g. divide()
2196         *  @param oper  An operand to this function if any
2197         *  @param def   The default return value if trap not enabled
2198         *  @param result    The result that is specified to be delivered per
2199         *                   IEEE 854, if any
2200         *  @return the value that should be return by the operation triggering the trap
2201         */
2202        protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
2203            return def;
2204        }
2205    
2206        /** Returns the type - one of FINITE, INFINITE, SNAN, QNAN.
2207         * @return type of the number
2208         */
2209        public int classify() {
2210            return nans;
2211        }
2212    
2213        /** Creates an instance that is the same as x except that it has the sign of y.
2214         * abs(x) = dfp.copysign(x, dfp.one)
2215         * @param x number to get the value from
2216         * @param y number to get the sign from
2217         * @return a number with the value of x and the sign of y
2218         */
2219        public static Dfp copysign(final Dfp x, final Dfp y) {
2220            Dfp result = x.newInstance(x);
2221            result.sign = y.sign;
2222            return result;
2223        }
2224    
2225        /** Returns the next number greater than this one in the direction of x.
2226         * If this==x then simply returns this.
2227         * @param x direction where to look at
2228         * @return closest number next to instance in the direction of x
2229         */
2230        public Dfp nextAfter(final Dfp x) {
2231    
2232            // make sure we don't mix number with different precision
2233            if (field.getRadixDigits() != x.field.getRadixDigits()) {
2234                field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
2235                final Dfp result = newInstance(getZero());
2236                result.nans = QNAN;
2237                return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
2238            }
2239    
2240            // if this is greater than x
2241            boolean up = false;
2242            if (this.lessThan(x)) {
2243                up = true;
2244            }
2245    
2246            if (compare(this, x) == 0) {
2247                return newInstance(x);
2248            }
2249    
2250            if (lessThan(getZero())) {
2251                up = !up;
2252            }
2253    
2254            final Dfp inc;
2255            Dfp result;
2256            if (up) {
2257                inc = newInstance(getOne());
2258                inc.exp = this.exp-mant.length+1;
2259                inc.sign = this.sign;
2260    
2261                if (this.equals(getZero())) {
2262                    inc.exp = MIN_EXP-mant.length;
2263                }
2264    
2265                result = add(inc);
2266            } else {
2267                inc = newInstance(getOne());
2268                inc.exp = this.exp;
2269                inc.sign = this.sign;
2270    
2271                if (this.equals(inc)) {
2272                    inc.exp = this.exp-mant.length;
2273                } else {
2274                    inc.exp = this.exp-mant.length+1;
2275                }
2276    
2277                if (this.equals(getZero())) {
2278                    inc.exp = MIN_EXP-mant.length;
2279                }
2280    
2281                result = this.subtract(inc);
2282            }
2283    
2284            if (result.classify() == INFINITE && this.classify() != INFINITE) {
2285                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2286                result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2287            }
2288    
2289            if (result.equals(getZero()) && this.equals(getZero()) == false) {
2290                field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2291                result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2292            }
2293    
2294            return result;
2295    
2296        }
2297    
2298        /** Convert the instance into a double.
2299         * @return a double approximating the instance
2300         * @see #toSplitDouble()
2301         */
2302        public double toDouble() {
2303    
2304            if (isInfinite()) {
2305                if (lessThan(getZero())) {
2306                    return Double.NEGATIVE_INFINITY;
2307                } else {
2308                    return Double.POSITIVE_INFINITY;
2309                }
2310            }
2311    
2312            if (isNaN()) {
2313                return Double.NaN;
2314            }
2315    
2316            Dfp y = this;
2317            boolean negate = false;
2318            if (lessThan(getZero())) {
2319                y = negate();
2320                negate = true;
2321            }
2322    
2323            /* Find the exponent, first estimate by integer log10, then adjust.
2324             Should be faster than doing a natural logarithm.  */
2325            int exponent = (int)(y.log10() * 3.32);
2326            if (exponent < 0) {
2327                exponent--;
2328            }
2329    
2330            Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
2331            while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
2332                tempDfp = tempDfp.multiply(2);
2333                exponent++;
2334            }
2335            exponent--;
2336    
2337            /* We have the exponent, now work on the mantissa */
2338    
2339            y = y.divide(DfpMath.pow(getTwo(), exponent));
2340            if (exponent > -1023) {
2341                y = y.subtract(getOne());
2342            }
2343    
2344            if (exponent < -1074) {
2345                return 0;
2346            }
2347    
2348            if (exponent > 1023) {
2349                return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
2350            }
2351    
2352    
2353            y = y.multiply(newInstance(4503599627370496l)).rint();
2354            String str = y.toString();
2355            str = str.substring(0, str.length()-1);
2356            long mantissa = Long.parseLong(str);
2357    
2358            if (mantissa == 4503599627370496L) {
2359                // Handle special case where we round up to next power of two
2360                mantissa = 0;
2361                exponent++;
2362            }
2363    
2364            /* Its going to be subnormal, so make adjustments */
2365            if (exponent <= -1023) {
2366                exponent--;
2367            }
2368    
2369            while (exponent < -1023) {
2370                exponent++;
2371                mantissa >>>= 1;
2372            }
2373    
2374            long bits = mantissa | ((exponent + 1023L) << 52);
2375            double x = Double.longBitsToDouble(bits);
2376    
2377            if (negate) {
2378                x = -x;
2379            }
2380    
2381            return x;
2382    
2383        }
2384    
2385        /** Convert the instance into a split double.
2386         * @return an array of two doubles which sum represent the instance
2387         * @see #toDouble()
2388         */
2389        public double[] toSplitDouble() {
2390            double split[] = new double[2];
2391            long mask = 0xffffffffc0000000L;
2392    
2393            split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
2394            split[1] = subtract(newInstance(split[0])).toDouble();
2395    
2396            return split;
2397        }
2398    
2399    }