gnu/usr.bin/as/flonum-mult.c

   1 /* flonum_mult.c - multiply two flonums
   2    Copyright (C) 1987, 1990, 1991, 1992 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 #ifndef lint
  22 static char rcsid[] = "$FreeBSD: src/gnu/usr.bin/as/flonum-mult.c,v 1.6 1999/08/27 23:34:15 peter Exp $";
  23 #endif
  24
  25 #include "flonum.h"
  26
  27 /*      plan for a . b => p(roduct)
  28
  29
  30         +-------+-------+-/   /-+-------+-------+
  31         | a     | a     |  ...  | a     | a     |
  32         |  A    |  A-1  |       |  1    |  0    |
  33         +-------+-------+-/   /-+-------+-------+
  34
  35
  36         +-------+-------+-/   /-+-------+-------+
  37         | b     | b     |  ...  | b     | b     |
  38         |  B    |  B-1  |       |  1    |  0    |
  39         +-------+-------+-/   /-+-------+-------+
  40
  41
  42         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  43         | p     | p     |  ...  | p     |  ...  | p     | p     |
  44         |  A+B+1|  A+B  |       |  N    |       |  1    |  0    |
  45         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  46
  47         /^\
  48         (carry) a .b       ...      |      ...   a .b    a .b
  49         A  B                |             0  1    0  0
  50         |
  51         ...         |      ...   a .b
  52         |                 1  0
  53         |
  54         |          ...
  55         |
  56         |
  57         |
  58         |                 ___
  59         |                 \
  60         +-----  P  =   >  a .b
  61         N         /__  i  j
  62
  63         N = 0 ... A+B
  64
  65         for all i,j where i+j=N
  66         [i,j integers > 0]
  67
  68         a[], b[], p[] may not intersect.
  69         Zero length factors signify 0 significant bits: treat as 0.0.
  70         0.0 factors do the right thing.
  71         Zero length product OK.
  72
  73         I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
  74         because I felt the ForTran way was more intuitive. The C way would
  75         probably yield better code on most C compilers. Dean Elsner.
  76         (C style also gives deeper insight [to me] ... oh well ...)
  77         */
  78 \f
  79 void flonum_multip (a, b, product)
  80 const FLONUM_TYPE *a;
  81 const FLONUM_TYPE *b;
  82 FLONUM_TYPE *product;
  83 {
  84         int                     size_of_a;              /* 0 origin */
  85         int                     size_of_b;              /* 0 origin */
  86         int                     size_of_product;        /* 0 origin */
  87         int                     size_of_sum;            /* 0 origin */
  88         int                     extra_product_positions;/* 1 origin */
  89         unsigned long   work;
  90         unsigned long   carry;
  91         long            exponent;
  92         LITTLENUM_TYPE *        q;
  93         long            significant;            /* TRUE when we emit a non-0 littlenum  */
  94         /* ForTran accent follows. */
  95         int                     P;      /* Scan product low-order -> high. */
  96         int                     N;      /* As in sum above.  */
  97         int                     A;      /* Which [] of a? */
  98         int                     B;      /* Which [] of b? */
  99
 100         if ((a->sign != '-' && a->sign != '+') || (b->sign != '-' && b->sign != '+')) {
 101                 /* ...
 102                    Got to fail somehow.  Any suggestions? */
 103                 product->sign=0;
 104                 return;
 105         }
 106         product->sign = (a->sign == b->sign) ? '+' : '-';
 107         size_of_a               = a->leader     -  a->low;
 108         size_of_b               = b->leader     -  b->low;
 109         exponent                = a->exponent   +  b->exponent;
 110         size_of_product = product->high -  product->low;
 111         size_of_sum             = size_of_a             +  size_of_b;
 112         extra_product_positions  =  size_of_product  -  size_of_sum;
 113         if (extra_product_positions < 0)
 114             {
 115                     P = extra_product_positions; /* P < 0 */
 116                     exponent -= extra_product_positions; /* Increases exponent. */
 117             }
 118         else
 119             {
 120                     P = 0;
 121             }
 122         carry = 0;
 123         significant = 0;
 124         for (N = 0;
 125              N <= size_of_sum;
 126              N++)
 127             {
 128                     work = carry;
 129                     carry = 0;
 130                     for (A = 0;
 131                          A <= N;
 132                          A ++)
 133                         {
 134                                 B = N - A;
 135                                 if (A <= size_of_a   &&   B <= size_of_b  &&  B >= 0)
 136                                     {
 137 #ifdef TRACE
 138                                             printf("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n", A, a->low[A], B, b->low[B], work);
 139 #endif
 140                                             work += a->low[A]   *   b->low[B];
 141                                             carry += work >> LITTLENUM_NUMBER_OF_BITS;
 142                                             work &= LITTLENUM_MASK;
 143 #ifdef TRACE
 144                                             printf("work=%08x carry=%04x\n", work, carry);
 145 #endif
 146                                     }
 147                         }
 148                     significant |= work;
 149                     if (significant || P<0)
 150                         {
 151                                 if (P >= 0)
 152                                     {
 153                                             product->low[P] = work;
 154 #ifdef TRACE
 155                                             printf("P=%d. work[p]:=%04x\n", P, work);
 156 #endif
 157                                     }
 158                                 P ++;
 159                         }
 160                     else
 161                         {
 162                                 extra_product_positions ++;
 163                                 exponent ++;
 164                         }
 165             }
 166         /*
 167          * [P]->position # size_of_sum + 1.
 168          * This is where 'carry' should go.
 169          */
 170 #ifdef TRACE
 171         printf("final carry =%04x\n", carry);
 172 #endif
 173         if (carry)
 174             {
 175                     if (extra_product_positions > 0)
 176                         {
 177                                 product->low[P] = carry;
 178                         }
 179                     else
 180                         {
 181                                 /* No room at high order for carry littlenum. */
 182                                 /* Shift right 1 to make room for most significant littlenum. */
 183                                 exponent ++;
 184                                 P --;
 185                                 for (q  = product->low + P;
 186                                      q >= product->low;
 187                                      q --)
 188                                     {
 189                                             work = * q;
 190                                             * q = carry;
 191                                             carry = work;
 192                                     }
 193                         }
 194             }
 195         else
 196             {
 197                     P --;
 198             }
 199         product->leader = product->low + P;
 200         product->exponent       = exponent;
 201 }
 202
 203 /* end of flonum_mult.c */