112 lines
4.4 KiB
112 lines
4.4 KiB
// binary operator *
|
|
|
|
// General includes.
|
|
#include "cl_sysdep.h"
|
|
|
|
// Specification.
|
|
#include "cln/lfloat.h"
|
|
|
|
|
|
// Implementation.
|
|
|
|
#include "cl_LF.h"
|
|
#include "cl_LF_impl.h"
|
|
#include "cl_DS.h"
|
|
#include "cl_F.h"
|
|
|
|
namespace cln {
|
|
|
|
const cl_LF operator* (const cl_LF& x1, const cl_LF& x2)
|
|
{
|
|
// Methode:
|
|
// Falls x1=0.0 oder x2=0.0 -> Ergebnis 0.0
|
|
// Sonst: Ergebnis-Vorzeichen = VZ von x1 xor VZ von x2.
|
|
// Ergebnis-Exponent = Summe der Exponenten von x1 und x2.
|
|
// Produkt der Mantissen bilden (2n Digits).
|
|
// Falls das fhrende Bit =0 ist: Mantissenprodukt um 1 Bit nach links
|
|
// schieben (die vorderen n+1 Digits gengen)
|
|
// und Exponent decrementieren.
|
|
// Runden auf n Digits liefert die Ergebnis-Mantisse.
|
|
var uintC len1 = TheLfloat(x1)->len;
|
|
var uintC len2 = TheLfloat(x2)->len;
|
|
var uintC len = (len1 < len2 ? len1 : len2); // min. L�ge n von x1 und x2
|
|
var uintL uexp1 = TheLfloat(x1)->expo;
|
|
if (uexp1==0) // x1=0.0 -> Ergebnis 0.0
|
|
{ if (len < len1) return shorten(x1,len); else return x1; }
|
|
var uintL uexp2 = TheLfloat(x2)->expo;
|
|
if (uexp2==0) // x2=0.0 -> Ergebnis 0.0
|
|
{ if (len < len2) return shorten(x2,len); else return x2; }
|
|
// Exponenten addieren:
|
|
// (uexp1-LF_exp_mid) + (uexp2-LF_exp_mid) = (uexp1+uexp2-LF_exp_mid)-LF_exp_mid
|
|
uexp1 = uexp1 + uexp2;
|
|
if (uexp1 >= uexp2)
|
|
// kein Carry
|
|
{ if (uexp1 < LF_exp_mid+LF_exp_low)
|
|
{ if (underflow_allowed())
|
|
{ cl_error_floating_point_underflow(); }
|
|
else
|
|
{ return encode_LF0(len); } // Ergebnis 0.0
|
|
} }
|
|
else
|
|
// Carry
|
|
{ if (uexp1 > (uintL)(LF_exp_mid+LF_exp_high+1)) { cl_error_floating_point_overflow(); } }
|
|
uexp1 = uexp1 - LF_exp_mid;
|
|
// Nun ist LF_exp_low <= uexp1 <= LF_exp_high+1.
|
|
// neues Long-Float allozieren:
|
|
var Lfloat y = allocate_lfloat(len,uexp1,
|
|
TheLfloat(x1)->sign ^ TheLfloat(x2)->sign // Vorzeichen kombinieren
|
|
);
|
|
// Produkt bilden:
|
|
var const uintD* x1_LSDptr = arrayLSDptr(TheLfloat(x1)->data,len1);
|
|
var const uintD* x2_LSDptr = arrayLSDptr(TheLfloat(x2)->data,len2);
|
|
#ifndef CL_LF_PEDANTIC
|
|
if (len1 > len2)
|
|
{ x1_LSDptr = x1_LSDptr lspop (len1-(len2+1)); len1 = len2+1; }
|
|
else if (len1 < len2)
|
|
{ x2_LSDptr = x2_LSDptr lspop (len2-(len1+1)); len2 = len1+1; }
|
|
#endif
|
|
var uintD* MSDptr;
|
|
CL_ALLOCA_STACK;
|
|
UDS_UDS_mul_UDS(len1,x1_LSDptr,
|
|
len2,x2_LSDptr,
|
|
MSDptr=,,);
|
|
{var uintD* midptr = MSDptr mspop len; // Pointer in die Mitte der len1+len2 Digits
|
|
if ((sintD)mspref(MSDptr,0) >= 0) // fhrendes Bit abtesten
|
|
{ // erste n+1 Digits um 1 Bit nach links schieben:
|
|
shift1left_loop_lsp(midptr mspop 1,len+1);
|
|
// Exponenten decrementieren:
|
|
if (--(TheLfloat(y)->expo) == LF_exp_low-1)
|
|
{ if (underflow_allowed())
|
|
{ cl_error_floating_point_underflow(); }
|
|
else
|
|
{ return encode_LF0(len); } // Ergebnis 0.0
|
|
}
|
|
}
|
|
// erste H�fte des Mantissenprodukts bertragen:
|
|
{var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
|
|
var uintD* y_mantLSDptr = copy_loop_msp(MSDptr,y_mantMSDptr,len);
|
|
// Runden:
|
|
if ( ((sintD)mspref(midptr,0) >= 0) // n�hstes Bit =0 -> abrunden
|
|
|| ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
|
|
&& !test_loop_msp(midptr mspop 1,len1+len2-len-1)
|
|
// round-to-even
|
|
&& ((lspref(midptr,0) & bit(0)) ==0)
|
|
) )
|
|
// abrunden
|
|
{}
|
|
else
|
|
// aufrunden
|
|
{ if ( inc_loop_lsp(y_mantLSDptr,len) )
|
|
{ // �ertrag durchs Aufrunden (kann nur auftreten,
|
|
// wenn vorhin um 1 Bit nach links geschoben wurde)
|
|
mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
|
|
(TheLfloat(y)->expo)++; // Exponent wieder zurck-erhhen
|
|
} }
|
|
// LF_exp_low <= exp <= LF_exp_high sicherstellen:
|
|
if (TheLfloat(y)->expo == LF_exp_high+1) { cl_error_floating_point_overflow(); }
|
|
}}
|
|
return y;
|
|
}
|
|
// Bit complexity (N = max(length(x1),length(x2))): O(M(N)).
|
|
|
|
} // namespace cln
|