src/float/output/cl_F_dprint.cc

   1 // print_float().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/float_io.h"
   8
   9
  10 // Implementation.
  11
  12 // Michael Stoll 10.2.1990 - 26.3.1990
  13 // Bruno Haible 8.9.1990 - 10.9.1990
  14
  15 // Grundgedanken:
  16 // Jede Real-Zahl /= 0 repräsentiert ein (offenes) Intervall. Es wird die-
  17 // jenige Dezimalzahl mit möglichst wenig Stellen ausgegeben, die in diesem
  18 // Intervall liegt.
  19 // Um auch große Exponenten zu behandeln, werden Zweier- in Zehnerpotenzen
  20 // erst einmal näherungsweise umgerechnet. Nötigenfalls wird die Rechen-
  21 // genauigkeit erhöht. Hierbei wird von den Long-Floats beliebiger
  22 // Genauigkeit Gebrauch gemacht.
  23
  24 // Stützt sich auf:
  25 // cl_ln2(digits) liefert ln(2) mit mindestens digits Mantissenbits.
  26 // cl_ln10(digits) liefert ln(10) mit mindestens digits Mantissenbits.
  27 // cl_decimal_string(integer) liefert zu einem Integer >0
  28 //   einen String mit seiner Dezimaldarstellung.
  29 // (substring string start [end]) wie subseq, jedoch für Strings schneller.
  30
  31 CL_REQUIRE(cl_F_ln2_var)
  32 CL_REQUIRE(cl_F_ln10_var)
  33 #include <cstring>
  34 #include "cln/output.h"
  35 #include "cl_sstring.h"
  36 #include "cln/float.h"
  37 #include "cl_F.h"
  38 #include "cl_LF.h"
  39 #include "cl_F_tran.h"
  40 #include "cln/rational.h"
  41 #include "cln/integer.h"
  42 #include "cln/integer_io.h"
  43 #include "cl_I.h"
  44
  45 namespace cln {
  46
  47 // Hauptfunktion zur Umwandlung von Floats ins Dezimalsystem:
  48 // Zu einem Float x werden ein String as und drei Integers k,e,s
  49 // berechnet mit folgenden Eigenschaften:
  50 // s = sign(x).
  51 // Falls x/=0, betrachte |x| statt x. Also oBdA x>0.
  52 //   Seien x1 und x2 die nächstkleinere bzw. die nächstgrößere Zahl zu x
  53 //   vom selben Floating-Point-Format. Die Zahl x repräsentiert somit das
  54 //   offene Intervall von (x+x1)/2 bis (x+x2)/2.
  55 //   a ist ein Integer >0, mit genau k Dezimalstellen (k>=1), und es gilt
  56 //   (x+x1)/2 < a*10^(-k+e) < (x+x2)/2 .
  57 //   Dabei ist k minimal, also a nicht durch 10 teilbar.
  58 // Falls x=0: a=0, k=1, e=0.
  59 // as ist die Ziffernfolge von a, der Länge k.
  60
  61 // typedef
  62 struct cl_decimal_decoded_float {
  63         char * a;
  64         uintC k;
  65         cl_I e;
  66         cl_I s;
  67 // Constructor.
  68         cl_decimal_decoded_float (char * ap, uintC kp, const cl_I& ep, const cl_I& sp) : a(ap), k(kp), e(ep), s(sp) {}
  69 };
  70
  71
  72 static const cl_decimal_decoded_float decode_float_decimal (const cl_F& x)
  73 {
  74   var cl_idecoded_float x_idecoded = integer_decode_float(x);
  75   var cl_I& binmant = x_idecoded.mantissa;
  76   var cl_I& binexpo = x_idecoded.exponent;
  77   var cl_I& sign = x_idecoded.sign;
  78   if (eq(binmant,0)) // x=0 ?
  79     // a=0, k=1, e=0, s=0
  80     return cl_decimal_decoded_float(cl_sstring("0",1), 1, 0, 0);
  81   // x/=0, also ist sign das Vorzeichen von x und
  82   // |x| = 2^binexpo * float(binmant,x) . Ab jetzt oBdA x>0.
  83   // Also x = 2^binexpo * float(binmant,x) .
  84   var uintC l = integer_length(binmant); // Anzahl der Bits von binmant, >=3
  85   var cl_I binmant2 = ash(binmant,1); // 2*binmant
  86   var cl_I oben = plus1(binmant2); // obere Intervallgrenze ist
  87                                    // (x+x2)/2 = 2^(binexpo-1) * oben
  88   var cl_I unten = minus1(binmant2); // untere Intervallgrenze ist
  89   var uintL untenshift = 0;          // (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
  90   if (integer_length(unten) == l) {
  91     // Normalerweise integerlength(unten) = 1+integerlength(binmant).
  92     // Hier integerlength(unten) = l = integerlength(binmant),
  93     // also war binmant eine Zweierpotenz. In diesem Fall ist die
  94     // die Toleranz nach oben 1/2 Einheit, aber die Toleranz nach unten
  95     // nur 1/4 Einheit: (x+x1)/2 = 2^(binexpo-2) * (4*binmant-1)
  96     unten = minus1(ash(binmant2,1));
  97     untenshift = 1;
  98   }
  99   // Bestimme d (ganz) und a1,a2 (ganz, >0) so, daß
 100   // die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
 101   // die ganzen a mit a1 <= a <= a2 sind und 0 <= a2-a1 < 20 gilt.
 102   // Wandle dazu 2^e := 2^(binexpo-1) ins Dezimalsystem um.
 103   var cl_I e = binexpo - 1;
 104   var bool e_gross = (abs(e) > ash(l,1)); // Ist |e| recht groß, >2*l ?
 105   var uintC g;     // Hilfsvariablen für den Fall, daß |e| groß ist
 106   var cl_I f;      //
 107   var cl_I zehn_d; // Hilfsvariable 10^|d| für den Fall, daß |e| klein ist
 108   var cl_I d;  // Ergebnisvariablen
 109   var cl_I a1; //
 110   var cl_I a2; //
 111   if (e_gross) { // Ist |e| recht groß ?
 112     // Da 2^e nur näherungsweise gehen kann, braucht man Schutzbits.
 113     var uintL h = 16; // Anzahl der Schutzbits, muß >= 3 sein
 114     neue_schutzbits:
 115     // Ziel: 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10.
 116     g = l + h; // Anzahl der gültigen Bits von f
 117     // Schätze d = floor(e*lg(2))
 118     // mit Hilfe der Näherungsbrüche von lg(2):
 119     // (0 1/3 3/10 28/93 59/196 146/485 643/2136 4004/13301
 120     //  8651/28738 12655/42039 21306/70777 76573/254370 97879/325147
 121     //  1838395/6107016 1936274/6432163 13456039/44699994
 122     //  15392313/51132157 44240665/146964308 59632978/198096465
 123     //  103873643/345060773 475127550/1578339557 579001193/1923400330
 124     //  24793177656/82361153417 149338067129/496090320832
 125     //  174131244785/578451474249 845863046269/2809896217828
 126     //  1865857337323/6198243909905 6443435058238/21404627947543
 127     // )
 128     // e>=0 : wähle lg(2) < a/b < lg(2) + 1/e,
 129     //        dann ist d <= floor(e*a/b) <= d+1 .
 130     // e<0  : wähle lg(2) - 1/abs(e) < a/b < lg(2),
 131     //        dann ist d <= floor(e*a/b) <= d+1 .
 132     // Es ist bekannt, dass abs(e) <= 2^31 + 2^32*64, falls intEsize == 32,
 133     //            bzw. dass abs(e) <= 2^63 + 2^64*64, falls intEsize == 64.
 134     // (Hierbei steht 64 für die maximale intDsize und es wurde benutzt,
 135     // dass intEsize >= intCsize.)
 136     // Unser d sei := floor(e*a/b)-1. (d /= 0, da abs(e) >= 7.)
 137     d = minus1(minusp(e)
 138                ? (e >= -970
 139                   ? floor1(e*3,10) // Näherungsbruch 3/10
 140 #if (intEsize==32)
 141                   : floor1(e*97879,325147) // Näherungsbruch 97879/325147
 142 #else
 143                   : (e >= -1800000000LL
 144                      ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
 145                      : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
 146                     )
 147 #endif
 148                  )
 149                : (e <= 22000
 150                   ? floor1(e*28,93) // Näherungsbruch 28/93
 151 #if (intEsize==32)
 152                   : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
 153 #else
 154                   : (e <= 3300000000LL
 155                      ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
 156                      : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
 157                     )
 158 #endif
 159                  )
 160               );
 161     // Das wahre d wird durch diese Schätzung entweder getroffen
 162     // oder um 1 unterschätzt.
 163     // Anders ausgedrückt: 0 < e*log(2)-d*log(10) < 2*log(10).
 164     // Nun f/2^g als exp(e*log(2)-d*log(10)) berechnen.
 165     // Da f < 100*2^g < 2^(g+7), sind g+7 Bits relative Genauigkeit
 166     // des Ergebnisses, also g+7 Bits absolute Genauigkeit von
 167     // e*log(2)-d*log(10) nötig. Dazu mit l'=integerlength(e)
 168     // für log(2): g+7+l' Bits abs. Gen., g+7+l' Bits rel. Gen.,
 169     // für log(10): g+7+l' Bits abs. Gen., g+7+l'+2 Bist rel. Gen.
 170     var float_format_t gen = (float_format_t)(g + integer_length(e) + 9); // Genauigkeit
 171     var cl_F f2g = exp(The(cl_F)(e * cl_ln2(gen)) - The(cl_F)(d * cl_ln10(gen))); // f/2^g
 172     // Das so berechnete f/2^g ist >1, <100.
 173     // Mit 2^g multiplizieren und auf eine ganze Zahl runden:
 174     f = round1(scale_float(f2g,g)); // liefert f
 175     // Eventuell f und d korrigieren:
 176     if (f >= ash(10,g)) // f >= 10*2^g ?
 177       { f = floor1(f,10); d = d+1; }
 178     // Nun ist 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10 und
 179     // f ein Integer ist, der um höchstens 1 vom wahren Wert abweicht:
 180     // 10^d * (f-1)/2^g < 2^e < 10^d * (f+1)/2^g
 181     // Wir verkleinern nun das offene Intervall
 182     // von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
 183     // bis (x+x2)/2 = 2^(binexpo-1) * oben
 184     // zu einem abgeschlossenen Intervall
 185     // von 10^d * (f+1)/2^(g+untenshift) * unten
 186     // bis 10^d * (f-1)/2^g * oben
 187     // und suchen darin Zahlen der Form 10^d * a mit ganzem a.
 188     // Wegen  oben - unten/2^untenshift >= 3/2
 189     // und  oben + unten/2^untenshift <= 4*binmant+1 < 2^(l+2) <= 2^(g-1)
 190     // ist die Intervall-Länge
 191     // = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
 192     // = 10^d * ( f * (oben - unten/2^untenshift)
 193     //            - (oben + unten/2^untenshift) ) / 2^g
 194     // >= 10^d * (2^g * 3/2 - 2^(g-1)) / 2^g
 195     // = 10^d * (3/2 - 2^(-1)) = 10^d
 196     // und daher gibt es in dem Intervall mindestens eine Zahl
 197     // dieser Form.
 198     // Die Zahlen der Form 10^d * a in diesem Intervall sind die
 199     // mit a1 <= a <= a2, wobei a2 = floor((f-1)*oben/2^g) und
 200     // a1 = ceiling((f+1)*unten/2^(g+untenshift))
 201     //    = floor(((f+1)*unten-1)/2^(g+untenshift))+1 .
 202     // Wir haben eben gesehen, daß a1 <= a2 sein muß.
 203     a1 = plus1(ash(minus1((f+1)*unten),-(g+untenshift)));
 204     a2 = ash((f-1)*oben,-g);
 205     // Wir können auch das offene Intervall
 206     // von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
 207     // bis (x+x2)/2 = 2^(binexpo-1) * oben
 208     // in das (abgeschlossene) Intervall
 209     // von 10^d * (f-1)/2^(g+untenshift) * unten
 210     // bis 10^d * (f+1)/2^g * oben
 211     // einschachteln. Hierin sind die Zahlen der Form 10^d * a
 212     // die mit a1' <= a <= a2', wobei a1' <= a1 <= a2 <= a2' ist
 213     // und sich a1' und a2' analog zu a1 und a2 berechnen.
 214     // Da (f-1)*oben/2^g und (f+1)*oben/2^g sich um 2*oben/2^g
 215     // < 2^(l+2-g) < 1 unterscheiden, unterscheiden sich a2 und
 216     // a2' um höchstens 1.
 217     // Ebenso, wenn 'oben' durch 'unten/2^untenshift' ersetzt
 218     // wird: a1' und a1 unterscheiden sich um höchstens 1.
 219     // Ist nun a1' < a1 oder a2 < a2' , so ist die Zweierpotenz-
 220     // Näherung 10^d * f/2^g für 2^e nicht genau genug gewesen,
 221     // und man hat das Ganze mit erhöhtem h zu wiederholen.
 222     // Ausnahme (da hilft auch keine höhere Genauigkeit):
 223     //   Wenn die obere oder untere Intervallgrenze (x+x2)/2 bzw.
 224     //   (x+x1)/2 selbst die Gestalt 10^d * a mit ganzem a hat.
 225     //   Dies testet man so:
 226     //     (x+x2)/2 = 2^e * oben == 10^d * a  mit ganzem a, wenn
 227     //     - für e>=0, (dann 0 <= d <= e): 5^d | oben,
 228     //     - für e<0, (dann e <= d < 0): 2^(d-e) | oben, was
 229     //                nur für d-e=0 der Fall ist.
 230     //     (x+x1)/2 = 2^(e-untenshift) * unten == 10^d * a
 231     //     mit ganzem a, wenn
 232     //     - für e>0, (dann 0 <= d < e): 5^d | unten,
 233     //     - für e<=0, (dann e <= d <= 0): 2^(d-e+untenshift) | unten,
 234     //                 was nur für d-e+untenshift=0 der Fall ist.
 235     // Da wir es jedoch mit großem |e| zu tun haben, kann dieser
 236     // Ausnahmefall hier gar nicht eintreten!
 237     // Denn im Falle e>=0: Aus e>=2*l und l>=11 folgt
 238     //   e >= (l+2)*ln(10)/ln(5) + ln(10)/ln(2),
 239     //   d >= e*ln(2)/ln(10)-1 >= (l+2)*ln(2)/ln(5),
 240     //   5^d >= 2^(l+2),
 241     //   und wegen 0 < unten < 2^(l+2) und 0 < oben < 2^(l+1)
 242     //   sind unten und oben nicht durch 5^d teilbar.
 243     // Und im Falle e<=0: Aus -e>=2*l und l>=6 folgt
 244     //   -e >= (l+2)*ln(10)/ln(5),
 245     //   d-e >= e*ln(2)/ln(10)-1-e = (1-ln(2)/ln(10))*(-e)-1
 246     //          = (-e)*ln(5)/ln(10)-1 >= l+1,
 247     //   2^(d-e) >= 2^(l+1),
 248     //   und wegen 0 < unten < 2^(l+1+untenshift) ist unten nicht
 249     //   durch 2^(d-e+untenshift) teilbar, und wegen
 250     //   0 < oben < 2^(l+1) ist oben nicht durch 2^(d-e) teilbar.
 251     {
 252       var cl_I a1prime = plus1(ash(minus1((f-1)*unten),-(g+untenshift)));
 253       if (a1prime < a1)
 254         { h = 2*h; goto neue_schutzbits; } // h verdoppeln und alles wiederholen
 255       var cl_I a2prime = ash((f+1)*oben,-g);
 256       if (a2 < a2prime)
 257         { h = 2*h; goto neue_schutzbits; } // h verdoppeln und alles wiederholen
 258     }
 259     // Jetzt ist a1 der kleinste und a2 der größte Wert, der
 260     // für a möglich ist.
 261     // Wegen  oben - unten/2^untenshift <= 2
 262     // ist die obige Intervall-Länge
 263     // = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
 264     // < 10^d * ((f-1)*oben - (f-1)*unten/2^untenshift) / 2^g
 265     // = 10^d * (f-1)/2^g * (oben - unten/2^untenshift)
 266     // < 10^d * 10 * 2,
 267     // also gibt es höchstens 20 mögliche Werte für a.
 268   } else {
 269     // |e| ist recht klein -> man kann 2^e und 10^d exakt ausrechnen
 270     if (!minusp(e)) {
 271       // e >= 0. Schätze d = floor(e*lg(2)) wie oben.
 272       // Es ist e<=2*l<2^39, falls intCsize == 32,
 273       //   bzw. e<=2*l<2^71, falls intCsize == 64.
 274       d = (e <= 22000
 275            ? floor1(e*28,93) // Näherungsbruch 28/93
 276 #if (intCsize==32)
 277            : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
 278 #else
 279            : (e <= 3300000000LL
 280               ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
 281               : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
 282              )
 283 #endif
 284           );
 285       // Das wahre d wird durch diese Schätzung entweder getroffen
 286       // oder um 1 überschätzt, aber das können wir leicht feststellen.
 287       zehn_d = The(cl_I)(expt(10,d)); // zehn_d = 10^d
 288       if (ash(1,e) < zehn_d) // falls 2^e < 10^d,
 289         { d = d-1; zehn_d = exquo(zehn_d,10); } // Schätzung korrigieren
 290       // Nun ist 10^d <= 2^e < 10^(d+1) und zehn_d = 10^d.
 291       // a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
 292       // a2 sei das größte ganze a < 2^e * oben / 10^d.
 293       // a1 = 1+floor(unten*2^e/(2^untenshift*10^d)),
 294       // a2 = floor((oben*2^e-1)/10^d).
 295       a1 = plus1(floor1(ash(unten,e),ash(zehn_d,untenshift)));
 296       a2 = floor1(minus1(ash(oben,e)),zehn_d);
 297     } else {
 298       // e < 0. Schätze d = floor(e*lg(2)) wie oben.
 299       // Es ist |e|<=2*l<2^39, falls intCsize == 32,
 300       //   bzw. |e|<=2*l<2^71, falls intCsize == 64.
 301       d = (e >= -970
 302            ? floor1(e*3,10) // Näherungsbruch 3/10
 303 #if (intCsize==32)
 304            : floor1(e*97879,325147) // Näherungsbruch 97879/325147
 305 #else
 306            : (e >= -1800000000LL
 307               ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
 308               : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
 309              )
 310 #endif
 311           );
 312       // Das wahre d wird durch diese Schätzung entweder getroffen
 313       // oder um 1 überschätzt, aber das können wir leicht feststellen.
 314       zehn_d = The(cl_I)(expt(10,-d)); // zehn_d = 10^(-d)
 315       if (integer_length(zehn_d) <= -e) // falls 2^e < 10^d,
 316         { d = d-1; zehn_d = zehn_d*10; } // Schätzung korrigieren
 317       // Nun ist 10^d <= 2^e < 10^(d+1) und zehn_d = 10^(-d).
 318       // a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
 319       // a2 sei das größte ganze a < 2^e * oben / 10^d.
 320       // a1 = 1+floor(unten*10^(-d)/2^(-e+untenshift)),
 321       // a2 = floor((oben*10^(-d)-1)/2^(-e))
 322       a1 = plus1(ash(unten*zehn_d,e-untenshift));
 323       a2 = ash(minus1(oben*zehn_d),e);
 324     }
 325   }
 326   // Nun sind die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
 327   // die ganzen a mit a1 <= a <= a2. Deren gibt es höchstens 20.
 328   // Diese werden in drei Schritten auf einen einzigen reduziert:
 329   // 1. Enthält der Bereich eine durch 10 teilbare Zahl a ?
 330   //    ja -> setze a1:=ceiling(a1/10), a2:=floor(a2/10), d:=d+1.
 331   // Danach enthält der Bereich a1 <= a <= a2 höchstens 10
 332   // mögliche Werte für a.
 333   // 2. Falls jetzt einer der möglichen Werte durch 10 teilbar ist
 334   //    (es kann nur noch einen solchen geben),
 335   //    wird er gewählt, die anderen vergessen.
 336   // 3. Sonst wird unter allen noch möglichen Werten der zu x
 337   //    nächstgelegene gewählt.
 338   var bool d_shift = false; // Flag, ob im 1. Schritt d incrementiert wurde
 339   var cl_I a; // das ausgewählte a
 340   // 1.
 341   {
 342     var cl_I b1 = ceiling1(a1,10);
 343     var cl_I b2 = floor1(a2,10);
 344     if (b1 <= b2) // noch eine durch 10 teilbare Zahl a ?
 345       { a1 = b1; a2 = b2; d = d+1; d_shift = true; }
 346       else
 347       goto keine_10_mehr;
 348   }
 349   // 2.
 350   a = floor1(a2,10);
 351   if (10*a >= a1) {
 352     // Noch eine durch 10 teilbare Zahl -> durch 10 teilen.
 353     d = d+1; // noch d erhöhen, zehn-d wird nicht mehr gebraucht
 354     // Nun a in einen Dezimalstring umwandeln
 355     // und dann Nullen am Schluß streichen:
 356     var char* as = cl_decimal_string(a); // Ziffernfolge zu a>0
 357     var uintC las = ::strlen(as); // Länge der Ziffernfolge
 358     var uintC k = las; // Länge ohne die gestrichenen Nullen am Schluß
 359     var cl_I ee = k+d; // a * 10^d = a * 10^(-k+ee)
 360     while (as[k-1] == '0') // eine 0 am Schluß?
 361       { // ja -> a := a / 10 (wird aber nicht mehr gebraucht),
 362         // d := d+1 (wird aber nicht mehr gebraucht),
 363         k = k-1; as[k] = '\0';
 364       }
 365     return cl_decimal_decoded_float(as,k,ee,sign);
 366   }
 367   // 3.
 368   keine_10_mehr:
 369   if (a1 == a2) {
 370     // a1=a2 -> keine Frage der Auswahl mehr:
 371     a = a1;
 372   } else {
 373     // a1<a2 -> zu x nächstgelegenes 10^d * a wählen:
 374     if (e_gross) {
 375       // a = round(f*2*binmant/2^g/(1oder10)) (beliebige Rundung)
 376       //   = ceiling(floor(f*2*binmant/(1oder10)/2^(g-1))/2) wählen:
 377       var cl_I temp = f * binmant2;
 378       if (d_shift) { temp = floor1(temp,10); }
 379       a = ash(plus1(ash(temp,1-g)),-1);
 380     } else {
 381       // |e| klein -> analog wie oben a2 berechnet wurde
 382       if (!minusp(e)) {
 383         // e>=0: a = round(2^e*2*binmant/10^d)
 384         if (d_shift) { zehn_d = 10*zehn_d; }
 385         a = round1(ash(binmant2,e),zehn_d);
 386       } else {
 387         // e<0, also war d<0, jetzt (wegen Schritt 1) d<=0.
 388         // a = round(2*binmant*10^(-d)/2^(-e))
 389         if (d_shift) { zehn_d = floor1(zehn_d,10); }
 390         a = ash(plus1(ash(binmant2*zehn_d,e+1)),-1);
 391       }
 392     }
 393   }
 394   var char* as = cl_decimal_string(a); // Ziffernfolge zu a>0
 395   var uintC k = ::strlen(as);
 396   ASSERT(as[k-1] != '0');
 397   return cl_decimal_decoded_float(as,k,k+d,sign);
 398 }
 399
 400 // Ausgabefunktion:
 401 void print_float (std::ostream& stream, const cl_print_float_flags& flags, const cl_F& z)
 402 {
 403   var cl_decimal_decoded_float z_decoded = decode_float_decimal(z);
 404   var char * & mantstring = z_decoded.a;
 405   var uintC& mantlen = z_decoded.k;
 406   var cl_I& expo = z_decoded.e;
 407   var cl_I& sign = z_decoded.s;
 408   // arg in Dezimaldarstellung: +/- 0.mant * 10^expo, wobei
 409   //  mant die Mantisse: als Simple-String mantstring mit Länge mantlen,
 410   //  expo der Dezimal-Exponent,
 411   //  sign das Vorzeichen (-1 oder 0 oder 1).
 412   if (eq(sign,-1)) // z < 0 ?
 413     fprintchar(stream,'-');
 414   var bool flag = (expo >= -2) && (expo <= 7); // z=0 oder 10^-3 <= |z| < 10^7 ?
 415   // Was ist auszugeben? Fallunterscheidung:
 416   // flag gesetzt -> "fixed-point notation":
 417   //   expo <= 0 -> Null, Punkt, -expo Nullen, alle Ziffern
 418   //   0 < expo < mantlen ->
 419   //     die ersten expo Ziffern, Punkt, die restlichen Ziffern
 420   //   expo >= mantlen -> alle Ziffern, expo-mantlen Nullen, Punkt, Null
 421   //   Nach Möglichkeit kein Exponent// wenn nötig, Exponent 0.
 422   // flag gelöscht -> "scientific notation":
 423   //   erste Ziffer, Punkt, die restlichen Ziffern, bei mantlen=1 eine Null
 424   //   Exponent.
 425   if (flag && !plusp(expo)) {
 426     // "fixed-point notation" mit expo <= 0
 427     // erst Null und Punkt, dann -expo Nullen, dann alle Ziffern
 428     fprintchar(stream,'0');
 429     fprintchar(stream,'.');
 430     for (uintV i = -FN_to_V(expo); i > 0; i--)
 431       fprintchar(stream,'0');
 432     fprint(stream,mantstring);
 433     expo = 0; // auszugebender Exponent ist 0
 434   } else {
 435     // "fixed-point notation" mit expo > 0 oder "scientific notation"
 436     var uintV scale = (flag ? FN_to_V(expo) : 1);
 437     // Der Dezimalpunkt wird um scale Stellen nach rechts geschoben,
 438     // d.h. es gibt scale Vorkommastellen. scale > 0.
 439     if (scale < mantlen) {
 440       // erst scale Ziffern, dann Punkt, dann restliche Ziffern:
 441       { for (uintL i = 0; i < scale; i++)
 442           fprintchar(stream,mantstring[i]);
 443       }
 444       fprintchar(stream,'.');
 445       { for (uintC i = scale; i < mantlen; i++)
 446           fprintchar(stream,mantstring[i]);
 447       }
 448     } else {
 449       // scale>=mantlen -> es bleibt nichts für die Nachkommastellen.
 450       // alle Ziffern, dann scale-mantlen Nullen, dann Punkt und Null
 451       fprint(stream,mantstring);
 452       for (uintV i = mantlen; i < scale; i++)
 453         fprintchar(stream,'0');
 454       fprintchar(stream,'.');
 455       fprintchar(stream,'0');
 456     }
 457     expo = expo - scale; // der auszugebende Exponent ist um scale kleiner.
 458   }
 459   // Nun geht's zum Exponenten:
 460   var char exp_marker;
 461   floattypecase(z
 462   ,     exp_marker = 's';
 463   ,     exp_marker = 'f';
 464   ,     exp_marker = 'd';
 465   ,     exp_marker = 'L';
 466   );
 467   if (!flags.float_readably) {
 468     floatformatcase(flags.default_float_format
 469     ,   if (exp_marker=='s') { exp_marker = 'E'; }
 470     ,   if (exp_marker=='f') { exp_marker = 'E'; }
 471     ,   if (exp_marker=='d') { exp_marker = 'E'; }
 472     ,   if ((exp_marker=='L') && (len == TheLfloat(z)->len)) { exp_marker = 'E'; }
 473     );
 474   }
 475   if (!(flag && (exp_marker=='E'))) { // evtl. Exponent ganz weglassen
 476     fprintchar(stream,exp_marker);
 477     print_integer(stream,10,expo);
 478   }
 479   // Fertig. Aufräumen.
 480   free_hook(mantstring);
 481 }
 482
 483 }  // namespace cln