// Bruno Haible 8.9.1990 - 10.9.1990
// Grundgedanken:
-// Jede Real-Zahl /= 0 repräsentiert ein (offenes) Intervall. Es wird die-
-// jenige Dezimalzahl mit möglichst wenig Stellen ausgegeben, die in diesem
+// Jede Real-Zahl /= 0 repräsentiert ein (offenes) Intervall. Es wird die-
+// jenige Dezimalzahl mit möglichst wenig Stellen ausgegeben, die in diesem
// Intervall liegt.
-// Um auch große Exponenten zu behandeln, werden Zweier- in Zehnerpotenzen
-// erst einmal näherungsweise umgerechnet. Nötigenfalls wird die Rechen-
-// genauigkeit erhöht. Hierbei wird von den Long-Floats beliebiger
+// Um auch große Exponenten zu behandeln, werden Zweier- in Zehnerpotenzen
+// erst einmal näherungsweise umgerechnet. Nötigenfalls wird die Rechen-
+// genauigkeit erhöht. Hierbei wird von den Long-Floats beliebiger
// Genauigkeit Gebrauch gemacht.
-// Stützt sich auf:
+// Stützt sich auf:
// cl_ln2(digits) liefert ln(2) mit mindestens digits Mantissenbits.
// cl_ln10(digits) liefert ln(10) mit mindestens digits Mantissenbits.
// cl_decimal_string(integer) liefert zu einem Integer >0
// einen String mit seiner Dezimaldarstellung.
-// (substring string start [end]) wie subseq, jedoch für Strings schneller.
+// (substring string start [end]) wie subseq, jedoch für Strings schneller.
-CL_REQUIRE(cl_F_ln2_var)
-CL_REQUIRE(cl_F_ln10_var)
#include <cstring>
#include "cln/output.h"
#include "cl_sstring.h"
// berechnet mit folgenden Eigenschaften:
// s = sign(x).
// Falls x/=0, betrachte |x| statt x. Also oBdA x>0.
-// Seien x1 und x2 die nächstkleinere bzw. die nächstgrößere Zahl zu x
-// vom selben Floating-Point-Format. Die Zahl x repräsentiert somit das
+// Seien x1 und x2 die nächstkleinere bzw. die nächstgrößere Zahl zu x
+// vom selben Floating-Point-Format. Die Zahl x repräsentiert somit das
// offene Intervall von (x+x1)/2 bis (x+x2)/2.
// a ist ein Integer >0, mit genau k Dezimalstellen (k>=1), und es gilt
// (x+x1)/2 < a*10^(-k+e) < (x+x2)/2 .
// Dabei ist k minimal, also a nicht durch 10 teilbar.
// Falls x=0: a=0, k=1, e=0.
-// as ist die Ziffernfolge von a, der Länge k.
+// as ist die Ziffernfolge von a, der Länge k.
// typedef
struct cl_decimal_decoded_float {
unten = minus1(ash(binmant2,1));
untenshift = 1;
}
- // Bestimme d (ganz) und a1,a2 (ganz, >0) so, daß
+ // Bestimme d (ganz) und a1,a2 (ganz, >0) so, daß
// die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
// die ganzen a mit a1 <= a <= a2 sind und 0 <= a2-a1 < 20 gilt.
// Wandle dazu 2^e := 2^(binexpo-1) ins Dezimalsystem um.
var cl_I e = binexpo - 1;
- var bool e_gross = (abs(e) > ash(l,1)); // Ist |e| recht groß, >2*l ?
- var uintC g; // Hilfsvariablen für den Fall, daß |e| groß ist
+ var bool e_gross = (abs(e) > ash(l,1)); // Ist |e| recht groß, >2*l ?
+ var uintC g; // Hilfsvariablen für den Fall, daß |e| groß ist
var cl_I f; //
- var cl_I zehn_d; // Hilfsvariable 10^|d| für den Fall, daß |e| klein ist
+ var cl_I zehn_d; // Hilfsvariable 10^|d| für den Fall, daß |e| klein ist
var cl_I d; // Ergebnisvariablen
var cl_I a1; //
var cl_I a2; //
- if (e_gross) { // Ist |e| recht groß ?
- // Da 2^e nur näherungsweise gehen kann, braucht man Schutzbits.
- var uintL h = 16; // Anzahl der Schutzbits, muß >= 3 sein
+ if (e_gross) { // Ist |e| recht groß ?
+ // Da 2^e nur näherungsweise gehen kann, braucht man Schutzbits.
+ var uintL h = 16; // Anzahl der Schutzbits, muß >= 3 sein
neue_schutzbits:
// Ziel: 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10.
- g = l + h; // Anzahl der gültigen Bits von f
- // Schätze d = floor(e*lg(2))
- // mit Hilfe der Näherungsbrüche von lg(2):
+ g = l + h; // Anzahl der gültigen Bits von f
+ // Schätze d = floor(e*lg(2))
+ // mit Hilfe der Näherungsbrüche von lg(2):
// (0 1/3 3/10 28/93 59/196 146/485 643/2136 4004/13301
// 8651/28738 12655/42039 21306/70777 76573/254370 97879/325147
// 1838395/6107016 1936274/6432163 13456039/44699994
// 174131244785/578451474249 845863046269/2809896217828
// 1865857337323/6198243909905 6443435058238/21404627947543
// )
- // e>=0 : wähle lg(2) < a/b < lg(2) + 1/e,
+ // e>=0 : wähle lg(2) < a/b < lg(2) + 1/e,
// dann ist d <= floor(e*a/b) <= d+1 .
- // e<0 : wähle lg(2) - 1/abs(e) < a/b < lg(2),
+ // e<0 : wähle lg(2) - 1/abs(e) < a/b < lg(2),
// dann ist d <= floor(e*a/b) <= d+1 .
// Es ist bekannt, dass abs(e) <= 2^31 + 2^32*64, falls intEsize == 32,
// bzw. dass abs(e) <= 2^63 + 2^64*64, falls intEsize == 64.
- // (Hierbei steht 64 für die maximale intDsize und es wurde benutzt,
+ // (Hierbei steht 64 für die maximale intDsize und es wurde benutzt,
// dass intEsize >= intCsize.)
// Unser d sei := floor(e*a/b)-1. (d /= 0, da abs(e) >= 7.)
d = minus1(minusp(e)
? (e >= -970
- ? floor1(e*3,10) // Näherungsbruch 3/10
+ ? floor1(e*3,10) // Näherungsbruch 3/10
#if (intEsize==32)
- : floor1(e*97879,325147) // Näherungsbruch 97879/325147
+ : floor1(e*97879,325147) // Näherungsbruch 97879/325147
#else
: (e >= -1800000000LL
- ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
- : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
+ ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
+ : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
)
#endif
)
: (e <= 22000
- ? floor1(e*28,93) // Näherungsbruch 28/93
+ ? floor1(e*28,93) // Näherungsbruch 28/93
#if (intEsize==32)
- : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
+ : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
#else
: (e <= 3300000000LL
- ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
- : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
+ ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
+ : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
)
#endif
)
);
- // Das wahre d wird durch diese Schätzung entweder getroffen
- // oder um 1 unterschätzt.
- // Anders ausgedrückt: 0 < e*log(2)-d*log(10) < 2*log(10).
+ // Das wahre d wird durch diese Schätzung entweder getroffen
+ // oder um 1 unterschätzt.
+ // Anders ausgedrückt: 0 < e*log(2)-d*log(10) < 2*log(10).
// Nun f/2^g als exp(e*log(2)-d*log(10)) berechnen.
// Da f < 100*2^g < 2^(g+7), sind g+7 Bits relative Genauigkeit
// des Ergebnisses, also g+7 Bits absolute Genauigkeit von
- // e*log(2)-d*log(10) nötig. Dazu mit l'=integerlength(e)
- // für log(2): g+7+l' Bits abs. Gen., g+7+l' Bits rel. Gen.,
- // für log(10): g+7+l' Bits abs. Gen., g+7+l'+2 Bist rel. Gen.
+ // e*log(2)-d*log(10) nötig. Dazu mit l'=integerlength(e)
+ // für log(2): g+7+l' Bits abs. Gen., g+7+l' Bits rel. Gen.,
+ // für log(10): g+7+l' Bits abs. Gen., g+7+l'+2 Bist rel. Gen.
var float_format_t gen = (float_format_t)(g + integer_length(e) + 9); // Genauigkeit
var cl_F f2g = exp(The(cl_F)(e * cl_ln2(gen)) - The(cl_F)(d * cl_ln10(gen))); // f/2^g
// Das so berechnete f/2^g ist >1, <100.
if (f >= ash(10,g)) // f >= 10*2^g ?
{ f = floor1(f,10); d = d+1; }
// Nun ist 2^e ~= 10^d * f/2^g, wobei 1 <= f/2^g < 10 und
- // f ein Integer ist, der um höchstens 1 vom wahren Wert abweicht:
+ // f ein Integer ist, der um höchstens 1 vom wahren Wert abweicht:
// 10^d * (f-1)/2^g < 2^e < 10^d * (f+1)/2^g
// Wir verkleinern nun das offene Intervall
// von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
// und suchen darin Zahlen der Form 10^d * a mit ganzem a.
// Wegen oben - unten/2^untenshift >= 3/2
// und oben + unten/2^untenshift <= 4*binmant+1 < 2^(l+2) <= 2^(g-1)
- // ist die Intervall-Länge
+ // ist die Intervall-Länge
// = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
// = 10^d * ( f * (oben - unten/2^untenshift)
// - (oben + unten/2^untenshift) ) / 2^g
// mit a1 <= a <= a2, wobei a2 = floor((f-1)*oben/2^g) und
// a1 = ceiling((f+1)*unten/2^(g+untenshift))
// = floor(((f+1)*unten-1)/2^(g+untenshift))+1 .
- // Wir haben eben gesehen, daß a1 <= a2 sein muß.
+ // Wir haben eben gesehen, daß a1 <= a2 sein muß.
a1 = plus1(ash(minus1((f+1)*unten),-(g+untenshift)));
a2 = ash((f-1)*oben,-g);
- // Wir können auch das offene Intervall
+ // Wir können auch das offene Intervall
// von (x+x1)/2 = 2^(binexpo-1-untenshift) * unten
// bis (x+x2)/2 = 2^(binexpo-1) * oben
// in das (abgeschlossene) Intervall
// und sich a1' und a2' analog zu a1 und a2 berechnen.
// Da (f-1)*oben/2^g und (f+1)*oben/2^g sich um 2*oben/2^g
// < 2^(l+2-g) < 1 unterscheiden, unterscheiden sich a2 und
- // a2' um höchstens 1.
+ // a2' um höchstens 1.
// Ebenso, wenn 'oben' durch 'unten/2^untenshift' ersetzt
- // wird: a1' und a1 unterscheiden sich um höchstens 1.
+ // wird: a1' und a1 unterscheiden sich um höchstens 1.
// Ist nun a1' < a1 oder a2 < a2' , so ist die Zweierpotenz-
- // Näherung 10^d * f/2^g für 2^e nicht genau genug gewesen,
- // und man hat das Ganze mit erhöhtem h zu wiederholen.
- // Ausnahme (da hilft auch keine höhere Genauigkeit):
+ // Näherung 10^d * f/2^g für 2^e nicht genau genug gewesen,
+ // und man hat das Ganze mit erhöhtem h zu wiederholen.
+ // Ausnahme (da hilft auch keine höhere Genauigkeit):
// Wenn die obere oder untere Intervallgrenze (x+x2)/2 bzw.
// (x+x1)/2 selbst die Gestalt 10^d * a mit ganzem a hat.
// Dies testet man so:
// (x+x2)/2 = 2^e * oben == 10^d * a mit ganzem a, wenn
- // - für e>=0, (dann 0 <= d <= e): 5^d | oben,
- // - für e<0, (dann e <= d < 0): 2^(d-e) | oben, was
- // nur für d-e=0 der Fall ist.
+ // - für e>=0, (dann 0 <= d <= e): 5^d | oben,
+ // - für e<0, (dann e <= d < 0): 2^(d-e) | oben, was
+ // nur für d-e=0 der Fall ist.
// (x+x1)/2 = 2^(e-untenshift) * unten == 10^d * a
// mit ganzem a, wenn
- // - für e>0, (dann 0 <= d < e): 5^d | unten,
- // - für e<=0, (dann e <= d <= 0): 2^(d-e+untenshift) | unten,
- // was nur für d-e+untenshift=0 der Fall ist.
- // Da wir es jedoch mit großem |e| zu tun haben, kann dieser
+ // - für e>0, (dann 0 <= d < e): 5^d | unten,
+ // - für e<=0, (dann e <= d <= 0): 2^(d-e+untenshift) | unten,
+ // was nur für d-e+untenshift=0 der Fall ist.
+ // Da wir es jedoch mit großem |e| zu tun haben, kann dieser
// Ausnahmefall hier gar nicht eintreten!
// Denn im Falle e>=0: Aus e>=2*l und l>=11 folgt
// e >= (l+2)*ln(10)/ln(5) + ln(10)/ln(2),
if (a2 < a2prime)
{ h = 2*h; goto neue_schutzbits; } // h verdoppeln und alles wiederholen
}
- // Jetzt ist a1 der kleinste und a2 der größte Wert, der
- // für a möglich ist.
+ // Jetzt ist a1 der kleinste und a2 der größte Wert, der
+ // für a möglich ist.
// Wegen oben - unten/2^untenshift <= 2
- // ist die obige Intervall-Länge
+ // ist die obige Intervall-Länge
// = 10^d * ((f-1)*oben - (f+1)*unten/2^untenshift) / 2^g
// < 10^d * ((f-1)*oben - (f-1)*unten/2^untenshift) / 2^g
// = 10^d * (f-1)/2^g * (oben - unten/2^untenshift)
// < 10^d * 10 * 2,
- // also gibt es höchstens 20 mögliche Werte für a.
+ // also gibt es höchstens 20 mögliche Werte für a.
} else {
// |e| ist recht klein -> man kann 2^e und 10^d exakt ausrechnen
if (!minusp(e)) {
- // e >= 0. Schätze d = floor(e*lg(2)) wie oben.
+ // e >= 0. Schätze d = floor(e*lg(2)) wie oben.
// Es ist e<=2*l<2^39, falls intCsize == 32,
// bzw. e<=2*l<2^71, falls intCsize == 64.
d = (e <= 22000
- ? floor1(e*28,93) // Näherungsbruch 28/93
+ ? floor1(e*28,93) // Näherungsbruch 28/93
#if (intCsize==32)
- : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
+ : floor1(e*1838395,6107016) // Näherungsbruch 1838395/6107016
#else
: (e <= 3300000000LL
- ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
- : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
+ ? floor1(e*12655,42039) // Näherungsbruch 12655/42039
+ : floor1(e*149338067129LL,496090320832LL) // Näherungsbruch 149338067129/496090320832
)
#endif
);
- // Das wahre d wird durch diese Schätzung entweder getroffen
- // oder um 1 überschätzt, aber das können wir leicht feststellen.
+ // Das wahre d wird durch diese Schätzung entweder getroffen
+ // oder um 1 überschätzt, aber das können wir leicht feststellen.
zehn_d = The(cl_I)(expt(10,d)); // zehn_d = 10^d
if (ash(1,e) < zehn_d) // falls 2^e < 10^d,
- { d = d-1; zehn_d = exquo(zehn_d,10); } // Schätzung korrigieren
+ { d = d-1; zehn_d = exquo(zehn_d,10); } // Schätzung korrigieren
// Nun ist 10^d <= 2^e < 10^(d+1) und zehn_d = 10^d.
// a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
- // a2 sei das größte ganze a < 2^e * oben / 10^d.
+ // a2 sei das größte ganze a < 2^e * oben / 10^d.
// a1 = 1+floor(unten*2^e/(2^untenshift*10^d)),
// a2 = floor((oben*2^e-1)/10^d).
a1 = plus1(floor1(ash(unten,e),ash(zehn_d,untenshift)));
a2 = floor1(minus1(ash(oben,e)),zehn_d);
} else {
- // e < 0. Schätze d = floor(e*lg(2)) wie oben.
+ // e < 0. Schätze d = floor(e*lg(2)) wie oben.
// Es ist |e|<=2*l<2^39, falls intCsize == 32,
// bzw. |e|<=2*l<2^71, falls intCsize == 64.
d = (e >= -970
- ? floor1(e*3,10) // Näherungsbruch 3/10
+ ? floor1(e*3,10) // Näherungsbruch 3/10
#if (intCsize==32)
- : floor1(e*97879,325147) // Näherungsbruch 97879/325147
+ : floor1(e*97879,325147) // Näherungsbruch 97879/325147
#else
: (e >= -1800000000LL
- ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
- : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
+ ? floor1(e*8651,28738) // Näherungsbruch 8651/28738
+ : floor1(e*24793177656LL,82361153417LL) // Näherungsbruch 24793177656/82361153417
)
#endif
);
- // Das wahre d wird durch diese Schätzung entweder getroffen
- // oder um 1 überschätzt, aber das können wir leicht feststellen.
+ // Das wahre d wird durch diese Schätzung entweder getroffen
+ // oder um 1 überschätzt, aber das können wir leicht feststellen.
zehn_d = The(cl_I)(expt(10,-d)); // zehn_d = 10^(-d)
if (integer_length(zehn_d) <= -e) // falls 2^e < 10^d,
- { d = d-1; zehn_d = zehn_d*10; } // Schätzung korrigieren
+ { d = d-1; zehn_d = zehn_d*10; } // Schätzung korrigieren
// Nun ist 10^d <= 2^e < 10^(d+1) und zehn_d = 10^(-d).
// a1 sei das kleinste ganze a > 2^(e-untenshift) * unten / 10^d,
- // a2 sei das größte ganze a < 2^e * oben / 10^d.
+ // a2 sei das größte ganze a < 2^e * oben / 10^d.
// a1 = 1+floor(unten*10^(-d)/2^(-e+untenshift)),
// a2 = floor((oben*10^(-d)-1)/2^(-e))
a1 = plus1(ash(unten*zehn_d,e-untenshift));
}
}
// Nun sind die ganzen a mit (x+x1)/2 < 10^d * a < (x+x2)/2 genau
- // die ganzen a mit a1 <= a <= a2. Deren gibt es höchstens 20.
+ // die ganzen a mit a1 <= a <= a2. Deren gibt es höchstens 20.
// Diese werden in drei Schritten auf einen einzigen reduziert:
- // 1. Enthält der Bereich eine durch 10 teilbare Zahl a ?
+ // 1. Enthält der Bereich eine durch 10 teilbare Zahl a ?
// ja -> setze a1:=ceiling(a1/10), a2:=floor(a2/10), d:=d+1.
- // Danach enthält der Bereich a1 <= a <= a2 höchstens 10
- // mögliche Werte für a.
- // 2. Falls jetzt einer der möglichen Werte durch 10 teilbar ist
+ // Danach enthält der Bereich a1 <= a <= a2 höchstens 10
+ // mögliche Werte für a.
+ // 2. Falls jetzt einer der möglichen Werte durch 10 teilbar ist
// (es kann nur noch einen solchen geben),
- // wird er gewählt, die anderen vergessen.
- // 3. Sonst wird unter allen noch möglichen Werten der zu x
- // nächstgelegene gewählt.
+ // wird er gewählt, die anderen vergessen.
+ // 3. Sonst wird unter allen noch möglichen Werten der zu x
+ // nächstgelegene gewählt.
var bool d_shift = false; // Flag, ob im 1. Schritt d incrementiert wurde
- var cl_I a; // das ausgewählte a
+ var cl_I a; // das ausgewählte a
// 1.
{
var cl_I b1 = ceiling1(a1,10);
a = floor1(a2,10);
if (10*a >= a1) {
// Noch eine durch 10 teilbare Zahl -> durch 10 teilen.
- d = d+1; // noch d erhöhen, zehn-d wird nicht mehr gebraucht
+ d = d+1; // noch d erhöhen, zehn-d wird nicht mehr gebraucht
// Nun a in einen Dezimalstring umwandeln
- // und dann Nullen am Schluß streichen:
+ // und dann Nullen am Schluß streichen:
var char* as = cl_decimal_string(a); // Ziffernfolge zu a>0
- var uintC las = ::strlen(as); // Länge der Ziffernfolge
- var uintC k = las; // Länge ohne die gestrichenen Nullen am Schluß
+ var uintC las = ::strlen(as); // Länge der Ziffernfolge
+ var uintC k = las; // Länge ohne die gestrichenen Nullen am Schluß
var cl_I ee = k+d; // a * 10^d = a * 10^(-k+ee)
- while (as[k-1] == '0') // eine 0 am Schluß?
+ while (as[k-1] == '0') // eine 0 am Schluß?
{ // ja -> a := a / 10 (wird aber nicht mehr gebraucht),
// d := d+1 (wird aber nicht mehr gebraucht),
k = k-1; as[k] = '\0';
// a1=a2 -> keine Frage der Auswahl mehr:
a = a1;
} else {
- // a1<a2 -> zu x nächstgelegenes 10^d * a wählen:
+ // a1<a2 -> zu x nächstgelegenes 10^d * a wählen:
if (e_gross) {
// a = round(f*2*binmant/2^g/(1oder10)) (beliebige Rundung)
- // = ceiling(floor(f*2*binmant/(1oder10)/2^(g-1))/2) wählen:
+ // = ceiling(floor(f*2*binmant/(1oder10)/2^(g-1))/2) wählen:
var cl_I temp = f * binmant2;
if (d_shift) { temp = floor1(temp,10); }
a = ash(plus1(ash(temp,1-g)),-1);
var cl_I& expo = z_decoded.e;
var cl_I& sign = z_decoded.s;
// arg in Dezimaldarstellung: +/- 0.mant * 10^expo, wobei
- // mant die Mantisse: als Simple-String mantstring mit Länge mantlen,
+ // mant die Mantisse: als Simple-String mantstring mit Länge mantlen,
// expo der Dezimal-Exponent,
// sign das Vorzeichen (-1 oder 0 oder 1).
if (eq(sign,-1)) // z < 0 ?
// 0 < expo < mantlen ->
// die ersten expo Ziffern, Punkt, die restlichen Ziffern
// expo >= mantlen -> alle Ziffern, expo-mantlen Nullen, Punkt, Null
- // Nach Möglichkeit kein Exponent// wenn nötig, Exponent 0.
- // flag gelöscht -> "scientific notation":
+ // Nach Möglichkeit kein Exponent// wenn nötig, Exponent 0.
+ // flag gelöscht -> "scientific notation":
// erste Ziffer, Punkt, die restlichen Ziffern, bei mantlen=1 eine Null
// Exponent.
if (flag && !plusp(expo)) {
fprintchar(stream,mantstring[i]);
}
} else {
- // scale>=mantlen -> es bleibt nichts für die Nachkommastellen.
+ // scale>=mantlen -> es bleibt nichts für die Nachkommastellen.
// alle Ziffern, dann scale-mantlen Nullen, dann Punkt und Null
fprint(stream,mantstring);
for (uintV i = mantlen; i < scale; i++)
fprintchar(stream,exp_marker);
print_integer(stream,10,expo);
}
- // Fertig. Aufräumen.
+ // Fertig. Aufräumen.
free_hook(mantstring);
}
git://www.ginac.de / cln.git / blobdiff