#pragma once // // FILE: statHelpers.h // AUTHOR: Rob Tillaart // VERSION: 0.2.0 // 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.2.0")) /////////////////////////////////////////////////////////////////////////// // // 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(array, 100); template 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); // SEMI_FACTORIAL // 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); // SKIP_FACTORIAL (experimental) // note step must be larger than 0 // note when step == 1 ==> factorial // note when step == 2 ==> semiFactorial uint32_t skipFactorial(uint32_t n, uint32_t skip); uint64_t skipFactorial64(uint32_t n, uint32_t skip); double dSkipFactorial(uint32_t n, uint32_t skip); /////////////////////////////////////////////////////////////////////////// // // 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 --