flt_mul.c 1.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879
  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. #include "flt_misc.h"
  7. void
  8. flt_mul(e1,e2,e3)
  9. register flt_arith *e1,*e2,*e3;
  10. {
  11. /* Multiply the extended numbers e1 and e2, and put the
  12. result in e3.
  13. */
  14. register int i,j; /* loop control */
  15. unsigned short mp[4];
  16. unsigned short mc[4];
  17. unsigned short result[8]; /* result */
  18. register unsigned short *pres;
  19. flt_status = 0;
  20. /* first save the sign (XOR) */
  21. e3->flt_sign = e1->flt_sign ^ e2->flt_sign;
  22. /* compute new exponent */
  23. e3->flt_exp = e1->flt_exp + e2->flt_exp + 1;
  24. /* 128 bit multiply of mantissas */
  25. /* assign unknown long formats */
  26. /* to known unsigned word formats */
  27. flt_split(e1, mp);
  28. flt_split(e2, mc);
  29. for (i = 8; i--;) {
  30. result[i] = 0;
  31. }
  32. /*
  33. * fill registers with their components
  34. */
  35. for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
  36. unsigned short k = 0;
  37. long mpi = mp[i];
  38. for(j=4;j--;) {
  39. long tmp = (long)pres[j] + k;
  40. if (mc[j]) tmp += mpi * mc[j];
  41. pres[j] = tmp & 0xFFFF;
  42. k = (tmp >> 16) & 0xFFFF;
  43. }
  44. pres[-1] = k;
  45. }
  46. if (! (result[0] & 0x8000)) {
  47. e3->flt_exp--;
  48. for (i = 0; i <= 3; i++) {
  49. result[i] <<= 1;
  50. if (result[i+1]&0x8000) result[i] |= 1;
  51. }
  52. result[4] <<= 1;
  53. }
  54. /*
  55. * combine the registers to a total
  56. */
  57. e3->m1 = ((long)result[0] << 16) + result[1];
  58. e3->m2 = ((long)result[2] << 16) + result[3];
  59. if (result[4] & 0x8000) {
  60. if (++e3->m2 == 0 || (e3->m2 & ~ 0xFFFFFFFF)) {
  61. e3->m2 = 0;
  62. if (++e3->m1 == 0 || (e3->m1 & ~ 0xFFFFFFFF)) {
  63. e3->m1 = 0x80000000;
  64. e3->flt_exp++;
  65. }
  66. }
  67. }
  68. flt_chk(e3);
  69. }