feat: more tests

src/code/geo2d/area_poc.hpp

@@ -1,6 +1,7 @@
 #ifndef TIFALIBS_GEO2D_AREA_POC
 #define TIFALIBS_GEO2D_AREA_POC
 
+#include "../math/kahan.hpp"
 #include "circle.hpp"
 #include "polygon.hpp"
 #include "sarea_ct.hpp"
@@ -9,11 +10,11 @@ namespace tifa_libs::geo {
 
 template <class FP>
 CEXP FP area_PoC(polygon<FP> CR po, circle<FP> CR c) {
-  FP ans{};
+  math::kahan<FP> ans{};
   u32 sz = (u32)po.vs.size();
   if (sz < 3) return ans;
   flt_ (u32, i, 0, sz) ans += sarea_CT(c, po[i], po[po.next(i)]);
-  return abs(ans);
+  return abs<FP>(ans);
 }
 
 }  // namespace tifa_libs::geo

src/code/geo2d/coverage_rect_with_min_area.hpp

@@ -2,26 +2,26 @@
 #define TIFALIBS_GEO2D_COVERAGE_RECT_WITH_MIN_AREA
 
 #include "cvh.hpp"
+#include "dot.hpp"
 
 namespace tifa_libs::geo {
 
 // Coverage rectangle with min area
 template <class FP>
-CEXP polygon<FP> coverage_rect_with_min_area(cT_(cvh<FP>) ch) {
+CEXP polygon<FP> coverage_rect_with_min_area(cvh<FP> CR ch) {
   u32 n = (u32)ch.vs.size();
   if (n == 0) return ch;
   if (n == 1) return polygon{vec<point<FP>>{ch[0], ch[0], ch[0], ch[0]}};
   if (n == 2) return polygon{vec<point<FP>>{ch[0], ch[0], ch[1], ch[1]}};
   FP ans = std::numeric_limits<FP>::max();
   u32 r = 1, p = 1, q = 1;
-  u32 ans_i, ans_r, ans_p, ans_q;
+  u32 ans_i = 0, ans_r = 0, ans_p = 0, ans_q = 0;
   flt_ (u32, i, 0, n) {
     while (!is_neg(cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[ch.next(r)]) - cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[r]))) r = ch.next(r);
     while (!is_neg(dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[ch.next(p)]) - dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[p]))) p = ch.next(p);
     if (i == 0) q = p;
     while (!is_pos(dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[ch.next(q)]) - dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[q]))) q = ch.next(q);
-    FP tmp = cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[r]) * (dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[p]) - dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[q])) / (ch.vs[i] - ch.vs[ch.next(i)]).norm2();
-    if (ans > tmp) ans = tmp, ans_i = i, ans_r = r, ans_p = p, ans_q = q;
+    if (FP _ = cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[r]) * (dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[p]) - dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[q])) / (ch.vs[i] - ch.vs[ch.next(i)]).norm2(); ans > _) ans = _, ans_i = i, ans_r = r, ans_p = p, ans_q = q;
   }
   point dir = line{ch[ans_i], ch[ch.next(ans_i)]}.direction(), vdir = rot90(dir);
   line li{ch[ans_i], ch[ans_i] + dir}, lp{ch[ans_p], ch[ans_p] + vdir}, lr{ch[ans_r], ch[ans_r] + dir}, lq{ch[ans_q], ch[ans_q] + vdir};

src/code/geo2d/coverage_rect_with_min_circum.hpp

@@ -14,7 +14,7 @@ CEXP polygon<FP> coverage_rect_with_min_circum(cT_(cvh<FP>) ch) {
   if (n == 2) return polygon{vec<point<FP>>{ch[0], ch[0], ch[1], ch[1]}};
   FP ans = std::numeric_limits<FP>::max();
   u32 r = 1, p = 1, q = 1;
-  u32 ans_i = 0, ans_r, ans_p, ans_q;
+  u32 ans_i = 0, ans_r = 0, ans_p = 0, ans_q = 0;
   flt_ (u32, i, 0, n) {
     while (!is_neg(cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[ch.next(r)]) - cross(ch.vs[i], ch.vs[ch.next(i)], ch.vs[r]))) r = ch.next(r);
     while (!is_neg(dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[ch.next(p)]) - dot(ch.vs[i], ch.vs[ch.next(i)], ch.vs[p]))) p = ch.next(p);

src/code/geo2d/dist_ps.hpp

@@ -7,14 +7,18 @@
 
 namespace tifa_libs::geo {
 
-// min dist_PP from a point to another point which belongs to a segment
+// the point which is min dist_PP from a point to another point which belongs to a segment
 template <class FP>
-CEXP FP dist_PS(point<FP> CR p, line<FP> CR s) {
-  if (s.l == s.r) return dist_PP(s.l, p);
-  if (point h = proj(s, p); is_in_middle(s.l, h, s.r)) return dist_PP(p, h);
-  return min(dist_PP(s.l, p), dist_PP(s.r, p));
+CEXP point<FP> dist_PS_P(point<FP> CR p, line<FP> CR s) {
+  if (s.l == s.r) return s.l;
+  if (point<FP> h = proj(s, p); is_in_middle(s.l, h, s.r)) return h;
+  return is_lt(dist_PP(s.l, p), dist_PP(s.r, p)) ? s.l : s.r;
 }
 
+// min dist_PP from a point to another point which belongs to a segment
+template <class FP>
+CEXP FP dist_PS(point<FP> CR p, line<FP> CR s) { return dist_PP(dist_PS_P(p, s), p); }
+
 }  // namespace tifa_libs::geo
 
 #endif

src/code/geo2d/fermatp.hpp

@@ -1,6 +1,7 @@
 #ifndef TIFALIBS_GEO2D_FERMATP
 #define TIFALIBS_GEO2D_FERMATP
 
+#include "../math/kahan.hpp"
 #include "../opt/heuristic_sa.hpp"
 #include "make_p_polar.hpp"
 #include "massp.hpp"
@@ -12,8 +13,8 @@ template <class FP>
 point<FP> fermatp(polygon<FP> CR po, FP begin = 1e10, FP end = eps_v<FP>, FP delta = .999) {
   static rand::Gen<std::uniform_real_distribution<FP>> gen_angle(0, 2 * std::numbers::pi_v<FP>);
   auto gen = [](point<FP> CR pre, FP t) { return pre + make_P_polar(t, gen_angle()); };
-  auto fitness = [&po](point<FP> CR p) {
-    FP dis = 0;
+  auto fitness = [&po](point<FP> CR p) -> FP {
+    math::kahan<FP> dis = 0;
     for (u32 i = 0; i < po.vs.size(); ++i) dis += dist_PP(po.vs[i], p);
     return dis;
   };

src/code/geo2d/make_c_rcc_ex.hpp

@@ -11,8 +11,8 @@ namespace tifa_libs::geo {
 template <class FP>
 CEXP std::optional<ptt<circle<FP>>> make_C_rCC_ex(FP r, circle<FP> CR c1, circle<FP> CR c2) {
   if (relation_CC(c1, c2) == RELCC::lyingin_cc) return {};
-  if (auto ps = ins_CC({c1.o, c1.r + r}, {c2.o, c2.r + r}); !ps.has_value()) return {};
-  else return {{ps.value().first, r}, {ps.value().second, r}};
+  if (auto ps = ins_CC<FP>({c1.o, c1.r + r}, {c2.o, c2.r + r}); !ps.has_value()) return {};
+  else return ptt<circle<FP>>{{ps.value().first, r}, {ps.value().second, r}};
 }
 
 }  // namespace tifa_libs::geo

src/code/geo2d/make_c_rll.hpp

@@ -15,7 +15,7 @@ CEXP std::optional<pt4<circle<FP>>> make_C_rLL(FP r, line<FP> CR l1, line<FP> CR
   dir2 *= r / dir2.norm();
   point<FP> dirl1 = rot90(dir1), dirr1 = rot270(dir1), dirl2 = rot90(dir2), dirr2 = rot270(dir2);
   line<FP> u1 = {l1.l + dirl1, l1.r + dirl1}, u2 = {l1.l + dirr1, l1.r + dirr1}, v1 = {l2.l + dirl2, l2.r + dirl2}, v2 = {l2.l + dirr2, l2.r + dirr2};
-  return {{ins_LL(u1, v1), r}, {ins_LL(u1, v2), r}, {ins_LL(u2, v1), r}, {ins_LL(u2, v2), r}};
+  return pt4<circle<FP>>{{ins_LL(u1, v1), r}, {ins_LL(u1, v2), r}, {ins_LL(u2, v1), r}, {ins_LL(u2, v2), r}};
 }
 
 }  // namespace tifa_libs::geo

src/code/geo2d/make_c_rpl.hpp

@@ -1,6 +1,7 @@
 #ifndef TIFALIBS_GEO2D_MAKE_C_RPL
 #define TIFALIBS_GEO2D_MAKE_C_RPL
 
+#include "dist_pl.hpp"
 #include "ins_cl.hpp"
 
 namespace tifa_libs::geo {
@@ -14,10 +15,10 @@ CEXP std::optional<ptt<circle<FP>>> make_C_rPL(FP r, point<FP> CR p, line<FP> CR
   point dir = l.direction();
   dir *= r / dir.norm();
   point dirl = rot90(dir), dirr = rot270(dir);
-  if (is_zero(dis)) return {{p + dirl, r}, {p + dirr, r}};
+  if (is_zero(dis)) return ptt<circle<FP>>{{p + dirl, r}, {p + dirr, r}};
   circle c{p, r};
   if (auto ps = ins_CL(c, {l.l + dirl, l.r + dirl}); !ps.has_value() && !(ps = ins_CL(c, {l.l + dirr, l.r + dirr})).has_value()) return {};
-  else return {{ps.value().first, r}, {ps.value().second, r}};
+  else return ptt<circle<FP>>{{ps.value().first, r}, {ps.value().second, r}};
 }
 
 }  // namespace tifa_libs::geo

src/code/geo2d/massp.hpp

@@ -8,17 +8,17 @@ namespace tifa_libs::geo {
 template <class FP>
 CEXP point<FP> massp(polygon<FP> CR po) {
   point<FP> ret{};
-  FP area{};
+  math::kahan<FP> area{};
   u32 n = (u32)po.vs.size();
   if (n == 0) return ret;
   if (n == 1) return po.vs[0];
   for (u32 i = 1; i < n - 1; ++i) {
-    FP tmp = cross(po.vs[0], po.vs[i], po.vs[i + 1]);
-    if (is_zero(tmp)) continue;
-    area += tmp;
-    ret += (po.vs[0] + po.vs[i] + po.vs[i + 1]) * (tmp / 3);
+    FP _ = cross(po.vs[0], po.vs[i], po.vs[i + 1]);
+    if (is_zero(_)) continue;
+    area += _;
+    ret += (po.vs[0] + po.vs[i] + po.vs[i + 1]) * (_ / 3);
   }
-  if (!is_zero(area)) ret /= area;
+  if (!is_zero<FP>(area)) ret /= area;
   return ret;
 }
 

src/code/geo2d/polygon.hpp

@@ -2,6 +2,7 @@
 #define TIFALIBS_GEO2D_POLYGON
 
 #include "../edh/discretization.hpp"
+#include "../math/kahan.hpp"
 #include "cross.hpp"
 #include "dist_pp.hpp"
 
@@ -42,13 +43,13 @@ struct polygon {
   CEXP u32 next(u32 idx) const { return idx + 1 == (u32)vs.size() ? 0 : idx + 1; }
 
   CEXP FP circum() const {
-    FP ret = dist_PP(vs.back(), vs.front());
+    math::kahan<FP> ret = dist_PP(vs.back(), vs.front());
     for (u32 i = 0; i < (u32)vs.size() - 1; ++i) ret += dist_PP(vs[i], vs[i + 1]);
     return ret;
   }
   CEXP FP area() const {
     if (vs.size() < 3) return 0;
-    FP ret = vs.back() ^ vs.front();
+    math::kahan<FP> ret = vs.back() ^ vs.front();
     for (u32 i = 0; i < (u32)vs.size() - 1; ++i) ret += vs[i] ^ vs[i + 1];
     return ret / 2;
   }

src/code/geo2d/rel_pop.hpp

@@ -2,6 +2,7 @@
 #define TIFALIBS_GEO2D_REL_POP
 
 #include "cvh.hpp"
+#include "dot.hpp"
 #include "is_on_s.hpp"
 
 namespace tifa_libs::geo {

src/code/math/kahan.hpp

@@ -5,16 +5,21 @@
 
 namespace tifa_libs::math {
 
-template <class FP>
-CEXP FP kahan(vec<FP> CR vec) {
-  FP sum = 0, c = 0;
-  for (auto x : vec) {
-    FP y = x - c, t = sum + y;
-    c = (t - sum) - y;
-    sum = t;
+template <std::floating_point FP>
+class kahan {
+  FP sum, c;
+
+ public:
+  CEXP kahan(FP val = 0) : sum(val), c(0) {}
+
+  CEXP kahan& operator+=(FP x) {
+    FP y = x - c;
+    volatile FP t = sum + y, z = t - sum;
+    c = z - y, sum = t;
+    return *this;
   }
-  return sum;
-}
+  CEXP operator FP() const { return sum; }
+};
 
 }  // namespace tifa_libs::math
 

src/code/util/traits.hpp

@@ -12,7 +12,7 @@ concept iterable_c = requires(T v) {
 };
 
 template <class T>
-concept container_c = iterable_c<T> && !std::derived_from<T, std::basic_string<typename T::value_type>>;
+concept container_c = iterable_c<T> && !std::derived_from<T, std::basic_string<typename T::value_type>> && !std::derived_from<T, std::basic_string_view<typename T::value_type>>;
 
 template <class T>
 CEXP bool is_char_v = std::is_same_v<T, char> || std::is_same_v<T, signed char> || std::is_same_v<T, unsigned char>;

src/data/uva/10173/1.in

@@ -0,0 +1,10 @@
+3
+-3.000 5.000
+7.000 9.000
+17.000 5.000
+4
+10.000 10.000
+10.000 20.000
+20.000 20.000
+20.000 10.000
+0

src/data/uva/10173/1.out

@@ -0,0 +1,2 @@
+80.0000
+100.0000