3 #include "../dpdefs/csprogsdefs.qh"
6 #include "../dpdefs/dpextensions.qh"
7 #include "../dpdefs/progsdefs.qh"
10 int fpclassify(float x)
20 bool isfinite(float x)
22 return !(isnan(x) || isinf(x));
26 return (x != 0) && (x + x == x);
34 bool isnormal(float x)
45 return log(x + sqrt(x*x - 1));
49 return log(x + sqrt(x*x + 1));
53 return 0.5 * log((1+x) / (1-x));
57 return 0.5 * (exp(x) + exp(-x));
61 return 0.5 * (exp(x) - exp(-x));
65 return sinh(x) / cosh(x);
91 return floor(log2(fabs(x)));
93 float ldexp(float x, int e)
97 float logn(float x, float base)
99 return log(x) / log(base);
103 return log(x) * M_LOG10E;
111 return log(x) * M_LOG2E;
115 return floor(log2(fabs(x)));
119 return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
122 float scalbn(float x, int n)
124 return x * pow(2, n);
129 return copysign(pow(fabs(x), 1.0/3.0), x);
131 float hypot(float x, float y)
133 return sqrt(x*x + y*y);
138 // approximation taken from wikipedia
141 return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
147 vector lgamma(float x)
149 // TODO improve accuracy
151 return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0';
152 if(x < 1 && x == floor(x))
153 return nan("gamma") * '1 1 1';
158 // reflection formula:
159 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
160 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
161 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
163 v.x = log(M_PI) - log(fabs(v.z)) - v.x;
170 return lgamma(x + 1) - log(x) * '1 0 0';
172 return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
174 float tgamma(float x)
178 return exp(v.x) * v.y;
190 float pymod(float x, float y)
192 return x - y * floor(x / y);
195 float nearbyint(float x)
201 return (x>=0) ? floor(x) : ceil(x);
204 float fmod(float x, float y)
206 return x - y * trunc(x / y);
208 float remainder(float x, float y)
210 return x - y * rint(x / y);
212 vector remquo(float x, float y)
221 float copysign(float x, float y)
223 return fabs(x) * ((y>0) ? 1 : -1);
225 float nan(string tag)
229 float nextafter(float x, float y)
233 return nan("nextafter");
235 return -nextafter(-x, -y);
236 // now we know that x < y
237 // so we need the next number > x
239 d = max(fabs(x), 0.00000000000000000000001);
250 float nexttoward(float x, float y)
252 return nextafter(x, y);
255 float fdim(float x, float y)
259 float fmax(float x, float y)
263 float fmin(float x, float y)
267 float fma(float x, float y, float z)
272 int isgreater(float x, float y)
276 int isgreaterequal(float x, float y)
280 int isless(float x, float y)
284 int islessequal(float x, float y)
288 int islessgreater(float x, float y)
290 return x < y || x > y;
292 int isunordered(float x, float y)
294 return !(x < y || x == y || x > y);
297 vector cross(vector a, vector b)
300 '1 0 0' * (a.y * b.z - a.z * b.y)
301 + '0 1 0' * (a.z * b.x - a.x * b.z)
302 + '0 0 1' * (a.x * b.y - a.y * b.x);