gas/flonum-mult.c

   1 /* flonum_multip.c - multiply two flonums
   2    Copyright (C) 1987, 1990, 1991 Free Software Foundation, Inc.
   3
   4    This file is part of Gas, the GNU Assembler.
   5
   6    The GNU assembler is distributed in the hope that it will be
   7    useful, but WITHOUT ANY WARRANTY.  No author or distributor
   8    accepts responsibility to anyone for the consequences of using it
   9    or for whether it serves any particular purpose or works at all,
  10    unless he says so in writing.  Refer to the GNU Assembler General
  11    Public License for full details.
  12
  13    Everyone is granted permission to copy, modify and redistribute
  14    the GNU Assembler, but only under the conditions described in the
  15    GNU Assembler General Public License.  A copy of this license is
  16    supposed to have been given to you along with the GNU Assembler
  17    so you can know your rights and responsibilities.  It should be
  18    in a file named COPYING.  Among other things, the copyright
  19    notice and this notice must be preserved on all copies.  */
  20
  21 #include "flonum.h"
  22
  23 /*      plan for a . b => p(roduct)
  24
  25
  26         +-------+-------+-/   /-+-------+-------+
  27         | a     | a     |  ...  | a     | a     |
  28         |  A    |  A-1  |       |  1    |  0    |
  29         +-------+-------+-/   /-+-------+-------+
  30
  31
  32         +-------+-------+-/   /-+-------+-------+
  33         | b     | b     |  ...  | b     | b     |
  34         |  B    |  B-1  |       |  1    |  0    |
  35         +-------+-------+-/   /-+-------+-------+
  36
  37
  38         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  39         | p     | p     |  ...  | p     |  ...  | p     | p     |
  40         |  A+B+1|  A+B  |       |  N    |       |  1    |  0    |
  41         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  42
  43         /^\
  44         (carry) a .b       ...      |      ...   a .b    a .b
  45         A  B                |             0  1    0  0
  46         |
  47         ...         |      ...   a .b
  48         |                 1  0
  49         |
  50         |          ...
  51         |
  52         |
  53         |
  54         |                 ___
  55         |                 \
  56         +-----  P  =   >  a .b
  57         N         /__  i  j
  58
  59         N = 0 ... A+B
  60
  61         for all i,j where i+j=N
  62         [i,j integers > 0]
  63
  64         a[], b[], p[] may not intersect.
  65         Zero length factors signify 0 significant bits: treat as 0.0.
  66         0.0 factors do the right thing.
  67         Zero length product OK.
  68
  69         I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
  70         because I felt the ForTran way was more intuitive. The C way would
  71         probably yield better code on most C compilers. Dean Elsner.
  72         (C style also gives deeper insight [to me] ... oh well ...)
  73         */
  74 \f
  75 void flonum_multip (a, b, product)
  76 const FLONUM_TYPE *a;
  77 const FLONUM_TYPE *b;
  78 FLONUM_TYPE *product;
  79 {
  80         int                     size_of_a;              /* 0 origin */
  81         int                     size_of_b;              /* 0 origin */
  82         int                     size_of_product;        /* 0 origin */
  83         int                     size_of_sum;            /* 0 origin */
  84         int                     extra_product_positions;/* 1 origin */
  85         unsigned long   work;
  86         unsigned long   carry;
  87         long            exponent;
  88         LITTLENUM_TYPE *        q;
  89         long            significant;            /* TRUE when we emit a non-0 littlenum  */
  90         /* ForTran accent follows. */
  91         int                     P;      /* Scan product low-order -> high. */
  92         int                     N;      /* As in sum above.  */
  93         int                     A;      /* Which [] of a? */
  94         int                     B;      /* Which [] of b? */
  95
  96         if((a->sign!='-' && a->sign!='+') || (b->sign!='-' && b->sign!='+')) {
  97                 /* ...
  98                    Got to fail somehow.  Any suggestions? */
  99                 product->sign=0;
 100                 return;
 101         }
 102         product -> sign = (a->sign == b->sign) ? '+' : '-';
 103         size_of_a               = a       -> leader     -  a       -> low;
 104         size_of_b               = b       -> leader     -  b       -> low;
 105         exponent                = a       -> exponent   +  b       -> exponent;
 106         size_of_product = product -> high       -  product -> low;
 107         size_of_sum             = size_of_a             +  size_of_b;
 108         extra_product_positions  =  size_of_product  -  size_of_sum;
 109         if (extra_product_positions < 0)
 110             {
 111                     P = extra_product_positions; /* P < 0 */
 112                     exponent -= extra_product_positions; /* Increases exponent. */
 113             }
 114         else
 115             {
 116                     P = 0;
 117             }
 118         carry = 0;
 119         significant = 0;
 120         for (N = 0;
 121              N <= size_of_sum;
 122              N++)
 123             {
 124                     work = carry;
 125                     carry = 0;
 126                     for (A = 0;
 127                          A <= N;
 128                          A ++)
 129                         {
 130                                 B = N - A;
 131                                 if (A <= size_of_a   &&   B <= size_of_b  &&  B >= 0)
 132                                     {
 133 #ifdef TRACE
 134                                             printf("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n", A, a->low[A], B, b->low[B], work);
 135 #endif
 136                                             work += a -> low [A]   *   b -> low [B];
 137                                             carry += work >> LITTLENUM_NUMBER_OF_BITS;
 138                                             work &= LITTLENUM_MASK;
 139 #ifdef TRACE
 140                                             printf("work=%08x carry=%04x\n", work, carry);
 141 #endif
 142                                     }
 143                         }
 144                     significant |= work;
 145                     if (significant || P<0)
 146                         {
 147                                 if (P >= 0)
 148                                     {
 149                                             product -> low [P] = work;
 150 #ifdef TRACE
 151                                             printf("P=%d. work[p]:=%04x\n", P, work);
 152 #endif
 153                                     }
 154                                 P ++;
 155                         }
 156                     else
 157                         {
 158                                 extra_product_positions ++;
 159                                 exponent ++;
 160                         }
 161             }
 162         /*
 163          * [P]-> position # size_of_sum + 1.
 164          * This is where 'carry' should go.
 165          */
 166 #ifdef TRACE
 167         printf("final carry =%04x\n", carry);
 168 #endif
 169         if (carry)
 170             {
 171                     if (extra_product_positions > 0)
 172                         {
 173                                 product -> low [P] = carry;
 174                         }
 175                     else
 176                         {
 177                                 /* No room at high order for carry littlenum. */
 178                                 /* Shift right 1 to make room for most significant littlenum. */
 179                                 exponent ++;
 180                                 P --;
 181                                 for (q  = product -> low + P;
 182                                      q >= product -> low;
 183                                      q --)
 184                                     {
 185                                             work = * q;
 186                                             * q = carry;
 187                                             carry = work;
 188                                     }
 189                         }
 190             }
 191         else
 192             {
 193                     P --;
 194             }
 195         product -> leader       = product -> low + P;
 196         product -> exponent     = exponent;
 197 }
 198
 199 /* end: flonum_multip.c */