Back to index

cell-binutils  2.17cvs20070401
flonum-mult.c
Go to the documentation of this file.
00001 /* flonum_mult.c - multiply two flonums
00002    Copyright 1987, 1990, 1991, 1992, 1995, 2000, 2002, 2003
00003    Free Software Foundation, Inc.
00004 
00005    This file is part of Gas, the GNU Assembler.
00006 
00007    The GNU assembler is distributed in the hope that it will be
00008    useful, but WITHOUT ANY WARRANTY.  No author or distributor
00009    accepts responsibility to anyone for the consequences of using it
00010    or for whether it serves any particular purpose or works at all,
00011    unless he says so in writing.  Refer to the GNU Assembler General
00012    Public License for full details.
00013 
00014    Everyone is granted permission to copy, modify and redistribute
00015    the GNU Assembler, but only under the conditions described in the
00016    GNU Assembler General Public License.  A copy of this license is
00017    supposed to have been given to you along with the GNU Assembler
00018    so you can know your rights and responsibilities.  It should be
00019    in a file named COPYING.  Among other things, the copyright
00020    notice and this notice must be preserved on all copies.  */
00021 
00022 #include "ansidecl.h"
00023 #include "flonum.h"
00024 
00025 /*     plan for a . b => p(roduct)
00026 
00027        +-------+-------+-/   /-+-------+-------+
00028        | a    | a    |  ... | a    | a    |
00029        |  A   |  A-1 |      |  1   |  0   |
00030        +-------+-------+-/   /-+-------+-------+
00031 
00032        +-------+-------+-/   /-+-------+-------+
00033        | b    | b    |  ... | b    | b    |
00034        |  B   |  B-1 |      |  1   |  0   |
00035        +-------+-------+-/   /-+-------+-------+
00036 
00037        +-------+-------+-/   /-+-------+-/   /-+-------+-------+
00038        | p    | p    |  ... | p    |  ... | p    | p    |
00039        |  A+B+1|  A+B       |      |  N   |      |  1   |  0   |
00040        +-------+-------+-/   /-+-------+-/   /-+-------+-------+
00041 
00042        /^\
00043        (carry) a .b     ...     |     ...  a .b   a .b
00044        A  B              |           0  1   0  0
00045        |
00046        ...        |     ...  a .b
00047        |               1  0
00048        |
00049        |         ...
00050        |
00051        |
00052        |
00053        |               ___
00054        |               \
00055        +-----  P  =   >  a .b
00056        N        /__  i  j
00057 
00058        N = 0 ... A+B
00059 
00060        for all i,j where i+j=N
00061        [i,j integers > 0]
00062 
00063        a[], b[], p[] may not intersect.
00064        Zero length factors signify 0 significant bits: treat as 0.0.
00065        0.0 factors do the right thing.
00066        Zero length product OK.
00067 
00068        I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
00069        because I felt the ForTran way was more intuitive. The C way would
00070        probably yield better code on most C compilers. Dean Elsner.
00071        (C style also gives deeper insight [to me] ... oh well ...)  */
00072 
00073 void
00074 flonum_multip (const FLONUM_TYPE *a, const FLONUM_TYPE *b,
00075               FLONUM_TYPE *product)
00076 {
00077   int size_of_a;            /* 0 origin  */
00078   int size_of_b;            /* 0 origin  */
00079   int size_of_product;             /* 0 origin  */
00080   int size_of_sum;          /* 0 origin  */
00081   int extra_product_positions;     /* 1 origin  */
00082   unsigned long work;
00083   unsigned long carry;
00084   long exponent;
00085   LITTLENUM_TYPE *q;
00086   long significant;         /* TRUE when we emit a non-0 littlenum  */
00087   /* ForTran accent follows.  */
00088   int P;                    /* Scan product low-order -> high.  */
00089   int N;                    /* As in sum above.  */
00090   int A;                    /* Which [] of a?  */
00091   int B;                    /* Which [] of b?  */
00092 
00093   if ((a->sign != '-' && a->sign != '+')
00094       || (b->sign != '-' && b->sign != '+'))
00095     {
00096       /* Got to fail somehow.  Any suggestions?  */
00097       product->sign = 0;
00098       return;
00099     }
00100   product->sign = (a->sign == b->sign) ? '+' : '-';
00101   size_of_a = a->leader - a->low;
00102   size_of_b = b->leader - b->low;
00103   exponent = a->exponent + b->exponent;
00104   size_of_product = product->high - product->low;
00105   size_of_sum = size_of_a + size_of_b;
00106   extra_product_positions = size_of_product - size_of_sum;
00107   if (extra_product_positions < 0)
00108     {
00109       P = extra_product_positions; /* P < 0  */
00110       exponent -= extra_product_positions;       /* Increases exponent.  */
00111     }
00112   else
00113     {
00114       P = 0;
00115     }
00116   carry = 0;
00117   significant = 0;
00118   for (N = 0; N <= size_of_sum; N++)
00119     {
00120       work = carry;
00121       carry = 0;
00122       for (A = 0; A <= N; A++)
00123        {
00124          B = N - A;
00125          if (A <= size_of_a && B <= size_of_b && B >= 0)
00126            {
00127 #ifdef TRACE
00128              printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n",
00129                     A, a->low[A], B, b->low[B], work);
00130 #endif
00131              /* Watch out for sign extension!  Without the casts, on
00132                the DEC Alpha, the multiplication result is *signed*
00133                int, which gets sign-extended to convert to the
00134                unsigned long!  */
00135              work += (unsigned long) a->low[A] * (unsigned long) b->low[B];
00136              carry += work >> LITTLENUM_NUMBER_OF_BITS;
00137              work &= LITTLENUM_MASK;
00138 #ifdef TRACE
00139              printf ("work=%08x carry=%04x\n", work, carry);
00140 #endif
00141            }
00142        }
00143       significant |= work;
00144       if (significant || P < 0)
00145        {
00146          if (P >= 0)
00147            {
00148              product->low[P] = work;
00149 #ifdef TRACE
00150              printf ("P=%d. work[p]:=%04x\n", P, work);
00151 #endif
00152            }
00153          P++;
00154        }
00155       else
00156        {
00157          extra_product_positions++;
00158          exponent++;
00159        }
00160     }
00161   /* [P]-> position # size_of_sum + 1.
00162      This is where 'carry' should go.  */
00163 #ifdef TRACE
00164   printf ("final carry =%04x\n", carry);
00165 #endif
00166   if (carry)
00167     {
00168       if (extra_product_positions > 0)
00169        product->low[P] = carry;
00170       else
00171        {
00172          /* No room at high order for carry littlenum.  */
00173          /* Shift right 1 to make room for most significant littlenum.  */
00174          exponent++;
00175          P--;
00176          for (q = product->low + P; q >= product->low; q--)
00177            {
00178              work = *q;
00179              *q = carry;
00180              carry = work;
00181            }
00182        }
00183     }
00184   else
00185     P--;
00186   product->leader = product->low + P;
00187   product->exponent = exponent;
00188 }