flt_div.c 2.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129
  1. /*
  2. (c) copyright 1989 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. #include "flt_misc.h"
  7. void
  8. flt_div(e1,e2,e3)
  9. register flt_arith *e1,*e2,*e3;
  10. {
  11. long result[2];
  12. register long *lp;
  13. unsigned short u[9], v[5];
  14. register int j;
  15. register unsigned short *u_p = u;
  16. int maxv = 4;
  17. flt_arith cpe1 = *e1;
  18. flt_status = 0;
  19. e3->flt_sign = e1->flt_sign ^ e2->flt_sign;
  20. if ((e2->m1 | e2->m2) == 0) {
  21. flt_status = FLT_DIV0;
  22. e3->flt_exp = EXT_MAX;
  23. e3->m1 = 0x80000000;
  24. e3->m2 = 0L;
  25. return;
  26. }
  27. if ((e1->m1 | e1->m2) == 0) { /* 0 / anything == 0 */
  28. e3->flt_exp = 0; /* make sure */
  29. e3->m1 = e3->m2 = 0L;
  30. e3->flt_sign = 0;
  31. return;
  32. }
  33. e3->flt_exp = e1->flt_exp - e2->flt_exp;
  34. u[4] = (e1->m2 & 1) << 15;
  35. flt_b64_sft(&(cpe1.flt_mantissa),1);
  36. flt_split(&cpe1, u);
  37. u[5] = 0; u[6] = 0; u[7] = 0;
  38. flt_split(e2, &v[1]);
  39. while (! v[maxv]) maxv--;
  40. result[0] = 0;
  41. result[1] = 0;
  42. lp = result;
  43. /*
  44. * Use an algorithm of Knuth (The art of programming, Seminumerical
  45. * algorithms), to divide u by v. u and v are both seen as numbers
  46. * with base 65536.
  47. */
  48. for (j = 0; j <= 3; j++, u_p++) {
  49. long q_est, temp;
  50. long v1 = v[1];
  51. if (j == 2) lp++;
  52. if (u_p[0] == 0 && u_p[1] < v[1]) continue;
  53. temp = ((long)u_p[0] << 16) + u_p[1];
  54. if (u_p[0] >= v[1]) {
  55. q_est = 0x0000FFFFL;
  56. }
  57. else if (v[1] == 1) {
  58. q_est = temp;
  59. }
  60. else if (temp >= 0) {
  61. q_est = temp / v1;
  62. }
  63. else {
  64. long rem;
  65. q_est = (0x7FFFFFFF/v1)+((temp&0x7FFFFFFF)/v1);
  66. rem = (0x7FFFFFFF%v1)+((temp&0x7FFFFFFF)%v1)+1;
  67. while (rem >= v1) {
  68. q_est++;
  69. rem -= v1;
  70. }
  71. }
  72. temp -= q_est * v1;
  73. while (!(temp&0xFFFF0000) &&
  74. ucmp((long)v[2]*q_est,(temp<<16)+(long)u_p[2]) > 0) {
  75. q_est--;
  76. temp += v1;
  77. }
  78. /* Now, according to Knuth, we have an estimate of the
  79. quotient, that is either correct or one too big, but
  80. almost always correct.
  81. */
  82. if (q_est != 0) {
  83. int i;
  84. long k = 0;
  85. int borrow = 0;
  86. for (i = maxv; i > 0; i--) {
  87. long tmp = q_est * (long)v[i] + k + borrow;
  88. unsigned short md = tmp & 0xFFFF;
  89. borrow = (md > u_p[i]);
  90. u_p[i] -= md;
  91. k = (tmp >> 16) & 0xFFFF;
  92. }
  93. k += borrow;
  94. borrow = (long)u_p[0] < k;
  95. u_p[0] -= k;
  96. if (borrow) {
  97. /* So, this does not happen often; the estimate
  98. was one too big; correct this
  99. */
  100. q_est--;
  101. borrow = 0;
  102. for (i = maxv; i > 0; i--) {
  103. long tmp
  104. = v[i]+(long)u_p[i]+borrow;
  105. u_p[i] = tmp & 0xFFFF;
  106. borrow = (tmp >> 16) & 0xFFFF;
  107. }
  108. u_p[0] += borrow;
  109. }
  110. *lp |= (j & 1) ? q_est : (q_est<<16);
  111. }
  112. }
  113. e3->m1 = result[0];
  114. e3->m2 = result[1];
  115. flt_nrm(e3);
  116. flt_chk(e3);
  117. }