math.c

/*
    This file is part of duckOS.
    
    duckOS is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
    
    duckOS is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
    
    You should have received a copy of the GNU General Public License
    along with duckOS.  If not, see <https://www.gnu.org/licenses/>.
    
    Copyright (c) Byteduck 2016-2020. All rights reserved.
*/

#include "math.h"
#include "stdlib.h"

#define X87_FUNC(name, instr) \
    float name##f(float arg) { float res; asm(#instr : "=t"(res) : "0"(arg)); return res; } \
	double name(double arg) { double res; asm(#instr : "=t"(res) : "0"(arg)); return res; } \
	long double name##l(long double arg) { long double res; asm(#instr : "=t"(res) : "0"(arg)); return res; }

#define X87_FUNC2(name, instr) \
    float name##f(float arga, float argb) { float res; asm(#instr : "=t"(res) : "0"(arga), "u"(argb) : "st(1)"); return res; } \
	double name(double arga, double argb) { double res; asm(#instr : "=t"(res) : "0"(arga), "u"(argb) : "st(1)"); return res; } \
	long double name##l(long double arga, long double argb) { long double res; asm(#instr : "=t"(res) : "0"(arga), "u"(argb) : "st(1)"); return res; }

X87_FUNC(sqrt, fsqrt);
X87_FUNC(sin, fsin);
X87_FUNC(cos, fcos);
X87_FUNC(atan, fpatan);
X87_FUNC2(atan2, fpatan);
X87_FUNC(tan, fptan);

// Does this work? No clue.
#define asin_func(name, T) \
	T name(T x) { \
		if (x > 1 || x < -1) \
			return NAN; \
		T a = x, sum = x; \
		for (int i = 0; i < 5; i++) { \
			a = a * (((2 * i - 1) * (2 * i - 1)) / ((2 * i) * (2 * i + 1) * x * x)); \
			sum += a; \
		} \
		return sum; \
	}
asin_func(asin, double);
asin_func(asinf, float);
asin_func(asinl, long double);

double acos(double x) { return 0.5 * M_PI - asin(x); }
float acosf(float x) { return 0.5f * (float) M_PI - asinf(x);}
long double acosl(long double x) { return 0.5f * (long double) M_PI - asinl(x);}

double ceil(double x) { return __builtin_ceil(x); }
float ceilf(float x) { return __builtin_ceilf(x); }
long double ceill(long double x) { return __builtin_ceill(x); }
double floor(double x) { return __builtin_floor(x); }
float floorf(float x) { return __builtin_floorf(x); }
long double floorl(long double x) { return __builtin_floorl(x); }
double trunc(double x) { return (double) ((int) x); }
float truncf(float x) { return (float) ((int) x); }
long double truncl(long double x) { return (long double) ((long) x); }

double copysign(double dest, double src) {
	union {
		double val;
		uint64_t bits;
	} *destu = (void*) &dest, *srcu = (void*) &src;
	destu->bits = (destu->bits & ~(UINT64_C(1) << 63)) | (srcu->bits & (UINT64_C(1) << 63));
	return destu->val;
}

float copysignf(float dest, float src) {
	union {
		float val;
		uint32_t bits;
	} *destu = (void*) &dest, *srcu = (void*) &src;
	destu->bits = (destu->bits & ~(UINT32_C(1) << 31)) | (srcu->bits & (UINT32_C(1) << 31));
	return destu->val;
}

#define expfunc(name, preamble) \
	double name(double arg) { double res; asm(preamble"\nfld1\nfld %%st(1)\nfprem\nf2xm1\nfaddp\nfscale\nfstp %%st(1)" : "=t"(res) : "0"(arg)); return res; } \
	float name##f(float arg) { float res; asm(preamble"\nfld1\nfld %%st(1)\nfprem\nf2xm1\nfaddp\nfscale\nfstp %%st(1)" : "=t"(res) : "0"(arg)); return res; } \
	long double name##l(long double arg) { long double res; asm(preamble"\nfld1\nfld %%st(1)\nfprem\nf2xm1\nfaddp\nfscale\nfstp %%st(1)" : "=t"(res) : "0"(arg)); return res; }
expfunc(exp, "fldl2e\nfmulp");
expfunc(exp2, "");

double ldexp(double x, int exp) { return x * (0x1u >> exp); }
float ldexpf(float x, int exp) { return x * (0x1u >> exp); }

double fmod(double x, double y) { return fmodf(x, y); }
float fmodf(float x, float y) {
	uint16_t fpu;
	do {
		asm("fprem\nfnstsw %%ax" : "+t"(x), "=a"(fpu) : "u"(y));
	} while (fpu & 0x400);
	return x;
}

#define logfunc(name, ldinstr) \
	double name(double arg) { double res; asm(#ldinstr"\nfxch %%st(1)\nfyl2x" : "=t"(res) : "0"(arg)); return res; } \
	float name##f(float arg) { float res; asm(#ldinstr"\nfxch %%st(1)\nfyl2x" : "=t"(res) : "0"(arg)); return res; } \
	long double name##l(long double arg) { long double res; asm(#ldinstr"\nfxch %%st(1)\nfyl2x" : "=t"(res) : "0"(arg)); return res; }
logfunc(log, fldln2);
logfunc(log10, fldlg2);
logfunc(log2, fld1);

#define powfunc(name, T, exp2func, log2func) \
	T name(T a, T b) { \
		if (b == 0) \
			return 1; \
		if (b == 1) \
			return a; \
		if (a == 0) \
			return 0; \
		if (__builtin_isnan(b)) \
			return b; \
		int bint = b; \
		if ((T) bint == b) { \
			T out = a; \
			for (int i = 0; i < bint - 1; i++) \
				out *= a; \
			if (a < 0) \
				a = 1.01 / a; \
			return a; \
		} \
		return exp2func(b * log2func(a)); \
	}
powfunc(pow, double, exp2, log2);
powfunc(powf, float, exp2f, log2f);

#define roundfunc(name, Ta, Tb, Tc) \
	Ta name(double arg) { return (Ta) (arg > 0 ? floor(arg + 0.5) : ceil(arg - 0.5)); } \
	Tb name##f(float arg) { return (Tb) (arg > 0 ? floorf(arg + 0.5) : ceilf(arg - 0.5)); } \
	Tc name##l(long double arg) { return (Tc) (arg > 0 ? floorl(arg + 0.5) : ceill(arg - 0.5)); }
roundfunc(round, double, float, long double);
roundfunc(lround, long int, long int, long int);
roundfunc(llround, long long int, long long int, long long int);

#define scalbnfunc(name, T) T name(T arg, int exp) { return arg * (1 << exp); }
scalbnfunc(scalbn, double);
scalbnfunc(scalbnf, float);
scalbnfunc(scalbnl, long double);

#define modffunc(name, T, IntT) T name(T arg, T* iptr) { \
		*iptr = (T) (IntT) arg; \
		float ret = arg - *iptr; \
		return arg < 0 ? ret - *iptr : ret; \
	}
modffunc(modf, double, long);
modffunc(modff, float, int);
modffunc(modfl, long double, long long);

int ilogb(double arg) {
	if (arg == 0)
		return FP_ILOGB0;
	else if (isinf(arg))
		return INT_MAX;
	else if (isnan(arg))
		return FP_ILOGBNAN;

	union {
		double val;
		uint64_t bits;
	} *argp = (void*) &arg;
	return (int) ((argp->bits >> 52) & 0x7ff) - 1023;
}

int ilogbf(float arg) {
	if (arg == 0)
		return FP_ILOGB0;
	else if (isinf(arg))
		return INT_MAX;
	else if (isnan(arg))
		return FP_ILOGBNAN;

	union {
		float val;
		uint32_t bits;
	} *argp = (void*) &arg;
	return (int) ((argp->bits >> 23) & 0xff) - 127;
}

int ilogbl(long double arg) {
	if (arg == 0)
		return FP_ILOGB0;
	else if (isinf(arg))
		return INT_MAX;
	else if (isnan(arg))
		return FP_ILOGBNAN;

	union {
		long double val;
		uint64_t bits[6];
	} __attribute__((packed)) *argp = (void*) &arg;
	return (int) (argp->bits[4] & 0x7FF) - 16383;
}

#define frexpfunc(suffix, T) T frexp##suffix(T arg, int* exp) { *exp = (arg == 0) ? 0 : (1 + ilogb##suffix(arg)); return scalbn##suffix(arg, -(*exp)); }
frexpfunc( , double);
frexpfunc(f, float);
frexpfunc(l, long double);

#define hypotfunc(suffix, T) T hypot##suffix(T x, T y) { return sqrt##suffix(x*x + y*y); }
hypotfunc( , double);
hypotfunc(f, float);
hypotfunc(l, long double);