div_ext.c 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
  1. /*
  2. (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
  3. See the copyright notice in the ACK home directory, in the file "Copyright".
  4. */
  5. /* $Id$ */
  6. /*
  7. DIVIDE EXTENDED FORMAT
  8. */
  9. #include "FP_bias.h"
  10. #include "FP_trap.h"
  11. #include "FP_types.h"
  12. /*
  13. November 15, 1984
  14. This is a routine to do the work.
  15. There are two versions:
  16. One is based on the partial products method
  17. and makes no use possible machine instructions
  18. to divide (hardware dividers).
  19. The other is used when USE_DIVIDE is defined. It is much faster on
  20. machines with fast 4 byte operations.
  21. */
  22. /********************************************************/
  23. void
  24. div_ext(e1,e2)
  25. EXTEND *e1,*e2;
  26. {
  27. short error = 0;
  28. B64 result;
  29. register unsigned long *lp;
  30. #ifndef USE_DIVIDE
  31. short count;
  32. #else
  33. unsigned short u[9], v[5];
  34. register int j;
  35. register unsigned short *u_p = u;
  36. int maxv = 4;
  37. #endif
  38. if ((e2->m1 | e2->m2) == 0) {
  39. /*
  40. * Exception 8.2 - Divide by zero
  41. */
  42. trap(EFDIVZ);
  43. e1->m1 = e1->m2 = 0L;
  44. e1->exp = EXT_MAX;
  45. return;
  46. }
  47. if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */
  48. e1->exp = 0; /* make sure */
  49. return;
  50. }
  51. #ifndef USE_DIVIDE
  52. /*
  53. * numbers are right shifted one bit to make sure
  54. * that m1 is quaranteed to be larger if its
  55. * maximum bit is set
  56. */
  57. b64_rsft(&e1->mantissa); /* 64 bit shift right */
  58. b64_rsft(&e2->mantissa); /* 64 bit shift right */
  59. e1->exp++;
  60. e2->exp++;
  61. #endif
  62. /* check for underflow, divide by zero, etc */
  63. e1->sign ^= e2->sign;
  64. e1->exp -= e2->exp;
  65. #ifndef USE_DIVIDE
  66. /* do division of mantissas */
  67. /* uses partial product method */
  68. /* init control variables */
  69. count = 64;
  70. result.h_32 = 0L;
  71. result.l_32 = 0L;
  72. /* partial product division loop */
  73. while (count--) {
  74. /* first left shift result 1 bit */
  75. /* this is ALWAYS done */
  76. b64_lsft(&result);
  77. /* compare dividend and divisor */
  78. /* if dividend >= divisor add a bit */
  79. /* and subtract divisior from dividend */
  80. if ( (e1->m1 < e2->m1) ||
  81. ((e1->m1 == e2->m1) && (e1->m2 < e2->m2) ))
  82. ; /* null statement */
  83. /* i.e., don't add or subtract */
  84. else {
  85. result.l_32++; /* ADD */
  86. if (e2->m2 > e1->m2)
  87. e1->m1 -= 1; /* carry in */
  88. e1->m1 -= e2->m1; /* do SUBTRACTION */
  89. e1->m2 -= e2->m2; /* SUBTRACTION */
  90. }
  91. /* shift dividend left one bit OR */
  92. /* IF it equals ZERO we can break out */
  93. /* of the loop, but still must shift */
  94. /* the quotient the remaining count bits */
  95. /* NB save the results of this test in error */
  96. /* if not zero, then the result is inexact. */
  97. /* this would be reported in IEEE standard */
  98. /* lp points to dividend */
  99. lp = &e1->m1;
  100. error = ((*lp | *(lp+1)) != 0L) ? 1 : 0;
  101. if (error) { /* more work */
  102. /* assume max bit == 0 (see above) */
  103. b64_lsft(&e1->mantissa);
  104. continue;
  105. }
  106. else
  107. break; /* leave loop */
  108. } /* end of divide by subtraction loop */
  109. if (count > 0) {
  110. lp = &result.h_32;
  111. if (count > 31) { /* move to higher word */
  112. *lp = *(lp+1);
  113. count -= 32;
  114. *(lp+1) = 0L; /* clear low word */
  115. }
  116. if (*lp)
  117. *lp <<= count; /* shift rest of way */
  118. lp++; /* == &result.l_32 */
  119. if (*lp) {
  120. result.h_32 |= (*lp >> 32-count);
  121. *lp <<= count;
  122. }
  123. }
  124. #else /* USE_DIVIDE */
  125. u[4] = (e1->m2 & 1) << 15;
  126. b64_rsft(&(e1->mantissa));
  127. u[0] = e1->m1 >> 16;
  128. u[1] = e1->m1;
  129. u[2] = e1->m2 >> 16;
  130. u[3] = e1->m2;
  131. u[5] = 0; u[6] = 0; u[7] = 0;
  132. v[1] = e2->m1 >> 16;
  133. v[2] = e2->m1;
  134. v[3] = e2->m2 >> 16;
  135. v[4] = e2->m2;
  136. while (! v[maxv]) maxv--;
  137. result.h_32 = 0;
  138. result.l_32 = 0;
  139. lp = &result.h_32;
  140. /*
  141. * Use an algorithm of Knuth (The art of programming, Seminumerical
  142. * algorithms), to divide u by v. u and v are both seen as numbers
  143. * with base 65536.
  144. */
  145. for (j = 0; j <= 3; j++, u_p++) {
  146. unsigned long q_est, temp;
  147. if (j == 2) lp++;
  148. if (u_p[0] == 0 && u_p[1] < v[1]) continue;
  149. temp = ((unsigned long)u_p[0] << 16) + u_p[1];
  150. if (u_p[0] >= v[1]) {
  151. q_est = 0x0000FFFFL;
  152. }
  153. else {
  154. q_est = temp / v[1];
  155. }
  156. temp -= q_est * v[1];
  157. while (temp < 0x10000 && v[2]*q_est > ((temp<<16)+u_p[2])) {
  158. q_est--;
  159. temp += v[1];
  160. }
  161. /* Now, according to Knuth, we have an estimate of the
  162. quotient, that is either correct or one too big, but
  163. almost always correct.
  164. */
  165. if (q_est != 0) {
  166. int i;
  167. unsigned long k = 0;
  168. int borrow = 0;
  169. for (i = maxv; i > 0; i--) {
  170. unsigned long tmp = q_est * v[i] + k + borrow;
  171. unsigned short md = tmp;
  172. borrow = (md > u_p[i]);
  173. u_p[i] -= md;
  174. k = tmp >> 16;
  175. }
  176. k += borrow;
  177. borrow = u_p[0] < k;
  178. u_p[0] -= k;
  179. if (borrow) {
  180. /* So, this does not happen often; the estimate
  181. was one too big; correct this
  182. */
  183. *lp |= (j & 1) ? (q_est - 1) : ((q_est-1)<<16);
  184. borrow = 0;
  185. for (i = maxv; i > 0; i--) {
  186. unsigned long tmp
  187. = v[i]+(unsigned long)u_p[i]+borrow;
  188. u_p[i] = tmp;
  189. borrow = tmp >> 16;
  190. }
  191. u_p[0] += borrow;
  192. }
  193. else *lp |= (j & 1) ? q_est : (q_est<<16);
  194. }
  195. }
  196. #ifdef EXCEPTION_INEXACT
  197. u_p = &u[0];
  198. for (j = 7; j >= 0; j--) {
  199. if (*u_p++) {
  200. error = 1;
  201. break;
  202. }
  203. }
  204. #endif
  205. #endif
  206. #ifdef EXCEPTION_INEXACT
  207. if (error) {
  208. /*
  209. * report here exception 8.5 - Inexact
  210. * from Draft 8.0 of IEEE P754:
  211. * In the absence of an invalid operation exception,
  212. * if the rounded result of an operation is not exact or if
  213. * it overflows without a trap, then the inexact exception
  214. * shall be assigned. The rounded or overflowed result
  215. * shall be delivered to the destination.
  216. */
  217. INEXACT();
  218. #endif
  219. e1->mantissa = result;
  220. nrm_ext(e1);
  221. if (e1->exp < EXT_MIN) {
  222. /*
  223. * Exception 8.4 - Underflow
  224. */
  225. trap(EFUNFL); /* underflow */
  226. e1->exp = EXT_MIN;
  227. e1->m1 = e1->m2 = 0L;
  228. return;
  229. }
  230. if (e1->exp >= EXT_MAX) {
  231. /*
  232. * Exception 8.3 - Overflow
  233. */
  234. trap(EFOVFL); /* overflow */
  235. e1->exp = EXT_MAX;
  236. e1->m1 = e1->m2 = 0L;
  237. return;
  238. }
  239. }