rational.c

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <inttypes.h>

struct Rational {
	ssize_t numerator;
	 size_t denominator;
};

enum {
	RAT_STR_LEN = 42,
};

struct Rational const rat_pi = {
	.numerator   = 355,
	.denominator = 113,
};

struct Rational const rat_e = {
	.numerator   = 2721,
	.denominator = 1001,
};


static inline ssize_t _abs(ssize_t const a) {
	return a<0? -a : a;
}

ssize_t fast_ipow(ssize_t x, ssize_t expo) {
	if( expo < 0 ) {
		return 0;
	} else if( expo==0 || x==1 ) {
		return 1;
	} else if( expo >= 64 ) {
		return 0x7fFFffFFFFffFFffL;
	}
	ssize_t b = 1;
	ssize_t acc = 1;
	while (b <= expo) {
		acc *= x * (expo & b);
		x *= x;
		b <<= 1;
	}
	return acc;
}

size_t gcd(size_t a, size_t b) {
	while( b != 0 ) {
		size_t const temp = a % b;
		a = b;
		b = temp;
	}
	return a;
}

size_t gcf(size_t const a, size_t const b) {
	return (a*b) / gcd(a,b);
}

struct Rational make_rational(ssize_t numerator, size_t denominator) {
	if( denominator==0 ) {
		return ( struct Rational ){ -1, 1 };
	}
	size_t const divisor = gcd(( size_t )(abs(numerator)), denominator);
	if( divisor != 0 ) {
		numerator   /= divisor;
		denominator /= divisor;
	}
	return ( struct Rational ){ numerator, denominator };
}

struct Rational rational_reciprocal(struct Rational const a) {
	return make_rational(a.denominator, a.numerator);
}

int rational_eq(struct Rational const a, struct Rational const b) {
	return a.denominator==b.denominator && a.numerator==b.numerator;
}

struct Rational rational_mult(struct Rational const a, struct Rational const b) {
	return make_rational(a.numerator * b.numerator, a.denominator * b.denominator);
}

struct Rational rational_div(struct Rational const a, struct Rational const b) {
	return rational_mult(a, rational_reciprocal(b));
}

struct Rational rational_add(struct Rational const a, struct Rational const b) {
	if( a.denominator==b.denominator ) {
		return make_rational(a.numerator + b.numerator, a.denominator);
	}
	size_t const common_factor_denom = gcf(a.denominator, b.denominator);
	ssize_t const sum = (a.numerator * (common_factor_denom / a.denominator)) + (b.numerator * (common_factor_denom / b.denominator));
	return make_rational(sum, common_factor_denom);
}

struct Rational rational_sub(struct Rational const a, struct Rational const b) {
	return rational_add(a, ( struct Rational ){ -b.numerator, b.denominator });
}

struct Rational rational_pow(struct Rational a, ssize_t pwr) {
	if( pwr < 0 ) {
		a = rational_reciprocal(a);
	}
	a.numerator = fast_ipow(a.numerator, _abs(pwr));
	a.denominator = ( size_t )(fast_ipow(( ssize_t )(a.denominator), _abs(pwr)));
	return a;
}

struct Rational rational_factorial(size_t const n) {
	struct Rational result = make_rational(1, 1);
	for( size_t i=1; i <= n; i++ ) {
		result = rational_mult(result, make_rational(i, 1));
	}
	return result;
}

struct Rational rational_sin(struct Rational const x) {
	/// TODO:
	return make_rational(0, 1);
}


double rational_to_float(struct Rational const r) {
	if( r.denominator==0 ) {
		union {
			uint64_t u64;
			double   f64;
		} const c = {-1ULL};
		return c.f64;
	}
	return ( double )(r.numerator) / ( double )(r.denominator);
}

struct Rational float_to_rational(double const f) {
	double const tolerance = 1e-6;
	double value = f;
	size_t denominator = 1;
	while( fabs(value - round(value)) > tolerance && denominator < 1000000 ) {
		value *= 10;
		denominator *= 10;
	}
	ssize_t numerator = ( ssize_t )(round(value));
	size_t const divisor = gcd(( size_t )(abs(numerator)), denominator);
	numerator   /= divisor;
	denominator /= divisor;
	return ( struct Rational ){numerator, denominator};
}

char *rational_to_string(struct Rational const r, char (*buf)[RAT_STR_LEN]) {
	snprintf(&(*buf[0]), sizeof *buf, "%zd/%zu", r.numerator, r.denominator);
	return *buf;
}

void rational_to_stream(struct Rational const r, FILE *const stream) {
	char r_str[RAT_STR_LEN] = {0};
	fprintf(stream, "%s", rational_to_string(r, &r_str));
}

int main(void) {
	#define RAT_STR    (char[RAT_STR_LEN]){0}
	double const f = 4.5;
	struct Rational const r1 = float_to_rational(f);
	printf("double: %f approximated as %s\n", f, rational_to_string(r1, &RAT_STR));
	
	struct Rational const r2 = make_rational(32, 128);
	printf("r2: %s\n", rational_to_string(r2, &RAT_STR));
	printf("test multiplying: %s * %s = %s\n", rational_to_string(r1, &RAT_STR), rational_to_string(r2, &RAT_STR), rational_to_string(rational_mult(r1, r2), &RAT_STR));
	printf("test dividing: %s / %s = %s\n", rational_to_string(r1, &RAT_STR), rational_to_string(r2, &RAT_STR), rational_to_string(rational_div(r1, r2), &RAT_STR));
	printf("test adding: %s + %s = %s\n", rational_to_string(r1, &RAT_STR), rational_to_string(r2, &RAT_STR), rational_to_string(rational_add(r1, r2), &RAT_STR));
	printf("test power: %s^%i = %s\n", rational_to_string(r1, &RAT_STR), 3, rational_to_string(rational_pow(r1, 3), &RAT_STR));
	printf("test negative power: %s^%i = %s\n", rational_to_string(r1, &RAT_STR), -3, rational_to_string(rational_pow(r1, -3), &RAT_STR));
	
	printf("test multiplying pi: %s * %s = %s\n", rational_to_string(r1, &RAT_STR), rational_to_string(rat_pi, &RAT_STR), rational_to_string(rational_mult(r1, rat_pi), &RAT_STR));
	printf("test multiplying pi out: %s * %s = %s\n", rational_to_string(r1, &RAT_STR), rational_to_string(rat_pi, &RAT_STR), rational_to_string(rational_div(rational_mult(r1, rat_pi), r1), &RAT_STR));
	
	struct Rational const one_over_two = make_rational(1,2);
	struct Rational const pi_over_two  = rational_mult(one_over_two, rat_pi);
	printf( "test sin(%s*%s) = sin(%s) = %s\n", rational_to_string(one_over_two, &RAT_STR), rational_to_string(rat_pi, &RAT_STR), rational_to_string(pi_over_two, &RAT_STR), rational_to_string(rational_sin(pi_over_two), &RAT_STR) );
}