Buffons Nadelwurf

[zeropage]

Ja, dieses Jupyter Notebook nutzt dieser Prof in seinen regulären Vorlesungen auch. Jedenfalls früher, er macht kaum noch welche.

[wob]

Ich habe mal aus Spaß einen Test gemacht:

import timeit
import numpy as np

n = 1000000
rng = np.random.default_rng()  # das ist ein PCG64, kein MT19937
val = rng.uniform(0, 1, n) + np.cos(rng.uniform(0, np.pi, n))
# val = val.astype(np.float32)


def v1():
    return np.sum((val < 0) | (val >= 1)) / n

def v2():
    return np.sum(np.floor(val) != 0) / n

def v3():
    return np.sum(np.abs(np.floor(val))) / n


print("v1", timeit.timeit(v1, number=1000))
print("v2", timeit.timeit(v2, number=1000))
print("v3", timeit.timeit(v3, number=1000))

Da kommt bei mir relativ konsistent (+-0.1s) raus (Einheit ist Sekunden für die je 1000 Calls, die je auf der 1000000 Elemente großen random+cos()-Summe arbeiten):

v1 1.6318554740282707
v2 2.3389774310053326
v3 3.911751442006789

D.h. Variante 1 mit 2 Vergleichen ist schneller als Variante 2 (das Original) und deine Variante ist die langsamste.

Allerdings wird hier mit double bzw np.float64 gearbeitet. Spannend ist es auch, oben val = val.astype(np.float32) einzukommentieren - plötzlich ist deine Variante dann schnell:

(mit float32)
v1 1.0741378079983406
v2 1.216420892975293
v3 1.0514918970293365

Aber hier ist der Unterschied zwischen v1 und v3 sehr klein und ich sollte ggf. mehr Iterationen machen. v2 ist aber konsistent am langsamsten mit float32.

Kannst ja auch mal mit C++ vergleichen. Jedenfalls sieht der Code mit numpy gleich dem in C++ ähnlicher, zumindest was den rng betrifft.

[zeropage]

Rein von der Messung ist bei mir die Version1 mit double am schnellsten. Aber ich bin mir nicht sicher, ob ich korrekt übersetzt habe.

Gemessen habe ich mit

#pragma once
#include <chrono>

// ist das sinnvoll, hier ein template zu nehmen, oder sollte man sich auf float oder double festlegen?
template <typename T> class Timer
{
public:
    using clock_t = std::chrono::high_resolution_clock;
    using duration_t = std::chrono::duration<T>;

private:
    clock_t::time_point startTime;
    duration_t elapsedTime = {};

public:
    Timer() : startTime(clock_t::now()) {}

    void reset() { startTime = clock_t::now(); }
    void stop() { elapsedTime = clock_t::now() - startTime; }

    T seconds() const { return elapsedTime.count(); }
    T milliSeconds() const { return elapsedTime.count() * 1000; }
    T microSeconds() const { return elapsedTime.count() * 1000000; }
};

Der sonstige Code ist dabei recht umfangreich geworden

#include <iostream>
#include <cmath>
#include <random>
#include "Timer.h"

const float PI = std::acosf(-1);

// Hilfsklasse um random-Werte zu generieren
template <typename T> class Random
{
	std::mt19937 randomGen;

public:
	Random(const std::mt19937& randomGen) : randomGen(randomGen) {}

	T uniform(const T begin, const T end)
	{
		std::uniform_real_distribution<T> random(begin, end);
		return random(randomGen);
	}
};

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

// Varianten für N Nadelwürfe
template <typename T> 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> 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> 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;
}

// Varianten in float oder double
template <typename T> void buffon(Random<T>& rnd, const int N, const int tries)
{
	Timer<float> timer;
	for (int i = 0; i < tries; ++i)
	{
		timer.reset();
		auto buffon = buffonV1(rnd, N);
		timer.stop();

		auto elapsedTimeMilliSec = timer.milliSeconds();
		auto elapsedTimeMicroSec = timer.microSeconds();
		std::cout << "V1: " << buffon << "  number: " << N << "  time: " << elapsedTimeMilliSec << " milliSec  time/number: " << elapsedTimeMicroSec / N << " microSec\n";
	}
	std::cout << '\n';
	for (int i = 0; i < tries; ++i)
	{
		timer.reset();
		auto buffon = buffonV2(rnd, N);
		timer.stop();

		auto elapsedTimeMilliSec = timer.milliSeconds();
		auto elapsedTimeMicroSec = timer.microSeconds();
		std::cout << "V2: " << buffon << "  number: " << N << "  time: " << elapsedTimeMilliSec << " milliSec  time/number: " << elapsedTimeMicroSec / N << " microSec\n";
	}
	std::cout << '\n';
	for (int i = 0; i < tries; ++i)
	{
		timer.reset();
		auto buffon = buffonV3(rnd, N);
		timer.stop();

		auto elapsedTimeMilliSec = timer.milliSeconds();
		auto elapsedTimeMicroSec = timer.microSeconds();
		std::cout << "V3: " << buffon << "  number: " << N << "  time: " << elapsedTimeMilliSec << " milliSec  time/number: " << elapsedTimeMicroSec / N << " microSec\n";
	}
}

int main()
{
	std::random_device rd;
	std::mt19937 randomGen(rd());

	const int N = 1000000;

	Random<float> rndf(randomGen);
	std::cout << "<float>\n";
	buffon(rndf, N, 3);

	std::cout << '\n';
	Random<double> rndd(randomGen);
	std::cout << "<double>\n";
	buffon(rndd, N, 3);
}

[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)