133965Sjdp/* flonum_mult.c - multiply two flonums 2218822Sdim Copyright 1987, 1990, 1991, 1992, 1995, 2000, 2002, 2003 377298Sobrien Free Software Foundation, Inc. 433965Sjdp 533965Sjdp This file is part of Gas, the GNU Assembler. 633965Sjdp 733965Sjdp The GNU assembler is distributed in the hope that it will be 833965Sjdp useful, but WITHOUT ANY WARRANTY. No author or distributor 933965Sjdp accepts responsibility to anyone for the consequences of using it 1033965Sjdp or for whether it serves any particular purpose or works at all, 1133965Sjdp unless he says so in writing. Refer to the GNU Assembler General 1233965Sjdp Public License for full details. 1333965Sjdp 1433965Sjdp Everyone is granted permission to copy, modify and redistribute 1533965Sjdp the GNU Assembler, but only under the conditions described in the 1633965Sjdp GNU Assembler General Public License. A copy of this license is 1733965Sjdp supposed to have been given to you along with the GNU Assembler 1833965Sjdp so you can know your rights and responsibilities. It should be 1933965Sjdp in a file named COPYING. Among other things, the copyright 2033965Sjdp notice and this notice must be preserved on all copies. */ 2133965Sjdp 22104834Sobrien#include "ansidecl.h" 2333965Sjdp#include "flonum.h" 2433965Sjdp 2533965Sjdp/* plan for a . b => p(roduct) 2677298Sobrien 2733965Sjdp +-------+-------+-/ /-+-------+-------+ 2833965Sjdp | a | a | ... | a | a | 2933965Sjdp | A | A-1 | | 1 | 0 | 3033965Sjdp +-------+-------+-/ /-+-------+-------+ 3177298Sobrien 3233965Sjdp +-------+-------+-/ /-+-------+-------+ 3333965Sjdp | b | b | ... | b | b | 3433965Sjdp | B | B-1 | | 1 | 0 | 3533965Sjdp +-------+-------+-/ /-+-------+-------+ 3677298Sobrien 3733965Sjdp +-------+-------+-/ /-+-------+-/ /-+-------+-------+ 3833965Sjdp | p | p | ... | p | ... | p | p | 3933965Sjdp | A+B+1| A+B | | N | | 1 | 0 | 4033965Sjdp +-------+-------+-/ /-+-------+-/ /-+-------+-------+ 4177298Sobrien 4233965Sjdp /^\ 4333965Sjdp (carry) a .b ... | ... a .b a .b 4433965Sjdp A B | 0 1 0 0 4533965Sjdp | 4633965Sjdp ... | ... a .b 4733965Sjdp | 1 0 4833965Sjdp | 4933965Sjdp | ... 5033965Sjdp | 5133965Sjdp | 5233965Sjdp | 5333965Sjdp | ___ 5433965Sjdp | \ 5533965Sjdp +----- P = > a .b 5633965Sjdp N /__ i j 5777298Sobrien 5833965Sjdp N = 0 ... A+B 5977298Sobrien 6033965Sjdp for all i,j where i+j=N 6133965Sjdp [i,j integers > 0] 6277298Sobrien 6333965Sjdp a[], b[], p[] may not intersect. 6433965Sjdp Zero length factors signify 0 significant bits: treat as 0.0. 6533965Sjdp 0.0 factors do the right thing. 6633965Sjdp Zero length product OK. 6777298Sobrien 6833965Sjdp I chose the ForTran accent "foo[bar]" instead of the C accent "*garply" 6933965Sjdp because I felt the ForTran way was more intuitive. The C way would 7033965Sjdp probably yield better code on most C compilers. Dean Elsner. 7177298Sobrien (C style also gives deeper insight [to me] ... oh well ...) */ 7233965Sjdp 7377298Sobrienvoid 74130561Sobrienflonum_multip (const FLONUM_TYPE *a, const FLONUM_TYPE *b, 75130561Sobrien FLONUM_TYPE *product) 7633965Sjdp{ 7777298Sobrien int size_of_a; /* 0 origin */ 7877298Sobrien int size_of_b; /* 0 origin */ 7977298Sobrien int size_of_product; /* 0 origin */ 8077298Sobrien int size_of_sum; /* 0 origin */ 8177298Sobrien int extra_product_positions; /* 1 origin */ 8233965Sjdp unsigned long work; 8333965Sjdp unsigned long carry; 8433965Sjdp long exponent; 8533965Sjdp LITTLENUM_TYPE *q; 8633965Sjdp long significant; /* TRUE when we emit a non-0 littlenum */ 8777298Sobrien /* ForTran accent follows. */ 8877298Sobrien int P; /* Scan product low-order -> high. */ 8933965Sjdp int N; /* As in sum above. */ 9077298Sobrien int A; /* Which [] of a? */ 9177298Sobrien int B; /* Which [] of b? */ 9233965Sjdp 9377298Sobrien if ((a->sign != '-' && a->sign != '+') 9477298Sobrien || (b->sign != '-' && b->sign != '+')) 9533965Sjdp { 9677298Sobrien /* Got to fail somehow. Any suggestions? */ 9733965Sjdp product->sign = 0; 9833965Sjdp return; 9933965Sjdp } 10033965Sjdp product->sign = (a->sign == b->sign) ? '+' : '-'; 10133965Sjdp size_of_a = a->leader - a->low; 10233965Sjdp size_of_b = b->leader - b->low; 10333965Sjdp exponent = a->exponent + b->exponent; 10433965Sjdp size_of_product = product->high - product->low; 10533965Sjdp size_of_sum = size_of_a + size_of_b; 10633965Sjdp extra_product_positions = size_of_product - size_of_sum; 10733965Sjdp if (extra_product_positions < 0) 10833965Sjdp { 10977298Sobrien P = extra_product_positions; /* P < 0 */ 11077298Sobrien exponent -= extra_product_positions; /* Increases exponent. */ 11133965Sjdp } 11233965Sjdp else 11333965Sjdp { 11433965Sjdp P = 0; 11533965Sjdp } 11633965Sjdp carry = 0; 11733965Sjdp significant = 0; 11833965Sjdp for (N = 0; N <= size_of_sum; N++) 11933965Sjdp { 12033965Sjdp work = carry; 12133965Sjdp carry = 0; 12233965Sjdp for (A = 0; A <= N; A++) 12333965Sjdp { 12433965Sjdp B = N - A; 12533965Sjdp if (A <= size_of_a && B <= size_of_b && B >= 0) 12633965Sjdp { 12733965Sjdp#ifdef TRACE 12877298Sobrien printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n", 12977298Sobrien A, a->low[A], B, b->low[B], work); 13033965Sjdp#endif 13133965Sjdp /* Watch out for sign extension! Without the casts, on 13233965Sjdp the DEC Alpha, the multiplication result is *signed* 13333965Sjdp int, which gets sign-extended to convert to the 13433965Sjdp unsigned long! */ 13533965Sjdp work += (unsigned long) a->low[A] * (unsigned long) b->low[B]; 13633965Sjdp carry += work >> LITTLENUM_NUMBER_OF_BITS; 13733965Sjdp work &= LITTLENUM_MASK; 13833965Sjdp#ifdef TRACE 13933965Sjdp printf ("work=%08x carry=%04x\n", work, carry); 14033965Sjdp#endif 14133965Sjdp } 14233965Sjdp } 14333965Sjdp significant |= work; 14433965Sjdp if (significant || P < 0) 14533965Sjdp { 14633965Sjdp if (P >= 0) 14733965Sjdp { 14833965Sjdp product->low[P] = work; 14933965Sjdp#ifdef TRACE 15033965Sjdp printf ("P=%d. work[p]:=%04x\n", P, work); 15133965Sjdp#endif 15233965Sjdp } 15333965Sjdp P++; 15433965Sjdp } 15533965Sjdp else 15633965Sjdp { 15733965Sjdp extra_product_positions++; 15833965Sjdp exponent++; 15933965Sjdp } 16033965Sjdp } 16177298Sobrien /* [P]-> position # size_of_sum + 1. 16277298Sobrien This is where 'carry' should go. */ 16333965Sjdp#ifdef TRACE 16433965Sjdp printf ("final carry =%04x\n", carry); 16533965Sjdp#endif 16633965Sjdp if (carry) 16733965Sjdp { 16833965Sjdp if (extra_product_positions > 0) 16977298Sobrien product->low[P] = carry; 17033965Sjdp else 17133965Sjdp { 17277298Sobrien /* No room at high order for carry littlenum. */ 17377298Sobrien /* Shift right 1 to make room for most significant littlenum. */ 17433965Sjdp exponent++; 17533965Sjdp P--; 17633965Sjdp for (q = product->low + P; q >= product->low; q--) 17733965Sjdp { 17833965Sjdp work = *q; 17933965Sjdp *q = carry; 18033965Sjdp carry = work; 18133965Sjdp } 18233965Sjdp } 18333965Sjdp } 18433965Sjdp else 18577298Sobrien P--; 18633965Sjdp product->leader = product->low + P; 18733965Sjdp product->exponent = exponent; 18833965Sjdp} 189