mul_ext.c 2.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798
  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. ROUTINE TO MULTIPLY TWO EXTENDED FORMAT NUMBERS
  8. */
  9. # include "FP_bias.h"
  10. # include "FP_trap.h"
  11. # include "FP_types.h"
  12. # include "FP_shift.h"
  13. void
  14. mul_ext(e1,e2)
  15. EXTEND *e1,*e2;
  16. {
  17. register int i,j; /* loop control */
  18. unsigned short mp[4]; /* multiplier */
  19. unsigned short mc[4]; /* multipcand */
  20. unsigned short result[8]; /* result */
  21. register unsigned short *pres;
  22. /* first save the sign (XOR) */
  23. e1->sign ^= e2->sign;
  24. /* compute new exponent */
  25. e1->exp += e2->exp + 1;
  26. /* 128 bit multiply of mantissas */
  27. /* assign unknown long formats */
  28. /* to known unsigned word formats */
  29. mp[0] = e1->m1 >> 16;
  30. mp[1] = (unsigned short) e1->m1;
  31. mp[2] = e1->m2 >> 16;
  32. mp[3] = (unsigned short) e1->m2;
  33. mc[0] = e2->m1 >> 16;
  34. mc[1] = (unsigned short) e2->m1;
  35. mc[2] = e2->m2 >> 16;
  36. mc[3] = (unsigned short) e2->m2;
  37. for (i = 8; i--;) {
  38. result[i] = 0;
  39. }
  40. /*
  41. * fill registers with their components
  42. */
  43. for(i=4, pres = &result[4];i--;pres--) if (mp[i]) {
  44. unsigned short k = 0;
  45. unsigned long mpi = mp[i];
  46. for(j=4;j--;) {
  47. unsigned long tmp = (unsigned long)pres[j] + k;
  48. if (mc[j]) tmp += mpi * mc[j];
  49. pres[j] = tmp;
  50. k = tmp >> 16;
  51. }
  52. pres[-1] = k;
  53. }
  54. if (! (result[0] & 0x8000)) {
  55. e1->exp--;
  56. for (i = 0; i <= 3; i++) {
  57. result[i] <<= 1;
  58. if (result[i+1]&0x8000) result[i] |= 1;
  59. }
  60. result[4] <<= 1;
  61. }
  62. /*
  63. * combine the registers to a total
  64. */
  65. e1->m1 = ((unsigned long)(result[0]) << 16) + result[1];
  66. e1->m2 = ((unsigned long)(result[2]) << 16) + result[3];
  67. if (result[4] & 0x8000) {
  68. if (++e1->m2 == 0)
  69. if (++e1->m1 == 0) {
  70. e1->m1 = NORMBIT;
  71. e1->exp++;
  72. }
  73. }
  74. /* check for overflow */
  75. if (e1->exp >= EXT_MAX) {
  76. trap(EFOVFL);
  77. /* if caught */
  78. /* return signed infinity */
  79. e1->exp = EXT_MAX;
  80. infinity: e1->m1 = e1->m2 =0L;
  81. return;
  82. }
  83. /* check for underflow */
  84. if (e1->exp < EXT_MIN) {
  85. trap(EFUNFL);
  86. e1->exp = EXT_MIN;
  87. goto infinity;
  88. }
  89. }