Buffons Nadelwurf

[zeropage]

PS:
Hier noch meine Ergebnisse, um nicht vermuten zu müssen (im Release mit (/O2))

<float>
V1: 63.6544  number: 1000000  time: 136.928 milliSec  time/number: 0.136928 microSec
V1: 63.6435  number: 1000000  time: 135.168 milliSec  time/number: 0.135168 microSec
V1: 63.6315  number: 1000000  time: 137.597 milliSec  time/number: 0.137597 microSec

V2: 63.6136  number: 1000000  time: 149.032 milliSec  time/number: 0.149032 microSec
V2: 63.5668  number: 1000000  time: 150.392 milliSec  time/number: 0.150392 microSec
V2: 63.5212  number: 1000000  time: 149.777 milliSec  time/number: 0.149777 microSec

V3: 63.6275  number: 1000000  time: 150.238 milliSec  time/number: 0.150238 microSec
V3: 63.6409  number: 1000000  time: 148.403 milliSec  time/number: 0.148403 microSec
V3: 63.5773  number: 1000000  time: 146.317 milliSec  time/number: 0.146317 microSec

<double>
V1: 63.6624  number: 1000000  time: 130.27 milliSec  time/number: 0.13027 microSec
V1: 63.5813  number: 1000000  time: 129.87 milliSec  time/number: 0.12987 microSec
V1: 63.5432  number: 1000000  time: 130.164 milliSec  time/number: 0.130164 microSec

V2: 63.5812  number: 1000000  time: 141.18 milliSec  time/number: 0.141181 microSec
V2: 63.6809  number: 1000000  time: 140.912 milliSec  time/number: 0.140912 microSec
V2: 63.644  number: 1000000  time: 140.453 milliSec  time/number: 0.140453 microSec

V3: 63.6482  number: 1000000  time: 139.01 milliSec  time/number: 0.13901 microSec
V3: 63.614  number: 1000000  time: 139.286 milliSec  time/number: 0.139287 microSec
V3: 63.7084  number: 1000000  time: 138.801 milliSec  time/number: 0.138801 microSec

[hustbaer]

Hm. Mit https://quick-bench.com bekomm ich da andere Ergebnisse: https://quick-bench.com/q/kk2E4TjABCANWBeAslK8gDowgS4
EDIT: V5 in dem Link ist falsch, es müsste

    auto const ival = long(val + 1);
    sum += (ival != 1);

heissen. /EDIT

Mit Clang 12 -O3 sind alle Varianten inetwa gleich schnell, und float ist konsistent ein gutes Stück schneller als double. Das entspricht ziemlich genau dem was ich erwartet hätte. Das ganze wird vermutlich von der cos Berechnung limitiert, wie man den Rest macht sollte keine grosse Rolle spielen.

Wobei ich zwei Dinge angepasst habe: ich hab die Definition von PI nach buffonVal reingezogen, den Typ zu T geändert und nen Literal statt acosf(-1) verwendet. Und ich hab N deutlich kleiner gemacht.

Mit GCC 10.3 (ebenfalls -O3) sieht man ein paar Unterschiede zwischen den Varianten: V2 ist ein wenig langsamer und V3 ist ein gutes Stück langsamer. Und V4 und V5 die ich dazugebastelt habe sind eine Spur schneller als V1. Ansonsten auch hier float deutlich schneller als double.

Der von GCC generierte Code ist insgesamt viel schneller als der von Clang generierte, was vermutlich der Grund ist dass hier die Unterschiede zwischen den Varianten grösser sind.

[hustbaer]

Nachtrag: wie man sich täuschen kann. Hab das ganze nochmal mit __attribute__((always_inline)) probiert und mir die Disassembly angesehen - wo man dann sieht wo die ganze Zeit draufgegangen ist. Und die geht bei diversen Divisionen und Konvertierungen drauf. cos scheint quasi keine Rolle zu spielen. uniform_real_distribution?

[hustbaer]

OK. Hab das ganze jetzt nochmal abgeändert: RNG durch eine Xoroshiro Variante ersetzt und die Distribution selbst implementiert. Das macht die Sache ca. 8x schneller (mit Clang). (Meine "V5" hab ich auch repariert.)

https://quick-bench.com/q/tS_6RYW_h4UWseO3Qlnl3skwCyM

Jetzt scheint auch endlich cos im Profiling auf, wenn auch nur mit ~8% (BM_Float_V4).

Code:

#include <cmath>
#include <random>

#define FORCE_INLINE __attribute__((always_inline))

template <typename T>
class Random {

  FORCE_INLINE uint64_t rotl(const uint64_t x, int k) {
    return (x << k) | (x >> (64 - k));
  }

  uint64_t s[2];

  FORCE_INLINE uint64_t next(void) {
    const uint64_t s0 = s[0];
    uint64_t s1 = s[1];
    const uint64_t result = s0 + s1;

    s1 ^= s0;
    s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b
    s[1] = rotl(s1, 37); // c

    return result;
  }

public:
  Random() {
    std::random_device rd;
    s[0] = uint64_t(rd());
    s[1] = uint64_t(rd());
  }

	FORCE_INLINE T uniform1() {
    static constexpr T scale = T(1) / 0x1p64;
		return T(next()) * scale;
	}

	FORCE_INLINE T uniformPi() {
    static constexpr T scale = T(3.14159265358979323846264) / 0x1p64;
		return T(next()) * scale;
	}
};

template <typename T>
FORCE_INLINE T buffonVal(Random<T>& rnd) {
	return rnd.uniform1() + std::cos(rnd.uniformPi());
}

template <typename T>
FORCE_INLINE T buffonV1(Random<T>& rnd, const int N) { 
	T sum = 0;
	for (int i = 0; i < N; ++i) {
		auto val = buffonVal(rnd);
		sum += (val < 0) | (val >= 1);
	}
	return sum / N * 100;
}

template <typename T>
FORCE_INLINE T buffonV2(Random<T>& rnd, const int N) {
	T sum = 0;
	for (int i = 0; i < N; ++i)
		sum += std::floor(buffonVal(rnd)) !=0;
	return sum / N * 100;
}

template <typename T>
FORCE_INLINE T buffonV3(Random<T>& rnd, const int N) {
	T sum = 0;
	for (int i = 0; i < N; ++i)
		sum += std::abs(std::floor(buffonVal(rnd)));
	return sum / N * 100;
}

template <typename T>
FORCE_INLINE T buffonV4(Random<T>& rnd, const int N) {
  long sum1 = 0;
  long sum2 = 0;
	for (int i = 0; i < N; ++i) {
    auto const val = buffonVal(rnd);
    if (val < 0)
      sum1++;
    if (val >= 1)
      sum2++;
  }
	return T(sum1 + sum2) / N * 100;
}

template <typename T>
FORCE_INLINE T buffonV5(Random<T>& rnd, const int N) {
  long sum = 0;
  for (int i = 0; i < N; ++i) {
    auto const val = buffonVal(rnd);
    auto const ival = long(val + 1);
    sum += (ival != 1);
  }
  return T(sum) / N * 100;
}

template <class T, class F>
FORCE_INLINE void bench(benchmark::State& state, F func) {
  Random<T> rng;
  for (auto _ : state) {
    benchmark::DoNotOptimize(func(rng));
  }
}

static constexpr int N = 10000;

void BM_Double_V1(benchmark::State& state) {
  bench<double>(state, [](auto& rng){ return buffonV1(rng, N); });
}
BENCHMARK(BM_Double_V1);

void BM_Double_V2(benchmark::State& state) {
  bench<double>(state, [](auto& rng){ return buffonV2(rng, N); });
}
BENCHMARK(BM_Double_V2);

void BM_Double_V3(benchmark::State& state) {
  bench<double>(state, [](auto& rng){ return buffonV3(rng, N); });
}
BENCHMARK(BM_Double_V3);

void BM_Double_V4(benchmark::State& state) {
  bench<double>(state, [](auto& rng){ return buffonV4(rng, N); });
}
BENCHMARK(BM_Double_V4);

void BM_Double_V5(benchmark::State& state) {
  bench<double>(state, [](auto& rng){ return buffonV5(rng, N); });
}
BENCHMARK(BM_Double_V5);

void BM_Float_V1(benchmark::State& state) {
  bench<float>(state, [](auto& rng){ return buffonV1(rng, N); });
}
BENCHMARK(BM_Float_V1);

void BM_Float_V2(benchmark::State& state) {
  bench<float>(state, [](auto& rng){ return buffonV2(rng, N); });
}
BENCHMARK(BM_Float_V2);

void BM_Float_V3(benchmark::State& state) {
  bench<float>(state, [](auto& rng){ return buffonV3(rng, N); });
}
BENCHMARK(BM_Float_V3);

void BM_Float_V4(benchmark::State& state) {
  bench<float>(state, [](auto& rng){ return buffonV4(rng, N); });
}
BENCHMARK(BM_Float_V4);

void BM_Float_V5(benchmark::State& state) {
  bench<float>(state, [](auto& rng){ return buffonV5(rng, N); });
}
BENCHMARK(BM_Float_V5);

[hustbaer]

Nochmal Nachtrag: Clang ist bloss mit libstdc++ so langsam. Mit libc++ ist Clang sogar schneller als GCC mit libstdc++.

Naja, egal. Was man hier auf jeden Fall mitnehmen kann ist dass Compiler und Standard Library hier viel mehr ausmachen als alles andere - zumindest wenn man RNG und Distribution der Standard Library verwendet.

[zeropage]

@hustbaer sagte in Buffons Nadelwurf:

Naja, egal.

Ich find das schon alles mächtig interessant

[hustbaer]

@zeropage
Ist es eh. Nur könnte man damit sicher noch Stunden verbringen wenn man das alles noch genauer analysieren wollte. Und die Zeit mag ich mir grad nicht nehmen

[wob]

In meinen "Benchmarks" in Python hatte ich den random-Teil extra vorberechnet und daher nicht mit in mein v1/v2/v3 reingenommen - ansonsten hätte ich auch größere Schwankungen gehabt und vor allem wollte ich nicht den rng testen, sondern die verschiedenen Berechnugsarten.

Sicherlich kann man auch Argumente für den anderen Weg finden. Allerdings muss man in Python immer so denken, dass man eine Berechnung für viele Elemente macht (mit numpy oder pandas-Funktionen) - wenn man manuelle Schleifen macht, wirds langsam.

[zeropage]

Ich habe noch eine Frage zu der Konstante PI.

Wenn ich die als globale Konstante haben möchte, kann ich dann einfach zB

namespace math
{
	template <typename T> static constexpr T PI = T(3.14159265358979323846264);
        //... evtl weitere eigene Funktionen...
}

übernehmen?

Aufgerufen würde das dann so:

// einzelner Nadelwurf
template <typename T> T buffonVal(Random<T>& rnd) 
{ 
	return rnd.uniform(0, 1) + std::cos(rnd.uniform(0, math::PI<T>)); 
}

Funktionieren tut es offenbar, aber ist es richtig?

[5cript]

@zeropage
pffff, das ist offensichtlich viel zu einfach für eine C++ Lösung für Konstanten.

<joke>

#include <iostream>
#include <iomanip>
#include <limits>

#define CAST(X) static_cast <T> (X)

template <typename T>
constexpr T sqrt (T of, T x = CAST(1), int rbi = 15) {
    return (rbi == 1) ? x : sqrt(of, x - (x * x - of) / (2. * x), rbi - 1);
}

namespace internal {
    constexpr long long factorial (long long index) {
        return (index != 1LL && index != 0LL) ? factorial(index - 1LL) * index : 1LL;
    }

    template <typename T>
    constexpr T euler_constant_estimator(int rbi = 20) {
        return (rbi == 1) ? CAST (1) :
            euler_constant_estimator <T> (rbi - 1) +
            CAST (rbi % 2 ? CAST (1LL) : CAST (-1LL)) /
            CAST (factorial(CAST (rbi-1)))
        ;
    }

    template <typename T>
    constexpr T sqrt_2_estimator(T previous, int rbi = 20) {
        return (rbi == 1) ? previous : sqrt_2_estimator((previous + (CAST(2) / previous))/CAST(2), rbi-1);
    }
}

constexpr long double e = 1.L / internal::euler_constant_estimator <long double> ();
constexpr long double sqrt_2 = internal::sqrt_2_estimator <long double> (1.);

namespace internal {
    template <typename T>
    constexpr T pi_estimator(T a = CAST(1), T b = CAST(1) / CAST(sqrt_2), T t = CAST(1)/CAST(4), T p = CAST(1), int rbi = 5) {
        return (rbi == 1) ? ((((a+b)/ CAST(2)) + sqrt(a*b)) * (((a+b)/ CAST(2)) + sqrt(a*b))) / (4 * (t - p * ( (a - (a + b) / CAST(2)) * (a - (a + b) / CAST(2)) )))
            : pi_estimator(/* a */ (a+b)/ CAST(2),
                            /* b */ sqrt(a*b),
                            /* t */ t - p * ( (a - (a + b) / CAST(2)) * (a - (a + b) / CAST(2)) ),
                            /* p */ CAST(2) * p,
                            rbi - 1)
        ;
    }
}

constexpr long double pi = internal::pi_estimator <long double>();

int main() {
    std::cout << std::numeric_limits<long double>::digits10 << '\n';
    std::cout << std::setprecision(21) << pi << '\n';
}

</joke>

[hustbaer]

@zeropage sagte in Buffons Nadelwurf:

Ich habe noch eine Frage zu der Konstante PI.

Wenn ich die als globale Konstante haben möchte, kann ich dann einfach zB

namespace math
{
	template <typename T> static constexpr T PI = T(3.14159265358979323846264);
        //... evtl weitere eigene Funktionen...
}

übernehmen?

Denke schon.
Wobei ich dazusagen sollte dass ich hier einfach die ersten paar Dezimalstellen von einer Webseite kopiert habe - ohne nachzusehen wie viele überhaupt Sinn machen oder auch nur zu zählen wie viele Stellen ich kopiert habe.

Und vermutlich wäre es sinnvoll ein L anzuhängen. Sonst könnte bei PI<long double> u.U. Unsinn rauskommen.

[zeropage]

Ich hatte mal gelesen, um den größten denkbaren Kreis (Radius des beobachtbaren Universum) mit der größten denkbaren Genauigkeit (Planck-Länge) zu berrechnen, benötigt man nur eine Handvoll (6 oder 7 oder so) Dezimal-Stellen von PI.

Find das aber nicht mehr.

Ich würde meinen, eben mindestens soviele wie in einem zB <long double> reinpassen?

[Quiche Lorraine]

@zeropage sagte in Buffons Nadelwurf:

Ich hatte mal gelesen, um den größten denkbaren Kreis (Radius des beobachtbaren Universum) mit der größten denkbaren Genauigkeit (Planck-Länge) zu berrechnen, benötigt man nur eine Handvoll (6 oder 7 oder so) Dezimal-Stellen von PI.

Wie soll das funktionieren wo der Dynamikbereich astronomisch groß ist?

Ich habe mal nachgerechnet. Der Unfang der Milchstraße dürfte cirka 1,189E19 m betragen. Die Planklänge beträgt 1,616E-35 m. Im Worse Case ist der Umfang 1,189E19 + 1.616E-35 groß. Alleine um die Zahl darzustellen benötigt man cirka 62 Stellen.

Da müsste man eigentlich mal die Fehlerfortpflanzung durchführen.

Edit: Habe mich beim Umfang verrechnet.

[SeppJ]

@zeropage sagte in Buffons Nadelwurf:

Ich hatte mal gelesen, um den größten denkbaren Kreis (Radius des beobachtbaren Universum) mit der größten denkbaren Genauigkeit (Planck-Länge) zu berrechnen, benötigt man nur eine Handvoll (6 oder 7 oder so) Dezimal-Stellen von PI.

Zweifel. Ich weiß zwar nicht, was du mit "berechnen" hier meinst, aber zwischen Radius des Universums und Plancklänge liegen gut 60 Größenordnungen. Hast du vielleicht eine 0 veregessen? 60-70 Stellen klingt realistischer und ist immer noch wenig gegenüber den Leuten, die als Hobby Pi auf 1000 Stellen auswendig lernen.

[zeropage]

Das wird wohl wirklich nicht hinhauen. Wenn ich nur wüsste, woher ich das habe...

[SeppJ]

@zeropage sagte in Buffons Nadelwurf:

Das wird wohl wirklich nicht hinhauen. Wenn ich nur wüsste, woher ich das habe...

Es könnte ungefähr hinhauen, wenn man die Unkenntnis über die Größen der Plancklänge und des beobachtbaren Universums mit betrachtet. Denn die Plancklänge ist tatsächlich nur auf ein paar Stellen genau bekannt, und das Universum ist nur ganz grob geschätzt vermessen (Pi mal Daumen, haha ). Das ist dann auch die alte Weisheit über double precision Zahlen: Die sind viel genauer als alles was die Menschheit je vermessen hat, und von daher ist es immer mit Vorsicht zu genießen, wenn jemand behauptet, doubles seien für irgendetwas nicht genau genug.

PS: Der Rekord für genaueste Messung ist, wenn meine Information noch aktuell ist, das magnetische Moment des Elektrons. Das steht bei einer Genauigkeit von knap 10^-12. Und dieser außergewöhnliche Rekord ist immer noch 2-3 Größenordnungen vom double entfernt.

[wob]

Nur zur allgemeinen Info: die aktuell besten bekannten Werte der Naturkonstanten findet man bei der Particle Data Group. Für @SeppJ schaut man unter Leptonen, ansonsten guckt man immer unter Shortcuts->Physical constants.

Magnetic moment anomaly
(g−2)/2 = (1159.65218091±0.00000026)×10^(-6)

[SeppJ]

Aber das ist doch weniger genau als das magnetische Moment des Elektrons!? Muss es ja auch, schließlich ist das anomale Moment die Abweichung des Messwerts von dem, was man (naiv) erwarten würde. Gewissermaßen kann man natürlich sagen, es sei die gleiche Genauigkeit, wenn man alles in Einheiten der naiven Erwartung angibt. Und das sehe ich so auch ein. Macht den Case für die mehr als ausreichende Genauigkeit des double um so besser, denn nun brauchen wir nur 10 Stellen, um den den größten Triumph der Mess- und Rechentechnik abzuspeichern.

[zeropage]

Spannend ist, das es insgesamt ^80 Größenordnungen sind, aber von unseren Maßen (ein Meter, eine Sekunde und so) in größer und kleiner Richtung ^40 Größenordnungen, also die Hälfte sind.

Das habe ich mal in einem Astronomie Buch von 1969 gelesen. Ohne nachzurechen denke ich, diese Werte gelten heute auch noch?

[hustbaer]

Hä?
Das beobachtbare Universum ist ca. 10^27m gross (Durchmesser).
Die Planck-Länge ist 1.6... * 10^-35m.
Das ergibt in Summe weder 80, noch sind es ähnliche viele Grössenordnungen.