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