[C++] hoge precisie berekeningen

[JanWillem32]

Vanwege het ontwerpen van letttertypen moet ik regelmatig cirkels en cirkelbogen maken. Omdat ik geen zin meer had om bijna-optimale cirkelbogen uit mijn hoofd te rekenen, besloot ik om een keer te programmeren. Ik heb dus een beetje knip- en plakwerk gedaan van voorbeeldscripts en -functies, en het werkt allemaal prima, behalve als ik een betere boog wil dan 295661502 groot. Uit mijn hoofd kwam ik zelf weliswaar niet verder dan de boog van 198 groot, maar ik weet niet hoe ik de hoeveelheid decimalen van een "double"-type vergroot.
R is de maximale radius van de boog, en daarna geeft het programma de coöridinaten van alle mogelijke bogen, met de afwijking ten opzichte van een werkelijke boog op het eind.

C++:

#include <iostream>
#include <cmath>
using namespace std;
int main() {
int h;
cout << "R=";
cin >> h;
double b = 1;
double e = 1;
int f = 1;
int g = 1;
  for (int a = 1; h >= 2*a; a++) {
    b = a*sqrt(2);
    int c = b;
    double d = c-b;
    if (d < 0) {
      d = d*-1; 
      }
    if (d < e) {
      e = d;
      f = b;
      g = a;
      cout << "(0,0) (0," << g*2 << ") (" << 2*(f-g) << "," << g*2 << ") (" << f << "," << f << ") (" << g*2 << "," << 2*(f-g) << ") (" << g*2 << ",0) D:" << e << endl;
      }
    a = a+1;
    }
return 0;
}

[NMe]

Gebruik eens een long double?

offtopic:
En ik hoop dat jij weet wat die code doet, want zo'n code met vaiabelenamen van één letter begrijp ik niet snel. En ik weet wel zeker dat ik daar niet alleen in ben.

[ Voor 75% gewijzigd door NMe op 03-06-2006 14:43 ]

[MSalters]

De hoeveelheid decimalen kun en hoef je niet te vergroten. Dit is typisch numerieke wiskunde, je lost zoiets op met goede algoritmes.

Ik heb echter geen enkel idee waar je mee bezig bent, dus met die algoritmes kan ik je niet helpen. Bogen zijn cirkelsegmenten, en straal is ook een term die met cirkels te maken heeft, maar voor de rest heb ik het gevoel dat je ook zelf niet helemaal weet wat je aan het doen bent. Het feit dat je blijkbaar je variabele 'h' gebruikt voor een straal zegt eigenlijk al genoeg. Als je al zonodig 1 letter wil gebruiken, gebruik dan R !

[Verwijderd]

Ik heb geen verstand van cirkelbogen of van C++ (kan 't wel lezen), maar ik word altijd helemaal kriegel wanneer iemand variabelen en parameters met 1 letter aanduidt. Kom op zeg, we leven niet meer in de tijd van de homecomputers waarbij de BASIC interpreter geen langere namen aankon.

Geef die dingen een zinvolle naam, dat leest en interpreteert een stuk prettiger.

Overigens lijkt me een boog van max ~ 300 miljoen wel aardig werkbaar. Of mis ik iets?

[MSalters]

Ja, 300 miljoen picometer is niet veel. Maar zoals gezegd is dat meer een probleem van de TS; die begrijpt vermoedelijk niet helemaal wat het idee van floating point is.

[JanWillem32]

Het gebruiken van long double heeft wel geholpen, Bij het opgeven van 10.000.000.000 is nu het maximum wat er als straal uitkomt 777767720. De functie van dit programma is eigenlijk een een cirkelboog van boven naar rechts te berekenen als alleen gehele getallen mgen opgegeven worden. Met behulp van de eenheidcirkel is het duidelijk is dat sin 45° = ½√2 . Oftewel ik zoek een geheel getal dat als ik ½√2 daarmee vermenigvuldig dat ik dan ongeveer een geheel getal krijg. Daarvoor gebruik ik de letter a, die van 1 tot en met Rmax optelt. Deviation staat voor de 1-dimensionale afwijking (horizontaal of verticaal) van een werkelijke cirkelboog tot de opgegeven cirkelboog. Uiteindelijk geeft het programma nu een set coördinaten om een cirkelboog van telkens toenemende precisie te maken. de tweede coördinaat vormt de hoogte van de cirkelboog , de zesde de breedte van de cirkelboog. De letter f is gebruikt om de x- en y-coördinaten in te vullen voor het middenpunt op de cirkelboog.

C++:

#include <iostream>
#include <cmath>
using namespace std;
int main() {
long double Rmax;
cout << "Rmax=";
cin >> Rmax;
long double b = 1;
long double Deviation = 1;
int f = 1;
int R = 1;
  for (int a = 1; Rmax >= a; a++) {
    b = a*0.5*sqrt(2);
    int c = b;
    long double d = c-b;
    if (d < 0) {
      d = d*-1; 
      }
    if (d < Deviation) {
      Deviation = d;
      f = b;
      R = a;
      cout << "(0,0) (0," << R << ") (" << 2*f-R << "," << R << ") (" << f << "," << f << ") (" << R << "," << 2*f-R << ") (" << R << ",0) Deviation:" << Deviation << endl;
      }
    }
return 0;
}

edit:
Is een long double nu 64-bit of 80-bit? On-line zie ik verschillen staan tussen de opgaves en in mijn compiler (Dev C++ v4.9.9.2 bèta) is het niet opgegeven. Ik wil wel graag weten in hoe verre "Deviation" zelf betrouwbaar is.

[ Voor 28% gewijzigd door JanWillem32 op 03-06-2006 17:42 ]

[NMe]

JanWillem32 schreef op zaterdag 03 juni 2006 @ 17:35:
Is een long double nu 64-bit of 80-bit? On-line zie ik verschillen staan tussen de opgaves en in mijn compiler (Dev C++ v4.9.9.2 bèta) is het niet opgegeven. Ik wil wel graag weten in hoe verre "Deviation" zelf betrouwbaar is.

Als ik het goed heb verschilt dat per besturingssysteem/compiler. Ik vermoed dat het meestal 64 bits is.

Overigens zijn je variabelenamen nog steeds voor verbetering vatbaar. De letters a-f gebruiken is heel erg vatbaar voor interpretatiefouten. Bekijk over twee weken je eigen code eens en je snapt er waarschijnlijk geen snars meer van.

[MSalters]

Verschilt per compiler, niet per OS. Uit m'n hoofd is MSVC bijvoorbeeld 64 bits, maar ICC 80 bits.

Terug naar het eigenlijk onderwerp: ik snap er nog steeds geen hout van. Ik begrijp nu dat b staat voor het punt (b,b) waar een cirkel met straal a en middelpunt (a,0) doorheen gaat.

Verder redenerend: c is een onzinnige variabele, d is std::floor(b), en je zoekt het minimum van std::abs(d) .

Dat suggereert verder dat je eens moet kijken welke functies er allemaal in <cmath> zitten. Blijft staan dat ik denk dat de numerieke implementatie van je algoritme ernstig suboptimaal is. Ten eerste vergelijk je 10.000.000.000,0 (10 miljard) met een int (terwijl je INT_MAX waarschijnlijk maar 4 miljard is). Ik denk dus dat het een oneindige loop is. Vervolgens assign je ook nog een waarde van mogelijk zo'n 7 miljard aan int c. Kortom, het resultaat is waarschijnlijk volstrekt onbetrouwbaar.

[JanWillem32]

Om maar weer eens een poging te wagen: ik heb weer eens gekeken hoe C++ eigenlijk ging. Aangezien ik al die weken nep-C+ van de scripts in TES4:Oblivion voor mijn kiezen heb gehad, is nette C++ dus niet zo goed meer blijven hangen. Met dank aan u allen , denk ik dat ik eindelijk iets heb dat ook met hoge precisie werkt voor de grotere cirkels. Deze is in ieder geval beter leesbaar.

C++:

#include <iostream>
#include <cmath>
using namespace std;
int main() {
long double Rmax;
cout << "Rmax=";
cin >> Rmax;
long double Deviation = 1;
for (long double R = 1; Rmax >= R; R++) {
  long double a = R*0.5*sqrt(2);
  long double d = a-floor(a);
  if (d < Deviation) {
    Deviation = d;
    cout << "(0,0) (0," << R << ") (" << 2*floor(a)-R << "," << R << ") (" << floor(a) << "," << floor(a) << ") (" << R << "," << 2*floor(a)-R << ") (" << R << ",0) Deviation:" << Deviation << endl;
    }
  d = ceil(a)-a;
  if (d < Deviation) {
    Deviation = d;
    cout << "(0,0) (0," << R << ") (" << 2*ceil(a)-R << "," << R << ") (" << ceil(a) << "," << ceil(a) << ") (" << R << "," << 2*ceil(a)-R << ") (" << R << ",0) Deviation:" << Deviation << endl;
    }
  }
return 0;
}

off-topic: Het derde en vijfde punt dat deze genereert is trouwens geen punt op de cirkelboog, maar een off-curve punt, die gebruikt wordt om de lijn die tussen het vorige en volgende punt zit, krom te trekken naar zichzelf toe. Als je deze twee weglaat, en de delen van de figuur in de (x,-y), (-x,y) en (-x,-y) kwadranten toevoegt, heb je een regelmatige achthoek. Als deze code nu wel goed werkt, zal ik eens kijken hoe ik andere ۞(ster) figuurtjes met deze methode kan maken.

[Verwijderd]

Ik begrijp er eigenlijk weinig van wat je doet. Je lijkt het getal (x) tussen de 1 en RMAX (=10.000.000) te zoeken waarvan het punt y = x * 1/2 *sqrt(2) zo dicht mogelijk bij een integer uitkomt. Maar uit je verhaal maak ik op dat je op zoek bent naar het grootste getal (????).

Wat je nu eigenlijk aan het doen bent is om de volgende formule in te kloppen

min(abs(y-round(y)) waarbij y = x * C. Dit is natuurlijk minimaal wanneer x = 1/C dus wanneer x = sqr(2) of x een veelvoud van sqr(2) is. 14142135623730950488016887242097 lijkt mij een leuke kandidaat. Maar je speurtocht naar het grootste getal lijkt mij dan ook zinloos. Aangezien sqrt(2) een oneindig aantal decimalen heeft (het is te bewijzen dat sqrt(2) geen breuk is.) kun je theoretisch altijd een betere oplossing vinden door je getal te vergroten.

Als je op zoek bent naar alle cirkelbogen met een maximale afwijking D lijkt het mij logisch je Deviation vast te prikken. Nu komen er op een gegeven moment geen grotere cirkelbogen meer bij omdat de minimale Deviation in je set van toegestaane stralen al is gevonden. Ook al hebben die cirkelbogen een kleinere afwijking dan sommige van je eerder gevonden cirkelbogen. Anderzijds, als je echt op zoek bent naar de cirkelboog met de minimale afwijking in je set, is het zinloos om al je tussenresultaten af te drukken. Daar wordt je algoritme namelijk onnodig traag van. Maar in beide gevallen snap ik niet wat je met een "betere boog dan 295661502 groot" bedoelt.

[JanWillem32]

Omdat het een beetje moeilijk is om voor te stellen wat dit programma eigenlijk produceert, heb ik even een voorbeeld gemaakt:

Rmax=1000000000
(0,0) (0,1) (-1,1) (0,0) (1,-1) (1,0) Deviation:0.707107
(0,0) (0,1) (1,1) (1,1) (1,1) (1,0) Deviation:0.292893
(0,0) (0,3) (1,3) (2,2) (3,1) (3,0) Deviation:0.12132
(0,0) (0,7) (3,7) (5,5) (7,3) (7,0) Deviation:0.0502525
(0,0) (0,17) (7,17) (12,12) (17,7) (17,0) Deviation:0.0208153
(0,0) (0,41) (17,41) (29,29) (41,17) (41,0) Deviation:0.00862197
(0,0) (0,99) (41,99) (70,70) (99,41) (99,0) Deviation:0.00357134
(0,0) (0,239) (99,239) (169,169) (239,99) (239,0) Deviation:0.0014793
(0,0) (0,577) (239,577) (408,408) (577,239) (577,0) Deviation:0.000612745
(0,0) (0,1393) (577,1393) (985,985) (1393,577) (1393,0) Deviation:0.000253807
(0,0) (0,3363) (1393,3363) (2378,2378) (3363,1393) (3363,0) Deviation:0.00010513
(0,0) (0,8119) (3363,8119) (5741,5741) (8119,3363) (8119,0) Deviation:4.35464e-005
(0,0) (0,19601) (8119,19601) (13860,13860) (19601,8119) (19601,0) Deviation:1.80375e-005
(0,0) (0,47321) (19601,47321) (33461,33461) (47321,19601) (47321,0) Deviation:7.47138e-006
(0,0) (0,114243) (47321,114243) (80782,80782) (114243,47321) (114243,0) Deviation:3.09475e-006
(0,0) (0,275807) (114243,275807) (195025,195025) (275807,114243) (275807,0) Deviation:1.28187e-006
(0,0) (0,665857) (275807,665857) (470832,470832) (665857,275807) (665857,0) Deviation:5.31007e-007
(0,0) (0,1.60752e+006) (665857,1.60752e+006) (1.13669e+006,1.13669e+006) (1.60752e+006,665857) (1.60752e+006,0) Deviation:2.19859e-007
(0,0) (0,3.8809e+006) (1.60752e+006,3.8809e+006) (2.74421e+006,2.74421e+006) (3.8809e+006,1.60752e+006) (3.8809e+006,0) Deviation:9.12885e-008
(0,0) (0,9.36932e+006) (3.8809e+006,9.36932e+006) (6.62511e+006,6.62511e+006) (9.36932e+006,3.8809e+006) (9.36932e+006,0) Deviation:3.72825e-008
(0,0) (0,2.26195e+007) (9.36932e+006,2.26195e+007) (1.59944e+007,1.59944e+007) (2.26195e+007,9.36932e+006) (2.26195e+007,0) Deviation:1.67238e-008
(0,0) (0,5.46084e+007) (2.26195e+007,5.46084e+007) (3.8614e+007,3.8614e+007) (5.46084e+007,2.26195e+007) (5.46084e+007,0) Deviation:3.83443e-009
(0,0) (0,2.41053e+008) (9.98475e+007,2.41053e+008) (1.7045e+008,1.7045e+008) (2.41053e+008,9.98475e+007) (2.41053e+008,0) Deviation:1.38243e-009
(0,0) (0,5.36715e+008) (2.22314e+008,5.36715e+008) (3.79515e+008,3.79515e+008) (5.36715e+008,2.22314e+008) (5.36715e+008,0) Deviation:1.07684e-009
(0,0) (0,7.77768e+008) (3.22162e+008,7.77768e+008) (5.49965e+008,5.49965e+008) (7.77768e+008,3.22162e+008) (7.77768e+008,0) Deviation:2.91038e-010

Dit een lijst net als de lijst met priemgetallen: zonder eind, maar dan met een veel spaarzamere spreiding. Bij Rmax=1000000000 geeft dat dus 24 antwoorden (die eerste is mijn fout; het commando "long double Deviation = 1;" hoort "long double Deviation = 0.5;" te zijn).
@Jan Klaasen: 14142135623730950488016887242097 is een aardige kandidaat, maar niet bijzonder waarschijnlijk. De decimale verschuivingen van ½√2 zijn in het voorbeeld enkel aanwezig in 1 en 1393 , terwijl de andere zeven decimale verschuivingen niet in deze lijst staan. En traag is het programma zeker niet: het voorbeeld genereren duurde slechts een minuutje. Daarbij, de verhouding antwoorden:loop cycli is 1:40000000, dus de tussenresultaten wegschijven is nou niet echt significant voor de totale berekeningstijd. Het is voor het maken van lettertypen van belang een basisvorm te hebben met zo min mogelijk afwijkingen. De uiteindelijke grootte op het scherm en op het papier kan toch verschaald worden van ongeveer 2 keer zo groot tot ongeveer 10000 keer zo klein (die verschalingen ondersteunen wel floating-point berekeningen). Daarom heb ik niets aan cirkels waarbij de precisie lager is dan bij die van een kleinere cirkel.

[MSalters]

Dit algoritme lijkt wel te kloppen, nu loop je inderdaad tegen de beperkte precisie van sqrt(2) aan. Het vervelende is dat als a=7E+10, dan is a-floor(a) niet heel erg precies meer.

Nu is het moment op Numerical Recipes of een dergelijk goed boek erbij te pakken. Ik vermoed dat de truc is om a te splitsen in a1+a2, waarbij geldt dat floor(a1)==0 en 0<=a2<=1. Het idee is dat je in elke stap 0.5*sqrt(2) bij a1+a2 optelt. Dat levert een cumulatieve fout op, als je niet corrigeert. Je weet echter dat (a1+a2)*(a1+a2) een geheel getal is. Daarom is ook (2*a1*a2+a2*a2) een gehele getal. Je hebt nu de grote term a1*a1 weggestreept. In de praktijk mis je elke keer 1 of 2 bits, maar door het uitrekenen van (2*a1*a2+a2*a2) kun je bepalen welke dat zijn.
Vervolgens is de bepaling van d triviaal, dat is namelijk òf a2 òf (1-a2).

De enige lastige truc is bepalen wat precies de aanpassing van a2 moet zijn, maar dat kan met Newton eenvoudig genoeg. De eerste schatting van a2 is goed genoeg, dus de eerstbepaalde waarde van (2*a1*a2+a2*a2) is waarschijnlijk minder dan 0.5 van de correcte waarde verwijderd. Noem die eerste schatting even Z, en zoek dan het nulpunt van (2*a1*a2+a2*a2)-Z

[JanWillem32]

@MSalters
de precisie van de functies met onder andere sqrt zijn gelimiteerd aan het gebruik van cmath. Daarom ben ik nu overschakeld op M_APM. De static library daarvan op een windows pc in elkaar zetten was een heel debakel, maar op een gegeven moment lukte dat toch.

C++:

#include <iostream>
#include "M_APM.H"
using namespace std;
int main() {
m_apm_cpp_precision(40);
char R[40];
char A[40];
char B[40];
char Deviation[40];
MAPM r, b;
MAPM c = .5*sqrt(2);
MAPM rmax = 1E40;
MAPM deviation=.5;
for (r = 1; rmax >= r; r+=2) {
  MAPM a = r*c;
  MAPM d = a-floor(a);
  if (d > .5) {
    d = 1-d;
    a = a+1;
  }
  if (d < deviation) {
    deviation = d;
    b = 2*floor(a)-r;
    r.toIntegerString(R);
    a.toIntegerString(A);
    b.toIntegerString(B);
    deviation.toString(Deviation, -1);
    cout << "(0,0) (0," << R << ") (" << B << "," << R << ") (" << A << "," << A << ") (" << R << "," << B << ") (" << R << ",0) Deviation=" << Deviation << endl;
  }
}
return 0;
}

Deze methode heeft een instelbare (hoge) precisie, maar is daardoor ook gelijk een stuk trager. Onderstaande code is behoorlijk sneller, en werkt tot Rmax=2,5×10⁸.

C++:

#include <iostream>
#include <cmath>
using namespace std;
int main() {
long double Rmax;
cout << "Rmax=";
cin >> Rmax;
long double Deviation = .5;
long double c = .5*sqrt(2);
for (long double R = 1; Rmax >= R; R+=2) {
  long double a = R*c;
  long double d = a-floor(a);
  if (d > .5) {
    d = 1-d;
    a = a+1;
  }
  if (d < Deviation) {
    Deviation = d;
    cout << "(0,0) (0," << R << ") (" << 2*floor(a)-R << "," << R << ") (" << floor(a) << "," << floor(a) << ") (" << R << "," << 2*floor(a)-R << ") (" << R << ",0) Deviation=" << Deviation << endl;
  }
}
return 0;
}

edit: @hvdberg - Goed idee! ; verwerkt

[ Voor 57% gewijzigd door JanWillem32 op 07-06-2006 18:24 ]

[Verwijderd]

Waarom reken je iedere iteratie .5*sqrt(2) uit? Als je die in een CONST zet rekent het waarschijnlijk nog sneller.

[MSalters]

JanWillem32 schreef op woensdag 07 juni 2006 @ 12:17:
@MSalters
de precisie van de functies met onder andere sqrt zijn gelimiteerd aan het gebruik van cmath.

Klopt. Daarom bereken ik per iteratie een foutterm. Aangezien dat kan met simpele vermenigvuldigingen en optellingen is de precisie van sqrt(2) alleen bepalend voor de snelheid van convergentie.

Daarom ben ik nu overschakeld op M_APM. De static library daarvan op een windows pc in elkaar zetten was een heel debakel, maar op een gegeven moment lukte dat toch.

C++:

#include <iostream>
#include "M_APM.H"
using namespace std;
int main() {
m_apm_cpp_precision(40);
char R[40];
char A[40];
char B[40];
char Deviation[40];
MAPM r, b;
MAPM c = .5*sqrt(2);
...
}

Deze methode heeft een instelbare (hoge) precisie, maar is daardoor ook gelijk een stuk trager. Onderstaande code is behoorlijk sneller, en werkt tot Rmax=2,5×10⁸.

Niet helemaal een hogere precisie, want je maakt in deze implementatie nog steeds dezelfde ongecorrigeerde fout met sqrt(2). Dat is namelijk geen MAPM functie, en de return waarde daarvan heeft dus waarschijnlijk maar 9 decimalen. De 40 decimalen in c zijn dus 9 echte en 31 random.

[MSalters]

Verwijderd schreef op woensdag 07 juni 2006 @ 16:18:
Waarom reken je iedere iteratie .5*sqrt(2) uit? Als je die in een CONST zet rekent het waarschijnlijk nog sneller.

Moderne compilers kunnen dat vermoedelijk optimaliseren. Meten is weten, maar is zou het niet als eerste veranderen.

[MSalters]

Op zch werkt het algoritme wel wat ik bedacht had om de precisie van sqrt(2) te compenseren. Alleen zit je daarin met een heel erg vergelijkbaar probleem van significante cijfers. Je blijkt dus daar vergelijkbare trucen uit te moeten halen met het splitsen van grote getallen in integer delen en fracties.

Overigens is het wel makkelijker om a2 in de range [-0.5, +0.5] te houden, dat is dan namelijk direct je Deviation term.

[MSalters]

Merkwaardig. Ik krijg (verschil vanaf Deviation < 3.09475e-006)

Rmax=1000000000
(0,0) (0,1) (-1,1) (0,0) (1,-1) (1,0) Deviation=0.292893
(0,0) (0,3) (1,3) (2,2) (3,1) (3,0) Deviation=0.12132
(0,0) (0,7) (1,7) (4,4) (7,1) (7,0) Deviation=0.0502525
(0,0) (0,17) (7,17) (12,12) (17,7) (17,0) Deviation=0.0208153
(0,0) (0,41) (15,41) (28,28) (41,15) (41,0) Deviation=0.00862197
(0,0) (0,99) (41,99) (70,70) (99,41) (99,0) Deviation=0.00357134
(0,0) (0,239) (97,239) (168,168) (239,97) (239,0) Deviation=0.0014793
(0,0) (0,577) (239,577) (408,408) (577,239) (577,0) Deviation=0.000612745
(0,0) (0,1393) (575,1393) (984,984) (1393,575) (1393,0) Deviation=0.000253807
(0,0) (0,3363) (1393,3363) (2378,2378) (3363,1393) (3363,0) Deviation=0.00010513

(0,0) (0,8119) (3361,8119) (5740,5740) (8119,3361) (8119,0) Deviation=4.35464e-0
05
(0,0) (0,19601) (8119,19601) (13860,13860) (19601,8119) (19601,0) Deviation=1.80
375e-005
(0,0) (0,47321) (19599,47321) (33460,33460) (47321,19599) (47321,0) Deviation=7.
47138e-006
(0,0) (0,114243) (47321,114243) (80782,80782) (114243,47321) (114243,0) Deviatio
n=3.09476e-006
(0,0) (0,275807) (114241,275807) (195024,195024) (275807,114241) (275807,0) Devi
ation=1.28187e-006
(0,0) (0,665857) (275807,665857) (470832,470832) (665857,275807) (665857,0) Devi
ation=5.3102e-007
(0,0) (0,1.60752e+006) (665855,1.60752e+006) (1.13669e+006,1.13669e+006) (1.6075
2e+006,665855) (1.60752e+006,0) Deviation=2.19829e-007
(0,0) (0,3.8809e+006) (1.60752e+006,3.8809e+006) (2.74421e+006,2.74421e+006) (3.
8809e+006,1.60752e+006) (3.8809e+006,0) Deviation=9.13624e-008
(0,0) (0,9.36932e+006) (3.8809e+006,9.36932e+006) (6.62511e+006,6.62511e+006) (9
.36932e+006,3.8809e+006) (9.36932e+006,0) Deviation=3.71039e-008
(0,0) (0,2.26195e+007) (9.36932e+006,2.26195e+007) (1.59944e+007,1.59944e+007) (
2.26195e+007,9.36932e+006) (2.26195e+007,0) Deviation=1.61547e-008
(0,0) (0,5.46084e+007) (2.26195e+007,5.46084e+007) (3.8614e+007,3.8614e+007) (5.
46084e+007,2.26195e+007) (5.46084e+007,0) Deviation=5.79456e-009
(0,0) (0,1.31836e+008) (5.46084e+007,1.31836e+008) (9.32224e+007,9.32224e+007) (
1.31836e+008,5.46084e+007) (1.31836e+008,0) Deviation=3.56553e-009
(0,0) (0,1.86445e+008) (7.72279e+007,1.86445e+008) (1.31836e+008,1.31836e+008) (
1.86445e+008,7.72279e+007) (1.86445e+008,0) Deviation=2.22903e-009
(0,0) (0,4.82106e+008) (1.99695e+008,4.82106e+008) (3.40901e+008,3.40901e+008) (
4.82106e+008,1.99695e+008) (4.82106e+008,0) Deviation=9.52823e-010
...

[JanWillem32]

Ik heb de MAPM bibliotheek expres er in gezet om sqrt uit te voeren. Kijk maar eens goed, er is geen cmath meer om anders sqrt uit te voeren. 'm_apm_cpp_precision(40);" bepaalt de minimumwaarde voor het aantal decimalen voor de uitkomst van elke berekinging, mits ondersteund door MAPM.
Daarnaast is deviation er alleen als indicatie voor de precisie, bij een radius van 5.46084×10⁷ wijkt bij het snelle programma de waarde voor 2.26195×10⁷ af (gezien vanuit een gegenereerd bestand om de antwoorden in te zetten).

[ Voor 32% gewijzigd door JanWillem32 op 08-06-2006 00:57 ]

[MSalters]

Voor de volledigheid:

C++:

void RootAdjust(long double a1, long double& a2);
void Overflow(long double& a1, long double& a2);

int main() { 
    long double Rmax; 
    cout << "Rmax="; 
    cin >> Rmax; 
    long double Deviation = 1.0; 
    long double root_1_2 = .5*sqrt(2.0L); 

    // Instead of keeping a in a single FP number, we use two where the first one
    // is really an integer and the second one is in the range [-.5, +.5). These
    // summed together form our real number. al is our "deviation from integer"
    long double ah = 0; 
    long double al = 0;

    for (long double R = 1; Rmax >= R; R+=1.0L) {
        // Increment (ah,al)
        al += root_1_2;
        Overflow(ah, al);
        // It should now be the root of R/2, but we might lose a LSB. Fix it.
        RootAdjust(ah, al);
        // Overflow after adjustment (quite unlikely, but better safe than sorry)
        Overflow(ah,al);

        if (abs(al) < Deviation) {
            Deviation = abs(al);
            long double a = ah + al; // loses precision only on output.
            cout << "(0,0) (0," << R << ") (" << 2*floor(a)-R << "," << R << ") (" << floor(a) << "," << floor(a) << ") (" << R << "," << 2*floor(a)-R << ") (" << R << ",0) Deviation=" << Deviation << endl; 
        } 
    } 
    return 0; 
}

void Overflow(long double& ah, long double& al)
{
    // We only add .7 at a time; -1 is enough.
    if (al >= 0.5L) {
        ah += 1.0L;
        al -= 1.0L;
    }
}

void RootAdjust(long double ah, long double& al)
{
    // 2 * (ah+al) * (ah+al) should be an integer. ah is integer, so (4*ah*al) + 2*al*al should be as well.
    // Error in al should be < 1e-9
    const long double error_goal = 1e-09;
    do {
        // Initial estimate - should be almost an integer. error1 is the difference, << 1
        // long double integer1 = (4.0L*ah*al) + (2.0L*al*al); // loses significant digts :(
        long double ih = (4.0L*ah*al); // If this is too big we'll drop significant digits here. too!
        long double il = (2.0L*al*al);
        // Obviously we can still save two bits from ih
        // ih+il is now the integer we're looking for, but adding them introduces additional error.
        // By splitting each term separately in integer+fraction we can add two fracton terms which
        // are >= 0 and < 1. This preserves significant digits. The result is in the range (-2,2) but
        // typically is 0.00x, 0.999x, 1.00x or 1.999x. Ignore that minor detail when determining the
        // real error (the absolute distance to the closest integer)
        long double error_i = abs(modf(ih, &ih) +  modf(il, &il)); // ih and il now hold integers, but they aren't needed.
        while (error_i > 0.5L) {
            error_i -= 1.0L;
        }
        long double delta_a = error_i / ((4.0L*ah) + (2.0L*al)); // See the formula for ih+il
        if (delta_a < error_goal) {
            return; // Wouldn't really change al
        }
        // Adjust al based on the size of the measured error.
        al -= delta_a;
    } while (true); // Breaks in loop.
}

Lijkt wel een 8086, ah+al

[JanWillem32]

Het is een tijdje geleden, maar ik heb toch nog een verbetering voor het programma in elkaar gezet. Deze is het allersnelste en wordt qua precisie enkel beperkt door de "1E-253". (Lagere exponenten worden niet meer geaccepteerd.) Ik heb namelijk een trucje uitgehaald, met een gegeven dat ik eerst niet zag: kijk maar naar regel 32.

C++:

#include <iostream>
#include "M_APM.H"
using namespace std;
int main() {
m_apm_cpp_precision(506);
char R[506];
char A[506];
char B[506];
char D[506];
MAPM a, b, d;
MAPM r = 1;
MAPM c = .5*sqrt(2);
MAPM deviation = .5;
MAPM f = sqrt(2)-1;
while (deviation > 1E-253) {
  a = r*c;
  d = a-floor(a);
  if (d > .5) {
    d = 1-d;
    a = ceil(a);
  }
  else {
    a = floor(a);
  }
  if (d < deviation) {
    deviation = d;
    r.toIntegerString(R);
    a.toIntegerString(A);
    b.toIntegerString(B);
    d.toString(D, 10);
    cout << "(0,0) (0," << R << ") (" << B << "," << R << ") (" << A << "," << A << ") (" << R << "," << B << ") (" << R << ",0) Deviation=" << D << endl;
    b = r;
    r = floor(r/f-1);
  }
  r += 1;
}
return 0;
}