flonum-mult.c revision 77298
133965Sjdp/* flonum_mult.c - multiply two flonums 277298Sobrien Copyright (C) 1987, 1990, 1991, 1992, 2000 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 2233965Sjdp#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 7433965Sjdpflonum_multip (a, b, product) 7533965Sjdp const FLONUM_TYPE *a; 7633965Sjdp const FLONUM_TYPE *b; 7733965Sjdp FLONUM_TYPE *product; 7833965Sjdp{ 7977298Sobrien int size_of_a; /* 0 origin */ 8077298Sobrien int size_of_b; /* 0 origin */ 8177298Sobrien int size_of_product; /* 0 origin */ 8277298Sobrien int size_of_sum; /* 0 origin */ 8377298Sobrien int extra_product_positions; /* 1 origin */ 8433965Sjdp unsigned long work; 8533965Sjdp unsigned long carry; 8633965Sjdp long exponent; 8733965Sjdp LITTLENUM_TYPE *q; 8833965Sjdp long significant; /* TRUE when we emit a non-0 littlenum */ 8977298Sobrien /* ForTran accent follows. */ 9077298Sobrien int P; /* Scan product low-order -> high. */ 9133965Sjdp int N; /* As in sum above. */ 9277298Sobrien int A; /* Which [] of a? */ 9377298Sobrien int B; /* Which [] of b? */ 9433965Sjdp 9577298Sobrien if ((a->sign != '-' && a->sign != '+') 9677298Sobrien || (b->sign != '-' && b->sign != '+')) 9733965Sjdp { 9877298Sobrien /* Got to fail somehow. Any suggestions? */ 9933965Sjdp product->sign = 0; 10033965Sjdp return; 10133965Sjdp } 10233965Sjdp product->sign = (a->sign == b->sign) ? '+' : '-'; 10333965Sjdp size_of_a = a->leader - a->low; 10433965Sjdp size_of_b = b->leader - b->low; 10533965Sjdp exponent = a->exponent + b->exponent; 10633965Sjdp size_of_product = product->high - product->low; 10733965Sjdp size_of_sum = size_of_a + size_of_b; 10833965Sjdp extra_product_positions = size_of_product - size_of_sum; 10933965Sjdp if (extra_product_positions < 0) 11033965Sjdp { 11177298Sobrien P = extra_product_positions; /* P < 0 */ 11277298Sobrien exponent -= extra_product_positions; /* Increases exponent. */ 11333965Sjdp } 11433965Sjdp else 11533965Sjdp { 11633965Sjdp P = 0; 11733965Sjdp } 11833965Sjdp carry = 0; 11933965Sjdp significant = 0; 12033965Sjdp for (N = 0; N <= size_of_sum; N++) 12133965Sjdp { 12233965Sjdp work = carry; 12333965Sjdp carry = 0; 12433965Sjdp for (A = 0; A <= N; A++) 12533965Sjdp { 12633965Sjdp B = N - A; 12733965Sjdp if (A <= size_of_a && B <= size_of_b && B >= 0) 12833965Sjdp { 12933965Sjdp#ifdef TRACE 13077298Sobrien printf ("a:low[%d.]=%04x b:low[%d.]=%04x work_before=%08x\n", 13177298Sobrien A, a->low[A], B, b->low[B], work); 13233965Sjdp#endif 13333965Sjdp /* Watch out for sign extension! Without the casts, on 13433965Sjdp the DEC Alpha, the multiplication result is *signed* 13533965Sjdp int, which gets sign-extended to convert to the 13633965Sjdp unsigned long! */ 13733965Sjdp work += (unsigned long) a->low[A] * (unsigned long) b->low[B]; 13833965Sjdp carry += work >> LITTLENUM_NUMBER_OF_BITS; 13933965Sjdp work &= LITTLENUM_MASK; 14033965Sjdp#ifdef TRACE 14133965Sjdp printf ("work=%08x carry=%04x\n", work, carry); 14233965Sjdp#endif 14333965Sjdp } 14433965Sjdp } 14533965Sjdp significant |= work; 14633965Sjdp if (significant || P < 0) 14733965Sjdp { 14833965Sjdp if (P >= 0) 14933965Sjdp { 15033965Sjdp product->low[P] = work; 15133965Sjdp#ifdef TRACE 15233965Sjdp printf ("P=%d. work[p]:=%04x\n", P, work); 15333965Sjdp#endif 15433965Sjdp } 15533965Sjdp P++; 15633965Sjdp } 15733965Sjdp else 15833965Sjdp { 15933965Sjdp extra_product_positions++; 16033965Sjdp exponent++; 16133965Sjdp } 16233965Sjdp } 16377298Sobrien /* [P]-> position # size_of_sum + 1. 16477298Sobrien This is where 'carry' should go. */ 16533965Sjdp#ifdef TRACE 16633965Sjdp printf ("final carry =%04x\n", carry); 16733965Sjdp#endif 16833965Sjdp if (carry) 16933965Sjdp { 17033965Sjdp if (extra_product_positions > 0) 17177298Sobrien product->low[P] = carry; 17233965Sjdp else 17333965Sjdp { 17477298Sobrien /* No room at high order for carry littlenum. */ 17577298Sobrien /* Shift right 1 to make room for most significant littlenum. */ 17633965Sjdp exponent++; 17733965Sjdp P--; 17833965Sjdp for (q = product->low + P; q >= product->low; q--) 17933965Sjdp { 18033965Sjdp work = *q; 18133965Sjdp *q = carry; 18233965Sjdp carry = work; 18333965Sjdp } 18433965Sjdp } 18533965Sjdp } 18633965Sjdp else 18777298Sobrien P--; 18833965Sjdp product->leader = product->low + P; 18933965Sjdp product->exponent = exponent; 19033965Sjdp} 191