math.h

View source code here on GitHub!

Includes

uintmax_t factorial(uint32_t n)

Warning

This function only works for numbers smaller than MAX_FACTORIAL_64 or MAX_FACTORIAL_128, depending on the size of uintmax_t.

uintmax_t n_choose_r(uint32_t n, uint32_t r)

Returns -1 if there is an overflow. Otherwise returns n choose r.

Warning

This function only works for numbers smaller than MAX_FACTORIAL_64 or MAX_FACTORIAL_128, depending on the size of uintmax_t.

Note

The following are only available or necessary for the pcc compiler, as it has a hard time including <stdlib.h>.

uint8_t imprecise_log10(uintmax_t x)
double sqrt(double S)
uintmax_t ceil(double x)
  1#ifndef EULER_MATH_H
  2#define EULER_MATH_H
  3
  4#include "macros.h"
  5
  6#if !PCC_COMPILER
  7    #include <stdlib.h>
  8#endif
  9
 10#include <stdint.h>
 11
 12uintmax_t factorial(uint32_t n);
 13inline uintmax_t factorial(uint32_t n) {
 14    // note that this function only works for numbers smaller than MAX_FACTORIAL_64
 15    if ((sizeof(uintmax_t) == 8 && n > MAX_FACTORIAL_64) || (sizeof(uintmax_t) == 16 && n > MAX_FACTORIAL_128))
 16        return -1;
 17    uintmax_t ret = 1;
 18    for (uint32_t i = 2; i <= n; ++i)
 19        ret *= i;
 20    return ret;
 21}
 22
 23uintmax_t n_choose_r(uint32_t n, uint32_t r) {
 24    // function returns -1 if it overflows
 25    if ((sizeof(uintmax_t) == 8 && n <= MAX_FACTORIAL_64) || (sizeof(uintmax_t) == 16 && n <= MAX_FACTORIAL_128))
 26        return factorial(n) / factorial(r) / factorial(n-r);  // fast path if small enough
 27    // slow path for larger numbers
 28    int *factors;
 29    uintmax_t answer, tmp;
 30    uint32_t i, j;
 31    factors = (int *) malloc(sizeof(int) * (n + 1));
 32    // collect factors of final number
 33    for (i = 2; i <= n; i++)
 34        factors[i] = 1;
 35    // negative factor values indicate need to divide
 36    for (i = 2; i <= r; i++)
 37        factors[i] -= 1;
 38    for (i = 2; i <= n - r; i++)
 39        factors[i] -= 1;
 40    // this loop reduces to prime factors only
 41    for (i = n; i > 1; i--)
 42        for (j = 2; j < i; j++)
 43            if (i % j == 0) {
 44                factors[j] += factors[i];
 45                factors[i / j] += factors[i];
 46                factors[i] = 0;
 47                break;
 48            }
 49    i = j = 2;
 50    answer = 1;
 51    while (i <= n) {
 52        while (factors[i] > 0) {
 53            tmp = answer;
 54            answer *= i;
 55            while (answer < tmp && j <= n) {
 56                while (factors[j] < 0) {
 57                    tmp /= j;
 58                    factors[j]++;
 59                }
 60                j++;
 61                answer = tmp * i;
 62            }
 63            if (answer < tmp)
 64                return -1;  // this indicates an overflow
 65            factors[i]--;
 66        }
 67        i++;
 68    }
 69    while (j <= n) {
 70        while (factors[j] < 0) {
 71            answer /= j;
 72            factors[j]++;
 73        }
 74        j++;
 75    }
 76    free(factors);
 77    return answer;
 78}
 79
 80#if PCC_COMPILER
 81    uint8_t imprecise_log10(uintmax_t x);
 82    inline uint8_t imprecise_log10(uintmax_t x) {
 83        uint8_t answer = 0;
 84        while (x) {
 85            x /= 10;
 86            ++answer;
 87        }
 88        return answer;
 89    }
 90
 91    double sqrt(double S);
 92    inline double sqrt(double S) {
 93        // implements the Bakhshali method of square root computation to fix a PCC error
 94        double a, x = S / 2;
 95        uint32_t i;
 96        for (i = 0; i < PCC_SQRT_ACCURACY; i++) {
 97            a = (S - x*x) / (2 * x);
 98            x = x + a - (a * a) / (2 * (x + a));
 99        }
100        return x;
101    }
102
103    uintmax_t ceil(double x);
104    inline uintmax_t ceil(double x) {
105        uintmax_t ret = (uintmax_t) x;
106        if (x == (double) ret)
107            return ret;
108        return ret + 1;
109    }
110#endif
111
112#endif