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)
99 return log(x) * M_LOG10E;
107 return log(x) * M_LOG2E;
111 return floor(log2(fabs(x)));
115 return '1 0 0' * (f - trunc(f)) + '0 1 0' * trunc(f);
118 float scalbn(float x, int n)
120 return x * pow(2, n);
125 return copysign(pow(fabs(x), 1.0/3.0), x);
127 float hypot(float x, float y)
129 return sqrt(x*x + y*y);
134 // approximation taken from wikipedia
137 return copysign(sqrt(1 - exp(-y * (1.273239544735163 + 0.14001228868667 * y) / (1 + 0.14001228868667 * y))), x);
143 vector lgamma(float x)
145 // TODO improve accuracy
147 return fabs(x) * '1 0 0' + copysign(1, x) * '0 1 0';
148 if(x < 1 && x == floor(x))
149 return nan("gamma") * '1 1 1';
154 // reflection formula:
155 // gamma(1-z) * gamma(z) = pi / sin(pi*z)
156 // lgamma(1-z) + lgamma(z) = log(pi) - log(sin(pi*z))
157 // sign of gamma(1-z) = sign of gamma(z) * sign of sin(pi*z)
159 v.x = log(M_PI) - log(fabs(v.z)) - v.x;
166 return lgamma(x + 1) - log(x) * '1 0 0';
168 return (0.5 * log(2 * M_PI * x) + x * (log(x) - 1)) * '1 0 0' + '0 1 0';
170 float tgamma(float x)
174 return exp(v.x) * v.y;
186 float pymod(float x, float y)
188 return x - y * floor(x / y);
191 float nearbyint(float x)
197 return (x>=0) ? floor(x) : ceil(x);
200 float fmod(float x, float y)
202 return x - y * trunc(x / y);
204 float remainder(float x, float y)
206 return x - y * rint(x / y);
208 vector remquo(float x, float y)
217 float copysign(float x, float y)
219 return fabs(x) * ((y>0) ? 1 : -1);
221 float nan(string tag)
225 float nextafter(float x, float y)
229 return nan("nextafter");
231 return -nextafter(-x, -y);
232 // now we know that x < y
233 // so we need the next number > x
235 d = max(fabs(x), 0.00000000000000000000001);
246 float nexttoward(float x, float y)
248 return nextafter(x, y);
251 float fdim(float x, float y)
255 float fmax(float x, float y)
259 float fmin(float x, float y)
263 float fma(float x, float y, float z)
268 int isgreater(float x, float y)
272 int isgreaterequal(float x, float y)
276 int isless(float x, float y)
280 int islessequal(float x, float y)
284 int islessgreater(float x, float y)
286 return x < y || x > y;
288 int isunordered(float x, float y)
290 return !(x < y || x == y || x > y);
293 vector cross(vector a, vector b)
296 '1 0 0' * (a.y * b.z - a.z * b.y)
297 + '0 1 0' * (a.z * b.x - a.x * b.z)
298 + '0 0 1' * (a.x * b.y - a.y * b.x);