C: the long and short of it

[PerranOak]

I am using MPLAB X IDE v5.25, XC8 v2.10 and a PIC16F1827.

I needed a large integer to calculate an area (an exercise) and so used:

unsigned long int area

I couldn't understand the result until I investigated and found that area was only 16-bit not, as I had expected, 32-bit.

Am I missing something?

Cheers

[SiliconWizard]

The XC8 manual defines "unsigned long" as a 32 bit unsigned integer, so unless there's a big bug with your compiler's version, your problem is not with the width of the 'area' variable.

It's very likely that it's a problem with the expression you're using to compute the area, related to how C handles integer promotion. Please post the line of code that computes it.

[andersm]

Remember that the data type that's on the left side of the assignment operator has absolutely no influence over how the right-hand side is calculated.

[brucehoult]

Remember that the data type that's on the left side of the assignment operator has absolutely no influence over how the right-hand side is calculated.

Yup. If HEIGHT and WIDTH are 16 bit variables then to get a 32 bit area you need to cast at least one of them to long. Otherwise all arithmetic (for char, short, int) is done at int precision.

Code: [Select]

area = (unsigned long)height * width;


[PerranOak]

I reduced the code to its essential elements:

       #include <xc.h>
       #include "stdio.h"

       void main(void)   
       {
         unsigned int radius;
         unsigned long int area;
         unsigned char outstring[20];


         while(1==1)
         {
           radius = 25;
           area = (radius * radius * 31416) / 10000;
           sprintf( outstring, "Area=%u", area);
           lcdputs( outstring );

         }
       }

Then, I changed "radius" to long int and, hey presto, it worked!

Thank you very much guys, I'd never have thought of that in a month of Sundays.

[SiliconWizard]

That's a very common trap related to how C deals with integers.
I recommend reading this: https://wiki.sei.cmu.edu/confluence/display/c/INT02-C.+Understand+integer+conversion+rules

And all this actually: https://wiki.sei.cmu.edu/confluence/display/c/SEI+CERT+C+Coding+Standard

Depending on what you do, you don't have to follow all of these rules, but they will certainly teach you a lot anyway.

[Jeroen3]

Yeah, C integer promotions are tricky. Please note that literals, like the 31416 are at least int unless otherwise specified.

Code: [Select]

integer-suffix:
  unsigned-suffix long-suffixopt
  unsigned-suffix long-long-suffix
  long-suffix unsigned-suffixopt
  long-long-suffix unsigned-suffixopt

unsigned-suffix: one of
  u U
long-suffix: one of
  l L
long-long-suffix: one of
  ll LLL
You may also be interested in the Q number format and libfixmath to stop scattering around those division by 10000.

[mrflibble]

Title fixed.

[Nominal Animal]

As an aside, 355/113 ≃ 3.14159292 gives you Pi with six/seven digits of precision, and it approximates Pi best among ratios of two sixteen bit (signed or unsigned) integers.  With 24-bit (signed or unsigned) integers, 5419351/1725033 ≃ 3.141592653589815 gives you 12/13 digits.  With eight (unsigned) bits, the best is 245/78 ≃ 3.141025641, just four digits of precision.  Below that, the best is 22/7 ≃ 3.142857, just three digits.

On an 8-bit or 16-bit microcontroller, Q15 or Q23 (QN) work well for angles.  As an integer, x = 2^N·degrees/180 or x = 2^N·radians/π with modulo arithmetic, so 2^N and -2^N describe the same angle.  Trigonometric functions can then be implemented via CORDIC, and you won't need π for anything (except when building the constant tables, but you do that at compile time on a host computer anyway).  You get the same or better overall precision as floats, but with much faster operations, and saving a byte per variable too. Things like versors (unit quaternions describing rotations, or orientations in 3D space; analogous to rotation matrices, but with much nicer features) can use four Q23 components, for a total of 12 bytes, nicely.  Then, the matrix you generate from the versor can use Q15.Q16 or even Q23.8, depending on the scale/range of your coordinates, without compounding rounding errors causing any grief.

[brucehoult]

As an aside, 355/113 ≃ 3.14159292 gives you Pi with six/seven digits of precision, and it approximates Pi best among ratios of two sixteen bit (signed or unsigned) integers.  With 24-bit (signed or unsigned) integers, 5419351/1725033 ≃ 3.141592653589815 gives you 12/13 digits.  With eight (unsigned) bits, the best is 245/78 ≃ 3.141025641, just four digits of precision.  Below that, the best is 22/7 ≃ 3.142857, just three digits.

And with 30 to 34 bit integers you won't do better than 817696623/260280919 = 3.141592653589793 gives three more digits -- basically full DP.

3.14159265358979311118 817696623/260280919
3.14159265358979323846 pi
3.14159265358981538324 5419351/1725033

[chriva]

This is way overkill if your compiler can emulate 32-bit operations and you don't need more but it's quite easy to extend the number of bits for multiplication, addition and subtraction if you know assembler.
as long as you can detect carry, you're basically only memory bound.

Division is ofc also possible but it becomes dirty, real fast. Not fun unless you really, really need it.

[westfw]

355/113 ≃ 3.14159292 gives you Pi with six/seven digits of precision

Almost on topic: Not in C (or C++), which will compute 355/113 as integer division and give you 3.  Better use "355.0/113", or M_PI that is defined off in math.h

[Siwastaja]

355/113 ≃ 3.14159292 gives you Pi with six/seven digits of precision

Almost on topic: Not in C (or C++), which will compute 355/113 as integer division and give you 3.  Better use "355.0/113", or M_PI that is defined off in math.h

Obviously, you would use 355/113 in a limited-resolution fixed point system (as explained by Nominal Animal); if you have floating point available and can use it, there is absolutely no point using 355.0/113, but use M_PI instead.

In fixed point math using 355/113, you need to use parenthesis to make sure the multiplication by 355 is done first with the input value, then result divided by 113; if the calculation order is wrong, 355/113 will indeed evaluate to 3. Further, you need to manually verify the intermediate multiplication fits the range of variables you are using, or provide necessary typecasting to make it fit.

[SiliconWizard]

Yes obviously, that's for fixed-point use.

The matter in general is related to how to handle integer multiply and divide properly. Not specific to this example, and important to know.

Example: you just need to compute: n * 2 / 3
In many languages, if you write it like this, it will yield zero, and using parentheses is required: (n * 2) / 3. This is very basic programming knowledge. And in this small example, you may not even quite realize you're dealing with fixed point (although it kind of is.)

[kjpye]

Example: you just need to compute: n * 2 / 3
In many languages, if you write it like this, it will yield zero, and using parentheses is required: (n * 2) / 3.

Really? Which languages are those? All of the (relatively few) languages I know use the usual rules of mathematics for such expressions, and the multiplications and divisions are performed left to right, so the multiplication is done before the division.

[ehughes]

One thing that I have found helpful is abandon the use of int, long, etc with *any* compiler.    Use stdint.h types.   uint8_t, in32_t, etc.


It makes tracking things down a bit easier.

[TK]

One thing that I have found helpful is abandon the use of int, long, etc with *any* compiler.    Use stdint.h types.   uint8_t, in32_t, etc.


It makes tracking things down a bit easier.

Especially when you are programming microcontrollers, it makes more sense to use uintN_t than the ambiguous int, long, short, etc.

[chriva]

More often than not, you don't have stdint, intdef and so on.
typedef that mofo and check default behaviour for that particular platforms and compiler. Some compiler flags will also change default sizes of things so be wary

[SiliconWizard]

Example: you just need to compute: n * 2 / 3
In many languages, if you write it like this, it will yield zero, and using parentheses is required: (n * 2) / 3.

Really? Which languages are those? All of the (relatively few) languages I know use the usual rules of mathematics for such expressions, and the multiplications and divisions are performed left to right, so the multiplication is done before the division.

I was thinking about C first, of course, but also in general languages in which * and / have the same precedence, and the order of evaluation is not guaranteed.

But you're making me question the relevance of my small example. I've just tested (with GCC -O3), and 'n * 2 / 3' yields the exact same code as '(n * 2) / 3'. Now this is of course no proof that this will always be the case with all compilers, all expressions, all contexts, etc. I've taken a look at the C standard (C99), and I haven't clearly seen stated that the order of operations for operators with the same precedence is guaranteed to be left to right. (I may have missed it, so if there's someone  - maybe like Nominal Animal - that knows the standard well and can point me to that, I'll take it!)

I've built a habit of explicitely using parentheses in order not to make any assumption about the order of evaluation (except of course when the order is certain from the precedence of operators).

Now a better example (more likely to be *really* problematic in C) would be when adding integer promotion rules in the mix.

[SiliconWizard]

One thing that I have found helpful is abandon the use of int, long, etc with *any* compiler.    Use stdint.h types.   uint8_t, in32_t, etc.

When you need an explicit fixed size, or at least a minimum guaranteed size that is not implementation-defined, yes, absolutely!

But keep in mind that using stdint types doesn't mitigate per se the 'integer promotion' issues you can encounter when dealing with integer arithmetic in C, so you still need to be careful with that.

[tggzzz]

Example: you just need to compute: n * 2 / 3
In many languages, if you write it like this, it will yield zero, and using parentheses is required: (n * 2) / 3.

Really? Which languages are those? All of the (relatively few) languages I know use the usual rules of mathematics for such expressions, and the multiplications and divisions are performed left to right, so the multiplication is done before the division.

The operator precedence is completely separate from the order of evaluation, especially with C/C++ compilation with optimisations turned on.

There is no concept of left-to-right or right-to-left evaluation in C, which is not to be confused with left-to-right and right-to-left associativity of operators: the expression f1() + f2() + f3() is parsed as (f1() + f2()) + f3() due to left-to-right associativity of operator+, but the function call to f3 may be evaluated first, last, or between f1() or f2() at run time.
https://en.cppreference.com/w/c/language/eval_order

In my experience most programmers that think they know C, don't. I know that I don't know C anymore

For amysement, you could also look at APL; at https://tryapl.org/ try typing in each of these expressions, and explain the results
      7 - 7
      7 - 7 - 7
      7 - 7 - 7 - 7
      7 - 7 - 7 - 7 - 7

[Nominal Animal]

To do fixed-point multiplication in C, you almost always need to add proper casting.  The one special case is when you know (from prior checking or casting) that the product cannot overflow.  The same applies to fixed-point division, because it starts with a multiplication (by the integer value corresponding to 1.0 in fixed-point).

As an example, here is a function that multiplies a 16-bit signed integer (assumed to be not greater than 10430 in magnitude) by Pi.

Code: [Select]

int16_t mul_pi16(const int16_t v) {
    return ((int32_t)v * 355) / 113;  /* Note: 355/113 is the best 16:16-bit approximation for Pi */
}

By the way, there are exactly 105 32-bit signed integer ratios (210 if you include negative numerator and denominators) matching Pi in IEEE 754-Binary64 ("double").  My test program is still comparing them against IEEE-754 Binary64 Pi (on operations over all 32-bit signed integers), and it is already 3.4% done; I'll report if it yields any practical results.


In C99 and later, casting to a numeric type (integer or floating-point type) has much stricter definition than in ANSI C/C89: the cast value is limited to the range and precision of the type it is cast to.  This means that a cast in an expression is equivalent to using a temporary variable of that type, no matter what optimizations are in use (except those that break the standard, of course).

I am not a language lawyer, and much prefer real world experience to the standard text, but since none of the standard arithmetic operators (+, -, *, /, bit shifts) are sequence points, the compiler is free to evaluate such in any order it chooses, as long as the mathematical rules of precedence are followed.  By that I mean that given a*b - c*d, the compiler is free to choose which product it computes first.  (In C99 n1256, this is stated in 6.5p3: "[Except for logical operators, conditional operator, and function call parameters], the order of evaluation of subexpressions and the order in which side effects take place are both unspecified.")

Parentheses in an expression (outside function call parameter list, which lexically is parenthesised) are not a sequence point either, so they are only used in their mathematical sense, and do not define the order in which the compiler must evaluate the subexpressions.  In particular, an expression such as (a*(b-c))*d does not mean that the compiler has to calculate (b-c) first, then multiply it by a, and finally by d; the compiler is free to choose the evaluation order, as long as the math rules are not violated.

[Nominal Animal]

The report came in.  Among 32-bit signed integer ratios, there are two particularly useful: 908679182/289241567 and 2025273829/644664682.  The former is better when rounding down or truncating, the latter is better when rounding towards nearest integer.  In hexadecimal, these are 0x3629580e/0x113d79df and 0x78b739e5/0x266ccd6a, respectively.

The following assume int32_t x;

Instead of (int32_t)((double)x * M_PI) , you can use
    ((int64_t)x * 908679182) / 289241567
It gives the same answer for almost all x – differs by 1 only in 109 cases among all nonnegative 32-bit signed integer x.

Instead of (int32_t)((double)x / M_PI) , you can use
    ((int64_t)x * 644664682) / 2025273829
It gives the same answer for almost all x – differs by 1 only in 13 cases among all nonnegative 32-bit signed integer x.

Instead of (int32_t)round((double)x * M_PI) , you can use
    ((int64_t)x * 908679182 + 144620783) / 289241567
It gives the same answer for almost all x – differs by 1 only in 104 cases among all nonnegative 32-bit signed integer x.
(For negative x, you must use -144620783 instead.)

Instead of (int32_t)round((double)x / M_PI) , you can use
    ((int64_t)x * 644664682 + 1012636914) / 2025273829
It gives the same answer for almost all x – differs by 1 only in 13 cases among all nonnegative 32-bit signed integer x.
(For negative x, you must use -1012636914 instead.)

I haven't bothered to check the 13 cases when rounding to see which one is actually better, the floating-point one or the integer approximation.  In all tests with rounding, I explicitly used half the divisor; should have checked with +1 since the divisors are odd.  Also, I forgot to check for ceil() cases.

[Nominal Animal]

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:

Code: [Select]

#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):

Code: [Select]

#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.

[westfw]

Use stdint.h types.   uint8_t, in32_t, etc.

Maybe.  You should probably be using "size_t", "uint_leastN_t", or "int_fastN_t", too...
The percentage of CPUs where "int" is an inefficient type for something like:

Code: [Select]

   for (int i=0; i < 100; i++) { foo(myarray[i]); /* stuff */ }is becoming pretty small.