/* flonum_mult.c - multiply two flonums
   Copyright (C) 1987-2014 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;
}