cgv
Loading...
Searching...
No Matches
functions.h
1#pragma once
2
3#include <cmath>
4#include <cassert>
5#include <limits>
6#include <vector>
7
8#include "constants.h"
9
10#include "lib_begin.h"
11
12namespace cgv {
13namespace math {
14
16template<typename T>
17T sgn(const T& v) { return (T(0) < v) - (v < T(0)); }
19template<typename T>
20T sign(const T&v) { return (v >= T(0)) ? T(1) : T(-1); }
22template<typename T>
23T sign(const T& a, const T& b) { return (b < 0) ? -std::abs(a) : std::abs(a); }
25template<typename T>
26T step(const T& a, const T& b) { return (b < a) ? T(0) : T(1); }
28template<typename T>
29T plus(const T& v) { return v >= T(0) ? v : T(0); }
31template<typename T>
32T clamp(const T& v, const T& a, const T& b) { return v > b ? b : (v < a ? a : v); }
34template<typename T>
35T saturate(const T& v) { return v > T(1) ? T(1) : (v < T(0) ? T(0) : v); }
37template <typename T>
38T lerp(const T& a, const T& b, T t) { return (1 - t) * a + t * b; }
40template<typename T>
41T map(const T& v, const T& a, const T& b, const T& c, const T& d) { return c + (d - c) * ((v - a) / (b - a)); }
42
43namespace detail {
44
46template <typename T>
47T erfccheb(const T z)
48{
49 static const int ncof = 28;
50 static const double cof[28] = {
51 -1.3026537197817094, 6.4196979235649026e-1,
52 1.9476473204185836e-2,-9.561514786808631e-3,-9.46595344482036e-4,
53 3.66839497852761e-4,4.2523324806907e-5,-2.0278578112534e-5,
54 -1.624290004647e-6,1.303655835580e-6,1.5626441722e-8,-8.5238095915e-8,
55 6.529054439e-9,5.059343495e-9,-9.91364156e-10,-2.27365122e-10,
56 9.6467911e-11, 2.394038e-12,-6.886027e-12,8.94487e-13, 3.13092e-13,
57 -1.12708e-13,3.81e-16,7.106e-15,-1.523e-15,-9.4e-17,1.21e-16,-2.8e-17
58 };
59
60 int j;
61 double t, ty, tmp, d = 0.0, dd = 0.0;
62 assert(z >= 0.0);
63
64 t = 2.0 / (2.0 + z);
65 ty = 4.0*t - 2.0;
66 for (j = ncof - 1; j > 0; j--)
67 {
68 tmp = d;
69 d = ty * d - dd + cof[j];
70 dd = tmp;
71 }
72 return static_cast<T>(t*exp(-z * z + 0.5*(cof[0] + ty * d) - dd));
73}
74
75} // namespace detail
76
78template <typename T>
79T erf(const T x)
80{
81 if (x >= 0.)
82 return T(1.0) - detail::erfccheb(x);
83 else
84 return detail::erfccheb(-x) - T(1.0);
85}
86
88template <typename T>
89T erfc(const T x)
90{
91 if (x >= 0.)
92 return detail::erfccheb(x);
93 else
94 return 2.0 - detail::erfccheb(-x);
95}
96
98template <typename T>
99T erfcx(const T x)
100{
101 return exp(x*x)* erfccheb(x);
102}
103
105template <typename T>
106T erfc_inv(const T p)
107{
108 double x, err, t, pp;
109 if (p >= 2.0) return -100.;
110 if (p <= 0.0) return 100.;
111 pp = (p < 1.0) ? p : 2. - p;
112 t = sqrt(-2.*log(pp / 2.));
113 x = -0.70711*((2.30753 + t * 0.27061) / (1. + t * (0.99229 + t * 0.04481)) - t);
114 for (int j = 0; j < 2; j++) {
115 err = erfc(x) - pp;
116 x += err / (1.12837916709551257*exp(-(x*x)) - x * err);
117 }
118 return (p < 1.0 ? x : -x);
119}
120
122template <typename T>
123T erf_inv(const T p)
124{
125 return erfc_inv(1. - p);
126}
127
128
130template <typename T>
131T Phi(const T& x)
132{
133 return (T)(0.5*erfc(-x*sqrt(0.5)));
134}
135
137template <typename T>
138T inv_Phi(const T& p)
139{
140 static T a[] = {(T)-3.969683028665376e+01,(T) 2.209460984245205e+02,(T) -2.759285104469687e+02,(T) 1.383577518672690e+02,(T) -3.066479806614716e+01,(T) 2.506628277459239e+00};
141 static T b[] = {(T)-5.447609879822406e+01,(T) 1.615858368580409e+02,(T) -1.556989798598866e+02,(T) 6.680131188771972e+01,(T) -1.328068155288572e+01};
142 static T c[] = {(T)-7.784894002430293e-03,(T)-3.223964580411365e-01,(T)-2.400758277161838e+00,(T)-2.549732539343734e+00,(T) 4.374664141464968e+00,(T) 2.938163982698783e+00};
143 static T d[] = {(T)7.784695709041462e-03,(T)3.224671290700398e-01,(T)2.445134137142996e+00,(T)3.754408661907416e+00};
144
145 if (p <= 0)
146 return -std::numeric_limits<T>::infinity();
147 if (p >= 1)
148 return std::numeric_limits<T>::infinity();
149 T x;
150 if (p < (T)0.02425) {
151 T q = sqrt(-(T)2.0*log(p));
152 x = (((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
153 ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1);
154 }
155 else if (p <= (T)0.97575) {
156 T q = p - (T)0.5;
157 T r = q*q;
158 x = (((((a(1)*r+a(2))*r+a(3))*r+a(4))*r+a(5))*r+a(6))*q /
159 (((((b(1)*r+b(2))*r+b(3))*r+b(4))*r+b(5))*r+1);
160 }
161 else {
162 T q = sqrt(-(T)2.0*log(1-p));
163 x = -(((((c(1)*q+c(2))*q+c(3))*q+c(4))*q+c(5))*q+c(6)) /
164 ((((d(1)*q+d(2))*q+d(3))*q+d(4))*q+1);
165 }
166 return x;
167}
169extern CGV_API double compute_unit_ball_volume(unsigned n);
171extern CGV_API double compute_ball_volume(unsigned n, double R);
173extern CGV_API double compute_unit_sphere_area(unsigned n);
175extern CGV_API double compute_sphere_area(unsigned n, double R);
176
177//returns v*v
178template <typename T>
179T sqr(const T& v)
180{
181 return v*v;
182}
183
184//return the maximum of a and b
185template <typename T>
186T maximum(const T& a, const T&b)
187{
188 return a > b ? a : b;
189}
190
191//return minimum of a and b
192template<typename T>
193T minimum(const T &a, const T&b)
194{
195 return a < b ? a : b;
196}
197
198//return the maximum of a, b and c
199template <typename T>
200T maximum(const T& a, const T&b, const T& c)
201{
202 if(a > b)
203 {
204 return a > c ? a : c;
205 }
206 else
207 {
208 return b > c ? b : c;
209 }
210}
211
212//return the minimum of a, b and c
213template <typename T>
214T minimum(const T& a, const T&b, const T& c)
215{
216 if(a < b)
217 {
218 return a < c ? a : c;
219 }
220 else
221 {
222 return b < c ? b : c;
223 }
224}
225
226//convert angles measured in degrees to angles measured in radians
227template <typename T>
228T deg2rad(T deg)
229{
230 return static_cast<T>(deg * static_cast<T>(cgv::math::constants::pi / 180.0));
231}
232
233//convert angles measured in radians to angles measured in degrees
234template <typename T>
235T rad2deg(T rad)
236{
237 return static_cast<T>(rad * static_cast<T>(180.0 * cgv::math::constants::inv_pi));
238}
239
240//ln(gamma(x))
241extern CGV_API double gamma_ln(const double xx);
242
243//factorial of n (n!)
244extern CGV_API double fac(const int n);
245
246//returns ln(n!)
247extern CGV_API double fac_ln(const int n);
248
249//beta function
250extern CGV_API double beta(const double z, const double w);
251
252//binomial coefficient n over k
253extern CGV_API double nchoosek(const int n, const int k);
254
255} // namespace math
256} // namespace cgv
257
258#include <cgv/config/lib_end.h>
cgv
the cgv namespace
Definition print.h:11