1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math4.legacy.core.dfp;
19
20 import java.util.Arrays;
21
22 import org.apache.commons.math4.legacy.core.RealFieldElement;
23 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96 public class Dfp implements RealFieldElement<Dfp> {
97
98
99 public static final int RADIX = 10000;
100
101
102
103 public static final int MIN_EXP = -32767;
104
105
106
107 public static final int MAX_EXP = 32768;
108
109
110 public static final int ERR_SCALE = 32760;
111
112
113 public static final byte FINITE = 0;
114
115
116 public static final byte INFINITE = 1;
117
118
119 public static final byte SNAN = 2;
120
121
122 public static final byte QNAN = 3;
123
124
125 private static final String NAN_STRING = "NaN";
126
127
128 private static final String POS_INFINITY_STRING = "Infinity";
129
130
131 private static final String NEG_INFINITY_STRING = "-Infinity";
132
133
134 private static final String ADD_TRAP = "add";
135
136
137 private static final String MULTIPLY_TRAP = "multiply";
138
139
140 private static final String DIVIDE_TRAP = "divide";
141
142
143 private static final String SQRT_TRAP = "sqrt";
144
145
146 private static final String ALIGN_TRAP = "align";
147
148
149 private static final String TRUNC_TRAP = "trunc";
150
151
152 private static final String NEXT_AFTER_TRAP = "nextAfter";
153
154
155 private static final String LESS_THAN_TRAP = "lessThan";
156
157
158 private static final String GREATER_THAN_TRAP = "greaterThan";
159
160
161 private static final String NEW_INSTANCE_TRAP = "newInstance";
162
163
164 protected int[] mant;
165
166
167 protected byte sign;
168
169
170 protected int exp;
171
172
173 protected byte nans;
174
175
176 private final DfpField field;
177
178
179
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
190
191
192
193 protected Dfp(final DfpField field, byte x) {
194 this(field, (long) x);
195 }
196
197
198
199
200
201 protected Dfp(final DfpField field, int x) {
202 this(field, (long) x);
203 }
204
205
206
207
208
209 protected Dfp(final DfpField field, long x) {
210
211
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
219
220 isLongMin = true;
221 ++x;
222 }
223
224
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
242
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
253
254
255
256 protected Dfp(final DfpField field, double x) {
257
258
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
271 if (x == 0) {
272
273 if ((bits & 0x8000000000000000L) != 0) {
274 sign = -1;
275 }
276 return;
277 }
278
279 exponent++;
280
281
282 while ((mantissa & 0x0010000000000000L) == 0) {
283 exponent--;
284 mantissa <<= 1;
285 }
286 mantissa &= 0x000fffffffffffffL;
287 }
288
289 if (exponent == 1024) {
290
291 if (Double.isNaN(x)) {
292 sign = (byte) 1;
293 nans = QNAN;
294 } else if (x < 0) {
295 sign = (byte) -1;
296 nans = INFINITE;
297 } else {
298 sign = (byte) 1;
299 nans = INFINITE;
300 }
301 return;
302 }
303
304 Dfp xdfp = new Dfp(field, mantissa);
305 xdfp = xdfp.divide(new Dfp(field, 4503599627370496L)).add(field.getOne());
306 xdfp = xdfp.multiply(DfpMath.pow(field.getTwo(), exponent));
307
308 if ((bits & 0x8000000000000000L) != 0) {
309 xdfp = xdfp.negate();
310 }
311
312 System.arraycopy(xdfp.mant, 0, mant, 0, mant.length);
313 sign = xdfp.sign;
314 exp = xdfp.exp;
315 nans = xdfp.nans;
316 }
317
318
319
320
321 public Dfp(final Dfp d) {
322 mant = d.mant.clone();
323 sign = d.sign;
324 exp = d.exp;
325 nans = d.nans;
326 field = d.field;
327 }
328
329
330
331
332
333 protected Dfp(final DfpField field, final String s) {
334
335
336 mant = new int[field.getRadixDigits()];
337 sign = 1;
338 exp = 0;
339 nans = FINITE;
340 this.field = field;
341
342 boolean decimalFound = false;
343 final int rsize = 4;
344 final int offset = 4;
345 final char[] striped = new char[getRadixDigits() * rsize + offset * 2];
346
347
348 if (s.equals(POS_INFINITY_STRING)) {
349 sign = (byte) 1;
350 nans = INFINITE;
351 return;
352 }
353
354 if (s.equals(NEG_INFINITY_STRING)) {
355 sign = (byte) -1;
356 nans = INFINITE;
357 return;
358 }
359
360 if (s.equals(NAN_STRING)) {
361 sign = (byte) 1;
362 nans = QNAN;
363 return;
364 }
365
366
367 int p = s.indexOf("e");
368 if (p == -1) {
369 p = s.indexOf("E");
370 }
371
372 final String fpdecimal;
373 int sciexp = 0;
374 if (p != -1) {
375
376 fpdecimal = s.substring(0, p);
377 String fpexp = s.substring(p + 1);
378 boolean negative = false;
379
380 for (int i = 0; i < fpexp.length(); i++) {
381 if (fpexp.charAt(i) == '-') {
382 negative = true;
383 continue;
384 }
385 if (fpexp.charAt(i) >= '0' && fpexp.charAt(i) <= '9') {
386 sciexp = sciexp * 10 + fpexp.charAt(i) - '0';
387 }
388 }
389
390 if (negative) {
391 sciexp = -sciexp;
392 }
393 } else {
394
395 fpdecimal = s;
396 }
397
398
399 if (fpdecimal.indexOf("-") != -1) {
400 sign = -1;
401 }
402
403
404 p = 0;
405
406
407 int decimalPos = 0;
408 for (;;) {
409 if (fpdecimal.charAt(p) >= '1' && fpdecimal.charAt(p) <= '9') {
410 break;
411 }
412
413 if (decimalFound && fpdecimal.charAt(p) == '0') {
414 decimalPos--;
415 }
416
417 if (fpdecimal.charAt(p) == '.') {
418 decimalFound = true;
419 }
420
421 p++;
422
423 if (p == fpdecimal.length()) {
424 break;
425 }
426 }
427
428
429 int q = offset;
430 striped[0] = '0';
431 striped[1] = '0';
432 striped[2] = '0';
433 striped[3] = '0';
434 int significantDigits = 0;
435 for (;;) {
436 if (p == (fpdecimal.length())) {
437 break;
438 }
439
440
441 if (q == mant.length * rsize + offset + 1) {
442 break;
443 }
444
445 if (fpdecimal.charAt(p) == '.') {
446 decimalFound = true;
447 decimalPos = significantDigits;
448 p++;
449 continue;
450 }
451
452 if (fpdecimal.charAt(p) < '0' || fpdecimal.charAt(p) > '9') {
453 p++;
454 continue;
455 }
456
457 striped[q] = fpdecimal.charAt(p);
458 q++;
459 p++;
460 significantDigits++;
461 }
462
463
464
465 if (decimalFound && q != offset) {
466 for (;;) {
467 q--;
468 if (q == offset) {
469 break;
470 }
471 if (striped[q] == '0') {
472 significantDigits--;
473 } else {
474 break;
475 }
476 }
477 }
478
479
480 if (decimalFound && significantDigits == 0) {
481 decimalPos = 0;
482 }
483
484
485 if (!decimalFound) {
486 decimalPos = q - offset;
487 }
488
489
490 q = offset;
491 p = significantDigits - 1 + offset;
492
493 while (p > q) {
494 if (striped[p] != '0') {
495 break;
496 }
497 p--;
498 }
499
500
501 int i = ((rsize * 100) - decimalPos - sciexp % rsize) % rsize;
502 q -= i;
503 decimalPos += i;
504
505
506 while ((p - q) < (mant.length * rsize)) {
507 for (i = 0; i < rsize; i++) {
508 striped[++p] = '0';
509 }
510 }
511
512
513
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
527 round((striped[q] - '0') * 1000);
528 }
529 }
530
531
532
533
534
535
536
537 protected Dfp(final DfpField field, final byte sign, final byte nans) {
538 this.field = field;
539 this.mant = new int[field.getRadixDigits()];
540 this.sign = sign;
541 this.exp = 0;
542 this.nans = nans;
543 }
544
545
546
547
548
549 public Dfp newInstance() {
550 return new Dfp(getField());
551 }
552
553
554
555
556
557 public Dfp newInstance(final byte x) {
558 return new Dfp(getField(), x);
559 }
560
561
562
563
564
565 public Dfp newInstance(final int x) {
566 return new Dfp(getField(), x);
567 }
568
569
570
571
572
573 public Dfp newInstance(final long x) {
574 return new Dfp(getField(), x);
575 }
576
577
578
579
580
581 public Dfp newInstance(final double x) {
582 return new Dfp(getField(), x);
583 }
584
585
586
587
588
589
590 public Dfp newInstance(final Dfp d) {
591
592
593 if (field.getRadixDigits() != d.field.getRadixDigits()) {
594 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
595 final Dfp result = newInstance(getZero());
596 result.nans = QNAN;
597 return dotrap(DfpField.FLAG_INVALID, NEW_INSTANCE_TRAP, d, result);
598 }
599
600 return new Dfp(d);
601 }
602
603
604
605
606
607
608 public Dfp newInstance(final String s) {
609 return new Dfp(field, s);
610 }
611
612
613
614
615
616
617
618 public Dfp newInstance(final byte sig, final byte code) {
619 return field.newDfp(sig, code);
620 }
621
622
623
624
625
626
627
628
629
630 @Override
631 public DfpField getField() {
632 return field;
633 }
634
635
636
637
638 public int getRadixDigits() {
639 return field.getRadixDigits();
640 }
641
642
643
644
645 public Dfp getZero() {
646 return field.getZero();
647 }
648
649
650
651
652 public Dfp getOne() {
653 return field.getOne();
654 }
655
656
657
658
659 public Dfp getTwo() {
660 return field.getTwo();
661 }
662
663
664
665 protected void shiftLeft() {
666 if (mant.length - 1 > 0) {
667 System.arraycopy(mant, 0, mant, 1, mant.length - 1);
668 }
669 mant[0] = 0;
670 exp--;
671 }
672
673
674
675
676
677 protected void shiftRight() {
678 if (mant.length - 1 > 0) {
679 System.arraycopy(mant, 1, mant, 0, mant.length - 1);
680 }
681 mant[mant.length - 1] = 0;
682 exp++;
683 }
684
685
686
687
688
689
690
691
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
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
722
723
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
747
748
749 public boolean lessThan(final Dfp x) {
750
751
752 if (field.getRadixDigits() != x.field.getRadixDigits()) {
753 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
754 final Dfp result = newInstance(getZero());
755 result.nans = QNAN;
756 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, result);
757 return false;
758 }
759
760
761 if (isNaN() || x.isNaN()) {
762 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
763 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, x, newInstance(getZero()));
764 return false;
765 }
766
767 return compare(this, x) < 0;
768 }
769
770
771
772
773
774 public boolean greaterThan(final Dfp x) {
775
776
777 if (field.getRadixDigits() != x.field.getRadixDigits()) {
778 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
779 final Dfp result = newInstance(getZero());
780 result.nans = QNAN;
781 dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, result);
782 return false;
783 }
784
785
786 if (isNaN() || x.isNaN()) {
787 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
788 dotrap(DfpField.FLAG_INVALID, GREATER_THAN_TRAP, x, newInstance(getZero()));
789 return false;
790 }
791
792 return compare(this, x) > 0;
793 }
794
795
796
797
798 public boolean negativeOrNull() {
799
800 if (isNaN()) {
801 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
802 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
803 return false;
804 }
805
806 return (sign < 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
807 }
808
809
810
811
812 public boolean strictlyNegative() {
813
814 if (isNaN()) {
815 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
816 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
817 return false;
818 }
819
820 return (sign < 0) && ((mant[mant.length - 1] != 0) || isInfinite());
821 }
822
823
824
825
826 public boolean positiveOrNull() {
827
828 if (isNaN()) {
829 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
830 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
831 return false;
832 }
833
834 return (sign > 0) || ((mant[mant.length - 1] == 0) && !isInfinite());
835 }
836
837
838
839
840 public boolean strictlyPositive() {
841
842 if (isNaN()) {
843 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
844 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
845 return false;
846 }
847
848 return (sign > 0) && ((mant[mant.length - 1] != 0) || isInfinite());
849 }
850
851
852
853
854
855 @Override
856 public Dfp abs() {
857 Dfp result = newInstance(this);
858 result.sign = 1;
859 return result;
860 }
861
862
863
864
865 public boolean isInfinite() {
866 return nans == INFINITE;
867 }
868
869
870
871
872 public boolean isNaN() {
873 return (nans == QNAN) || (nans == SNAN);
874 }
875
876
877
878
879 public boolean isZero() {
880
881 if (isNaN()) {
882 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
883 dotrap(DfpField.FLAG_INVALID, LESS_THAN_TRAP, this, newInstance(getZero()));
884 return false;
885 }
886
887 return (mant[mant.length - 1] == 0) && !isInfinite();
888 }
889
890
891
892
893
894 @Override
895 public boolean equals(final Object other) {
896
897 if (other instanceof Dfp) {
898 final Dfp x = (Dfp) other;
899 if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
900 return false;
901 }
902
903 return compare(this, x) == 0;
904 }
905
906 return false;
907 }
908
909
910
911
912
913 @Override
914 public int hashCode() {
915 return 17 + (isZero() ? 0 : (sign << 8)) + (nans << 16) + exp + Arrays.hashCode(mant);
916 }
917
918
919
920
921
922 public boolean unequal(final Dfp x) {
923 if (isNaN() || x.isNaN() || field.getRadixDigits() != x.field.getRadixDigits()) {
924 return false;
925 }
926
927 return greaterThan(x) || lessThan(x);
928 }
929
930
931
932
933
934
935
936 private static int compare(final Dfp a, final Dfp b) {
937
938 if (a.mant[a.mant.length - 1] == 0 && b.mant[b.mant.length - 1] == 0 &&
939 a.nans == FINITE && b.nans == FINITE) {
940 return 0;
941 }
942
943 if (a.sign != b.sign) {
944 if (a.sign == -1) {
945 return -1;
946 } else {
947 return 1;
948 }
949 }
950
951
952 if (a.nans == INFINITE && b.nans == FINITE) {
953 return a.sign;
954 }
955
956 if (a.nans == FINITE && b.nans == INFINITE) {
957 return -b.sign;
958 }
959
960 if (a.nans == INFINITE && b.nans == INFINITE) {
961 return 0;
962 }
963
964
965 if (b.mant[b.mant.length - 1] != 0 && a.mant[b.mant.length - 1] != 0) {
966 if (a.exp < b.exp) {
967 return -a.sign;
968 }
969
970 if (a.exp > b.exp) {
971 return a.sign;
972 }
973 }
974
975
976 for (int i = a.mant.length - 1; i >= 0; i--) {
977 if (a.mant[i] > b.mant[i]) {
978 return a.sign;
979 }
980
981 if (a.mant[i] < b.mant[i]) {
982 return -a.sign;
983 }
984 }
985
986 return 0;
987 }
988
989
990
991
992
993
994
995 @Override
996 public Dfp rint() {
997 return trunc(DfpField.RoundingMode.ROUND_HALF_EVEN);
998 }
999
1000
1001
1002
1003
1004
1005 @Override
1006 public Dfp floor() {
1007 return trunc(DfpField.RoundingMode.ROUND_FLOOR);
1008 }
1009
1010
1011
1012
1013
1014
1015 @Override
1016 public Dfp ceil() {
1017 return trunc(DfpField.RoundingMode.ROUND_CEIL);
1018 }
1019
1020
1021
1022
1023
1024
1025 @Override
1026 public Dfp remainder(final Dfp d) {
1027
1028 final Dfp result = this.subtract(this.divide(d).rint().multiply(d));
1029
1030
1031 if (result.mant[mant.length - 1] == 0) {
1032 result.sign = sign;
1033 }
1034
1035 return result;
1036 }
1037
1038
1039
1040
1041
1042 protected Dfp trunc(final DfpField.RoundingMode rmode) {
1043 boolean changed = false;
1044
1045 if (isNaN()) {
1046 return newInstance(this);
1047 }
1048
1049 if (nans == INFINITE) {
1050 return newInstance(this);
1051 }
1052
1053 if (mant[mant.length - 1] == 0) {
1054
1055 return newInstance(this);
1056 }
1057
1058
1059
1060 if (exp < 0) {
1061 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1062 Dfp result = newInstance(getZero());
1063 result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
1064 return result;
1065 }
1066
1067
1068
1069
1070
1071 if (exp >= mant.length) {
1072 return newInstance(this);
1073 }
1074
1075
1076
1077
1078 Dfp result = newInstance(this);
1079 for (int i = 0; i < mant.length - result.exp; i++) {
1080 changed |= result.mant[i] != 0;
1081 result.mant[i] = 0;
1082 }
1083
1084 if (changed) {
1085 switch (rmode) {
1086 case ROUND_FLOOR:
1087 if (result.sign == -1) {
1088
1089 result = result.add(newInstance(-1));
1090 }
1091 break;
1092
1093 case ROUND_CEIL:
1094 if (result.sign == 1) {
1095
1096 result = result.add(getOne());
1097 }
1098 break;
1099
1100 case ROUND_HALF_EVEN:
1101 default:
1102 final Dfp half = newInstance("0.5");
1103 Dfp a = subtract(result);
1104 a.sign = 1;
1105 if (a.greaterThan(half)) {
1106 a = newInstance(getOne());
1107 a.sign = sign;
1108 result = result.add(a);
1109 }
1110
1111
1112 if (a.equals(half) && result.exp > 0 && (result.mant[mant.length - result.exp] & 1) != 0) {
1113 a = newInstance(getOne());
1114 a.sign = sign;
1115 result = result.add(a);
1116 }
1117 break;
1118 }
1119
1120 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1121 result = dotrap(DfpField.FLAG_INEXACT, TRUNC_TRAP, this, result);
1122 return result;
1123 }
1124
1125 return result;
1126 }
1127
1128
1129
1130
1131
1132 public int intValue() {
1133 Dfp rounded;
1134 int result = 0;
1135
1136 rounded = rint();
1137
1138 if (rounded.greaterThan(newInstance(2147483647))) {
1139 return 2147483647;
1140 }
1141
1142 if (rounded.lessThan(newInstance(-2147483648))) {
1143 return -2147483648;
1144 }
1145
1146 for (int i = mant.length - 1; i >= mant.length - rounded.exp; i--) {
1147 result = result * RADIX + rounded.mant[i];
1148 }
1149
1150 if (rounded.sign == -1) {
1151 result = -result;
1152 }
1153
1154 return result;
1155 }
1156
1157
1158
1159
1160
1161
1162 public int log10K() {
1163 return exp - 1;
1164 }
1165
1166
1167
1168
1169
1170 public Dfp power10K(final int e) {
1171 Dfp d = newInstance(getOne());
1172 d.exp = e + 1;
1173 return d;
1174 }
1175
1176
1177
1178
1179
1180 public int intLog10() {
1181 if (mant[mant.length - 1] > 1000) {
1182 return exp * 4 - 1;
1183 }
1184 if (mant[mant.length - 1] > 100) {
1185 return exp * 4 - 2;
1186 }
1187 if (mant[mant.length - 1] > 10) {
1188 return exp * 4 - 3;
1189 }
1190 return exp * 4 - 4;
1191 }
1192
1193
1194
1195
1196
1197 public Dfp power10(final int e) {
1198 Dfp d = newInstance(getOne());
1199
1200 if (e >= 0) {
1201 d.exp = e / 4 + 1;
1202 } else {
1203 d.exp = (e + 1) / 4;
1204 }
1205
1206 switch ((e % 4 + 4) % 4) {
1207 case 0:
1208 break;
1209 case 1:
1210 d = d.multiply(10);
1211 break;
1212 case 2:
1213 d = d.multiply(100);
1214 break;
1215 default:
1216 d = d.multiply(1000);
1217 }
1218
1219 return d;
1220 }
1221
1222
1223
1224
1225
1226
1227
1228 protected int complement(int extra) {
1229
1230 extra = RADIX - extra;
1231 for (int i = 0; i < mant.length; i++) {
1232 mant[i] = RADIX - mant[i] - 1;
1233 }
1234
1235 int rh = extra / RADIX;
1236 extra -= rh * RADIX;
1237 for (int i = 0; i < mant.length; i++) {
1238 final int r = mant[i] + rh;
1239 rh = r / RADIX;
1240 mant[i] = r - rh * RADIX;
1241 }
1242
1243 return extra;
1244 }
1245
1246
1247
1248
1249
1250 @Override
1251 public Dfp add(final Dfp x) {
1252
1253
1254 if (field.getRadixDigits() != x.field.getRadixDigits()) {
1255 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1256 final Dfp result = newInstance(getZero());
1257 result.nans = QNAN;
1258 return dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1259 }
1260
1261
1262 if (nans != FINITE || x.nans != FINITE) {
1263 if (isNaN()) {
1264 return this;
1265 }
1266
1267 if (x.isNaN()) {
1268 return x;
1269 }
1270
1271 if (nans == INFINITE && x.nans == FINITE) {
1272 return this;
1273 }
1274
1275 if (x.nans == INFINITE && nans == FINITE) {
1276 return x;
1277 }
1278
1279 if (x.nans == INFINITE && nans == INFINITE && sign == x.sign) {
1280 return x;
1281 }
1282
1283 if (x.nans == INFINITE && nans == INFINITE && sign != x.sign) {
1284 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1285 Dfp result = newInstance(getZero());
1286 result.nans = QNAN;
1287 result = dotrap(DfpField.FLAG_INVALID, ADD_TRAP, x, result);
1288 return result;
1289 }
1290 }
1291
1292
1293 Dfp a = newInstance(this);
1294 Dfp b = newInstance(x);
1295
1296
1297 Dfp result = newInstance(getZero());
1298
1299
1300 final byte asign = a.sign;
1301 final byte bsign = b.sign;
1302
1303 a.sign = 1;
1304 b.sign = 1;
1305
1306
1307 byte rsign = bsign;
1308 if (compare(a, b) > 0) {
1309 rsign = asign;
1310 }
1311
1312
1313
1314
1315
1316
1317 if (b.mant[mant.length - 1] == 0) {
1318 b.exp = a.exp;
1319 }
1320
1321 if (a.mant[mant.length - 1] == 0) {
1322 a.exp = b.exp;
1323 }
1324
1325
1326 int aextradigit = 0;
1327 int bextradigit = 0;
1328 if (a.exp < b.exp) {
1329 aextradigit = a.align(b.exp);
1330 } else {
1331 bextradigit = b.align(a.exp);
1332 }
1333
1334
1335 if (asign != bsign) {
1336 if (asign == rsign) {
1337 bextradigit = b.complement(bextradigit);
1338 } else {
1339 aextradigit = a.complement(aextradigit);
1340 }
1341 }
1342
1343
1344 int rh = 0;
1345 for (int i = 0; i < mant.length; i++) {
1346 final int r = a.mant[i] + b.mant[i] + rh;
1347 rh = r / RADIX;
1348 result.mant[i] = r - rh * RADIX;
1349 }
1350 result.exp = a.exp;
1351 result.sign = rsign;
1352
1353
1354
1355
1356 if (rh != 0 && (asign == bsign)) {
1357 final int lostdigit = result.mant[0];
1358 result.shiftRight();
1359 result.mant[mant.length - 1] = rh;
1360 final int excp = result.round(lostdigit);
1361 if (excp != 0) {
1362 result = dotrap(excp, ADD_TRAP, x, result);
1363 }
1364 }
1365
1366
1367 for (int i = 0; i < mant.length; i++) {
1368 if (result.mant[mant.length - 1] != 0) {
1369 break;
1370 }
1371 result.shiftLeft();
1372 if (i == 0) {
1373 result.mant[0] = aextradigit + bextradigit;
1374 aextradigit = 0;
1375 bextradigit = 0;
1376 }
1377 }
1378
1379
1380 if (result.mant[mant.length - 1] == 0) {
1381 result.exp = 0;
1382
1383 if (asign != bsign) {
1384
1385 result.sign = 1;
1386 }
1387 }
1388
1389
1390 final int excp = result.round(aextradigit + bextradigit);
1391 if (excp != 0) {
1392 result = dotrap(excp, ADD_TRAP, x, result);
1393 }
1394
1395 return result;
1396 }
1397
1398
1399
1400
1401 @Override
1402 public Dfp negate() {
1403 Dfp result = newInstance(this);
1404 result.sign = (byte) -result.sign;
1405 return result;
1406 }
1407
1408
1409
1410
1411
1412 @Override
1413 public Dfp subtract(final Dfp x) {
1414 return add(x.negate());
1415 }
1416
1417
1418
1419
1420
1421 protected int round(int n) {
1422 boolean inc = false;
1423 switch (field.getRoundingMode()) {
1424 case ROUND_DOWN:
1425 inc = false;
1426 break;
1427
1428 case ROUND_UP:
1429 inc = n != 0;
1430 break;
1431
1432 case ROUND_HALF_UP:
1433 inc = n >= 5000;
1434 break;
1435
1436 case ROUND_HALF_DOWN:
1437 inc = n > 5000;
1438 break;
1439
1440 case ROUND_HALF_EVEN:
1441 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 1);
1442 break;
1443
1444 case ROUND_HALF_ODD:
1445 inc = n > 5000 || (n == 5000 && (mant[0] & 1) == 0);
1446 break;
1447
1448 case ROUND_CEIL:
1449 inc = sign == 1 && n != 0;
1450 break;
1451
1452 case ROUND_FLOOR:
1453 default:
1454 inc = sign == -1 && n != 0;
1455 break;
1456 }
1457
1458 if (inc) {
1459
1460 int rh = 1;
1461 for (int i = 0; i < mant.length; i++) {
1462 final int r = mant[i] + rh;
1463 rh = r / RADIX;
1464 mant[i] = r - rh * RADIX;
1465 }
1466
1467 if (rh != 0) {
1468 shiftRight();
1469 mant[mant.length - 1] = rh;
1470 }
1471 }
1472
1473
1474 if (exp < MIN_EXP) {
1475
1476 field.setIEEEFlagsBits(DfpField.FLAG_UNDERFLOW);
1477 return DfpField.FLAG_UNDERFLOW;
1478 }
1479
1480 if (exp > MAX_EXP) {
1481
1482 field.setIEEEFlagsBits(DfpField.FLAG_OVERFLOW);
1483 return DfpField.FLAG_OVERFLOW;
1484 }
1485
1486 if (n != 0) {
1487
1488 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
1489 return DfpField.FLAG_INEXACT;
1490 }
1491
1492 return 0;
1493 }
1494
1495
1496
1497
1498
1499 @Override
1500 public Dfp multiply(final Dfp x) {
1501
1502
1503 if (field.getRadixDigits() != x.field.getRadixDigits()) {
1504 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1505 final Dfp result = newInstance(getZero());
1506 result.nans = QNAN;
1507 return dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1508 }
1509
1510 Dfp result = newInstance(getZero());
1511
1512
1513 if (nans != FINITE || x.nans != FINITE) {
1514 if (isNaN()) {
1515 return this;
1516 }
1517
1518 if (x.isNaN()) {
1519 return x;
1520 }
1521
1522 if (nans == INFINITE && x.nans == FINITE && x.mant[mant.length - 1] != 0) {
1523 result = newInstance(this);
1524 result.sign = (byte) (sign * x.sign);
1525 return result;
1526 }
1527
1528 if (x.nans == INFINITE && nans == FINITE && mant[mant.length - 1] != 0) {
1529 result = newInstance(x);
1530 result.sign = (byte) (sign * x.sign);
1531 return result;
1532 }
1533
1534 if (x.nans == INFINITE && nans == INFINITE) {
1535 result = newInstance(this);
1536 result.sign = (byte) (sign * x.sign);
1537 return result;
1538 }
1539
1540 if ((x.nans == INFINITE && nans == FINITE && mant[mant.length - 1] == 0) ||
1541 (nans == INFINITE && x.nans == FINITE && x.mant[mant.length - 1] == 0)) {
1542 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1543 result = newInstance(getZero());
1544 result.nans = QNAN;
1545 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, x, result);
1546 return result;
1547 }
1548 }
1549
1550 int[] product = new int[mant.length * 2];
1551
1552 for (int i = 0; i < mant.length; i++) {
1553 int rh = 0;
1554 for (int j = 0; j < mant.length; j++) {
1555 int r = mant[i] * x.mant[j];
1556 r += product[i + j] + rh;
1557
1558 rh = r / RADIX;
1559 product[i + j] = r - rh * RADIX;
1560 }
1561 product[i + mant.length] = rh;
1562 }
1563
1564
1565 int md = mant.length * 2 - 1;
1566 for (int i = mant.length * 2 - 1; i >= 0; i--) {
1567 if (product[i] != 0) {
1568 md = i;
1569 break;
1570 }
1571 }
1572
1573
1574 for (int i = 0; i < mant.length; i++) {
1575 result.mant[mant.length - i - 1] = product[md - i];
1576 }
1577
1578
1579 result.exp = exp + x.exp + md - 2 * mant.length + 1;
1580 result.sign = (byte) ((sign == x.sign) ? 1 : -1);
1581
1582 if (result.mant[mant.length - 1] == 0) {
1583
1584 result.exp = 0;
1585 }
1586
1587 final int excp;
1588 if (md > (mant.length - 1)) {
1589 excp = result.round(product[md - mant.length]);
1590 } else {
1591 excp = result.round(0);
1592 }
1593
1594 if (excp != 0) {
1595 result = dotrap(excp, MULTIPLY_TRAP, x, result);
1596 }
1597
1598 return result;
1599 }
1600
1601
1602
1603
1604
1605 @Override
1606 public Dfp multiply(final int x) {
1607 if (x >= 0 && x < RADIX) {
1608 return multiplyFast(x);
1609 } else {
1610 return multiply(newInstance(x));
1611 }
1612 }
1613
1614
1615
1616
1617
1618
1619 private Dfp multiplyFast(final int x) {
1620 Dfp result = newInstance(this);
1621
1622
1623 if (nans != FINITE) {
1624 if (isNaN()) {
1625 return this;
1626 }
1627
1628 if (nans == INFINITE && x != 0) {
1629 result = newInstance(this);
1630 return result;
1631 }
1632
1633 if (nans == INFINITE && x == 0) {
1634 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1635 result = newInstance(getZero());
1636 result.nans = QNAN;
1637 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, newInstance(getZero()), result);
1638 return result;
1639 }
1640 }
1641
1642
1643 if (x < 0 || x >= RADIX) {
1644 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1645 result = newInstance(getZero());
1646 result.nans = QNAN;
1647 result = dotrap(DfpField.FLAG_INVALID, MULTIPLY_TRAP, result, result);
1648 return result;
1649 }
1650
1651 int rh = 0;
1652 for (int i = 0; i < mant.length; i++) {
1653 final int r = mant[i] * x + rh;
1654 rh = r / RADIX;
1655 result.mant[i] = r - rh * RADIX;
1656 }
1657
1658 int lostdigit = 0;
1659 if (rh != 0) {
1660 lostdigit = result.mant[0];
1661 result.shiftRight();
1662 result.mant[mant.length - 1] = rh;
1663 }
1664
1665 if (result.mant[mant.length - 1] == 0) {
1666 result.exp = 0;
1667 }
1668
1669 final int excp = result.round(lostdigit);
1670 if (excp != 0) {
1671 result = dotrap(excp, MULTIPLY_TRAP, result, result);
1672 }
1673
1674 return result;
1675 }
1676
1677
1678
1679
1680
1681 @Override
1682 public Dfp divide(Dfp divisor) {
1683 int[] dividend;
1684 int[] quotient;
1685 int[] remainder;
1686 int qd;
1687 int nsqd;
1688 int trial = 0;
1689 int minadj;
1690 boolean trialgood;
1691 int md = 0;
1692 int excp;
1693
1694
1695 if (field.getRadixDigits() != divisor.field.getRadixDigits()) {
1696 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1697 final Dfp result = newInstance(getZero());
1698 result.nans = QNAN;
1699 return dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1700 }
1701
1702 Dfp result = newInstance(getZero());
1703
1704
1705 if (nans != FINITE || divisor.nans != FINITE) {
1706 if (isNaN()) {
1707 return this;
1708 }
1709
1710 if (divisor.isNaN()) {
1711 return divisor;
1712 }
1713
1714 if (nans == INFINITE && divisor.nans == FINITE) {
1715 result = newInstance(this);
1716 result.sign = (byte) (sign * divisor.sign);
1717 return result;
1718 }
1719
1720 if (divisor.nans == INFINITE && nans == FINITE) {
1721 result = newInstance(getZero());
1722 result.sign = (byte) (sign * divisor.sign);
1723 return result;
1724 }
1725
1726 if (divisor.nans == INFINITE && nans == INFINITE) {
1727 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1728 result = newInstance(getZero());
1729 result.nans = QNAN;
1730 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, divisor, result);
1731 return result;
1732 }
1733 }
1734
1735
1736 if (divisor.mant[mant.length - 1] == 0) {
1737 field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1738 result = newInstance(getZero());
1739 result.sign = (byte) (sign * divisor.sign);
1740 result.nans = INFINITE;
1741 result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, divisor, result);
1742 return result;
1743 }
1744
1745 dividend = new int[mant.length + 1];
1746 quotient = new int[mant.length + 2];
1747 remainder = new int[mant.length + 1];
1748
1749
1750
1751 dividend[mant.length] = 0;
1752 quotient[mant.length] = 0;
1753 quotient[mant.length + 1] = 0;
1754 remainder[mant.length] = 0;
1755
1756
1757
1758 for (int i = 0; i < mant.length; i++) {
1759 dividend[i] = mant[i];
1760 quotient[i] = 0;
1761 remainder[i] = 0;
1762 }
1763
1764
1765 nsqd = 0;
1766 for (qd = mant.length + 1; qd >= 0; qd--) {
1767
1768
1769
1770 final int divMsb = dividend[mant.length] * RADIX + dividend[mant.length - 1];
1771 int min = divMsb / (divisor.mant[mant.length - 1] + 1);
1772 int max = (divMsb + 1) / divisor.mant[mant.length - 1];
1773
1774 trialgood = false;
1775 while (!trialgood) {
1776
1777 trial = (min + max) / 2;
1778
1779
1780 int rh = 0;
1781 for (int i = 0; i < mant.length + 1; i++) {
1782 int dm = (i < mant.length) ? divisor.mant[i] : 0;
1783 final int r = (dm * trial) + rh;
1784 rh = r / RADIX;
1785 remainder[i] = r - rh * RADIX;
1786 }
1787
1788
1789 rh = 1;
1790 for (int i = 0; i < mant.length + 1; i++) {
1791 final int r = ((RADIX - 1) - remainder[i]) + dividend[i] + rh;
1792 rh = r / RADIX;
1793 remainder[i] = r - rh * RADIX;
1794 }
1795
1796
1797 if (rh == 0) {
1798
1799 max = trial - 1;
1800 continue;
1801 }
1802
1803
1804 minadj = (remainder[mant.length] * RADIX) + remainder[mant.length - 1];
1805 minadj /= divisor.mant[mant.length - 1] + 1;
1806
1807 if (minadj >= 2) {
1808 min = trial + minadj;
1809 continue;
1810 }
1811
1812
1813
1814 trialgood = false;
1815 for (int i = mant.length - 1; i >= 0; i--) {
1816 if (divisor.mant[i] > remainder[i]) {
1817 trialgood = true;
1818 break;
1819 }
1820 if (divisor.mant[i] < remainder[i]) {
1821 break;
1822 }
1823 }
1824
1825 if (remainder[mant.length] != 0) {
1826 trialgood = false;
1827 }
1828
1829 if (!trialgood) {
1830 min = trial + 1;
1831 }
1832 }
1833
1834
1835 quotient[qd] = trial;
1836 if (trial != 0 || nsqd != 0) {
1837 nsqd++;
1838 }
1839
1840 if (field.getRoundingMode() == DfpField.RoundingMode.ROUND_DOWN && nsqd == mant.length) {
1841
1842 break;
1843 }
1844
1845 if (nsqd > mant.length) {
1846
1847 break;
1848 }
1849
1850
1851 dividend[0] = 0;
1852 System.arraycopy(remainder, 0, dividend, 1, mant.length);
1853 }
1854
1855
1856 md = mant.length;
1857 for (int i = mant.length + 1; i >= 0; i--) {
1858 if (quotient[i] != 0) {
1859 md = i;
1860 break;
1861 }
1862 }
1863
1864
1865 for (int i = 0; i < mant.length; i++) {
1866 result.mant[mant.length - i - 1] = quotient[md - i];
1867 }
1868
1869
1870 result.exp = exp - divisor.exp + md - mant.length;
1871 result.sign = (byte) ((sign == divisor.sign) ? 1 : -1);
1872
1873 if (result.mant[mant.length - 1] == 0) {
1874 result.exp = 0;
1875 }
1876
1877 if (md > (mant.length - 1)) {
1878 excp = result.round(quotient[md - mant.length]);
1879 } else {
1880 excp = result.round(0);
1881 }
1882
1883 if (excp != 0) {
1884 result = dotrap(excp, DIVIDE_TRAP, divisor, result);
1885 }
1886
1887 return result;
1888 }
1889
1890
1891
1892
1893
1894
1895 public Dfp divide(int divisor) {
1896
1897
1898 if (nans != FINITE) {
1899 if (isNaN()) {
1900 return this;
1901 }
1902
1903 if (nans == INFINITE) {
1904 return newInstance(this);
1905 }
1906 }
1907
1908
1909 if (divisor == 0) {
1910 field.setIEEEFlagsBits(DfpField.FLAG_DIV_ZERO);
1911 Dfp result = newInstance(getZero());
1912 result.sign = sign;
1913 result.nans = INFINITE;
1914 result = dotrap(DfpField.FLAG_DIV_ZERO, DIVIDE_TRAP, getZero(), result);
1915 return result;
1916 }
1917
1918
1919 if (divisor < 0 || divisor >= RADIX) {
1920 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1921 Dfp result = newInstance(getZero());
1922 result.nans = QNAN;
1923 result = dotrap(DfpField.FLAG_INVALID, DIVIDE_TRAP, result, result);
1924 return result;
1925 }
1926
1927 Dfp result = newInstance(this);
1928
1929 int rl = 0;
1930 for (int i = mant.length - 1; i >= 0; i--) {
1931 final int r = rl * RADIX + result.mant[i];
1932 final int rh = r / divisor;
1933 rl = r - rh * divisor;
1934 result.mant[i] = rh;
1935 }
1936
1937 if (result.mant[mant.length - 1] == 0) {
1938
1939 result.shiftLeft();
1940 final int r = rl * RADIX;
1941 final int rh = r / divisor;
1942 rl = r - rh * divisor;
1943 result.mant[0] = rh;
1944 }
1945
1946 final int excp = result.round(rl * RADIX / divisor);
1947 if (excp != 0) {
1948 result = dotrap(excp, DIVIDE_TRAP, result, result);
1949 }
1950
1951 return result;
1952 }
1953
1954
1955 @Override
1956 public Dfp reciprocal() {
1957 return field.getOne().divide(this);
1958 }
1959
1960
1961
1962
1963
1964 @Override
1965 public Dfp sqrt() {
1966
1967
1968 if (nans == FINITE && mant[mant.length - 1] == 0) {
1969
1970 return newInstance(this);
1971 }
1972
1973 if (nans != FINITE) {
1974 if (nans == INFINITE && sign == 1) {
1975
1976 return newInstance(this);
1977 }
1978
1979 if (nans == QNAN) {
1980 return newInstance(this);
1981 }
1982
1983 if (nans == SNAN) {
1984 Dfp result;
1985
1986 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1987 result = newInstance(this);
1988 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
1989 return result;
1990 }
1991 }
1992
1993 if (sign == -1) {
1994
1995 Dfp result;
1996
1997 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
1998 result = newInstance(this);
1999 result.nans = QNAN;
2000 result = dotrap(DfpField.FLAG_INVALID, SQRT_TRAP, null, result);
2001 return result;
2002 }
2003
2004 Dfp x = newInstance(this);
2005
2006
2007 if (x.exp < -1 || x.exp > 1) {
2008 x.exp = this.exp / 2;
2009 }
2010
2011
2012 switch (x.mant[mant.length - 1] / 2000) {
2013 case 0:
2014 x.mant[mant.length - 1] = x.mant[mant.length - 1] / 2 + 1;
2015 break;
2016 case 2:
2017 x.mant[mant.length - 1] = 1500;
2018 break;
2019 case 3:
2020 x.mant[mant.length - 1] = 2200;
2021 break;
2022 default:
2023 x.mant[mant.length - 1] = 3000;
2024 }
2025
2026 Dfp dx = newInstance(x);
2027
2028
2029
2030
2031 Dfp px = getZero();
2032 Dfp ppx = getZero();
2033 while (x.unequal(px)) {
2034 dx = newInstance(x);
2035 dx.sign = -1;
2036 dx = dx.add(this.divide(x));
2037 dx = dx.divide(2);
2038 ppx = px;
2039 px = x;
2040 x = x.add(dx);
2041
2042 if (x.equals(ppx)) {
2043
2044 break;
2045 }
2046
2047
2048
2049 if (dx.mant[mant.length - 1] == 0) {
2050 break;
2051 }
2052 }
2053
2054 return x;
2055 }
2056
2057
2058
2059
2060 @Override
2061 public String toString() {
2062 if (nans != FINITE) {
2063
2064 if (nans == INFINITE) {
2065 return (sign < 0) ? NEG_INFINITY_STRING : POS_INFINITY_STRING;
2066 } else {
2067 return NAN_STRING;
2068 }
2069 }
2070
2071 if (exp > mant.length || exp < -1) {
2072 return dfp2sci();
2073 }
2074
2075 return dfp2string();
2076 }
2077
2078
2079
2080
2081 protected String dfp2sci() {
2082 char[] rawdigits = new char[mant.length * 4];
2083 char[] outputbuffer = new char[mant.length * 4 + 20];
2084 int p;
2085 int q;
2086 int e;
2087 int ae;
2088 int shf;
2089
2090
2091 p = 0;
2092 for (int i = mant.length - 1; i >= 0; i--) {
2093 rawdigits[p++] = (char) ((mant[i] / 1000) + '0');
2094 rawdigits[p++] = (char) (((mant[i] / 100) % 10) + '0');
2095 rawdigits[p++] = (char) (((mant[i] / 10) % 10) + '0');
2096 rawdigits[p++] = (char) (((mant[i]) % 10) + '0');
2097 }
2098
2099
2100 for (p = 0; p < rawdigits.length; p++) {
2101 if (rawdigits[p] != '0') {
2102 break;
2103 }
2104 }
2105 shf = p;
2106
2107
2108 q = 0;
2109 if (sign == -1) {
2110 outputbuffer[q++] = '-';
2111 }
2112
2113 if (p != rawdigits.length) {
2114
2115 outputbuffer[q++] = rawdigits[p++];
2116 outputbuffer[q++] = '.';
2117
2118 while (p < rawdigits.length) {
2119 outputbuffer[q++] = rawdigits[p++];
2120 }
2121 } else {
2122 outputbuffer[q++] = '0';
2123 outputbuffer[q++] = '.';
2124 outputbuffer[q++] = '0';
2125 outputbuffer[q++] = 'e';
2126 outputbuffer[q++] = '0';
2127 return String.valueOf(outputbuffer, 0, 5);
2128 }
2129
2130 outputbuffer[q++] = 'e';
2131
2132
2133
2134 e = exp * 4 - shf - 1;
2135 ae = e;
2136 if (e < 0) {
2137 ae = -e;
2138 }
2139
2140
2141 p = 1000000000;
2142 while (p > ae) {
2143 p /= 10;
2144 }
2145
2146 if (e < 0) {
2147 outputbuffer[q++] = '-';
2148 }
2149
2150 while (p > 0) {
2151 outputbuffer[q++] = (char)(ae / p + '0');
2152 ae %= p;
2153 p /= 10;
2154 }
2155
2156 return String.valueOf(outputbuffer, 0, q);
2157 }
2158
2159
2160
2161
2162 protected String dfp2string() {
2163 char[] buffer = new char[mant.length * 4 + 20];
2164 int p = 1;
2165 int q;
2166 int e = exp;
2167 boolean pointInserted = false;
2168
2169 buffer[0] = ' ';
2170
2171 if (e <= 0) {
2172 buffer[p++] = '0';
2173 buffer[p++] = '.';
2174 pointInserted = true;
2175 }
2176
2177 while (e < 0) {
2178 buffer[p++] = '0';
2179 buffer[p++] = '0';
2180 buffer[p++] = '0';
2181 buffer[p++] = '0';
2182 e++;
2183 }
2184
2185 for (int i = mant.length - 1; i >= 0; i--) {
2186 buffer[p++] = (char) ((mant[i] / 1000) + '0');
2187 buffer[p++] = (char) (((mant[i] / 100) % 10) + '0');
2188 buffer[p++] = (char) (((mant[i] / 10) % 10) + '0');
2189 buffer[p++] = (char) (((mant[i]) % 10) + '0');
2190 if (--e == 0) {
2191 buffer[p++] = '.';
2192 pointInserted = true;
2193 }
2194 }
2195
2196 while (e > 0) {
2197 buffer[p++] = '0';
2198 buffer[p++] = '0';
2199 buffer[p++] = '0';
2200 buffer[p++] = '0';
2201 e--;
2202 }
2203
2204 if (!pointInserted) {
2205
2206 buffer[p++] = '.';
2207 }
2208
2209
2210 q = 1;
2211 while (buffer[q] == '0') {
2212 q++;
2213 }
2214 if (buffer[q] == '.') {
2215 q--;
2216 }
2217
2218
2219 while (buffer[p - 1] == '0') {
2220 p--;
2221 }
2222
2223
2224 if (sign < 0) {
2225 buffer[--q] = '-';
2226 }
2227
2228 return String.valueOf(buffer, q, p - q);
2229 }
2230
2231
2232
2233
2234
2235
2236
2237
2238 public Dfp dotrap(int type, String what, Dfp oper, Dfp result) {
2239 Dfp def = result;
2240
2241 switch (type) {
2242 case DfpField.FLAG_INVALID:
2243 def = newInstance(getZero());
2244 def.sign = result.sign;
2245 def.nans = QNAN;
2246 break;
2247
2248 case DfpField.FLAG_DIV_ZERO:
2249 if (nans == FINITE && mant[mant.length - 1] != 0) {
2250
2251 def = newInstance(getZero());
2252 def.sign = (byte) (sign * oper.sign);
2253 def.nans = INFINITE;
2254 }
2255
2256 if (nans == FINITE && mant[mant.length - 1] == 0) {
2257
2258 def = newInstance(getZero());
2259 def.nans = QNAN;
2260 }
2261
2262 if (nans == INFINITE || nans == QNAN) {
2263 def = newInstance(getZero());
2264 def.nans = QNAN;
2265 }
2266
2267 if (nans == INFINITE || nans == SNAN) {
2268 def = newInstance(getZero());
2269 def.nans = QNAN;
2270 }
2271 break;
2272
2273 case DfpField.FLAG_UNDERFLOW:
2274 if ((result.exp + mant.length) < MIN_EXP) {
2275 def = newInstance(getZero());
2276 def.sign = result.sign;
2277 } else {
2278 def = newInstance(result);
2279 }
2280 result.exp += ERR_SCALE;
2281 break;
2282
2283 case DfpField.FLAG_OVERFLOW:
2284 result.exp -= ERR_SCALE;
2285 def = newInstance(getZero());
2286 def.sign = result.sign;
2287 def.nans = INFINITE;
2288 break;
2289
2290 default: def = result; break;
2291 }
2292
2293 return trap(type, what, oper, def, result);
2294 }
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307 protected Dfp trap(int type, String what, Dfp oper, Dfp def, Dfp result) {
2308 return def;
2309 }
2310
2311
2312
2313
2314 public int classify() {
2315 return nans;
2316 }
2317
2318
2319
2320
2321
2322
2323
2324 public static Dfp copySign(final Dfp x, final Dfp y) {
2325 Dfp result = x.newInstance(x);
2326 result.sign = y.sign;
2327 return result;
2328 }
2329
2330
2331
2332
2333
2334
2335 public Dfp nextAfter(final Dfp x) {
2336
2337
2338 if (field.getRadixDigits() != x.field.getRadixDigits()) {
2339 field.setIEEEFlagsBits(DfpField.FLAG_INVALID);
2340 final Dfp result = newInstance(getZero());
2341 result.nans = QNAN;
2342 return dotrap(DfpField.FLAG_INVALID, NEXT_AFTER_TRAP, x, result);
2343 }
2344
2345
2346 boolean up = this.lessThan(x);
2347
2348 if (compare(this, x) == 0) {
2349 return newInstance(x);
2350 }
2351
2352 if (lessThan(getZero())) {
2353 up = !up;
2354 }
2355
2356 final Dfp inc;
2357 Dfp result;
2358 if (up) {
2359 inc = newInstance(getOne());
2360 inc.exp = this.exp - mant.length + 1;
2361 inc.sign = this.sign;
2362
2363 if (this.equals(getZero())) {
2364 inc.exp = MIN_EXP - mant.length;
2365 }
2366
2367 result = add(inc);
2368 } else {
2369 inc = newInstance(getOne());
2370 inc.exp = this.exp;
2371 inc.sign = this.sign;
2372
2373 if (this.equals(inc)) {
2374 inc.exp = this.exp - mant.length;
2375 } else {
2376 inc.exp = this.exp - mant.length + 1;
2377 }
2378
2379 if (this.equals(getZero())) {
2380 inc.exp = MIN_EXP - mant.length;
2381 }
2382
2383 result = this.subtract(inc);
2384 }
2385
2386 if (result.classify() == INFINITE && this.classify() != INFINITE) {
2387 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2388 result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2389 }
2390
2391 if (result.equals(getZero()) && !this.equals(getZero())) {
2392 field.setIEEEFlagsBits(DfpField.FLAG_INEXACT);
2393 result = dotrap(DfpField.FLAG_INEXACT, NEXT_AFTER_TRAP, x, result);
2394 }
2395
2396 return result;
2397 }
2398
2399
2400
2401
2402
2403 public double toDouble() {
2404
2405 if (isInfinite()) {
2406 if (lessThan(getZero())) {
2407 return Double.NEGATIVE_INFINITY;
2408 } else {
2409 return Double.POSITIVE_INFINITY;
2410 }
2411 }
2412
2413 if (isNaN()) {
2414 return Double.NaN;
2415 }
2416
2417 Dfp y = this;
2418 boolean negate = false;
2419 int cmp0 = compare(this, getZero());
2420 if (cmp0 == 0) {
2421 return sign < 0 ? -0.0 : +0.0;
2422 } else if (cmp0 < 0) {
2423 y = negate();
2424 negate = true;
2425 }
2426
2427
2428
2429 int exponent = (int)(y.intLog10() * 3.32);
2430 if (exponent < 0) {
2431 exponent--;
2432 }
2433
2434 Dfp tempDfp = DfpMath.pow(getTwo(), exponent);
2435 while (tempDfp.lessThan(y) || tempDfp.equals(y)) {
2436 tempDfp = tempDfp.multiply(2);
2437 exponent++;
2438 }
2439 exponent--;
2440
2441
2442
2443 y = y.divide(DfpMath.pow(getTwo(), exponent));
2444 if (exponent > -1023) {
2445 y = y.subtract(getOne());
2446 }
2447
2448 if (exponent < -1074) {
2449 return 0;
2450 }
2451
2452 if (exponent > 1023) {
2453 return negate ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
2454 }
2455
2456
2457 y = y.multiply(newInstance(4503599627370496L)).rint();
2458 String str = y.toString();
2459 str = str.substring(0, str.length() - 1);
2460 long mantissa = Long.parseLong(str);
2461
2462 if (mantissa == 4503599627370496L) {
2463
2464 mantissa = 0;
2465 exponent++;
2466 }
2467
2468
2469 if (exponent <= -1023) {
2470 exponent--;
2471 }
2472
2473 while (exponent < -1023) {
2474 exponent++;
2475 mantissa >>>= 1;
2476 }
2477
2478 long bits = mantissa | ((exponent + 1023L) << 52);
2479 double x = Double.longBitsToDouble(bits);
2480
2481 if (negate) {
2482 x = -x;
2483 }
2484
2485 return x;
2486 }
2487
2488
2489
2490
2491
2492 public double[] toSplitDouble() {
2493 double[] split = new double[2];
2494 long mask = 0xffffffffc0000000L;
2495
2496 split[0] = Double.longBitsToDouble(Double.doubleToLongBits(toDouble()) & mask);
2497 split[1] = subtract(newInstance(split[0])).toDouble();
2498
2499 return split;
2500 }
2501
2502
2503
2504
2505 @Override
2506 public double getReal() {
2507 return toDouble();
2508 }
2509
2510
2511
2512
2513 @Override
2514 public Dfp add(final double a) {
2515 return add(newInstance(a));
2516 }
2517
2518
2519
2520
2521 @Override
2522 public Dfp subtract(final double a) {
2523 return subtract(newInstance(a));
2524 }
2525
2526
2527
2528
2529 @Override
2530 public Dfp multiply(final double a) {
2531 return multiply(newInstance(a));
2532 }
2533
2534
2535
2536
2537 @Override
2538 public Dfp divide(final double a) {
2539 return divide(newInstance(a));
2540 }
2541
2542
2543
2544
2545 @Override
2546 public Dfp remainder(final double a) {
2547 return remainder(newInstance(a));
2548 }
2549
2550
2551
2552
2553 @Override
2554 public long round() {
2555 return Math.round(toDouble());
2556 }
2557
2558
2559
2560
2561 @Override
2562 public Dfp signum() {
2563 if (isNaN() || isZero()) {
2564 return this;
2565 } else {
2566 return newInstance(sign > 0 ? +1 : -1);
2567 }
2568 }
2569
2570
2571
2572
2573 @Override
2574 public Dfp copySign(final Dfp s) {
2575 if ((sign >= 0 && s.sign >= 0) || (sign < 0 && s.sign < 0)) {
2576 return this;
2577 }
2578 return negate();
2579 }
2580
2581
2582
2583
2584 @Override
2585 public Dfp copySign(final double s) {
2586 long sb = Double.doubleToLongBits(s);
2587 if ((sign >= 0 && sb >= 0) || (sign < 0 && sb < 0)) {
2588 return this;
2589 }
2590 return negate();
2591 }
2592
2593
2594
2595
2596 @Override
2597 public Dfp scalb(final int n) {
2598 return multiply(DfpMath.pow(getTwo(), n));
2599 }
2600
2601
2602
2603
2604 @Override
2605 public Dfp hypot(final Dfp y) {
2606 return multiply(this).add(y.multiply(y)).sqrt();
2607 }
2608
2609
2610
2611
2612 @Override
2613 public Dfp cbrt() {
2614 return rootN(3);
2615 }
2616
2617
2618
2619
2620 @Override
2621 public Dfp rootN(final int n) {
2622 return (sign >= 0) ?
2623 DfpMath.pow(this, getOne().divide(n)) :
2624 DfpMath.pow(negate(), getOne().divide(n)).negate();
2625 }
2626
2627
2628
2629
2630 @Override
2631 public Dfp pow(final double p) {
2632 return DfpMath.pow(this, newInstance(p));
2633 }
2634
2635
2636
2637
2638 @Override
2639 public Dfp pow(final int n) {
2640 return DfpMath.pow(this, n);
2641 }
2642
2643
2644
2645
2646 @Override
2647 public Dfp pow(final Dfp e) {
2648 return DfpMath.pow(this, e);
2649 }
2650
2651
2652
2653
2654 @Override
2655 public Dfp exp() {
2656 return DfpMath.exp(this);
2657 }
2658
2659
2660
2661
2662 @Override
2663 public Dfp expm1() {
2664 return DfpMath.exp(this).subtract(getOne());
2665 }
2666
2667
2668
2669
2670 @Override
2671 public Dfp log() {
2672 return DfpMath.log(this);
2673 }
2674
2675
2676
2677
2678 @Override
2679 public Dfp log1p() {
2680 return DfpMath.log(this.add(getOne()));
2681 }
2682
2683
2684
2685
2686 @Override
2687 public Dfp log10() {
2688 return DfpMath.log(this).divide(DfpMath.log(newInstance(10)));
2689 }
2690
2691
2692
2693
2694 @Override
2695 public Dfp cos() {
2696 return DfpMath.cos(this);
2697 }
2698
2699
2700
2701
2702 @Override
2703 public Dfp sin() {
2704 return DfpMath.sin(this);
2705 }
2706
2707
2708
2709
2710 @Override
2711 public Dfp tan() {
2712 return DfpMath.tan(this);
2713 }
2714
2715
2716
2717
2718 @Override
2719 public Dfp acos() {
2720 return DfpMath.acos(this);
2721 }
2722
2723
2724
2725
2726 @Override
2727 public Dfp asin() {
2728 return DfpMath.asin(this);
2729 }
2730
2731
2732
2733
2734 @Override
2735 public Dfp atan() {
2736 return DfpMath.atan(this);
2737 }
2738
2739
2740
2741
2742 @Override
2743 public Dfp atan2(final Dfp x)
2744 throws DimensionMismatchException {
2745
2746
2747 final Dfp r = x.multiply(x).add(multiply(this)).sqrt();
2748
2749 if (x.sign >= 0) {
2750
2751
2752 return getTwo().multiply(divide(r.add(x)).atan());
2753 } else {
2754
2755
2756 final Dfp tmp = getTwo().multiply(divide(r.subtract(x)).atan());
2757 final Dfp pmPi = newInstance((tmp.sign <= 0) ? -Math.PI : Math.PI);
2758 return pmPi.subtract(tmp);
2759 }
2760 }
2761
2762
2763
2764
2765 @Override
2766 public Dfp cosh() {
2767 return DfpMath.exp(this).add(DfpMath.exp(negate())).divide(2);
2768 }
2769
2770
2771
2772
2773 @Override
2774 public Dfp sinh() {
2775 return DfpMath.exp(this).subtract(DfpMath.exp(negate())).divide(2);
2776 }
2777
2778
2779
2780
2781 @Override
2782 public Dfp tanh() {
2783 final Dfp ePlus = DfpMath.exp(this);
2784 final Dfp eMinus = DfpMath.exp(negate());
2785 return ePlus.subtract(eMinus).divide(ePlus.add(eMinus));
2786 }
2787
2788
2789
2790
2791 @Override
2792 public Dfp acosh() {
2793 return multiply(this).subtract(getOne()).sqrt().add(this).log();
2794 }
2795
2796
2797
2798
2799 @Override
2800 public Dfp asinh() {
2801 return multiply(this).add(getOne()).sqrt().add(this).log();
2802 }
2803
2804
2805
2806
2807 @Override
2808 public Dfp atanh() {
2809 return getOne().add(this).divide(getOne().subtract(this)).log().divide(2);
2810 }
2811
2812
2813
2814
2815 @Override
2816 public Dfp linearCombination(final Dfp[] a, final Dfp[] b)
2817 throws DimensionMismatchException {
2818 if (a.length != b.length) {
2819 throw new DimensionMismatchException(a.length, b.length);
2820 }
2821 Dfp r = getZero();
2822 for (int i = 0; i < a.length; ++i) {
2823 r = r.add(a[i].multiply(b[i]));
2824 }
2825 return r;
2826 }
2827
2828
2829
2830
2831 @Override
2832 public Dfp linearCombination(final double[] a, final Dfp[] b)
2833 throws DimensionMismatchException {
2834 if (a.length != b.length) {
2835 throw new DimensionMismatchException(a.length, b.length);
2836 }
2837 Dfp r = getZero();
2838 for (int i = 0; i < a.length; ++i) {
2839 r = r.add(b[i].multiply(a[i]));
2840 }
2841 return r;
2842 }
2843
2844
2845
2846
2847 @Override
2848 public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2) {
2849 return a1.multiply(b1).add(a2.multiply(b2));
2850 }
2851
2852
2853
2854
2855 @Override
2856 public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2) {
2857 return b1.multiply(a1).add(b2.multiply(a2));
2858 }
2859
2860
2861
2862
2863 @Override
2864 public Dfp linearCombination(final Dfp a1, final Dfp b1,
2865 final Dfp a2, final Dfp b2,
2866 final Dfp a3, final Dfp b3) {
2867 return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
2868 }
2869
2870
2871
2872
2873 @Override
2874 public Dfp linearCombination(final double a1, final Dfp b1,
2875 final double a2, final Dfp b2,
2876 final double a3, final Dfp b3) {
2877 return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
2878 }
2879
2880
2881
2882
2883 @Override
2884 public Dfp linearCombination(final Dfp a1, final Dfp b1, final Dfp a2, final Dfp b2,
2885 final Dfp a3, final Dfp b3, final Dfp a4, final Dfp b4) {
2886 return a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
2887 }
2888
2889
2890
2891
2892 @Override
2893 public Dfp linearCombination(final double a1, final Dfp b1, final double a2, final Dfp b2,
2894 final double a3, final Dfp b3, final double a4, final Dfp b4) {
2895 return b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
2896 }
2897 }