Buffons Nadelwurf

[zeropage]

Hallo,
in einer Mathevorlesung kam zum Thema "Buffonsches Nadelproblem" dieser Python-Code:

from random import uniform
from math import pi, cos, floor

def button(n):
	return sum(floor(uniform(0,1)+cos(uniform(0, pi))) != 0
		for _ in range(n) / n


buffon(1000000)

Den habe ich nach C++ übersetzt, raus kam dabei das:

#include <iostream>
#include <cmath>
#include <random>

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

float buffon(std::mt19937& rndGen, const int N)
{
	float sum = 0;
	for (int i = 0; i < N; ++i)
	{
		std::uniform_real_distribution<float> rnd1(0, 1);
		std::uniform_real_distribution<float> rnd2(0, PI);
		sum += std::abs(std::floorf(rnd1(rndGen) + std::cosf(rnd2(rndGen))));
	}
	return sum / N * 100;
}

int main()
{
	std::random_device rd;
	std::mt19937 rndGen(rd());
	
	for (int i = 0; i < 10; ++i)
	{
		std::cout << buffon(rndGen, 1000000) << '\n';
	}

}

Die Ergebnisse sind die zu erwartenden ~ 63,7 %. Aber habe ich das richtig übersetzt, oder eher zu wortwörtlich?
Und muss ich wirklich den RandomGenerator der Funktion übergeben?

[wob]

Dein Python-Code ist kein korrektes Python (invalid Syntax). Dein C++-Code gibt aber was anderes aus...

In C++ rechnest du einen Absolutbetrag aus, in Python vergleichst du irgendwas Abgerundetes mit 0 in der List comprehension, die aber keine Klammer zu hat. Außerdem wird die Python-Funktion button, wenn sie denn korrekt wäre, gar nicht aufgerufen, sondern eine völlig fehlende Funktion buffon.

Ok, wenn ich ff durch tt ersetze und die Klammer einfüge, sind das trotzdem verschiedene Dinge:
Der eine Code macht Abrunden + Vergleich mit 0 und aufsummierten der Fälle, wo es != 0 ist, der andere Code hat eine float-Summe und summiert float-Absolutbeträge.

[zeropage]

Aha. Den Python Code habe ich von hier:

Mathematische Eleganz hoch drei (Weihnachtsvideo 2021) bei Zeit 12:33.

Ist natürlich keine echte Vorlesung, nur ne Weihnachtsvorlesung ohne echte Zuschauer. Ich dachte, der wird das schon richtig machen, kenne mich mit Python sonst gar nicht aus.

Durch den Absolutbetrag wollte ich um den Vergleich != 0 herumkommen.

edit: das button im Python-Code war mein Abtippfehler...

[zeropage]

@zeropage sagte in Buffons Nadelwurf:

Ich dachte, der wird das schon richtig machen, kenne mich mit Python sonst gar nicht aus.

Oder es war Absicht, denn er hatte öfters extra erwähnt, das er gerne mal Fehler macht?

[wob]

Also, ich habe ja nicht gesagt, dass die Codes nicht zum selben Ergebnis kommen, aber der Rechenweg ist eben anders bzw. du hast Umformungen vorgenommen.

[zeropage]

Ja, wegen diesen Umformungen habe ich gefragt. Ich hab das mehr oder weniger "mechanisch" runter getippt, ohne "kreativität" wie es bei den Videos von oben immer betont wird.
Wollte ich fragen, ob das trotzdem so ok ist.

bzw: war dann meine Fragestellung falsch. Wie würde man diesen Python-Code ohne Umformungen übersetzen?

[wob]

Ich denke, sollte ok sein.

Ich gehe immer gerne daher und plotte den Kram.

Also z.B. so (in einem Jupyter Notebook):
pd.Series([uniform(0, 1) + cos(uniform(0, pi)) for _ in range(1000000)]).hist(bins=50)
Dann kann man sich leichter überlegen, warum dein Algorithmus auch geht.

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