lib/ack/float/mul_ext.fc

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