]> www.ginac.de Git - cln.git/blob - src/float/lfloat/elem/cl_LF_square.cc
Remove exception hooks in favor of real C++ exceptions:
[cln.git] / src / float / lfloat / elem / cl_LF_square.cc
1 // square().
2
3 // General includes.
4 #include "cl_sysdep.h"
5
6 // Specification.
7 #include "cln/lfloat.h"
8
9
10 // Implementation.
11
12 #include "cl_LF.h"
13 #include "cl_LF_impl.h"
14 #include "cl_DS.h"
15 #include "cl_F.h"
16
17 namespace cln {
18
19 const cl_LF square (const cl_LF& x)
20 {
21 // Methode: wie operator*(x,x).
22       var uintC len = TheLfloat(x)->len;
23       var uintE uexp = TheLfloat(x)->expo;
24       if (uexp==0) // x=0.0 -> Ergebnis 0.0
25         { return x; }
26       // Exponenten addieren:
27       // (uexp-LF_exp_mid) + (uexp-LF_exp_mid) = (2*uexp-LF_exp_mid)-LF_exp_mid
28       if ((sintE)uexp >= 0)
29         // kein Carry
30         { uexp = 2*uexp;
31           if (uexp < LF_exp_mid+LF_exp_low)
32             { if (underflow_allowed())
33                 { throw floating_point_underflow_exception(); }
34                 else
35                 { return encode_LF0(len); } // Ergebnis 0.0
36         }   }
37         else
38         // Carry
39         { uexp = 2*uexp;
40           if (uexp > (uintE)(LF_exp_mid+LF_exp_high+1)) { throw floating_point_overflow_exception(); }
41         }
42       uexp = uexp - LF_exp_mid;
43       // Nun ist LF_exp_low <= uexp <= LF_exp_high+1.
44       // neues Long-Float allozieren:
45       var Lfloat y = allocate_lfloat(len,uexp,0);
46       // Produkt bilden:
47       var const uintD* x_LSDptr = arrayLSDptr(TheLfloat(x)->data,len);
48       var uintD* MSDptr;
49       var uintD* LSDptr;
50       CL_ALLOCA_STACK;
51       num_stack_alloc(2*len,MSDptr=,LSDptr=);
52       cl_UDS_mul_square(x_LSDptr,len,LSDptr);
53       {var uintD* midptr = MSDptr mspop len; // Pointer in die Mitte der 2*len Digits
54        if ((sintD)mspref(MSDptr,0) >= 0) // führendes Bit abtesten
55          { // erste n+1 Digits um 1 Bit nach links schieben:
56            shift1left_loop_lsp(midptr mspop 1,len+1);
57            // Exponenten decrementieren:
58            if ((TheLfloat(y)->expo)-- == LF_exp_low-1)
59              { if (underflow_allowed())
60                  { throw floating_point_underflow_exception(); }
61                  else
62                  { return encode_LF0(len); } // Ergebnis 0.0
63              }
64          }
65        // erste Hälfte des Mantissenprodukts übertragen:
66        {var uintD* y_mantMSDptr = arrayMSDptr(TheLfloat(y)->data,len);
67         var uintD* y_mantLSDptr = copy_loop_msp(MSDptr,y_mantMSDptr,len);
68         // Runden:
69         if ( ((sintD)mspref(midptr,0) >= 0) // nächstes Bit =0 -> abrunden
70              || ( ((mspref(midptr,0) & ((uintD)bit(intDsize-1)-1)) ==0) // Bit =1, weitere Bits >0 -> aufrunden
71                   && !test_loop_msp(midptr mspop 1,len-1)
72                   // round-to-even
73                   && ((lspref(midptr,0) & bit(0)) ==0)
74            )    )
75           // abrunden
76           {}
77           else
78           // aufrunden
79           { if ( inc_loop_lsp(y_mantLSDptr,len) )
80               { // Übertrag durchs Aufrunden (kann nur auftreten,
81                 // wenn vorhin um 1 Bit nach links geschoben wurde)
82                 mspref(y_mantMSDptr,0) = bit(intDsize-1); // Mantisse := 10...0
83                 (TheLfloat(y)->expo)++; // Exponent wieder zurück-erhöhen
84           }   }
85         // LF_exp_low <= exp <= LF_exp_high sicherstellen:
86         if (TheLfloat(y)->expo == LF_exp_high+1) { throw floating_point_overflow_exception(); }
87       }}
88       return y;
89 }
90 // Bit complexity (N = length(x)): O(M(N)).
91
92 }  // namespace cln