added 118 prak

16 files changed, 485 insertions, 14 deletions

diff --git a/106/J.hpp b/106/J.hpp
new file mode 100644
index 0000000..a77b1d3
--- /dev/null
+++ b/106/J.hpp
@@ -0,0 +1,47 @@
+#pragma once
+#include "prakcommon.hpp"
+#include <array>
+
+inline constexpr const f64 PI = 3.141592657; // ;)
+
+inline constexpr const f64 rho_stl = 7800; // g/dm^3 (dm = decimeter)
+// delta d1 d2 d3 d4_1_ d1_ d2_ d3_ delta_
+inline constexpr const std::array d_stl = 
+   {0.0102, 0.0053, 0.0072, 0.0103, 0.0162, 0.0053, 0.0072, 0.0103, 0.0102};
+// lambd l1 l2 l3 l4_1_ l1_ l2_ l3_ lambd_
+inline constexpr const std::array l_stl = 
+   {0.0024, 0.0334, 0.03515, 0.071, 0.0053, 0.033, 0.0347, 0.07345, 0.0022};
+
+inline constexpr const f64 rho_alm = 2800; // g/dm^3 
+// d4_1 D1 d4_2_ D2
+inline constexpr const std::array d_alm =
+   {0.01575, 0.0806, 0.01675, 0.00555};
+// l4_1 L_ l4_2_ L
+inline constexpr const std::array l_alm =
+   {0.00635, 0.0111, 0.0053, 0.0205};
+
+
+constexpr inline f64 getJ(f64 m) {
+	f64 V1 = 0, V2 = 0;
+	f64 m1 = 0, m2 = 0;
+	f64 dm;
+	f64 J1 = 0, J2 = 0;
+	f64 D = d_alm.back() + d_alm[1];
+	for (size_t i = 0; i < d_stl.size(); ++i)
+		V1 += d_stl[i] * d_stl[i] * l_stl[i];
+	for (size_t i = 0; i < d_alm.size() - 1; ++i)
+		V2 += d_alm[i] * d_alm[i] * l_alm[i];
+	
+	m1 = PI * V1 * rho_stl / 4;
+	m2 = PI * rho_alm * (V2 + (D*D - d_alm[1]*d_alm[1]) * l_alm.back()) / 4;
+	dm = m - m1 - m2;
+	for (size_t i = 0; i < d_stl.size(); ++i) 
+		J1 += d_stl[i] * d_stl[i];
+	J1 *= m1/8;
+
+	for (size_t i = 0; i < d_alm.size(); ++i) 
+		J2 += d_alm[i] * d_alm[i];
+	J2 += D*D + d_alm[1]*d_alm[1];
+	J2 *= m2/8;
+	return J1+J2;
+}
diff --git a/106/Makefile b/106/Makefile
index c3dda94..610abc1 100644
--- a/106/Makefile
+++ b/106/Makefile
@@ -3,6 +3,11 @@ CFLAGS = -std=c++2c -mavx -Iinclude -ggdb
 
 .PHONY: all run_main clean gnuplot 
 
+all: gnuplot
+
+gnuplot:  run_main plots.gp
+	gnuplot plots2.gp
+
 run_main: main
 	./main
 
diff --git a/106/data b/106/data
new file mode 100644
index 0000000..9d985d2
--- /dev/null
+++ b/106/data
@@ -0,0 +1,8 @@
+N1	N2	t1	t2	t3	h1	h2	h3	2T1	2T2	2T3
+0	1	2.616	2.491	2.601	31.8	31.9	31.1	3.76	3.78 	3.71
+1	2	4.722	4.698	4.715	28.2	27.9	27.8	3.77	3.78	3.73
+2	3	4.511	4.488	4.467	25.1	24.3	25.0	3.73 	3.87	3.73
+3	4	4.301	4.284	4.216	22.7	21.9	22.1	3.76 	3.75 	3.75
+4	5	4.030	3.972	4.027	19.7	19.6	20.1	3.74	3.75 	3.75
+5	6	3.819	3.793	3.796	18.0	18.0	18.3	3.73 	3.74 	3.76
+6	7	3.669	3.619	3.629	15.9	16.3	16.2	3.74 	3.74 	3.77
diff --git a/106/main.cpp b/106/main.cpp
index c83e5b5..e4c619e 100644
--- a/106/main.cpp
+++ b/106/main.cpp
@@ -1,13 +1,146 @@
 #include "include/praktable.hpp"
 
-using table = prak::table<double>;
+using table = prak::table<f64>;
+using vecarg = const std::vector<f64> &;
+using f64v = std::vector<f64>;
+using f64p = prak::pvalue<f64>;
 
-void ex1(void) {
-	table t;
-	t.print();
+table prepare(std::string f) {
+	return table(f)
+		.add_column("t", {})
+		.add_column("st", {})
+		.add_column("h", {})
+		.add_column("sh", {})
+		.add_column("2T", {})
+		.add_column("s2T", {})
+		.apply(prak::avg<f64>,		{"t1", "t2", "t3"},	"t")
+		.apply(prak::avg<f64>,		{"h1", "h2", "h3"},	"h")
+		.apply(prak::avg<f64>,		{"2T1", "2T2", "2T3"},	"2T")
+		.apply([] (vecarg a) -> f64 { return std::sqrt(std::pow(prak::stddev<f64>(a), 2) + 0.001*0.001); }, {"t1", "t2", "t3"},	"st")
+		.apply([] (vecarg a) -> f64 { return std::sqrt(std::pow(prak::stddev<f64>(a), 2) + 0.100*0.100); }, {"h1", "h2", "h3"},	"sh")
+		.apply([] (vecarg a) -> f64 { return std::sqrt(std::pow(prak::stddev<f64>(a), 2) + 0.050*0.050); }, {"2T1", "2T2", "2T3"},"s2T")
+		.delete_cols({"t1", "t2", "t3", "h1", "h2", "h3", "2T1", "2T2", "2T3"});
+}
+
+// ex1
+f64 a_ht(vecarg a) { return 8 * a[0]/100 / std::pow(a[1], 2); }
+f64 v_ht(vecarg a) { return 4 * a[0] / 100 / a[1]; }
+// ex2
+f64 tau_rvv(vecarg a) { return 2 * prak::PI * a[0] / (a[1] + a[2]); }
+f64 D2T_mvvr(vecarg a) { return a[0] * (a[1] + a[2]) * (a[1] + a[2]) / 4 / a[3]; }
+f64 Pcal(vecarg a) { return a[0] / a[1] + 1; }
+// ex3
+// [0] = m, [1] = r, [2] = g, [3] = a
+f64 j_exp(vecarg a) { return a[0] * a[1] * a[1] * (a[2] / a[3] - 1); }
+/* f64 getJ(f64 m) */ #include "J.hpp"
+// ex4
+f64 Wkin_mvjr(vecarg a) { return a[0] * a[1] * a[1] / 2 * (1 + a[2] / (a[0] * a[3] * a[3])); }
+f64 Wpot_mgh(vecarg a) { return a[0] * a[1] * a[2]; }
+
+table ex1(table &t) {
+	table result = table({"N", "a", "sa", "v", "sv", "K", "sK"}, t.rows, NAN);
+	result	.apply_function_n(
+			a_ht, 
+			{t.begin("h"), t.begin("t")},
+			{t.begin("sh"), t.begin("st")},
+			t.rows, "a", "sa"
+		).apply_function_n(
+			v_ht,
+			{t.begin("h"), t.begin("t")},
+			{t.begin("sh"), t.begin("st")},
+			t.rows, "v", "sv"
+		).apply_function_n(
+			[](vecarg a){ return a[0] / a[1]; },
+			{t.begin("h") + 1, t.begin("h")},
+			{t.begin("sh") + 1, t.begin("sh")},
+			t.rows - 1, "K", "sK");
+	result["a", 0] *= 0.25;
+	result["sa", 0] *= 0.25;
+	result["v", 0] *= 0.5;
+	result["sv", 0] *= 0.5;
+	for (size_t i = 0; i < result.rows; ++i)
+		result["N", i] = i+1;
+
+	std::cout << "Предупреждение: Колонки K, sK НАДО сдвинуть на 1 позицию вниз в тетради!\n" << result << std::endl;
+	std::cout << "среднее ускорение = " << result.col_avg("a") << std::endl;
+	result.write_plot("ex1_v.plot", "N", "v", "sv");
+	return result;
+}
+
+void ex2(table &t, table &e1) {
+	table result({"N", "tau", "stau", "D2Tmax", "sD2Tmax", "P_cal", "P_exp", "mg"}, e1.rows - 1, NAN);
+	e1.add_column("r", std::vector(e1.rows, 0.0036));
+	e1.add_column("sr", std::vector(e1.rows, 0.00005));
+	e1.add_column("m", std::vector(e1.rows, 0.385));
+	e1.add_column("sm", std::vector(e1.rows, 0.004));	
+	result.fill_column("mg", 3.78);
+	result	.apply_function_n(tau_rvv, 
+			{e1.begin("r"), e1.begin("v"), e1.begin("v")+1},
+			{e1.begin("sr"), e1.begin("sv"), e1.begin("sv")+1},
+			e1.rows - 1, "tau", "stau"
+		).apply_function_n(D2T_mvvr, 
+			{e1.begin("m"), e1.begin("v"), e1.begin("v")+1, e1.begin("r")},
+			{e1.begin("sm"), e1.begin("sv"), e1.begin("sv")+1, e1.begin("sr")},
+			e1.rows - 1, "D2Tmax", "sD2Tmax"
+		).apply(
+			Pcal, {"D2Tmax", "mg"}, "P_cal"
+		);
+	for (size_t i = 0; i < result.rows; ++i)
+		result["N", i] = i+1;
+	result.write_plot("ex2_t.plot", "N", "tau", "stau");
+	result.write_plot("ex2_T.plot", "N", "D2Tmax", "sD2Tmax");
+	std::cout << result << std::endl;
+}
+
+f64p ex3(table &ex1) {
+	f64p g = {9.81, 0.01};
+	f64p a = {ex1.col_avg("a"), ex1.col_stddev("a")};
+	f64p r = {ex1.col_avg("r"), ex1.col_avg("sr")};
+	f64p m = {0.385, 0.004};
+	std::vector args = {m.val, r.val, g.val, a.val};
+	std::vector sgms = {m.err, r.err, g.err, a.err};
+	f64p J_exp = {j_exp(args), prak::sigma<f64>(j_exp, args, sgms)};
+	f64 J_cal = getJ(m.val);
+	std::cout << "Experimental J = " << J_exp << "\nCalculated J = " << J_cal << std::endl;
+	return J_exp;
+}
+
+void ex4(table &t, table &e1, f64p &J) {
+	table res({"N", "W_kin", "sW_kin", "W_pot", "sW_pot", "sumW", "dW", "dW(%)"}, e1.rows, NAN);
+	e1.add_column("J", std::vector(e1.rows, J.val));
+	e1.add_column("g", std::vector(e1.rows, 9.81));
+	e1.add_column("sJ", std::vector(e1.rows, J.err));
+	e1.add_column("sg", std::vector(e1.rows, 0.01));
+	t.multiply_column("h", 0.01).multiply_column("sh", 0.01);
+
+	for (size_t i = 0; i < res.rows; ++i)
+		res["N", i] = i+1;
+
+
+	res.apply_function_n(Wkin_mvjr, 
+			{e1.begin("m"), e1.begin("v"), e1.begin("J"), e1.begin("r")},
+			{e1.begin("sm"), e1.begin("sv"), e1.begin("sJ"), e1.begin("sr")},
+			e1.rows, "W_kin", "sW_kin"
+		).apply_function_n(Wpot_mgh, 
+			{e1.begin("m"), e1.begin("g"), t.begin("h")},
+			{e1.begin("sm"), e1.begin("sg"), t.begin("sh")},
+			e1.rows, "W_pot", "sW_pot"
+		).apply([](vecarg a){return std::abs(a[0]-a[1]);}, {"W_kin", "W_pot"}, "dW")
+		.apply([](vecarg a){return std::abs(a[0]-a[1])/a[0]*100;}, {"W_kin", "W_pot"}, "dW(%)")
+		.apply([](vecarg a){return (a[0]+a[1])/2;}, {"W_kin", "W_pot"}, "sumW");
+	res.write_plot("ex4_WK.plot", "N", "W_kin", "sW_kin");
+	res.write_plot("ex4_WP.plot", "N", "W_pot", "sW_pot");
+	res.write_plot("ex4_Wfull.plot", "N", "sumW", std::nullopt);
+
+	std::cout << res;
 }
 
 int main() {
-	ex1();	
+	table t = prepare("data");
+	t.print();
+	table e1 = ex1(t);
+	ex2(t, e1);
+	f64p J = ex3(e1);
+	ex4(t, e1, J);
 	return 0;
 }
diff --git a/106/plots.gp b/106/plots.gp
new file mode 100644
index 0000000..aecb7d3
--- /dev/null
+++ b/106/plots.gp
@@ -0,0 +1,39 @@
+# Set common properties
+set terminal pngcairo size 600,800 enhanced
+set style line 1 lw 2 pt 7 lc rgb "blue"  # Line style 1
+set style line 2 lw 2 pt 5 lc rgb "red"   # Line style 2
+set style line 3 lw 2 pt 9 lc rgb "green" # Line style 3
+
+# Plot 1: ex1_v.plot with yerrorbars and lines
+set output "ex1_velocity.png"
+set xlabel "N, количество"
+set ylabel "Скорость (м/с)"
+set key top left
+plot "ex1_v.plot" using 1:2:3 with yerrorbars ls 1 title "Скорость", \
+     "ex1_v.plot" using 1:2 with lines ls 1 notitle
+
+# Plot 2: ex2_t.plot and ex2_T.plot with two Y axes and lines
+set output "ex2_time_force.png"
+set xlabel "N, количество"
+set ylabel "Время (с)"
+set y2label "Сила (Н)"
+set y2tics
+set ytics nomirror
+set key top left
+plot "ex2_t.plot" using 1:2:3 with yerrorbars ls 1 title "Время", \
+     "ex2_t.plot" using 1:2 with lines ls 1 notitle, \
+     "ex2_T.plot" using 1:2:3 axes x1y2 with yerrorbars ls 2 title "Сила", \
+     "ex2_T.plot" using 1:2 axes x1y2 with lines ls 2 notitle
+
+# Plot 3: ex4_WK.plot, ex4_WP.plot (yerrorbars and lines) and ex4_Wfull.plot (lines only)
+set output "ex4_energy.png"
+set xlabel "N, количество"
+set ylabel "Энергия (Дж)"
+set key top left
+plot "ex4_WK.plot" using 1:2:3 with yerrorbars ls 1 title "Энергия Wk", \
+     "ex4_WK.plot" using 1:2 with lines ls 1 notitle, \
+     "ex4_WP.plot" using 1:2:3 with yerrorbars ls 2 title "Энергия Wp", \
+     "ex4_WP.plot" using 1:2 with lines ls 2 notitle, \
+     "ex4_Wfull.plot" using 1:2 with lines ls 3 title "Полная энергия"
+
+
diff --git a/106/plots2.gp b/106/plots2.gp
new file mode 100644
index 0000000..5dbdc9e
--- /dev/null
+++ b/106/plots2.gp
@@ -0,0 +1,37 @@
+
+# Common settings
+set terminal pngcairo size 1000, 800 enhanced
+set grid
+set datafile separator whitespace
+set key top left
+set xlabel "N, количество" font ",14"
+
+# 1. Plot ex1_v.plot
+set output "plot1_velocity.png"
+set ylabel "Скорость (м/с)" font ",14"
+plot \
+    "ex1_v.plot" using 1:2:3 with yerrorbars linecolor rgb "blue" pointtype 7 title "Скорость", \
+    "ex1_v.plot" using 1:2 smooth mcspline with lines linetype 2 linecolor rgb "blue" title ""
+
+# 2. Plot ex2_t.plot and ex2_T.plot
+set output "plot2_time_force.png"
+set ylabel "t (с)" font ",14"
+set y2label "Сила (Н)" font ",14"
+set y2tics
+set ytics nomirror
+plot \
+    "ex2_t.plot" using 1:2:3 with yerrorbars linecolor rgb "red" pointtype 5 title "Время", \
+    "ex2_t.plot" using 1:2 smooth mcspline with lines linetype 2 linecolor rgb "red" title "", \
+    "ex2_T.plot" using 1:2:3 axes x1y2 with yerrorbars linecolor rgb "green" pointtype 7 title "Сила", \
+    "ex2_T.plot" using 1:2 smooth mcspline axes x1y2 with lines linetype 2 linecolor rgb "green" title ""
+
+# 3. Plot ex4_WK.plot, ex4_WP.plot, ex4_Wfull.plot
+set output "plot3_energy.png"
+set ylabel "Энергия (единицы энергии)" font ",14"
+plot \
+    "ex4_WK.plot" using 1:2:3 with yerrorbars linecolor rgb "magenta" pointtype 7 title "Wk", \
+    "ex4_WK.plot" using 1:2 smooth mcspline with lines linetype 2 linecolor rgb "magenta" title "", \
+    "ex4_WP.plot" using 1:2:3 with yerrorbars linecolor rgb "cyan" pointtype 5 title "Wp", \
+    "ex4_WP.plot" using 1:2 smooth mcspline with lines linetype 2 linecolor rgb "cyan" title "", \
+    "ex4_Wfull.plot" using 1:2 with lines linecolor rgb "black" linetype 3 linewidth 1 title "Wfull"
+
diff --git a/106/plots3.gp b/106/plots3.gp
new file mode 100644
index 0000000..c763292
--- /dev/null
+++ b/106/plots3.gp
@@ -0,0 +1,51 @@
+# Common settings
+set terminal pngcairo size 1000, 800 enhanced
+set grid
+set datafile separator whitespace
+set xlabel "N, ед." font ",14"
+
+# 1. Plot ex1_v.plot
+set output "plot1_velocity.png"
+set ylabel "v (м/с)" font ",14"
+f1(x) = a1 * x + b1
+fit f1(x) "ex1_v.plot" using 1:2 via a1, b1
+plot \
+    "ex1_v.plot" using 1:2:3 with yerrorbars linecolor rgb "blue" pointtype 7 title "v", \
+    f1(x) with lines linetype 2 linecolor rgb "blue" notitle
+
+# 2. Plot ex2_t.plot and ex2_T.plot
+set output "plot2_time_force.png"
+set ylabel "t (с)" font ",14"
+set y2label "F (Н)" font ",14"
+set y2tics
+set ytics nomirror
+f2(x) = a2 * x + b2
+fit f2(x) "ex2_t.plot" using 1:2 via a2, b2
+f3(x) = a3 * x + b3
+fit f3(x) "ex2_T.plot" using 1:2 via a3, b3
+plot \
+    "ex2_t.plot" using 1:2:3 with yerrorbars linecolor rgb "red" pointtype 5 title "t, c", \
+    f2(x) with lines linetype 2 linecolor rgb "red" notitle,\
+    "ex2_T.plot" using 1:2:3 axes x1y2 with yerrorbars linecolor rgb "green" pointtype 7 title "F, Н", \
+    f3(x) axes x1y2 with lines linetype 2 linecolor rgb "green" notitle
+
+unset y2tics
+unset y2label
+# 3. Plot ex4_WK.plot, ex4_WP.plot, ex4_Wfull.plot
+set output "plot3_energy.png"
+set ylabel "Энергия (единицы энергии)" font ",14"
+f4(x) = a4 * x + b4
+fit f4(x) "ex4_WK.plot" using 1:2 via a4, b4
+f5(x) = a5 * x + b5
+fit f5(x) "ex4_WP.plot" using 1:2 via a5, b5
+f6(x) = a6 * x + b6
+fit f6(x) "ex4_Wfull.plot" using 1:2 via a6, b6
+plot \
+    "ex4_WK.plot" using 1:2:3 with yerrorbars lc rgb "red" pointtype 5 title "W_k", \
+    f4(x) with lines linetype 2 linecolor rgb "green" notitle, \
+    "ex4_WP.plot" using 1:2:3 with yerrorbars lc rgb "green" pointtype 5 title "W_p", \
+    f5(x) with lines linetype 2 linecolor rgb "blue" notitle, \
+    "ex4_Wfull.plot" using 1:2 pt 5 linecolor rgb "blue" title "W_{full}", \
+    f6(x) with lines linetype 2 linecolor rgb "red" notitle
+
+
diff --git a/118/Makefile b/118/Makefile
new file mode 100644
index 0000000..eb309ff
--- /dev/null
+++ b/118/Makefile
@@ -0,0 +1,18 @@
+
+CFLAGS = -std=c++2c -mavx -Iinclude -ggdb
+
+.PHONY: all run_main clean gnuplot 
+
+run: gnuplot
+
+gnuplot: plots.gp run_main
+	gnuplot $<
+
+run_main: main
+	./main
+
+main: main.cpp include/*
+	$(CXX) -o $@ $< $(CFLAGS)
+
+clean: 
+	rm -fr main *.png *.plot
diff --git a/118/README.md b/118/README.md
new file mode 100644
index 0000000..96fd716
--- /dev/null
+++ b/118/README.md
@@ -0,0 +1,5 @@
+<!-- Шаблон для прака -->
+<!-- файлы, заканчивающиеся на .plot считаются генерируемыми и удаляются через make clean -->
+# Обработка <> прака
+
+
diff --git a/118/compile_flags.txt b/118/compile_flags.txt
new file mode 100644
index 0000000..34ae930
--- /dev/null
+++ b/118/compile_flags.txt
@@ -0,0 +1,4 @@
+-Iinclude
+-std=c++2c
+-mavx2
+
diff --git a/118/include b/118/include
new file mode 120000
index 0000000..2225752
--- /dev/null
+++ b/118/include
@@ -0,0 +1 @@
+../libprakpp/include/
\ No newline at end of file
diff --git a/118/main.cpp b/118/main.cpp
new file mode 100644
index 0000000..5e2c0a2
--- /dev/null
+++ b/118/main.cpp
@@ -0,0 +1,13 @@
+#include <iostream>
+
+#include "include/prakmath.hpp"
+
+int main() {
+	prak::vector<f64> dx = {0.003, 0.005, 0.013, 0.007, 0.015, 0.020, 0.022, 0.028, 0.030, 0.032, 0.036, 0.040, 0.042, 0.046, 0.047, 0.065, 0.052};
+	prak::vector<f64> ms = {0.010, 0.020, 0.050, 0.030, 0.060, 0.070, 0.080, 0.100, 0.110, 0.120, 0.130, 0.150, 0.160, 0.170, 0.180, 0.243, 0.193};
+	prak::vector<f64> ss = {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};
+	prak::pvalue<f64> A, B;
+	prak::least_squares_linear<f64>(ms, dx, ss, A, B);
+	std::cout << A << '\n' << B << std::endl;
+	return 0;
+}
diff --git a/118/plots.gp b/118/plots.gp
new file mode 100644
index 0000000..702408a
--- /dev/null
+++ b/118/plots.gp
@@ -0,0 +1,14 @@
+set term pngcairo size 1000, 800
+set tmargin at screen 0.95
+
+f1(x) = a1*x+b1
+fit f1(x) '.plot' using 1:2:3 yerr via a1, b1
+
+set output ''
+set label "" at graph 0.5, graph 1.025 center
+set xlabel ""
+set ylabel ""
+
+plot 	'.plot' using 1:2:3 with yerrorbars notitle lc 0 pt 1 lw 2,	\
+	f1(x) title "" lc rgb "red",				\
+
diff --git a/libprakpp/include/prakcommon.hpp b/libprakpp/include/prakcommon.hpp
index c58040e..3035ff8 100644
--- a/libprakpp/include/prakcommon.hpp
+++ b/libprakpp/include/prakcommon.hpp
@@ -23,7 +23,6 @@ typedef long double	f128;
 #else
 typedef long double	_f64;
 #endif
-
 #if defined(_MSC_VER) || !defined(__cpp_multidimensional_subscript) || __cplusplus < 202110L
 #warning "can not use multidimentional subscript operator: falling back to `operator()`"
 #undef MDSUBSCRIPT
diff --git a/libprakpp/include/prakmath.hpp b/libprakpp/include/prakmath.hpp
index a5eb0b8..f8024c9 100644
--- a/libprakpp/include/prakmath.hpp
+++ b/libprakpp/include/prakmath.hpp
@@ -2,6 +2,7 @@
 
 #include <cmath>
 #include <concepts>
+#include <iomanip>
 #include <numeric>
 
 #ifdef __x86_64__
diff --git a/libprakpp/include/praktable.hpp b/libprakpp/include/praktable.hpp
index 70f6d49..c2778e4 100644
--- a/libprakpp/include/praktable.hpp
+++ b/libprakpp/include/praktable.hpp
@@ -9,6 +9,7 @@
 #include <fstream>
 #include <iomanip>
 #include <ostream>
+#include <set>
 #include <string>
 #include <vector>
 #include <cmath>
@@ -62,8 +63,10 @@ class table {
 	using function_type = function_t<dtype>;
 	using stringvec = std::vector<std::string>;
 
-	size_t index(const std::string &str) const {
-		return std::distance(names.cbegin(), std::find(names.cbegin(), names.cend(), str));
+	size_t index(const std::string &str) const noexcept(false) {
+		auto ret = std::distance(names.cbegin(), std::find(names.cbegin(), names.cend(), str));
+		if (ret == columns) throw std::out_of_range("Column " + str + " does not exist");
+		return ret;
 	}
 
 	FILE* open_gnuplot() {
@@ -95,6 +98,7 @@ public:
 		}
 		iterator &operator++() { data_index += columns; return *this; }
 		iterator operator++(int) { iterator ret = *this; ++(*this); return ret; }
+		iterator &operator+(int x) { data_index += columns * x; return *this; }
 		bool operator==(iterator other) { return data_index == other.data_index && parent == other.parent && col_index == other.col_index; }
 		bool operator!=(iterator other) { return data_index != other.data_index || parent != other.parent || col_index != other.col_index; }
 		value_type &operator*() { return parent->data[data_index]; };
@@ -154,7 +158,6 @@ public:
 
 	dtype & SUBSCR_OPRTR (const std::string &column, size_t row) noexcept(false) {
 		size_t i = index(column);
-		if (i == columns) throw std::out_of_range("Column " + column + " does not exist");
 		return data.at(names.size() * row + index(column));
 	}
 
@@ -211,12 +214,47 @@ public:
 			if (result.has_value()) data[columns * i + result_index] = function(v);
 			else (void)function(v);
 		}
-		/* print(std::cerr); */
 		return *this;
 	}
-	
+
+	table &apply_n(function_type function, std::vector<iterator> cols, size_t n, std::optional<std::string> result) {
+		size_t result_index = result.has_value() ? index(*result) : 0;
+		for (size_t i = 0; i < n; ++i) {
+			std::vector<dtype> v(cols.size());
+			for (size_t j = 0; j < cols.size(); ++j) 
+				v[j] = *cols[j]++;
+			if (result.has_value()) data[columns * i + result_index] = function(v);
+			else (void)function(v);
+		}
+		return *this;
+	}
+
+	table &apply_function_n(
+		function_type function, 
+		std::vector<iterator> cols, 
+		std::vector<iterator> sigma_cols, 
+		size_t n, 
+		const std::string &resval, 
+		const std::string &ressigma) 
+	{
+		if (cols.size() != sigma_cols.size()) 
+			throw dimension_error("cols.size() is not equal to sigma_cols.size()");
+		size_t	val_index = index(resval),
+			sgm_index = index(ressigma);
+		for (size_t i = 0; i < n; ++i) {
+			std::vector<dtype> v(cols.size());
+			std::vector<dtype> s(cols.size());
+			for (size_t j = 0; j < cols.size(); ++j) {
+				v[j] = *cols[j]++;
+				s[j] = *sigma_cols[j]++;
+			}
+			data[columns * i + val_index] = function(v);
+			data[columns * i + sgm_index] = sigma(function, v, s);
+		}
+		return *this;
+	}
 	/// adds a column with name `name` and data `column_data`
-	void add_column(std::string name, std::vector<dtype> column_data) {
+	table &add_column(std::string name, std::vector<dtype> column_data) {
 		if (column_data.size() == 0) column_data = std::vector<dtype>(rows, dtype{});
 		std::vector<dtype> data_new(rows * (++columns));
 
@@ -227,9 +265,53 @@ public:
 		}
 		data = std::move(data_new);
 		names.push_back(name);
+		return *this;
 	}
 
-	/// Appends a column to the table. if name is set, appends it to `opt_rownames`
+	/// Deletes a column from a table.
+	table &delete_col(const std::string &colname) {
+		std::vector<dtype> data_new(rows * (--columns));
+		size_t idx = index(colname);
+		for (size_t column = 0; column < names.size(); ++column) {
+			if (column == idx) continue;
+			size_t _col = column - (column > idx);
+			for (size_t row = 0; row < rows; ++row)
+				data_new[row * columns + _col] = data[row * names.size() + column];
+		}
+		data = std::move(data_new);
+		names.erase(names.begin() + idx);
+		return *this;
+	}
+
+	/// Deletes several columns 
+	table &delete_cols(const stringvec &cols) noexcept(false) {
+		size_t columns_new = columns - cols.size();
+		std::vector<dtype> data_new(rows * columns_new);
+		std::set<size_t> idxs;
+		for (const std::string &col : cols)
+			idxs.insert(index(col));
+
+		size_t skipped = 0;
+		for (size_t column = 0; column < columns; ++column) {
+			if (idxs.count(column) != 0) {
+				++skipped;
+				continue;
+			}
+			size_t _col = column - skipped;
+			for (size_t row = 0; row < rows; ++row)
+				data_new[row * columns_new + _col] = data[row * columns + column];
+		}
+		stringvec names_new = stringvec(columns_new);
+		for (size_t i = 0; const std::string &name : names)
+			if (idxs.count(index(name)) == 0) names_new[i++] = name;
+		names = std::move(names_new);
+		data = std::move(data_new);
+		columns = columns_new;
+		return *this;
+	}
+
+
+	/// Appends a row to the table. if name is set, appends it to `opt_rownames`
 	void add_row(std::vector<dtype> values, std::optional<std::string> name = std::nullopt) {
 		if (values.size() == 0) values = std::vector<dtype>(columns, dtype{});
 		data.resize(columns * (++rows));
@@ -237,6 +319,7 @@ public:
 		if (name.has_value()) opt_rownames.push_back(*name);
 	}
 
+
 	friend std::ostream& operator<<(std::ostream &os, table<dtype> &t) {
 		t.print(os);
 		return os;
@@ -300,6 +383,15 @@ public:
 	void fill_column(const std::string &column, dtype v) {
 		apply([&v](const std::vector<dtype>& _) -> dtype { return v; }, {}, column);
 	}
+	
+	/// UNTESTED!
+	template <typename mult>
+	table &multiply_column(const std::string &column, mult s) {
+		size_t i;
+		for (i = index(column); i < rows * columns; i += columns)
+			data[i] *= s;
+		return *this;
+	}
 
 	/// returns an std::pair with coefficients A and B in that order
 	std::pair<prak::pvalue<dtype>, prak::pvalue<dtype>> 
@@ -331,7 +423,11 @@ public:
 
 	/// calculate standard deviation of the column
 	dtype col_stddev(const std::string &column) {
-		assert(0);
+		dtype accum = dtype{};
+		dtype avg = col_avg(column);
+		for (auto it = begin(column); it != end(column); ++it)
+			accum += (*it - avg)*(*it - avg);
+		return std::sqrt(accum);
 	}
 
 	/// Serialize data in format `data[args[0]][i] data[args[1]][i] data[args[2]][i]...`