mirror of
https://github.com/RobTillaart/Arduino.git
synced 2024-10-03 18:09:02 -04:00
268 lines
6.9 KiB
C++
268 lines
6.9 KiB
C++
#pragma once
|
||
//
|
||
// FILE: statHelpers.h
|
||
// AUTHOR: Rob Tillaart
|
||
// VERSION: 0.1.8
|
||
// PURPOSE: Arduino library with a number of statistic helper functions.
|
||
// DATE: 2020-07-01
|
||
// URL: https://github.com/RobTillaart/statHelpers
|
||
|
||
|
||
#include "Arduino.h"
|
||
|
||
|
||
#define STATHELPERS_LIB_VERSION (F("0.1.8"))
|
||
|
||
|
||
///////////////////////////////////////////////////////////////////////////
|
||
//
|
||
// PERMUTATIONS
|
||
//
|
||
uint32_t permutations(uint8_t n, uint8_t k);
|
||
|
||
uint64_t permutations64(uint8_t n, uint8_t k);
|
||
|
||
// can be optimized similar to dfactorial
|
||
double dpermutations(uint8_t n, uint8_t k);
|
||
|
||
|
||
/*
|
||
http://wordaligned.org/articles/next-permutation snippet
|
||
|
||
As an example consider finding the next permutation of:
|
||
|
||
8342666411
|
||
The longest monotonically decreasing tail is 666411, and the corresponding head is 8342.
|
||
|
||
8342 666411
|
||
666411 is, by definition, reverse-ordered, and cannot be increased by permuting its elements.
|
||
To find the next permutation, we must increase the head; a matter of finding the smallest tail
|
||
element larger than the head’s final 2.
|
||
|
||
8342 666411
|
||
Walking back from the end of tail, the first element greater than 2 is 4.
|
||
|
||
8342 666411
|
||
Swap the 2 and the 4
|
||
|
||
8344 666211
|
||
Since head has increased, we now have a greater permutation. To reduce to the next permutation,
|
||
we reverse tail, putting it into increasing order.
|
||
|
||
8344 112666
|
||
Join the head and tail back together. The permutation one greater than 8342666411 is 8344112666.
|
||
*/
|
||
|
||
|
||
// http://www.nayuki.io/page/next-lexicographical-permutation-algorithm
|
||
|
||
|
||
// b = nextPermutation<char>(array, 100);
|
||
template <typename T>
|
||
bool nextPermutation(T * array, uint16_t size)
|
||
{
|
||
// Find longest non-increasing suffix
|
||
int i = size - 1;
|
||
while (i > 0 && array[i - 1] >= array[i]) i--;
|
||
// Now i is the head index of the suffix
|
||
|
||
// Are we at the last permutation already?
|
||
if (i <= 0) return false;
|
||
|
||
// Let array[i - 1] be the pivot
|
||
// Find rightmost element that exceeds the pivot
|
||
int j = size - 1;
|
||
while (array[j] <= array[i - 1])
|
||
j--;
|
||
// Now the value array[j] will become the new pivot
|
||
// Assertion: j >= i
|
||
|
||
// Swap the pivot with j
|
||
T temp = array[i - 1];
|
||
array[i - 1] = array[j];
|
||
array[j] = temp;
|
||
|
||
// Reverse the suffix
|
||
j = size - 1;
|
||
while (i < j)
|
||
{
|
||
temp = array[i];
|
||
array[i] = array[j];
|
||
array[j] = temp;
|
||
i++;
|
||
j--;
|
||
}
|
||
return true;
|
||
}
|
||
|
||
|
||
|
||
///////////////////////////////////////////////////////////////////////////
|
||
//
|
||
// FACTORIAL
|
||
//
|
||
|
||
// exact ==> 12!
|
||
uint32_t factorial(uint8_t n);
|
||
|
||
// exact ==> 20!
|
||
uint64_t factorial64(uint8_t n);
|
||
|
||
// float (4 byte) => 34!
|
||
// double (8 byte) => 170!
|
||
double dfactorialReference(uint8_t n);
|
||
|
||
// FASTER VERSION
|
||
// does part of the math with integers.
|
||
// tested on UNO and ESP32, roughly 3x faster
|
||
// numbers differ slightly in the order of IEEE754 precision => acceptable.
|
||
// 10e-7 for 4 bit float
|
||
// 10e-16 for 8 bit double
|
||
double dfactorial(uint8_t n);
|
||
|
||
// stirling is an approximation function for factorial(n).
|
||
// it is slower but constant in time.
|
||
// float => 26!
|
||
// double => 143!
|
||
double stirling(uint8_t n);
|
||
|
||
|
||
// SEMIFACTORIAL
|
||
|
||
// exact ==> 20!!
|
||
uint32_t semiFactorial(uint8_t n);
|
||
|
||
// exact ==> 33!!
|
||
uint64_t semiFactorial64(uint8_t n);
|
||
|
||
// float (4 byte) => 56!!
|
||
// double (8 byte) => 300!!
|
||
double dSemiFactorial(uint16_t n);
|
||
|
||
|
||
|
||
///////////////////////////////////////////////////////////////////////////
|
||
//
|
||
// COMBINATIONS
|
||
//
|
||
|
||
// works for n = 0..30 for all k
|
||
uint32_t combinations(uint16_t n, uint16_t k);
|
||
|
||
// works for n = 0..61 for all k
|
||
uint64_t combinations64(uint16_t n, uint16_t k);
|
||
|
||
// experimental - not exact but allows large values
|
||
// float (4 bits) works till n = 125 for all k
|
||
// double (8 bits) works till n = 1020 for all k
|
||
double dcombinations(uint16_t n, uint16_t k);
|
||
|
||
// recursive (mind your stack and time)
|
||
// works for n = 0..30 for all k
|
||
// educational purpose
|
||
uint32_t rcombinations(uint16_t n, uint16_t k);
|
||
|
||
// recursive
|
||
// works for n = 0..61 for all k
|
||
// educational purpose
|
||
uint64_t rcombinations64(uint16_t n, uint16_t k);
|
||
|
||
// very slow double recursive way by means of Pascals triangle.
|
||
// works for n = 0..30 for all k (but takes a lot of time)
|
||
// educational purpose
|
||
uint32_t combPascal(uint16_t n, uint16_t k);
|
||
|
||
|
||
/////////////////////////////////////////////////
|
||
//
|
||
// EXPERIMENTAL
|
||
//
|
||
// BIG SECTION
|
||
//
|
||
// keep track of exponent myself in 32 bit unsigned integer
|
||
// - can be extended to a 64 bit integer
|
||
// however it already takes hours to calculate with 32 bits
|
||
|
||
/*
|
||
|
||
UNO
|
||
n n! millis() from earlier version
|
||
-----------------------------------------
|
||
1 1.00000e0 0
|
||
10 3.62880e6 1
|
||
100 9.33262e157 7
|
||
1000 4.02386e2567 105
|
||
2000 3.31627e5735 231
|
||
4000 1.82880e12673 504
|
||
8000 5.18416e27752 1086
|
||
10000 2.84625e35659 1389
|
||
16000 5.11880e60319 2330
|
||
32000 1.06550e130270 4978
|
||
100000 2.82428e456573 17201
|
||
1000000 8.26379e5565708 206421 (3.5 minutes)
|
||
10000000 an hour? too long...
|
||
|
||
ESP32 240MHz
|
||
n n! millis() from earlier version
|
||
-----------------------------------------
|
||
1 1.00000e0 0
|
||
10 3.62880e6 0
|
||
100 9.33262e157 0
|
||
1e3 4.02387e2567 8
|
||
1e4 2.84626e35659 110
|
||
1e5 2.82423e456573 1390
|
||
1e6 8.26393e5565708 16781
|
||
1e7 1.20242e65657059 196573 (3++ minutes)
|
||
1e8 1.61720e756570556 2253211
|
||
1e9 9.90463e4270738226 25410726 (7++ hrs to detect overflow! :(
|
||
|
||
|
||
1 0
|
||
10 0.6
|
||
100 1.57
|
||
1000 2.567
|
||
10000 3.5659
|
||
100000 4.56573
|
||
1000000 5.565708
|
||
10000000 6.5657059
|
||
100000000 7.56570556
|
||
|
||
// largest found - exponent is approaching max_uint32_t - 4294967296
|
||
518678058! 4.1873547e4294967283 // break condition was hit...
|
||
next one should fit too
|
||
518678059! 2.1718890e4294967292
|
||
|
||
*/
|
||
|
||
void bigFactorial(uint32_t n, double &mantissa, uint32_t &exponent);
|
||
|
||
|
||
// Should work full range for n = 518678059, k = { 0..n }
|
||
// and n > 518678059 => max k is definitely smaller ==> expect n+1 ==> k-15 at start
|
||
// performance: low. depends on k.
|
||
//
|
||
//
|
||
// for relative small k and large n (multiple orders of magnitude) one can get an estimate
|
||
// P(n,k) ~ raw = log10(n - k/2) * k;
|
||
// exponent = int(raw);
|
||
// mantissa = pow(10, raw - int(raw));
|
||
void bigPermutations(uint32_t n, uint32_t k, double &mantissa, uint32_t &exponent);
|
||
|
||
void bigCombinations(uint32_t n, uint32_t k, double &mantissa, uint32_t &exponent);
|
||
|
||
|
||
////////////////////////////////////////////////////////////
|
||
//
|
||
// EXPERIMENTAL 64 BIT
|
||
//
|
||
|
||
void bigFactorial64(uint64_t n, double &mantissa, uint64_t &exponent);
|
||
|
||
void bigPermutations64(uint64_t n, uint64_t k, double &mantissa, uint64_t &exponent);
|
||
|
||
void bigCombinations64(uint64_t n, uint64_t k, double &mantissa, uint64_t &exponent);
|
||
|
||
|
||
// -- END OF FILE --
|
||
|