Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- /* File: fapprox.c */
- /* Purpose: Fraction approximator */
- /* Standard headers */
- #include <errno.h>
- #include <float.h>
- #include <stdarg.h>
- #include <stdbool.h>
- #include <stdint.h>
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include <tgmath.h>
- /*** Attributes ***/
- #ifdef __GNUC__
- /* Nonvalent attribute syntax */
- # define ATTR(attr) \
- __attribute__((attr))
- /* Multivalent attribute syntax */
- # define ATTRS(attr, ...) \
- __attribute__((attr(__VA_ARGS__)))
- /* Variable is not used */
- # define unused \
- ATTR(unused)
- /* Function never returns */
- # define EXIT \
- ATTR(noreturn)
- /* Function with printf() semantics */
- # define PRINTF(fmt, args) \
- ATTRS(format, printf, fmt, args)
- /* Function has only nontemporal global state */
- # define PURE \
- ATTR(pure)
- /* Function has only nontemporal local state */
- # define CONST \
- ATTR(const)
- #else /* __GNUC__ */
- /* Remove attributes when not compiling with gcc */
- # define ATTR(...)
- # define ATTRS(...)
- #endif /* __GNUC__ */
- /*** Constants ***/
- /* Two pi (360°) */
- #define TWO_PI 6.283185307179586476925286766559005768394338798750211641949889185
- /* pi (180°) */
- #define PI 3.141592653589793238462643383279502884197169399375105820974944592
- /* One half pi (90°) */
- #define HALF_PI 1.570796326794896619231321691639751442098584699687552910487472296
- /* One third pi (60°) */
- #define THIRD_PI 1.047197551196597746154214461093167628065723133125035273658314864
- /* One quarter pi (45°) */
- #define QTR_PI 0.785398163397448309615660845819875721049292349843776455243736148
- /* One sixth pi (30°) */
- #define SIXTH_PI 0.523598775598298873077107230546583814032861566562517636829157432
- /* pi / 180 */
- #define RAD_PER_DEG 0.017453292519943295769236907684886127134428718885417254560971914
- /* 180 / pi */
- #define DEG_PER_RAD 57.29577951308232087679815481410517033240547246656432154916024386
- /* e */
- #define E 2.718281828459045235360287471352662497757247093699959574966967628
- /* No error */
- #define ENOERR 0
- /* Minimum difference between two long doubles */
- #define epsi LDBL_EPSILON
- /*** Macros ***/
- /*
- * Return true if X is greater than Y
- */
- #define GT(x, y) \
- ((x) - (y) >= epsi)
- /*
- * Return true if X is less than Y
- */
- #define LT(x, y) \
- GT(y, x)
- /*
- * Return true if X is greater or equal to Y
- */
- #define GE(x, y) \
- (!LT(x, y))
- /*
- * Return true if X is less or equal to Y
- */
- #define LE(x, y) \
- (!GT(x, y))
- /*
- * Return true if X is equal to Y
- */
- #define EQ(x, y) \
- (!GT(x, y) && !LT(x, y))
- /*
- * Return true if X is not equal to Y
- */
- #define NE(x, y) \
- (!EQ(x, y))
- /*
- * Return true if X is equal to zero
- */
- #define Z(x) \
- EQ(x, 0)
- /*
- * Return true if X is not equal to zero
- */
- #define NZ(x) \
- (!Z(x))
- /*
- * Return true if X is negative
- */
- #define NEG(x) \
- LT(x, 0)
- /*
- * Return true if X is positive
- */
- #define POS(x) \
- GT(x, 0)
- /*
- * Return the sign of X
- */
- #define SGN(x) \
- (NEG(x) ? (-1) : (POS(x) ? (+1) : (0)))
- /*
- * Return the sign of X as a string
- * - "-" if X is negative
- * - "" if X is equal to zero
- * - "+" if X is positive
- */
- #define SGNSTR(x) \
- (NEG(x) ? "-" : (POS(x) ? "+" : ""))
- /*
- * Return true if theta is in the first quadrant
- */
- #define NE_QUAD(theta) \
- (GT(theta, 0) && LT(theta, HALF_PI))
- /*
- * Return the maximum of X and Y
- */
- #define MAX(x, y) \
- (GT(x, y) ? (x) : (y))
- /*
- * Return the minimum of X and Y
- */
- #define MIN(x, y) \
- (LT(x, y) ? (x) : (y))
- /*
- * Return the absolute value of X
- */
- #define ABS(x) \
- (NEG(x) ? (-x) : (+x))
- /*
- * Return the maximum of X and Y by absolute value
- */
- #define MAXABS(x, y) \
- (GT(ABS(x), ABS(y)) ? (x) : (y))
- /*
- * Return the minimum of X and Y by absolute value
- */
- #define MINABS(x, y) \
- (LT(ABS(x), ABS(y)) ? (x) : (y))
- /*
- * Return the maximum magnitude of X and Y
- */
- #define MAXMAG(x, y) \
- MAX(ABS(x), ABS(y))
- /*
- * Return the minimum magnitude of X and Y
- */
- #define MINMAG(x, y) \
- MIN(ABS(x), ABS(y))
- /*** Types ***/
- typedef intmax_t smax; /* Signed integer of maximum size */
- typedef uintmax_t umax; /* Unsigned integer of maximum size */
- typedef long double real; /* Floating point */
- /*** Variables ***/
- /* Program name */
- static const char *prog;
- /*** Constant functions ***/
- /* Declare constant */
- CONST static umax gcd(const umax x, const umax y);
- CONST static real reldiff_signed(const real x, const real y);
- /*
- * Return the greatest common divisor of X and Y
- */
- static umax gcd(const umax x, const umax y)
- {
- /* Evenly divisible */
- if (!y) return x;
- /* Recur */
- return gcd(y, x % y);
- }
- /*
- * Return the signed relative difference from Y to X
- */
- static real reldiff_signed(const real x, const real y)
- {
- return (x - y) / MAXMAG(x, y);
- }
- /*** Formatting functions ***/
- /* Declare printf() semantics */
- PRINTF(1, 2) static void warn(const char *const fmt, ...);
- /*
- * Display a formatted message on standard error
- *
- * Arguments:
- * - format string or NULL
- * - format arguments
- *
- * Output:
- * - program name
- * - formatted arguments
- * - library error
- */
- static void warn(const char *const fmt, ...)
- {
- /* Remember the error */
- const int err = errno;
- /* Display the program name */
- fprintf(stderr, "%s: ", prog);
- /* Format arguments if necessary */
- if (fmt)
- {
- va_list arg;
- /* Prepare arguments */
- va_start(arg, fmt);
- /* Warning */
- vfprintf(stderr, fmt, arg);
- /* Clean up */
- va_end(arg);
- }
- /* Note errors */
- if (err) fprintf(stderr, ": %s", strerror(err));
- /* Terminate */
- fprintf(stderr, "\n");
- /* Restore the error */
- errno = err;
- }
- /*** Terminal functions ***/
- /* Declare terminal */
- EXIT static void usage(void);
- /*
- * Display the program usage and exit
- */
- static void usage(void)
- {
- /* Forget errors */
- errno = ENOERR;
- /* Warning */
- warn("usage: <denominator bound> <fraction>...");
- /* Error */
- exit(1);
- }
- /*** Normal functions ***/
- /*
- * Empty buffers to standard output
- */
- static void flush()
- {
- /* Just call fflush() */
- fflush(stdout);
- }
- /*
- * Program entry
- *
- * Arguments:
- * - denominator bound (positive integer)
- * - list of real numbers
- */
- int main(const int argc, const char *const *argv)
- {
- /* Save the program name */
- prog = argv[0];
- /* Require valid commandline */
- if (argc < 3) usage();
- char *end;
- /* Extract the denominator bound */
- const smax bound = strtoll(argv[1], &end, 0);
- /* Notice invalid bound */
- if (*end) usage();
- /* Require positive values */
- if (bound < 1) usage();
- /* Obtain the maximum fully representable value */
- const real max = pow(FLT_RADIX, LDBL_MANT_DIG) - 1;
- /* Notice excessive bound */
- if (bound > max)
- {
- /* Warning */
- warn(
- "warning: %jd exceeds floating point precision (%u digits in base %u: %.*Lg)",
- bound,
- LDBL_MANT_DIG, FLT_RADIX, LDBL_DIG, max
- );
- }
- /* Determine maximum integer length */
- const size_t int_len = snprintf(NULL, 0, "%jd", bound);
- /* Determine maximum possible floating point length */
- const size_t flt_len = snprintf(NULL, 0, "%+.*Le", LDBL_DIG, LDBL_MAX);
- /* Examine arguments */
- for (argv += 2; *argv; argv++)
- {
- const char *const arg = *argv;
- /* Extract a real number */
- real frac = strtold(arg, &end);
- /* Notice invalid number */
- if (*end) usage();
- /* Classify this fraction */
- switch (fpclassify(frac))
- {
- /* Invalid */
- default:
- {
- /* Display usage */
- usage();
- }
- /* Simple case */
- case FP_ZERO:
- {
- /* Display */
- printf("0: 0\n");
- /* Next... */
- continue;
- }
- /* Not representable */
- case FP_SUBNORMAL:
- {
- /* Warning */
- warn("%s: not representable (magnitude is less than %Le)", arg, LDBL_MIN);
- /* Next... */
- continue;
- }
- /* Normal fractional value */
- case FP_NORMAL:
- {
- break;
- }
- }
- /* Remember the sign */
- const int sgn = SGN(frac);
- /* Begin */
- printf("%s: ", arg);
- /* Display the sign if necessary */
- if (sgn < 0) printf("-");
- /* Forget the sign */
- frac = ABS(frac);
- real ip;
- /* Extract the integer part */
- const real fp = modfl(frac, &ip);
- /* Display the integer part */
- printf("%.*Lg", LDBL_DIG, ip);
- /* Notice completion */
- if (Z(fp))
- {
- /* Terminate */
- printf("\n");
- /* Next... */
- continue;
- }
- /* Display the sign */
- printf(" %s\n", SGNSTR(sgn));
- /* Locate endpoint */
- const real x = bound;
- const real y = fp * bound;
- /* Extract radius */
- unused const real r = hypot(x, y);
- /* Extract central angle */
- const real theta = atan(fp);
- /* Obtain normal */
- const real norm = theta + HALF_PI;
- /* Extract angular cross section */
- const real csect = 2 * atan2(1, bound);
- /* Obtain sector bounds */
- const real min = NE_QUAD(theta - csect) ? (theta - csect) : nextafter(0, HALF_PI);
- const real max = NE_QUAD(theta + csect) ? (theta + csect) : nextafter(HALF_PI, 0);
- /* Obtain boundary slopes */
- const real slope_min = tan(min);
- unused const real slope_max = tan(max);
- /* Assume no exact solutions */
- bool no_exact = true;
- /* Traverse the sector */
- for (smax x1 = 1; x1 <= bound; x1++)
- {
- /* Locate the sector beginning */
- const real y1 = x1 * slope_min;
- /* Locate lattice points */
- for (real tx = x1, ty = y1; true; tx += cos(norm), ty += sin(norm))
- {
- /* Obtain the nearest point */
- const smax num = llrint(ty);
- const smax den = llrint(tx);
- /* Display valid points */
- if ((num > 0) && (den > 0))
- {
- /* Obtain the approximate fraction */
- const real approx = (real)(num) / (real)(den);
- /* Extract the relative difference */
- const real rd = reldiff_signed(approx, fp);
- /* Notice exact results */
- const bool exact = Z(rd);
- /* Examine exact results */
- if (exact)
- {
- /* Ignore duplicates */
- if (!no_exact) continue;
- /* Remember exactness */
- no_exact = false;
- }
- /* Extract the greatest common divisor */
- const umax g = gcd(num, den);
- char num_fmt[256];
- char den_fmt[256];
- /* Format the approximation */
- sprintf(num_fmt, "%s%jd", (exact ? "<" : " "), num / g);
- sprintf(den_fmt, "%jd%s", den / g, (exact ? ">" : " "));
- /* Display the approximation */
- printf("\t%*s / %*s", +(int_len + 1), num_fmt, -(int_len + 1), den_fmt);
- /* Display the corresponding fraction */
- printf("\t%-*.*Lg", flt_len, LDBL_DIG, approx);
- /* Indicate exactness */
- if (exact)
- {
- printf("\t==");
- }
- /* Compare */
- else
- {
- /* Display the arithmetic difference */
- printf("\t%+-*.*Lg", flt_len, LDBL_DIG, approx - fp);
- /* Display the logarithmic difference */
- printf("\tlog %+-*.*Le", flt_len, LDBL_DIG, log10(approx) - log10(fp));
- /* Display the relative difference */
- printf("\trel %+-*.*Le", flt_len, LDBL_DIG, rd);
- }
- /* Terminate */
- printf("\n");
- /* Flush output */
- flush();
- }
- /* Notice completion */
- if (atan2(ty, tx) > max) break;
- }
- }
- }
- }
- /* End: fapprox.c */
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement