Subversion Repositories openrisc_2011-10-31

/* flonum_mult.c - multiply two flonums

   Copyright 1987, 1990, 1991, 1992, 1995, 2000, 2002, 2003, 2007

   Free Software Foundation, Inc.

   This file is part of GAS, the GNU Assembler.

   GAS is free software; you can redistribute it and/or modify

   it under the terms of the GNU General Public License as published by

   the Free Software Foundation; either version 3, or (at your option)

   any later version.

   GAS is distributed in the hope that it will be useful, but WITHOUT

   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public

   License for more details.

   You should have received a copy of the GNU General Public License

   along with GAS; see the file COPYING.  If not, write to the Free

   Software Foundation, 51 Franklin Street - Fifth Floor, Boston, MA

   02110-1301, USA.  */

#include "ansidecl.h"

#include "flonum.h"

/*      plan for a . b => p(roduct)

        +-------+-------+-/   /-+-------+-------+

        | a     | a     |  ...  | a     | a     |

        |  A    |  A-1  |       |  1    |  0    |

        +-------+-------+-/   /-+-------+-------+

        +-------+-------+-/   /-+-------+-------+

        | b     | b     |  ...  | b     | b     |

        |  B    |  B-1  |       |  1    |  0    |

        +-------+-------+-/   /-+-------+-------+

        +-------+-------+-/   /-+-------+-/   /-+-------+-------+

        | p     | p     |  ...  | p     |  ...  | p     | p     |

        |  A+B+1|  A+B  |       |  N    |       |  1    |  0    |

        +-------+-------+-/   /-+-------+-/   /-+-------+-------+

        /^\

        (carry) a .b       ...      |      ...   a .b    a .b

        A  B                |             0  1    0  0

        |

        ...         |      ...   a .b

        |                 1  0

        |

        |          ...

        |

        |

        |

        |                 ___

        |                 \

        +-----  P  =   >  a .b

        N         /__  i  j

        N = 0 ... A+B

        for all i,j where i+j=N

        [i,j integers > 0]

        a[], b[], p[] may not intersect.

        Zero length factors signify 0 significant bits: treat as 0.0.

        0.0 factors do the right thing.

        Zero length product OK.

        I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"

        because I felt the ForTran way was more intuitive. The C way would

        probably yield better code on most C compilers. Dean Elsner.

        (C style also gives deeper insight [to me] ... oh well ...)  */

void

flonum_multip (const FLONUM_TYPE *a, const FLONUM_TYPE *b,

               FLONUM_TYPE *product)

{

  int size_of_a;                /* 0 origin  */

  int size_of_b;                /* 0 origin  */

  int size_of_product;          /* 0 origin  */

  int size_of_sum;              /* 0 origin  */

  int extra_product_positions;  /* 1 origin  */

  unsigned long work;

  unsigned long carry;

  long exponent;

  LITTLENUM_TYPE *q;

  long significant;             /* TRUE when we emit a non-0 littlenum  */

  /* ForTran accent follows.  */

  int P;                        /* Scan product low-order -> high.  */

  int N;                        /* As in sum above.  */

  int A;                        /* Which [] of a?  */

  int B;                        /* Which [] of b?  */

  if ((a->sign != '-' && a->sign != '+')

      || (b->sign != '-' && b->sign != '+'))

    {

      /* Got to fail somehow.  Any suggestions?  */

      product->sign = 0;

      return;

    }

  product->sign = (a->sign == b->sign) ? '+' : '-';

  size_of_a = a->leader - a->low;

  size_of_b = b->leader - b->low;

  exponent = a->exponent + b->exponent;

  size_of_product = product->high - product->low;

  size_of_sum = size_of_a + size_of_b;

  extra_product_positions = size_of_product - size_of_sum;

  if (extra_product_positions < 0)

    {

      P = extra_product_positions;      /* P < 0  */

      exponent -= extra_product_positions;      /* Increases exponent.  */

    }

  else

    {

      P = 0;

    }

  carry = 0;

  significant = 0;

  for (N = 0; N <= size_of_sum; N++)

    {

      work = carry;

      carry = 0;

      for (A = 0; A <= N; A++)

        {

          B = N - A;

          if (A <= size_of_a && B <= size_of_b && B >= 0)

            {

#ifdef TRACE

              printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n",

                      A, a->low[A], B, b->low[B], work);

#endif

              /* Watch out for sign extension!  Without the casts, on

                 the DEC Alpha, the multiplication result is *signed*

                 int, which gets sign-extended to convert to the

                 unsigned long!  */

              work += (unsigned long) a->low[A] * (unsigned long) b->low[B];

              carry += work >> LITTLENUM_NUMBER_OF_BITS;

              work &= LITTLENUM_MASK;

#ifdef TRACE

              printf ("work=%08x carry=%04x\n", work, carry);

#endif

            }

        }

      significant |= work;

      if (significant || P < 0)

        {

          if (P >= 0)

            {

              product->low[P] = work;

#ifdef TRACE

              printf ("P=%d. work[p]:=%04x\n", P, work);

#endif

            }

          P++;

        }

      else

        {

          extra_product_positions++;

          exponent++;

        }

    }

  /* [P]-> position # size_of_sum + 1.

     This is where 'carry' should go.  */

#ifdef TRACE

  printf ("final carry =%04x\n", carry);

#endif

  if (carry)

    {

      if (extra_product_positions > 0)

        product->low[P] = carry;

      else

        {

          /* No room at high order for carry littlenum.  */

          /* Shift right 1 to make room for most significant littlenum.  */

          exponent++;

          P--;

          for (q = product->low + P; q >= product->low; q--)

            {

              work = *q;

              *q = carry;

              carry = work;

            }

        }

    }

  else

    P--;

  product->leader = product->low + P;

  product->exponent = exponent;

}