modules/src/flt_arith/flt_mul.c

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