ccalc.c for testing math functions

[Rich Felker <dalias@aerifal.cx>]

[-- Attachment #1: Type: text/plain, Size: 530 bytes --]

Here is a quick test program I wrote for playing around with musl's
(or any libc/libm's) math functions. It's an RPN calculator that
accepts input in argv[1] or stdin. It's by no means mature, but it's
very simple and usable for the purpose for which it was written
(testing), and the concepts in it could easily be extended to build a
nice infix-notation direct-libm-based floating point expression
evaluator library or utility.

Example usage:

$ ./ccalc '-1 csqrt 1 atan 2 * 3 / * cexp p'
0.866025403784438596588 + i*0.5

Rich

[-- Attachment #2: ccalc.c --]
[-- Type: text/plain, Size: 5058 bytes --]

/*
 * ccalc.c - C calculator (RPN)
 * by Rich Felker, placed in the public domain
 *
 * This program implements a simplistic RPN calculator that directly
 * uses the system's C floating point operators (on complex long
 * double) and functions (on their respective argument types) for the
 * purpose of comparing and testing implementation behavior. In
 * addition the the standard functions, the following commands are
 * available:
 *
 * p - print the top of the stack to DECIMAL_DIG significant figures
 * x - print the current floating point exception flags' status
 *
 * No error checking for stack overflow or underflow is performed in
 * this version, so it is not safe for use with untrusted input.
 */

#include <math.h>
#include <complex.h>
#include <float.h>
#include <fenv.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define COMPLEX 1
#define LDBL 2
#define FLOAT 4

#define FUNCR(x) \
{ #x, (void (*)())x, 0 }, \
{ #x "f", (void (*)())x##f, FLOAT }, \
{ #x "l", (void (*)())x##l, LDBL }

#define FUNCC(x) \
{ #x, (void (*)())x, COMPLEX }, \
{ #x "f", (void (*)())x##f, COMPLEX|FLOAT }, \
{ #x "l", (void (*)())x##l, COMPLEX|LDBL }

#define FUNC(x) FUNCR(x), FUNCC(c##x)

struct func {
	char name[16];
	void (*addr)();
	unsigned flags;
};

/* missing: fma, frexp, ilogb, ldexp, llrint, llround, lrint, lround,
 * nexttoward, remquo, scalbln, scalbn */

static const struct func unary_funcs[] = {
	FUNC(acos), FUNC(acosh), FUNC(asin), FUNC(asinh),
	FUNC(atan), FUNC(atanh), FUNC(cos), FUNC(cosh),
	FUNC(sin), FUNC(sinh), FUNC(tan), FUNC(tanh),
	FUNC(exp), FUNC(log), FUNC(sqrt),

	FUNCR(cbrt), FUNCR(fabs),
	FUNCR(exp2), FUNCR(expm1),
	FUNCR(erf), FUNCR(erfc),
	FUNCR(log10), FUNCR(log1p), FUNCR(log2), FUNCR(logb),
	FUNCR(ceil), FUNCR(floor), FUNCR(nearbyint),
	FUNCR(rint), FUNCR(round), FUNCR(trunc),
	FUNCR(lgamma), FUNCR(tgamma),

	FUNCC(carg), FUNCC(cimag), FUNCC(creal), FUNCC(cproj), FUNCC(conj),
 
	{0}
};

static const struct func binary_funcs[] = {
	FUNCR(atan2), FUNCR(copysign), FUNCR(fdim),
	FUNCR(fmax), FUNCR(fmin), FUNCR(fmod), FUNCR(hypot),
	FUNC(pow), FUNCR(remainder), FUNCR(nextafter),
	{0}
};

int main(int argc, char **argv)
{
	char token[20000], *end;
	complex long double stack[100];
	int sp = 0;
	const struct func *f;
	FILE *in = stdin;
	int except;

	if (argc>1) in = fmemopen(argv[1], strlen(argv[1]), "r");
	setvbuf(in, 0, _IONBF, 0);

	while (fscanf(in, "%19999s", token)==1) {
		end = 0;
		stack[sp] = strtod(token, &end);
		if (end != token) switch (*end) {
		case 'f': case 'F':
			stack[sp++] = strtof(token, 0);
			continue;
		case 'l': case 'L':
			stack[sp++] = strtold(token, 0);
			continue;
		case 0:
			sp++;
			continue;
		}
		sp--;
		if (!token[1]) switch (token[0]) {
		case '+':
			stack[sp-1] += stack[sp];
			continue;
		case '-':
			stack[sp-1] -= stack[sp];
			continue;
		case '*':
			stack[sp-1] *= stack[sp];
			continue;
		case '/':
			stack[sp-1] /= stack[sp];
			continue;
		case 'p':
			printf("%.*Lg + i*%.*Lg\n",
				DECIMAL_DIG, creall(stack[sp]),
				DECIMAL_DIG, cimagl(stack[sp]));
			sp++;
			continue;
		case 'x':
			except = fetestexcept(FE_ALL_EXCEPT);
			if (!except) printf("NONE");
			if (except & FE_INVALID) printf("INVALID ");
			if (except & FE_OVERFLOW) printf("OVERFLOW ");
			if (except & FE_UNDERFLOW) printf("UNDERFLOW ");
			if (except & FE_INEXACT) printf("INEXACT ");
			putchar('\n');
			sp++;
			continue;
		}
		sp++;
		f = unary_funcs;
		for (; f->addr; f++) {
			if (strcmp(f->name, token)) continue;
			switch(f->flags) {
			case 0:
				stack[sp-1] = ((double (*)(double))f->addr)(stack[sp-1]);
				break;
			case FLOAT:
				stack[sp-1] = ((float (*)(float))f->addr)(stack[sp-1]);
				break;
			case LDBL:
				stack[sp-1] = ((long double (*)(long double))f->addr)(stack[sp-1]);
				break;
			case COMPLEX:
				stack[sp-1] = ((complex double (*)(complex double))f->addr)(stack[sp-1]);
				break;
			case COMPLEX|FLOAT:
				stack[sp-1] = ((complex float (*)(complex float))f->addr)(stack[sp-1]);
				break;
			case COMPLEX|LDBL:
				stack[sp-1] = ((complex long double (*)(complex long double))f->addr)(stack[sp-1]);
				break;
			}
			break;
		}
		f = binary_funcs;
		for (; f->addr; f++) {
			if (strcmp(f->name, token)) continue;
			sp--;
			switch(f->flags) {
			case 0:
				stack[sp-1] = ((double (*)(double,double))f->addr)(stack[sp-1], stack[sp]);
				break;
			case FLOAT:
				stack[sp-1] = ((float (*)(float,float))f->addr)(stack[sp-1], stack[sp]);
				break;
			case LDBL:
				stack[sp-1] = ((long double (*)(long double,long double))f->addr)(stack[sp-1], stack[sp]);
				break;
			case COMPLEX:
				stack[sp-1] = ((complex double (*)(complex double,complex double))f->addr)(stack[sp-1], stack[sp]);
				break;
			case COMPLEX|FLOAT:
				stack[sp-1] = ((complex float (*)(complex float,complex float))f->addr)(stack[sp-1], stack[sp]);
				break;
			case COMPLEX|LDBL:
				stack[sp-1] = ((complex long double (*)(complex long double,complex long double))f->addr)(stack[sp-1], stack[sp]);
				break;
			}
			break;
		}
	}
	return 0;
}