I wrote a couple of programs that can be used to find and evaluate integer ratios.
The first one is a POSIXy one that uses Python3 to evaluate an expression (so that you can run it with e.g. 'cos(22.5*pi/180)' as a parameter, then finds the best 32-bit signed integer ratio to approximate it:
#define POSIX_C_SOURCE 200809L
#include <stdlib.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/wait.h>
#include <fcntl.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <errno.h>
#define MAX_RATIOS 1000
int safefd(int fd)
{
unsigned int closemask = 0U;
int saved_errno = errno;
while (fd == STDIN_FILENO || fd == STDOUT_FILENO || fd == STDERR_FILENO) {
closemask |= 1 << fd;
fd = dup(fd);
if (fd == -1) {
saved_errno = errno;
break;
}
}
if (closemask & (1 << STDIN_FILENO)) close(STDIN_FILENO);
if (closemask & (1 << STDOUT_FILENO)) close(STDOUT_FILENO);
if (closemask & (1 << STDERR_FILENO)) close(STDERR_FILENO);
errno = saved_errno;
return fd;
}
int devnull(const int fd, const int flags)
{
int newfd;
/* Sanity check. */
if (fd == -1)
return errno = EBADF;
if (flags & O_CREAT)
return errno = EINVAL;
newfd = open("/dev/null", flags);
if (newfd == -1)
return errno;
if (newfd != fd) {
if (dup2(fd, newfd) == -1) {
const int saved_errno = errno;
close(newfd);
return saved_errno;
}
if (close(newfd) == -1) {
const int saved_errno = errno;
close(fd);
return saved_errno;
}
}
return 0;
}
static inline int writeall(const int fd, const void *const buf, const size_t len)
{
const char *ptr = buf;
const char *const end = (const char *)buf + len;
ssize_t n;
if (fd == -1)
return errno = EBADF;
while (ptr < end) {
n = write(fd, ptr, (size_t)(end - ptr));
if (n > 0) {
ptr += n;
} else
if (n != -1) {
return errno = EIO;
} else
if (errno != EINTR) {
return errno;
}
}
return 0;
}
int execute(const char *name, char *const argv[], const char *override_locale, pid_t *pid_to, int *output_fd)
{
int status_pipe[2] = { -1, -1 };
int output_pipe[2] = { -1, -1 };
int error;
pid_t child, p;
/* Sanity checks. */
if (!name || !*name || !argv || !argv[0])
return errno = EINVAL;
/* Clear pid_to and output_fd, if specified. */
if (pid_to)
*pid_to = (pid_t)-1;
if (output_fd)
*output_fd = -1;
/* Construct output pipe. */
if (pipe(output_pipe) == -1) {
error = errno;
return errno = (error ? error : EIO);
}
/* Construct status pipe. */
if (pipe(status_pipe) == -1) {
error = errno;
close(output_pipe[0]);
close(output_pipe[1]);
return errno = (error ? error : EIO);
}
child = fork();
if (!child) {
/* Child process. */
/* Close parent (read) ends of the pipes. */
close(status_pipe[0]);
close(output_pipe[0]);
/* Make sure status pipe write end does not overlap a standard descriptor. */
status_pipe[1] = safefd(status_pipe[1]);
if (status_pipe[1] == -1) {
close(output_pipe[1]);
exit(99);
}
if (output_fd) {
/* Move write end of output pipe to standard output. */
if (dup2(output_pipe[1], STDOUT_FILENO) == -1) {
error = errno;
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
if (close(output_pipe[1]) == -1) {
error = errno;
close(STDOUT_FILENO);
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
} else {
/* Discard all output by redirecting it to /dev/null. */
close(output_pipe[1]);
if (devnull(STDOUT_FILENO, O_WRONLY)) {
error = errno;
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
}
/* Mark status pipe write end close-on-exec. */
if (fcntl(status_pipe[1], F_SETFD, FD_CLOEXEC) == -1) {
error = errno;
close(STDOUT_FILENO);
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
/* Redirect standard input and error to /dev/null. */
if (devnull(STDIN_FILENO, O_RDONLY) ||
devnull(STDERR_FILENO, O_WRONLY)) {
error = errno;
close(STDOUT_FILENO);
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
if (override_locale) {
setenv("LANG", override_locale, 1);
setenv("LC_ALL", override_locale, 1);
}
/* Execute. */
execvp(name, argv);
error = errno;
close(STDOUT_FILENO);
writeall(status_pipe[1], &error, sizeof error);
close(status_pipe[1]);
exit(99);
}
/* Parent process. */
/* Close child ends of pipes. */
close(status_pipe[1]);
close(output_pipe[1]);
/* Read status pipe. */
{
char buffer[2 * sizeof (int)];
size_t have = 0;
ssize_t n;
error = 0;
while (have < sizeof buffer) {
n = read(status_pipe[0], buffer + have, sizeof buffer - have);
if (n > 0) {
have += n;
} else
if (n == 0) {
break;
} else
if (n != -1) {
error = EIO;
break;
} else
if (errno != EINTR) {
error = errno;
break;
}
}
if (!error && have == sizeof error)
memcpy(&error, buffer, sizeof error);
close(status_pipe[0]);
if (have > 0) {
/* Failed. Make sure we have an error code. */
if (!error)
error = ENOTSUP;
/* Reap child process. */
do {
p = waitpid(child, NULL, 0);
} while (p == -1 && errno == EINTR);
return errno = error;
}
}
/* Exec successful. */
if (output_fd) {
*output_fd = output_pipe[0];
} else {
close(output_pipe[0]);
}
if (pid_to) {
*pid_to = child;
}
return errno = 0;
}
const char expr_prefix[] =
"from math import *\n"
"print(\"%.32g\\n\" % (";
#define EXPR_PREFIX_LEN (sizeof expr_prefix - 1)
const char expr_suffix[] =
"))\n";
#define EXPR_SUFFIX_LEN (sizeof expr_suffix - 1)
int evaluate(const char *expr, double *to)
{
const size_t expr_len = (expr) ? strlen(expr) : 0;
char *expression, *ptr;
char *argv[4];
int result, error, fd = -1;
pid_t child = -1;
if (expr_len < 1)
return errno = EINVAL;
expression = malloc(EXPR_PREFIX_LEN + expr_len + EXPR_SUFFIX_LEN + 1);
if (!expression)
return errno = ENOMEM;
ptr = expression;
memcpy(ptr, expr_prefix, EXPR_PREFIX_LEN);
ptr += EXPR_PREFIX_LEN;
memcpy(ptr, expr, expr_len);
ptr += expr_len;
memcpy(ptr, expr_suffix, EXPR_SUFFIX_LEN);
ptr += EXPR_SUFFIX_LEN;
*ptr = '\0';
argv[0] = "python3";
argv[1] = "-c";
argv[2] = expression;
argv[3] = NULL;
result = execute("python3", argv, "C", &child, &fd);
if (result)
return errno = result;
free(expression);
error = 0;
{
char buffer[128];
size_t have = 0;
ssize_t n;
while (have < sizeof buffer) {
n = read(fd, buffer + have, sizeof buffer - have);
if (n > 0)
have += n;
else
if (n == 0)
break;
else
if (n != -1) {
error = EIO;
break;
} else
if (errno != EINTR) {
error = errno;
break;
}
}
if (!error && have > 0 && have < sizeof buffer) {
double value;
char *end = buffer;
buffer[have] = '\0';
errno = 0;
value = strtod(buffer, &end);
if (errno) {
error = errno;
} else
if (end == buffer) {
error = EINVAL;
} else
if (*end != '\t' && *end != '\n' && *end != '\v' &&
*end != '\f' && *end != '\r' && *end != ' ') {
error = EINVAL;
}
if (!error && !isfinite(value))
error = EDOM;
if (!error && to)
*to = value;
}
}
/* Reap child process. */
{
int status = 0;
pid_t p;
do {
p = waitpid(child, &status, 0);
} while (p == -1 && errno == EINTR);
if (p == -1) {
if (!error)
error = errno;
}
if (!error) {
if (!WIFEXITED(status) || WEXITSTATUS(status))
error = EINVAL;
}
}
/* Done. */
return errno = error;
}
struct ratio {
int n;
int d;
};
static double ratio_delta = HUGE_VAL;
static size_t ratios_max = 0;
static size_t ratios = 0;
struct ratio *ratio = NULL;
static inline void clear_ratios(void)
{
ratios = 0;
}
static inline void add_ratio(const int n, const int d)
{
if (ratios >= ratios_max) {
ratios_max = (ratios | 32767) + 32769;
ratio = realloc(ratio, ratios_max * sizeof ratio[0]);
if (!ratio) {
fprintf(stderr, "\rOut of memory.\n");
exit(EXIT_FAILURE);
}
}
ratio[ratios].n = n;
ratio[ratios].d = d;
ratios++;
}
static double target = 0.0;
static inline void check(const int n, const int d)
{
if (d != 0) {
const double delta = fabs(target - (double)((double)n / (double)d));
if (delta < ratio_delta) {
clear_ratios();
ratio_delta = delta;
add_ratio(n, d);
} else
if (delta == ratio_delta) {
size_t i;
/* Check if a duplicate ratio (common factor) */
for (i = 0; i < ratios; i++)
if (n > ratio[i].n) {
const unsigned int factor = n / ratio[i].n;
if (ratio[i].n * factor == n && ratio[i].d * factor == d)
return;
} else
if (n < ratio[i].n) {
const unsigned int factor = ratio[i].n / n;
if (n * factor == ratio[i].n && d * factor == ratio[i].d)
return;
} else {
if (d == ratio[i].d)
return;
}
/* No, new one; add. */
add_ratio(n, d);
}
}
}
static int parse_uint(const char *s, unsigned int *to)
{
const char *z = s;
unsigned long ul;
unsigned int u;
if (!s || !*s)
return -1;
errno = 0;
ul = strtoul(s, (char **)&z, 0);
if (errno)
return -1;
if (z == s)
return -1;
while (*z == '\t' || *z == '\n' || *z == '\v' ||
*z == '\f' || *z == '\r' || *z == ' ')
z++;
if (*z)
return -1;
u = (unsigned int)ul;
if ((unsigned long)u != ul)
return -1;
if (to)
*to = u;
return 0;
}
int main(int argc, char *argv[])
{
unsigned int nmax = 2147483647;
unsigned int n;
size_t i;
if (argc < 1) {
fprintf(stderr, "No.\n");
exit(EXIT_FAILURE);
}
if (argc < 2 || argc > 3 || (argc > 1 && (!strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")))) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
fprintf(stderr, " %s MATH-EXPRESSION [ MAX-NUMERATOR ]\n", argv[0]);
fprintf(stderr, "\n");
fprintf(stderr, "This program uses Python3 to evaluate MATH-EXPRESSION\n");
fprintf(stderr, "(with all of math imported to local namespace),\n");
fprintf(stderr, "and then finds the best 32-bit ratio to describe the value.\n");
fprintf(stderr, "\n");
exit(EXIT_FAILURE);
}
if (argc > 2) {
if (parse_uint(argv[2], &nmax) || nmax < 1) {
fprintf(stderr, "%s: Invalid maximum numerator.\n", argv[2]);
return EXIT_FAILURE;
}
}
if (evaluate(argv[1], &target)) {
fprintf(stderr, "%s: %s.\n", argv[1], strerror(errno));
exit(EXIT_FAILURE);
}
fprintf(stderr, "Finding best approximation for value %.24g:\n", target);
fflush(stderr);
for (n = 0; n <= nmax; n++) {
int d = (int)((double)n / target);
check(n, d - 1);
check(n, d);
check(n, d + 1);
if (!(n & 16777215)) {
fprintf(stderr, "\r%5.1f%% done", 100.0 * (double)n / (double)nmax);
fflush(stderr);
}
}
if (ratios < 1) {
fprintf(stderr, "\rNo approximations found.\n");
return EXIT_FAILURE;
} else {
size_t maxi;
if (ratio_delta > 0.0)
fprintf(stderr, "\rDone. Found %zu approximations:\n", ratios);
else
fprintf(stderr, "\rDone. Found %zu approximations to within double precision:\n", ratios);
fflush(stderr);
maxi = ratios;
if (maxi > MAX_RATIOS)
maxi = MAX_RATIOS;
for (i = 0; i < maxi; i++)
printf("%zu. %10d / %10d (error %+.24g)\n", i + 1, ratio[i].n, ratio[i].d, (double)ratio[i].n / (double)ratio[i].d - target);
}
return EXIT_SUCCESS;
}
The second is a standard C one that examines a specified ratio among all 32-bit integer multiplicands (where the result does not overflow):
#include <stdlib.h>
#include <inttypes.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include <math.h>
static inline const char *skip_spaces(const char *s)
{
if (s)
while (*s == '\t' || *s == '\n' || *s == '\v' ||
*s == '\f' || *s == '\r' || *s == ' ')
s++;
return s;
}
static int parse_double(const char *s, double *to)
{
const char *end = NULL;
double val;
if (!s || !*s)
return -1;
errno = 0;
val = strtod(s, (char **)&end);
if (errno)
return -1;
if (!isfinite(val))
return -1;
if (!end || end == s)
return -1;
end = skip_spaces(end);
if (*end)
return -1;
if (to)
*to = val;
return 0;
}
static int parse_uint(const char *s, unsigned int *to)
{
const char *end = NULL;
unsigned long ul;
unsigned int u;
if (!s || !*s)
return -1;
errno = 0;
ul = strtoul(s, (char **)&end, 0);
if (errno)
return -1;
if (!end || end == s)
return -1;
end = skip_spaces(end);
if (*end)
return -1;
u = (unsigned int)ul;
if ((unsigned long)u != ul)
return -1;
if (to)
*to = u;
return 0;
}
int main(int argc, char *argv[])
{
unsigned int nall = 0;
unsigned int ntrunc = 0, overtrunc = 0, undertrunc = 0;
unsigned int nround = 0, overround = 0, underround = 0;
unsigned int nfloor = 0, overfloor = 0, underfloor = 0;
unsigned int nceil = 0, overceil = 0, underceil = 0;
unsigned int imax = 2147483647;
unsigned int o = 0;
unsigned int i, n, d, lmax;
const char *valstr;
double value;
if (argc < 1) {
fprintf(stderr, "No.\n");
return EXIT_FAILURE;
}
if (argc < 4 || argc > 6 || !strcmp(argv[1], "-h") || !strcmp(argv[1], "--help")) {
fprintf(stderr, "\n");
fprintf(stderr, "Usage: %s [ -h | --help ]\n", argv[0]);
fprintf(stderr, " %s TARGET N D [ OFFSET [ MAX ] ]\n", argv[0]);
fprintf(stderr, "\n");
fprintf(stderr, "This program calculates\n");
fprintf(stderr, " (int)( ((int64_t)x * N + OFFSET) / d)\n");
fprintf(stderr, "for all nonnegative 32-bit signed integers, or up to MAX,\n");
fprintf(stderr, "and reports its differences to\n");
fprintf(stderr, " (int)(TARGET * x)\n");
fprintf(stderr, " (int)floor(TARGET * x)\n");
fprintf(stderr, " (int)ceil(TARGET * x)\n");
fprintf(stderr, " (int)round(TARGET * x)\n");
fprintf(stderr, "\n");
return EXIT_FAILURE;
}
if (parse_double(argv[1], &value) || value <= 0.0) {
fprintf(stderr, "%s: Invalid target floating-point value. Must be positive.\n", argv[1]);
return EXIT_FAILURE;
}
valstr = argv[1];
if (parse_uint(argv[2], &n) || n < 1) {
fprintf(stderr, "%s: Invalid numerator (N). Must be a nonnegative integer.\n", argv[2]);
return EXIT_FAILURE;
}
if (parse_uint(argv[3], &d) || d < 1) {
fprintf(stderr, "%s: Invalid denominator (D). Must be a nonnegative integer.\n", argv[3]);
return EXIT_FAILURE;
}
if (argc > 4) {
if (parse_uint(argv[4], &o)) {
fprintf(stderr, "%s: Invalid offset.\n", argv[4]);
return EXIT_FAILURE;
}
}
if (argc > 5) {
if (parse_uint(argv[5], &imax) || imax < 1) {
fprintf(stderr, "%s: Invalid maximum integer to examine.\n", argv[5]);
return EXIT_FAILURE;
}
}
if (value > 1.0) {
lmax = (unsigned int)ceil(2147483648.0 / value) - 1;
if (imax > lmax)
imax = lmax;
}
for (i = 0; i <= imax; i++) {
const int target_trunc = (int)(value * (double)i);
const int target_floor = (int)floor(value * (double)i);
const int target_ceil = (int)ceil(value * (double)i);
const int target_round = (int)round(value * (double)i);
const int approx = (int)(((int64_t)i * n + o) / d);
nall++;
ntrunc += (approx == target_trunc);
overtrunc += (approx > target_trunc);
undertrunc += (approx < target_trunc);
nround += (approx == target_round);
overround += (approx > target_round);
underround += (approx < target_round);
nfloor += (approx == target_floor);
overfloor += (approx > target_floor);
underfloor += (approx < target_floor);
nceil += (approx == target_ceil);
overceil += (approx > target_ceil);
underceil += (approx < target_ceil);
if (!(i & 16777215)) {
fprintf(stderr, "\r%5.1f%% done", 100.0 * (double)i / (double)imax);
fflush(stderr);
}
}
fprintf(stderr, "\rDone. \n");
fflush(stderr);
printf("(int)(((int64_t)x * %u + %u) / %u) == (int)(%s * x):\n", n, o, d, valstr);
printf(" %u correct (%.3f%%), %u differences (%u under, %u over)\n", ntrunc, 100.0*(double)ntrunc/(double)nall, nall-ntrunc, undertrunc, overtrunc);
printf("\n");
printf("(int)(((int64_t)x * %u + %u) / %u) == (int)floor(%s * x):\n", n, o, d, valstr);
printf(" %u correct (%.3f%%), %u differences (%u under, %u over)\n", nfloor, 100.0*(double)nfloor/(double)nall, nall-nfloor, underfloor, overfloor);
printf("\n");
printf("(int)(((int64_t)x * %u + %u) / %u) == (int)ceil(%s * x):\n", n, o, d, valstr);
printf(" %u correct (%.3f%%), %u differences (%u under, %u over)\n", nceil, 100.0*(double)nceil/(double)nall, nall-nceil, underceil, overceil);
printf("\n");
printf("(int)(((int64_t)x * %u + %u) / %u) == (int)round(%s * x):\n", n, o, d, valstr);
printf(" %u correct (%.3f%%), %u differences (%u under, %u over)\n", nround, 100.0*(double)nround/(double)nall, nall-nround, underround, overround);
printf("\n");
return EXIT_SUCCESS;
}
With these, in just a couple of minutes, I can say that to approximate cos(22.5*pi/180) ≃ 0.923879532511286738483136, there are 156 possible ratios among 32-bit integers. Of these, the first one (smallest numerator),
(int)(((int64_t)x * 126426688 + 68421637) / 136843261) == (int)round(0.923879532511286738483136 * x)
is almost always true for positive 32-bit signed integer x. It is not true for only 58 different x; of which 27 times the left side is one less, and 31 times one greater than the right side.