/cln-1.3.2/src/float/lfloat/elem/cl_LF_I_mul.cc

 1// cl_LF_I_mul().
 2
 3// General includes.
 4#include "base/cl_sysdep.h"
 5
 6// Specification.
 7#include "float/lfloat/cl_LF.h"
 8
 9
10// Implementation.
11
12#include "float/lfloat/cl_LF_impl.h"
13#include "cln/integer.h"
14#include "integer/cl_I.h"
15#include "base/digitseq/cl_DS.h"
16#include "float/cl_F.h"
17
18namespace cln {
19
20const cl_R cl_LF_I_mul (const cl_LF& x, const cl_I& y)
21{
22// Method:
23// If y=0, return 0.
24// If x=0.0, return x.
25// If y is longer than x, convert y to a float and multiply.
26// Else multiply the mantissa of x with the absolute value of y, then round.
27	if (eq(y,0)) { return 0; }
28	if (TheLfloat(x)->expo == 0) { return x; }
29	var cl_signean sign = -(cl_signean)minusp(y); // Vorzeichen von y
30	var cl_I abs_y = (sign==0 ? y : -y);
31	var uintC y_exp = integer_length(abs_y);
32	var uintC len = TheLfloat(x)->len;
33#ifndef CL_LF_PEDANTIC
34	if (ceiling(y_exp,intDsize) > len)
35		return x * cl_I_to_LF(y,len);
36#endif
37	// x länger als y, direkt multiplizieren.
38	CL_ALLOCA_STACK;
39	var const uintD* y_MSDptr;
40	var uintC y_len;
41	var const uintD* y_LSDptr;
42	I_to_NDS_nocopy(abs_y, y_MSDptr=,y_len=,y_LSDptr=,false,); // NDS zu y bilden, y_len>0
43	if (mspref(y_MSDptr,0)==0) y_len--; // NUDS zu y bilden, y_len>0
44	// Multiplizieren.
45	var uintD* prodMSDptr;
46	var uintC prodlen;
47	UDS_UDS_mul_UDS(len,arrayLSDptr(TheLfloat(x)->data,len),
48	                y_len,y_LSDptr,
49	                prodMSDptr=,prodlen=,);
50	// x fing mit 0 Nullbits an, y mit maximal intDsize-1 Nullbits,
51	// daher fängt das Produkt mit maximal intDsize Nullbits an.
52	var uintL shiftcount;
53	if (mspref(prodMSDptr,0)==0) {
54		shiftcount = intDsize;
55		msshrink(prodMSDptr); prodlen--;
56	} else {
57		integerlengthD(mspref(prodMSDptr,0), shiftcount = intDsize -);
58		if (shiftcount > 0)
59			shiftleft_loop_lsp(prodMSDptr mspop (len+1),len+1,shiftcount,0);
60	}
61	// Produkt ist nun normalisiert: höchstes Bit =1.
62	// exponent := exponent(x) + intDsize*y_len - shiftcount
63	var uintE uexp = TheLfloat(x)->expo;
64	var uintE iexp = intDsize*y_len - shiftcount; // >= 0 !
65	uexp = uexp + iexp;
66	if ((uexp < iexp) || (uexp > LF_exp_high))
67		throw floating_point_overflow_exception();
68	// Runden:
69	var uintD* midptr = prodMSDptr mspop len;
70	var uintC restlen = prodlen - len;
71	if ( (restlen==0)
72	     || ((sintD)mspref(midptr,0) >= 0) // nächstes Bit =0 -> abrunden
73	     || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
74	          && !test_loop_msp(midptr mspop 1,restlen-1)
75	          // round-to-even
76	          && ((lspref(midptr,0) & bit(0)) ==0)
77	   )    )
78	  // abrunden
79	  {}
80	  else
81	  // aufrunden
82	  { if ( inc_loop_lsp(midptr,len) )
83	      // Übertrag durchs Aufrunden
84	      { mspref(prodMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
85	        if (++uexp == LF_exp_high+1) { throw floating_point_overflow_exception(); }
86	  }   }
87	return encode_LFu(TheLfloat(x)->sign ^ sign, uexp, prodMSDptr, len);
88}
89// Bit complexity (N = max(length(x),length(y))): O(M(N)).
90
91}  // namespace cln