Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. /* flonum_mult.c - multiply two flonums
  2.    Copyright 1987, 1990, 1991, 1992, 1995, 2000, 2002, 2003, 2007
  3.    Free Software Foundation, Inc.
  4.  
  5.    This file is part of GAS, the GNU Assembler.
  6.  
  7.    GAS is free software; you can redistribute it and/or modify
  8.    it under the terms of the GNU General Public License as published by
  9.    the Free Software Foundation; either version 3, or (at your option)
  10.    any later version.
  11.  
  12.    GAS is distributed in the hope that it will be useful, but WITHOUT
  13.    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14.    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
  15.    License for more details.
  16.  
  17.    You should have received a copy of the GNU General Public License
  18.    along with GAS; see the file COPYING.  If not, write to the Free
  19.    Software Foundation, 51 Franklin Street - Fifth Floor, Boston, MA
  20.    02110-1301, USA.  */
  21.  
  22. #include "ansidecl.h"
  23. #include "flonum.h"
  24.  
  25. /*      plan for a . b => p(roduct)
  26.  
  27.         +-------+-------+-/   /-+-------+-------+
  28.         | a     | a     |  ...  | a     | a     |
  29.         |  A    |  A-1  |       |  1    |  0    |
  30.         +-------+-------+-/   /-+-------+-------+
  31.  
  32.         +-------+-------+-/   /-+-------+-------+
  33.         | b     | b     |  ...  | b     | b     |
  34.         |  B    |  B-1  |       |  1    |  0    |
  35.         +-------+-------+-/   /-+-------+-------+
  36.  
  37.         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  38.         | p     | p     |  ...  | p     |  ...  | p     | p     |
  39.         |  A+B+1|  A+B  |       |  N    |       |  1    |  0    |
  40.         +-------+-------+-/   /-+-------+-/   /-+-------+-------+
  41.  
  42.         /^\
  43.         (carry) a .b       ...      |      ...   a .b    a .b
  44.         A  B                |             0  1    0  0
  45.         |
  46.         ...         |      ...   a .b
  47.         |                 1  0
  48.         |
  49.         |          ...
  50.         |
  51.         |
  52.         |
  53.         |                 ___
  54.         |                 \
  55.         +-----  P  =   >  a .b
  56.         N         /__  i  j
  57.  
  58.         N = 0 ... A+B
  59.  
  60.         for all i,j where i+j=N
  61.         [i,j integers > 0]
  62.  
  63.         a[], b[], p[] may not intersect.
  64.         Zero length factors signify 0 significant bits: treat as 0.0.
  65.         0.0 factors do the right thing.
  66.         Zero length product OK.
  67.  
  68.         I chose the ForTran accent "foo[bar]" instead of the C accent "*garply"
  69.         because I felt the ForTran way was more intuitive. The C way would
  70.         probably yield better code on most C compilers. Dean Elsner.
  71.         (C style also gives deeper insight [to me] ... oh well ...)  */
  73. void
  74. flonum_multip (const FLONUM_TYPE *a, const FLONUM_TYPE *b,
  75.                FLONUM_TYPE *product)
  76. {
  77.   int size_of_a;                /* 0 origin  */
  78.   int size_of_b;                /* 0 origin  */
  79.   int size_of_product;          /* 0 origin  */
  80.   int size_of_sum;              /* 0 origin  */
  81.   int extra_product_positions;  /* 1 origin  */
  82.   unsigned long work;
  83.   unsigned long carry;
  84.   long exponent;
  85.   LITTLENUM_TYPE *q;
  86.   long significant;             /* TRUE when we emit a non-0 littlenum  */
  87.   /* ForTran accent follows.  */
  88.   int P;                        /* Scan product low-order -> high.  */
  89.   int N;                        /* As in sum above.  */
  90.   int A;                        /* Which [] of a?  */
  91.   int B;                        /* Which [] of b?  */
  92.  
  93.   if ((a->sign != '-' && a->sign != '+')
  94.       || (b->sign != '-' && b->sign != '+'))
  95.     {
  96.       /* Got to fail somehow.  Any suggestions?  */
  97.       product->sign = 0;
  98.       return;
  99.     }
  100.   product->sign = (a->sign == b->sign) ? '+' : '-';
  101.   size_of_a = a->leader - a->low;
  102.   size_of_b = b->leader - b->low;
  103.   exponent = a->exponent + b->exponent;
  104.   size_of_product = product->high - product->low;
  105.   size_of_sum = size_of_a + size_of_b;
  106.   extra_product_positions = size_of_product - size_of_sum;
  107.   if (extra_product_positions < 0)
  108.     {
  109.       P = extra_product_positions;      /* P < 0  */
  110.       exponent -= extra_product_positions;      /* Increases exponent.  */
  111.     }
  112.   else
  113.     {
  114.       P = 0;
  115.     }
  116.   carry = 0;
  117.   significant = 0;
  118.   for (N = 0; N <= size_of_sum; N++)
  119.     {
  120.       work = carry;
  121.       carry = 0;
  122.       for (A = 0; A <= N; A++)
  123.         {
  124.           B = N - A;
  125.           if (A <= size_of_a && B <= size_of_b && B >= 0)
  126.             {
  127. #ifdef TRACE
  128.               printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n",
  129.                       A, a->low[A], B, b->low[B], work);
  130. #endif
  131.               /* Watch out for sign extension!  Without the casts, on
  132.                  the DEC Alpha, the multiplication result is *signed*
  133.                  int, which gets sign-extended to convert to the
  134.                  unsigned long!  */
  135.               work += (unsigned long) a->low[A] * (unsigned long) b->low[B];
  136.               carry += work >> LITTLENUM_NUMBER_OF_BITS;
  137.               work &= LITTLENUM_MASK;
  138. #ifdef TRACE
  139.               printf ("work=%08x carry=%04x\n", work, carry);
  140. #endif
  141.             }
  142.         }
  143.       significant |= work;
  144.       if (significant || P < 0)
  145.         {
  146.           if (P >= 0)
  147.             {
  148.               product->low[P] = work;
  149. #ifdef TRACE
  150.               printf ("P=%d. work[p]:=%04x\n", P, work);
  151. #endif
  152.             }
  153.           P++;
  154.         }
  155.       else
  156.         {
  157.           extra_product_positions++;
  158.           exponent++;
  159.         }
  160.     }
  161.   /* [P]-> position # size_of_sum + 1.
  162.      This is where 'carry' should go.  */
  163. #ifdef TRACE
  164.   printf ("final carry =%04x\n", carry);
  165. #endif
  166.   if (carry)
  167.     {
  168.       if (extra_product_positions > 0)
  169.         product->low[P] = carry;
  170.       else
  171.         {
  172.           /* No room at high order for carry littlenum.  */
  173.           /* Shift right 1 to make room for most significant littlenum.  */
  174.           exponent++;
  175.           P--;
  176.           for (q = product->low + P; q >= product->low; q--)
  177.             {
  178.               work = *q;
  179.               *q = carry;
  180.               carry = work;
  181.             }
  182.         }
  183.     }
  184.   else
  185.     P--;
  186.   product->leader = product->low + P;
  187.   product->exponent = exponent;
  188. }
  189.