]> www.ginac.de Git - cln.git/blob - src/float/ffloat/elem/cl_FF_from_I.cc
2006-04-25 Bruno Haible <bruno@clisp.org>
[cln.git] / src / float / ffloat / elem / cl_FF_from_I.cc
1 // cl_I_to_FF().
2
3 // General includes.
4 #include "cl_sysdep.h"
5
6 // Specification.
7 #include "cl_FF.h"
8
9
10 // Implementation.
11
12 #include "cln/integer.h"
13 #include "cl_I.h"
14 #include "cl_DS.h"
15 #include "cl_F.h"
16
17 namespace cln {
18
19 const cl_FF cl_I_to_FF (const cl_I& x)
20 {
21 // Methode:
22 // x=0 -> Ergebnis 0.0
23 // Merke Vorzeichen von x.
24 // x:=(abs x)
25 // Exponent:=(integer-length x)
26 //   Greife die 25 höchstwertigen Bits heraus (angeführt von einer 1).
27 //   Runde das letzte Bit weg:
28 //     Bit 0 = 0 -> abrunden,
29 //     Bit 0 = 1 und Rest =0 -> round-to-even,
30 //     Bit 0 = 1 und Rest >0 -> aufrunden.
31 //   Dabei um ein Bit nach rechts schieben.
32 //   Bei Aufrundung auf 2^24 (rounding overflow) Mantisse um 1 Bit nach rechts
33 //     schieben und Exponent incrementieren.
34       if (eq(x,0)) { return cl_FF_0; }
35       var cl_signean sign = -(cl_signean)minusp(x); // Vorzeichen
36       var cl_I abs_x = (sign==0 ? x : -x);
37       var uintC exp = integer_length(abs_x); // (integer-length x)
38       // NDS zu |x|>0 bilden:
39       var const uintD* MSDptr;
40       var uintC len;
41       I_to_NDS_nocopy(abs_x, MSDptr=,len=,,cl_false,);
42       // MSDptr/len/LSDptr ist die NDS zu x, len>0.
43       // Führende Digits holen: Brauche FF_mant_len+1 Bits, dazu intDsize
44       // Bits (die NDS kann mit bis zu intDsize Nullbits anfangen).
45       // Dann werden diese Bits um (exp mod intDsize) nach rechts geschoben.
46       var uintD msd = msprefnext(MSDptr); // erstes Digit
47       #if (intDsize==64)
48       var uintD msdd = 0; // weiteres Digit
49       if (--len == 0) goto ok;
50       msdd = msprefnext(MSDptr);
51       #else // (intDsize<=32)
52       var uint32 msdd = 0; // weitere min(len-1,32/intDsize) Digits
53       #define NEXT_DIGIT(i)  \
54         { if (--len == 0) goto ok;                                   \
55           msdd |= (uint32)msprefnext(MSDptr) << (32-(i+1)*intDsize); \
56         }
57       DOCONSTTIMES(32/intDsize,NEXT_DIGIT);
58       #undef NEXT_DIGIT
59       #endif
60       --len; ok:
61       #if (intDsize==64)
62       // Die NDS besteht aus msd, msdd, und len weiteren Digits.
63       // Das höchste in 2^intDsize*msd+msdd gesetzte Bit ist Bit Nummer
64       // intDsize-1 + (exp mod intDsize).
65       var uintL shiftcount = exp % intDsize;
66       var uint64 mant = // führende 64 Bits
67         (shiftcount==0
68          ? msdd
69          : ((msd << (64-shiftcount)) | (msdd >> shiftcount))
70         );
71       // Das höchste in mant gesetzte Bit ist Bit Nummer 63.
72       if ( ((mant & bit(62-FF_mant_len)) ==0) // Bit 39 =0 -> abrunden
73            || ( ((mant & (bit(62-FF_mant_len)-1)) ==0) // Bit 39 =1 und Bits 38..0 =0
74                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msdd =0
75                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
76                 // round-to-even, je nach Bit 40 :
77                 && ((mant & bit(63-FF_mant_len)) ==0)
78          )    )
79         // abrunden
80         { mant = mant >> (63-FF_mant_len); }
81         else
82         // aufrunden
83         { mant = mant >> (63-FF_mant_len);
84           mant += 1;
85           if (mant >= bit(FF_mant_len+1)) // rounding overflow?
86             { mant = mant>>1; exp = exp+1; }
87         }
88       #else
89       // Die NDS besteht aus msd, msdd, und len weiteren Digits.
90       // Das höchste in 2^32*msd+msdd gesetzte Bit ist Bit Nummer
91       // 31 + (exp mod intDsize).
92       var uintL shiftcount = exp % intDsize;
93       var uint32 mant = // führende 32 Bits
94         (shiftcount==0
95          ? msdd
96          : (((uint32)msd << (32-shiftcount)) | (msdd >> shiftcount))
97         );
98       // Das höchste in mant gesetzte Bit ist Bit Nummer 31.
99       if ( ((mant & bit(30-FF_mant_len)) ==0) // Bit 7 =0 -> abrunden
100            || ( ((mant & (bit(30-FF_mant_len)-1)) ==0) // Bit 7 =1 und Bits 6..0 =0
101                 && ((msdd & (bit(shiftcount)-1)) ==0) // und weitere Bits aus msdd =0
102                 && (!test_loop_msp(MSDptr,len)) // und alle weiteren Digits =0
103                 // round-to-even, je nach Bit 8 :
104                 && ((mant & bit(31-FF_mant_len)) ==0)
105          )    )
106         // abrunden
107         { mant = mant >> (31-FF_mant_len); }
108         else
109         // aufrunden
110         { mant = mant >> (31-FF_mant_len);
111           mant += 1;
112           if (mant >= bit(FF_mant_len+1)) // rounding overflow?
113             { mant = mant>>1; exp = exp+1; }
114         }
115       #endif
116       return encode_FF(sign,(sintL)exp,mant);
117 }
118
119 }  // namespace cln