dp_add.c 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165
  1. // SPDX-License-Identifier: GPL-2.0-only
  2. /* IEEE754 floating point arithmetic
  3. * double precision: common utilities
  4. */
  5. /*
  6. * MIPS floating point support
  7. * Copyright (C) 1994-2000 Algorithmics Ltd.
  8. */
  9. #include "ieee754dp.h"
  10. union ieee754dp ieee754dp_add(union ieee754dp x, union ieee754dp y)
  11. {
  12. int s;
  13. COMPXDP;
  14. COMPYDP;
  15. EXPLODEXDP;
  16. EXPLODEYDP;
  17. ieee754_clearcx();
  18. FLUSHXDP;
  19. FLUSHYDP;
  20. switch (CLPAIR(xc, yc)) {
  21. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
  22. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
  23. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
  24. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
  25. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
  26. return ieee754dp_nanxcpt(y);
  27. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
  28. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
  29. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
  30. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
  31. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
  32. case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
  33. return ieee754dp_nanxcpt(x);
  34. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
  35. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
  36. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
  37. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
  38. return y;
  39. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
  40. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
  41. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
  42. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
  43. case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
  44. return x;
  45. /*
  46. * Infinity handling
  47. */
  48. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
  49. if (xs == ys)
  50. return x;
  51. ieee754_setcx(IEEE754_INVALID_OPERATION);
  52. return ieee754dp_indef();
  53. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
  54. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
  55. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
  56. return y;
  57. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
  58. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
  59. case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
  60. return x;
  61. /*
  62. * Zero handling
  63. */
  64. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
  65. if (xs == ys)
  66. return x;
  67. else
  68. return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
  69. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
  70. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
  71. return x;
  72. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
  73. case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
  74. return y;
  75. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
  76. DPDNORMX;
  77. fallthrough;
  78. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
  79. DPDNORMY;
  80. break;
  81. case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
  82. DPDNORMX;
  83. break;
  84. case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
  85. break;
  86. }
  87. assert(xm & DP_HIDDEN_BIT);
  88. assert(ym & DP_HIDDEN_BIT);
  89. /*
  90. * Provide guard,round and stick bit space.
  91. */
  92. xm <<= 3;
  93. ym <<= 3;
  94. if (xe > ye) {
  95. /*
  96. * Have to shift y fraction right to align.
  97. */
  98. s = xe - ye;
  99. ym = XDPSRS(ym, s);
  100. ye += s;
  101. } else if (ye > xe) {
  102. /*
  103. * Have to shift x fraction right to align.
  104. */
  105. s = ye - xe;
  106. xm = XDPSRS(xm, s);
  107. xe += s;
  108. }
  109. assert(xe == ye);
  110. assert(xe <= DP_EMAX);
  111. if (xs == ys) {
  112. /*
  113. * Generate 28 bit result of adding two 27 bit numbers
  114. * leaving result in xm, xs and xe.
  115. */
  116. xm = xm + ym;
  117. if (xm >> (DP_FBITS + 1 + 3)) { /* carry out */
  118. xm = XDPSRS1(xm);
  119. xe++;
  120. }
  121. } else {
  122. if (xm >= ym) {
  123. xm = xm - ym;
  124. } else {
  125. xm = ym - xm;
  126. xs = ys;
  127. }
  128. if (xm == 0)
  129. return ieee754dp_zero(ieee754_csr.rm == FPU_CSR_RD);
  130. /*
  131. * Normalize to rounding precision.
  132. */
  133. while ((xm >> (DP_FBITS + 3)) == 0) {
  134. xm <<= 1;
  135. xe--;
  136. }
  137. }
  138. return ieee754dp_format(xs, xe, xm);
  139. }