autoppl  v0.8
A C++ template library for probabilistic programming
density.hpp
Go to the documentation of this file.
1 #pragma once
2 #include <cmath>
3 #include <autoppl/math/math.hpp>
5 #include <fastad_bits/util/type_traits.hpp>
6 #include <Eigen/Dense>
7 
8 // MSVC does not seem to support M_PI
9 #ifndef M_PI
10 #define M_PI 3.14159265358979323846
11 #endif
12 
13 namespace ppl {
14 namespace math {
15 
17 
19 // Compile-time Constants
21 
22 inline constexpr double SQRT_TWO_PI =
23  2.506628274631000502415765284811045;
24 inline constexpr double LOG_SQRT_TWO_PI =
25  0.918938533204672741780329736405617;
26 
28 // Normal Density
30 
31 // sss
32 template <class XType
33  , class MeanType
34  , class SigmaType
35  , class = std::enable_if_t<
36  std::is_arithmetic_v<XType> &&
37  std::is_arithmetic_v<MeanType> &&
38  std::is_arithmetic_v<SigmaType>
39  > >
40 inline dist_value_t normal_pdf(const XType& x,
41  const MeanType& mean,
42  const SigmaType& sigma)
43 {
44  if (sigma <= 0) return math::neg_inf<dist_value_t>;
45  dist_value_t z = (x - mean) / sigma;
46  return std::exp(-0.5 * z * z) / (sigma * SQRT_TWO_PI);
47 }
48 
49 // vss
50 template <class XType
51  , class MeanType
52  , class SigmaType
53  , class = std::enable_if_t<
54  std::is_arithmetic_v<MeanType> &&
55  std::is_arithmetic_v<SigmaType>
56  > >
57 inline dist_value_t normal_pdf(const Eigen::MatrixBase<XType>& x,
58  const MeanType& mean,
59  const SigmaType& sigma)
60 {
61  if (sigma <= 0) return math::neg_inf<dist_value_t>;
62  dist_value_t z_sq = (x.array() - mean).matrix().squaredNorm() / (sigma * sigma);
63  return std::exp(-0.5 * z_sq) / std::pow(sigma * SQRT_TWO_PI, x.size());
64 }
65 
66 // vvs
67 template <class XType
68  , class MeanType
69  , class SigmaType
70  , class = std::enable_if_t<
71  std::is_arithmetic_v<SigmaType>
72  > >
73 inline dist_value_t normal_pdf(const Eigen::MatrixBase<XType>& x,
74  const Eigen::MatrixBase<MeanType>& mean,
75  const SigmaType& sigma)
76 {
77  static_assert(ad::util::is_eigen_vector_v<MeanType>);
78  assert(x.size() == mean.size());
79  if (sigma <= 0) return math::neg_inf<dist_value_t>;
80  dist_value_t z_sq = (x.array() - mean.array()).matrix().squaredNorm() / (sigma * sigma);
81  return std::exp(-0.5 * z_sq) / std::pow(sigma * SQRT_TWO_PI, x.size());
82 }
83 
84 // vsv, vsm
85 template <class XType
86  , class MeanType
87  , class SigmaType
88  , class = std::enable_if_t<
89  std::is_arithmetic_v<MeanType>
90  > >
91 inline dist_value_t normal_pdf(const Eigen::MatrixBase<XType>& x,
92  const MeanType& mean,
93  const Eigen::MatrixBase<SigmaType>& sigma)
94 {
95  if constexpr (ad::util::is_eigen_vector_v<SigmaType>) {
96  assert(x.size() == sigma.size());
97  if ((sigma.array() <= 0).any()) return math::neg_inf<dist_value_t>;
98  dist_value_t z_sq = ((x.array() - mean)/sigma.array()).matrix().squaredNorm();
99  return std::exp(-0.5 * z_sq) / (std::pow(SQRT_TWO_PI, x.size()) * sigma.array().prod());
100 
101  } else if constexpr (ad::util::is_eigen_matrix_v<SigmaType>) {
102  assert(x.size() == sigma.rows());
103  Eigen::LLT<Eigen::MatrixXd> llt(sigma);
104  if (llt.info() != Eigen::Success) return math::neg_inf<dist_value_t>;
105  dist_value_t z_sq = llt.matrixL().solve((x.array() - mean).matrix()).squaredNorm();
106  dist_value_t det = llt.matrixL().determinant();
107  return std::exp(-0.5 * z_sq) / (std::pow(SQRT_TWO_PI, x.size()) * det);
108  }
109 }
110 
111 // vvv, vvm
112 template <class XType
113  , class MeanType
114  , class SigmaType>
115 inline dist_value_t normal_pdf(const Eigen::MatrixBase<XType>& x,
116  const Eigen::MatrixBase<MeanType>& mean,
117  const Eigen::MatrixBase<SigmaType>& sigma)
118 {
119  static_assert(ad::util::is_eigen_vector_v<MeanType>);
120 
121  if constexpr (ad::util::is_eigen_vector_v<SigmaType>) {
122  assert(x.size() == sigma.size());
123  assert(x.size() == mean.size());
124  if ((sigma.array() <= 0).any()) return math::neg_inf<dist_value_t>;
125  dist_value_t z_sq = ((x.array() - mean.array())/sigma.array()).matrix().squaredNorm();
126  return std::exp(-0.5 * z_sq) / (std::pow(SQRT_TWO_PI, x.size()) * sigma.array().prod());
127  }
128 
129  else if constexpr (ad::util::is_eigen_matrix_v<SigmaType>) {
130  assert(x.size() == sigma.rows());
131  assert(x.size() == mean.size());
132  Eigen::LLT<Eigen::MatrixXd> llt(sigma);
133  if (llt.info() != Eigen::Success) return math::neg_inf<dist_value_t>;
134  dist_value_t z_sq = llt.matrixL().solve((x.array() - mean.array()).matrix()).squaredNorm();
135  dist_value_t det = llt.matrixL().determinant();
136  return std::exp(-0.5 * z_sq) / (std::pow(SQRT_TWO_PI, x.size()) * det);
137  }
138 }
139 
140 // sss
141 template <class XType
142  , class MeanType
143  , class SigmaType
144  , class = std::enable_if_t<
145  std::is_arithmetic_v<XType> &&
146  std::is_arithmetic_v<MeanType> &&
147  std::is_arithmetic_v<SigmaType>
148  > >
149 inline dist_value_t normal_log_pdf(const XType& x,
150  const MeanType& mean,
151  const SigmaType& sigma)
152 {
153  if (sigma <= 0) return math::neg_inf<dist_value_t>;
154  dist_value_t z = (x - mean) / sigma;
155  return -0.5 * z * z - std::log(sigma) - LOG_SQRT_TWO_PI;
156 }
157 
158 // vss
159 template <class XType
160  , class MeanType
161  , class SigmaType
162  , class = std::enable_if_t<
163  std::is_arithmetic_v<MeanType> &&
164  std::is_arithmetic_v<SigmaType>
165  > >
166 inline dist_value_t normal_log_pdf(const Eigen::MatrixBase<XType>& x,
167  const MeanType& mean,
168  const SigmaType& sigma)
169 {
170  if (sigma <= 0) return math::neg_inf<dist_value_t>;
171  dist_value_t z_sq = (x.array() - mean).matrix().squaredNorm() / (sigma * sigma);
172  return -0.5 * z_sq - (x.size() * (std::log(sigma) + LOG_SQRT_TWO_PI));
173 }
174 
175 // vvs
176 template <class XType
177  , class MeanType
178  , class SigmaType
179  , class = std::enable_if_t<
180  std::is_arithmetic_v<SigmaType>
181  > >
182 inline dist_value_t normal_log_pdf(const Eigen::MatrixBase<XType>& x,
183  const Eigen::MatrixBase<MeanType>& mean,
184  const SigmaType& sigma)
185 {
186  static_assert(ad::util::is_eigen_vector_v<MeanType>);
187  assert(x.size() == mean.size());
188  if (sigma <= 0) return math::neg_inf<dist_value_t>;
189  dist_value_t z_sq = (x.array() - mean.array()).matrix().squaredNorm() / (sigma * sigma);
190  return -0.5 * z_sq - (x.size() * (std::log(sigma) + LOG_SQRT_TWO_PI));
191 }
192 
193 // vsv, vsm
194 template <class XType
195  , class MeanType
196  , class SigmaType
197  , class = std::enable_if_t<
198  std::is_arithmetic_v<MeanType>
199  > >
200 inline dist_value_t normal_log_pdf(const Eigen::MatrixBase<XType>& x,
201  const MeanType& mean,
202  const Eigen::MatrixBase<SigmaType>& sigma)
203 {
204  if constexpr (ad::util::is_eigen_vector_v<SigmaType>) {
205  assert(x.size() == sigma.size());
206  if ((sigma.array() <= 0).any()) return math::neg_inf<dist_value_t>;
207  dist_value_t z_sq = ((x.array() - mean)/sigma.array()).matrix().squaredNorm();
208  return -0.5 * z_sq - (x.size() * LOG_SQRT_TWO_PI) - std::log(sigma.array().prod());
209  }
210 
211  else if constexpr (ad::util::is_eigen_matrix_v<SigmaType>) {
212  assert(x.size() == sigma.rows());
213  Eigen::LLT<Eigen::MatrixXd> llt(sigma);
214  if (llt.info() != Eigen::Success) return math::neg_inf<dist_value_t>;
215  dist_value_t z_sq = llt.matrixL().solve((x.array() - mean).matrix()).squaredNorm();
216  dist_value_t det = llt.matrixL().determinant();
217  return -0.5 * z_sq - (x.size() * LOG_SQRT_TWO_PI) - std::log(det);
218  }
219 }
220 
221 // vvv, vvm
222 template <class XType
223  , class MeanType
224  , class SigmaType>
225 inline dist_value_t normal_log_pdf(const Eigen::MatrixBase<XType>& x,
226  const Eigen::MatrixBase<MeanType>& mean,
227  const Eigen::MatrixBase<SigmaType>& sigma)
228 {
229  if constexpr (ad::util::is_eigen_vector_v<SigmaType>) {
230  assert(x.size() == sigma.size());
231  assert(x.size() == mean.size());
232  if ((sigma.array() <= 0).any()) return math::neg_inf<dist_value_t>;
233  dist_value_t z_sq = ((x.array() - mean.array())/sigma.array()).matrix().squaredNorm();
234  return -0.5 * z_sq - (x.size() * LOG_SQRT_TWO_PI) - std::log(sigma.array().prod());
235  }
236 
237  else if constexpr (ad::util::is_eigen_matrix_v<SigmaType>) {
238  assert(x.size() == sigma.rows());
239  assert(x.size() == mean.size());
240  Eigen::LLT<Eigen::MatrixXd> llt(sigma);
241  if (llt.info() != Eigen::Success) return math::neg_inf<dist_value_t>;
242  dist_value_t z_sq = llt.matrixL().solve((x.array() - mean.array()).matrix()).squaredNorm();
243  dist_value_t det = llt.matrixL().determinant();
244  return -0.5 * z_sq - (x.size() * LOG_SQRT_TWO_PI) - std::log(det);
245  }
246 }
247 
249 // Cauchy Density
251 
252 // sss
253 template <class XType
254  , class LocType
255  , class ScaleType
256  , class = std::enable_if_t<
257  std::is_arithmetic_v<XType> &&
258  std::is_arithmetic_v<LocType> &&
259  std::is_arithmetic_v<ScaleType>
260  > >
261 inline dist_value_t cauchy_log_pdf(const XType& x,
262  const LocType& loc,
263  const ScaleType& scale)
264 {
265  auto diff = x - loc;
266  return (scale > 0) ? -std::log(scale + diff * diff / scale) : neg_inf<double>;
267 }
268 
269 // vss
270 template <class XType
271  , class LocType
272  , class ScaleType
273  , class = std::enable_if_t<
274  std::is_arithmetic_v<LocType> &&
275  std::is_arithmetic_v<ScaleType>
276  > >
277 inline dist_value_t cauchy_log_pdf(const Eigen::MatrixBase<XType>& x,
278  const LocType& loc,
279  const ScaleType& scale)
280 {
281  bool cond = scale > 0.;
282  auto diff = x.array() - loc;
283  return cond ? -(scale + (1./scale) * diff * diff).log().sum() :
284  neg_inf<double>;
285 }
286 
287 // vvs
288 template <class XType
289  , class LocType
290  , class ScaleType
291  , class = std::enable_if_t<
292  std::is_arithmetic_v<ScaleType>
293  > >
294 inline dist_value_t cauchy_log_pdf(const Eigen::MatrixBase<XType>& x,
295  const Eigen::MatrixBase<LocType>& loc,
296  const ScaleType& scale)
297 {
298  bool cond = scale > 0.;
299  auto diff = x.array() - loc.array();
300  return cond ? -(scale + (1./scale) * diff * diff).log().sum() : neg_inf<double>;
301 }
302 
303 // vsv
304 template <class XType
305  , class LocType
306  , class ScaleType
307  , class = std::enable_if_t<
308  std::is_arithmetic_v<LocType>
309  > >
310 inline dist_value_t cauchy_log_pdf(const Eigen::MatrixBase<XType>& x,
311  const LocType& loc,
312  const Eigen::MatrixBase<ScaleType>& scale)
313 {
314  bool cond = (scale.array() > 0.).all();
315  auto diff = x.array() - loc;
316  auto gamma = scale.array();
317  return cond ? -(gamma + (1./gamma) * diff * diff).log().sum() : neg_inf<double>;
318 }
319 
320 // vvv
321 template <class XType
322  , class LocType
323  , class ScaleType>
324 inline dist_value_t cauchy_log_pdf(const Eigen::MatrixBase<XType>& x,
325  const Eigen::MatrixBase<LocType>& loc,
326  const Eigen::MatrixBase<ScaleType>& scale)
327 {
328  bool cond = (scale.array() > 0.).all();
329  auto diff = x.array() - loc.array();
330  auto gamma = scale.array();
331  return cond ? -(gamma + (1./gamma) * diff * diff).log().sum() : neg_inf<double>;
332 }
333 
335 // Uniform Density
337 
338 // sss
339 template <class XType
340  , class MinType
341  , class MaxType
342  , class = std::enable_if_t<
343  std::is_arithmetic_v<XType> &&
344  std::is_arithmetic_v<MinType> &&
345  std::is_arithmetic_v<MaxType>
346  > >
347 inline dist_value_t uniform_pdf(const XType& x,
348  const MinType& min,
349  const MaxType& max)
350 {
351  return (min < x && x < max) ? 1. / (max - min) : 0;
352 }
353 
354 // vss
355 template <class XType
356  , class MinType
357  , class MaxType
358  , class = std::enable_if_t<
359  std::is_arithmetic_v<MinType> &&
360  std::is_arithmetic_v<MaxType>
361  > >
362 inline dist_value_t uniform_pdf(const Eigen::MatrixBase<XType>& x,
363  const MinType& min,
364  const MaxType& max)
365 {
366  bool cond = (min < x.array()).all() && (x.array() < max).all();
367  return cond ? std::pow(1./(max-min), x.size()) : 0;
368 }
369 
370 // vvs
371 template <class XType
372  , class MinType
373  , class MaxType
374  , class = std::enable_if_t<
375  std::is_arithmetic_v<MaxType>
376  > >
377 inline dist_value_t uniform_pdf(const Eigen::MatrixBase<XType>& x,
378  const Eigen::MatrixBase<MinType>& min,
379  const MaxType& max)
380 {
381  bool cond = (min.array() < x.array()).all() && (x.array() < max).all();
382  return cond ? (1./(max-min.array())).prod() : 0;
383 }
384 
385 // vsv
386 template <class XType
387  , class MinType
388  , class MaxType
389  , class = std::enable_if_t<
390  std::is_arithmetic_v<MinType>
391  > >
392 inline dist_value_t uniform_pdf(const Eigen::MatrixBase<XType>& x,
393  const MinType& min,
394  const Eigen::MatrixBase<MaxType>& max)
395 {
396  bool cond = (min < x.array()).all() && (x.array() < max.array()).all();
397  return cond ? (1./(max.array()-min)).prod() : 0;
398 }
399 
400 // vvv
401 template <class XType
402  , class MinType
403  , class MaxType>
404 inline dist_value_t uniform_pdf(const Eigen::MatrixBase<XType>& x,
405  const Eigen::MatrixBase<MinType>& min,
406  const Eigen::MatrixBase<MaxType>& max)
407 {
408  bool cond = (min.array() < x.array()).all() && (x.array() < max.array()).all();
409  return cond ? (1./(max.array()-min.array())).prod() : 0;
410 }
411 
412 // sss
413 template <class XType
414  , class MinType
415  , class MaxType
416  , class = std::enable_if_t<
417  std::is_arithmetic_v<XType> &&
418  std::is_arithmetic_v<MinType> &&
419  std::is_arithmetic_v<MaxType>
420  > >
421 inline dist_value_t uniform_log_pdf(const XType& x,
422  const MinType& min,
423  const MaxType& max)
424 {
425  return (min < x && x < max) ? -std::log(max - min) : neg_inf<double>;
426 }
427 
428 // vss
429 template <class XType
430  , class MinType
431  , class MaxType
432  , class = std::enable_if_t<
433  std::is_arithmetic_v<MinType> &&
434  std::is_arithmetic_v<MaxType>
435  > >
436 inline dist_value_t uniform_log_pdf(const Eigen::MatrixBase<XType>& x,
437  const MinType& min,
438  const MaxType& max)
439 {
440  bool cond = (min < x.array()).all() && (x.array() < max).all();
441  return cond ?
442  static_cast<dist_value_t>(x.size())*(-std::log(max-min)) :
443  neg_inf<double>;
444 }
445 
446 // vvs
447 template <class XType
448  , class MinType
449  , class MaxType
450  , class = std::enable_if_t<
451  std::is_arithmetic_v<MaxType>
452  > >
453 inline dist_value_t uniform_log_pdf(const Eigen::MatrixBase<XType>& x,
454  const Eigen::MatrixBase<MinType>& min,
455  const MaxType& max)
456 {
457  bool cond = (min.array() < x.array()).all() && (x.array() < max).all();
458  return cond ? -(max-min.array()).log().sum() : neg_inf<double>;
459 }
460 
461 // vsv
462 template <class XType
463  , class MinType
464  , class MaxType
465  , class = std::enable_if_t<
466  std::is_arithmetic_v<MinType>
467  > >
468 inline dist_value_t uniform_log_pdf(const Eigen::MatrixBase<XType>& x,
469  const MinType& min,
470  const Eigen::MatrixBase<MaxType>& max)
471 {
472  bool cond = (min < x.array()).all() && (x.array() < max.array()).all();
473  return cond ? -(max.array()-min).log().sum() : neg_inf<double>;
474 }
475 
476 // vvv
477 template <class XType
478  , class MinType
479  , class MaxType>
480 inline dist_value_t uniform_log_pdf(const Eigen::MatrixBase<XType>& x,
481  const Eigen::MatrixBase<MinType>& min,
482  const Eigen::MatrixBase<MaxType>& max)
483 {
484  bool cond = (min.array() < x.array()).all() && (x.array() < max.array()).all();
485  return cond ? -(max.array()-min.array()).log().sum() : neg_inf<double>;
486 }
487 
489 // Bernoulli Density
491 
499 // ss
500 template <class XType
501  , class PType
502  , class = std::enable_if_t<
503  std::is_arithmetic_v<XType> &&
504  std::is_arithmetic_v<PType>
505  > >
506 inline dist_value_t bernoulli_pdf(const XType& x,
507  const PType& p)
508 {
509  if (p <= 0) return (x == 0) + 0.;
510  else if (p >= 1) return (x == 1) + 0.;
511 
512  if (x == 1) return p;
513  else if (x == 0) return 1. - p;
514  else return 0.0;
515 }
516 
517 // vs
518 template <class XType
519  , class PType
520  , class = std::enable_if_t<
521  std::is_arithmetic_v<PType>
522  > >
523 inline dist_value_t bernoulli_pdf(const Eigen::MatrixBase<XType>& x,
524  const PType& p)
525 {
526  double pdf = 1.;
527  for (int i = 0; i < x.size(); ++i) {
528  pdf *= bernoulli_pdf(x(i), p);
529  }
530  return pdf;
531 }
532 
533 // vv
534 template <class XType
535  , class PType>
536 inline dist_value_t bernoulli_pdf(const Eigen::MatrixBase<XType>& x,
537  const Eigen::MatrixBase<PType>& p)
538 {
539  assert(x.size() == p.size());
540  double pdf = 1.;
541  for (int i = 0; i < x.size(); ++i) {
542  pdf *= bernoulli_pdf(x(i), p(i));
543  }
544  return pdf;
545 }
546 
547 // ss
548 template <class XType
549  , class PType
550  , class = std::enable_if_t<
551  std::is_arithmetic_v<XType> &&
552  std::is_arithmetic_v<PType>
553  > >
554 inline dist_value_t bernoulli_log_pdf(const XType& x,
555  const PType& p)
556 {
557  if (p <= 0) {
558  if (x == 0) return 0.;
559  else return neg_inf<PType>;
560  }
561  else if (p >= 1) {
562  if (x == 1) return 0.;
563  else return neg_inf<PType>;
564  }
565 
566  if (x == 1) return std::log(p);
567  else if (x == 0) return std::log(1. - p);
568  else return neg_inf<PType>;
569 }
570 
571 // vs
572 template <class XType
573  , class PType
574  , class = std::enable_if_t<
575  std::is_arithmetic_v<PType>
576  > >
577 inline dist_value_t bernoulli_log_pdf(const Eigen::MatrixBase<XType>& x,
578  const PType& p)
579 {
580  double logpdf = 0.;
581  for (int i = 0; i < x.size(); ++i) {
582  logpdf += bernoulli_log_pdf(x(i), p);
583  }
584  return logpdf;
585 }
586 
587 // vv
588 template <class XType
589  , class PType>
590 inline dist_value_t bernoulli_log_pdf(const Eigen::MatrixBase<XType>& x,
591  const Eigen::MatrixBase<PType>& p)
592 {
593  assert(x.size() == p.size());
594  double logpdf = 0.;
595  for (int i = 0; i < x.size(); ++i) {
596  logpdf += bernoulli_log_pdf(x(i), p(i));
597  }
598  return logpdf;
599 }
600 
602 // Wishart Density
603 // Note: drops unnecessary constant terms
605 template <class XType
606  , class VType
607  , class NType>
608 inline dist_value_t wishart_log_pdf(const Eigen::MatrixBase<XType>& x,
609  const Eigen::MatrixBase<VType>& v,
610  const NType& n)
611 {
612  Eigen::LLT<Eigen::MatrixXd> x_llt(x);
613  Eigen::LLT<Eigen::MatrixXd> v_llt(v);
614  auto log_det_x = std::log(x_llt.matrixL().determinant());
615  auto log_det_v = std::log(v_llt.matrixL().determinant());
616  auto tr = v_llt.solve(x).trace();
617  dist_value_t p = v.rows();
618  return (n - p - 1.) * log_det_x - 0.5 * tr - n * log_det_v;
619 }
620 
621 } // namespace math
622 } // namespace ppl
ppl::math::normal_pdf
dist_value_t normal_pdf(const XType &x, const MeanType &mean, const SigmaType &sigma)
Definition: density.hpp:40
ppl::math::bernoulli_pdf
dist_value_t bernoulli_pdf(const XType &x, const PType &p)
Definition: density.hpp:506
dist_expr_traits.hpp
ppl::math::uniform_log_pdf
dist_value_t uniform_log_pdf(const XType &x, const MinType &min, const MaxType &max)
Definition: density.hpp:421
ppl::math::dist_value_t
util::dist_value_t dist_value_t
Definition: density.hpp:16
ppl::math::bernoulli_log_pdf
dist_value_t bernoulli_log_pdf(const XType &x, const PType &p)
Definition: density.hpp:554
ppl::math::wishart_log_pdf
dist_value_t wishart_log_pdf(const Eigen::MatrixBase< XType > &x, const Eigen::MatrixBase< VType > &v, const NType &n)
Definition: density.hpp:608
ppl::math::LOG_SQRT_TWO_PI
constexpr double LOG_SQRT_TWO_PI
Definition: density.hpp:24
ppl::math::normal_log_pdf
dist_value_t normal_log_pdf(const XType &x, const MeanType &mean, const SigmaType &sigma)
Definition: density.hpp:149
ppl::math::SQRT_TWO_PI
constexpr double SQRT_TWO_PI
Definition: density.hpp:22
ppl::util::dist_value_t
double dist_value_t
Definition: dist_expr_traits.hpp:13
ppl::math::cauchy_log_pdf
dist_value_t cauchy_log_pdf(const XType &x, const LocType &loc, const ScaleType &scale)
Definition: density.hpp:261
ppl
Definition: bounded.hpp:11
ppl::math::uniform_pdf
dist_value_t uniform_pdf(const XType &x, const MinType &min, const MaxType &max)
Definition: density.hpp:347
math.hpp