Saving... I was in the process of writing my physics sim when I noticed that Ozon delivered my new headphones. Which would also imply that I am already too late

src/l1/anne/codegen.c

@@ -14,6 +14,7 @@
 #include "margaret/png_pixel_masses.h"
 #include "lucy.h"
 #include "alice.h"
+#include "precise_integers.h"
 
 int main() {
     mkdir_nofail("l1");
@@ -31,6 +32,7 @@ int main() {
     generate_margaret_png_pixel_masses_header();
     generate_l1_lucy_headers();
     generate_code_for_alice_on_l1();
+    generate_l1_header_for_precise_integers();
     finish_layer(cstr("l1"));
     return 0;
 }

src/l1/anne/geom.h

@@ -2,6 +2,8 @@
 
 #include "../codegen/codegen.h"
 
+// todo: rewrite all of this
+
 SpanU8 vec_field_name(int ci) {
     assert(0 <= ci && ci < 4);
     return ci == 0 ? cstr("x") : (ci == 1 ? cstr("y") : (ci == 2 ? cstr("z") : cstr("w")));
@@ -55,6 +57,17 @@ NODISCARD VecU8 generate_xvecn_struct_and_base_methods(SpanU8 xvec, SpanU8 memb,
         VecU8_append_vec(&res, VecU8_fmt("-A.%s", vec_field_name(ci)));
     }
     VecU8_append_span(&res, cstr(" };\n}\n\n"));
+    /* xvecn_mul_scal method */
+    VecU8_append_vec(&res, VecU8_fmt(
+    "%s %s_mul_scal(%s A, %s B) {\n"
+    SPACE "return (%s) { ",
+        xvecn, xvecn, xvecn, memb, xvecn));
+    for (int ci = 0; ci < n; ci++) {
+        if (ci)
+            VecU8_append_span(&res, cstr(", "));
+        VecU8_append_vec(&res, VecU8_fmt("A.%s * B", vec_field_name(ci)));
+    }
+    VecU8_append_span(&res, cstr(" };\n}\n\n"));
 
     VecU8_drop(g_xvecn);
     return res;
@@ -75,17 +88,6 @@ NODISCARD VecU8 generate_xvecn_struct_and_cool_methods(SpanU8 xvec, SpanU8 memb,
         VecU8_append_vec(&res, VecU8_fmt("A.%s * A.%s", vec_field_name(i), vec_field_name(i)));
     }
     VecU8_append_span(&res, cstr(");\n}\n\n"));
-    /* xvecn_mul_scal method */
-    VecU8_append_vec(&res, VecU8_fmt(
-    "%s %s_mul_scal(%s A, %s B) {\n"
-    SPACE "return (%s) { ",
-        xvecn, xvecn, xvecn, memb, xvecn));
-    for (int ci = 0; ci < n; ci++) {
-        if (ci)
-            VecU8_append_span(&res, cstr(", "));
-        VecU8_append_vec(&res, VecU8_fmt("A.%s * B", vec_field_name(ci)));
-    }
-    VecU8_append_span(&res, cstr(" };\n}\n\n"));
     /* xvecn_div_by_scal method */
     VecU8_append_vec(&res, VecU8_fmt(
     "%s %s_div_by_scal(%s A, %s B) {\n"
@@ -312,28 +314,93 @@ NODISCARD VecU8 generate_xmatnm_transpose_method(SpanU8 xmat, SpanU8 xvec, SpanU
     return res;
 }
 
-NODISCARD VecU8 generate_square_xmatn_methods(SpanU8 xmat, SpanU8 xvec, SpanU8 memb, int n) {
-    VecU8 g_xmatn = codegen_name_xmatnm(xmat, n, n);
-    SpanU8 xmatn = VecU8_to_span(&g_xmatn);
+NODISCARD VecU8 generate_xmat234_det_method(SpanU8 xmat, SpanU8 xvec, SpanU8 memb){
     VecU8 res = VecU8_new();
+    VecU8_append_vec(&res, VecU8_fmt(
+        "%s %s2_det(%s2 a){\n"
+        SPACE "return a.x.x * a.y.y - a.x.y * a.y.x;\n"
+        "}\n\n", memb, xmat, xmat));
+    VecU8_append_vec(&res, VecU8_fmt(
+        "%s %s3_det(%s3 a){\n"
+        SPACE "return", memb, xmat, xmat));
+    for (int i = 0; i < 3; i++) {
+        if (i)
+            VecU8_append_span(&res, cstr("\n" SPACE SPACE));
+        for (int j = 0; j < 3; j++) {
+            if (j == i)
+                continue;
+            int k = 3 - i - j;
+            assert(k != i && k != j);
+            bool minused = ((i > j) + (j > k) + (i > k)) % 2;
+            VecU8_append_vec(&res, VecU8_fmt(" %c a.x.%s * a.y.%s * a.z.%s",
+                minused ? '-' : '+', vec_field_name(i), vec_field_name(j), vec_field_name(k)));
+        }
+    }
+    VecU8_append_span(&res, cstr(";\n}\n\n"));
 
     VecU8_append_vec(&res, VecU8_fmt(
+        "%s %s4_det(%s4 a){\n"
+        SPACE "return", memb, xmat, xmat));
+    for (int i = 0; i < 4; i++) {
+        if (i)
+            VecU8_append_span(&res, cstr("\n" SPACE SPACE));
+        for (int j = 0; j < 4; j++) {
+            if (j == i)
+                continue;
+            for (int k = 0; k < 4; k++) {
+                if (i == k || j == k)
+                    continue;
+                int u = 6 - i - j - k;
+                assert(u != i && u != j && u != k);
+                bool minused = ((i > j) + (j > k) + (i > k) + (i > u) + (j > u) + (k > u)) % 2;
+                VecU8_append_vec(&res, VecU8_fmt(" %c a.x.%s * a.y.%s * a.z.%s * a.w.%s",
+                    minused ? '-' : '+', vec_field_name(i), vec_field_name(j), vec_field_name(k), vec_field_name(u)));
+            }
+        }
+    }
+
+    VecU8_append_span(&res, cstr(";\n}\n\n"));
+    return res;
+}
+
+void codegen_append_square_xmatn_methods(VecU8* res, SpanU8 xmat, SpanU8 xvec, SpanU8 memb, int n) {
+    VecU8 g_xmatn = codegen_name_xmatnm(xmat, n, n);
+    SpanU8 xmatn = VecU8_to_span(&g_xmatn);
+    VecU8 g_xvecn = codegen_name_xvecn(xvec, n);
+    SpanU8 xvecn = VecU8_to_span(&g_xvecn);
+
+    VecU8_append_vec(res, VecU8_fmt(
         "const %s %s_E = { ", xmatn, xmatn));
     for (int x = 0; x < n; x++) {
         if (x)
-            VecU8_append_span(&res, cstr(", "));
-        VecU8_append_vec(&res, VecU8_fmt(".%s = {", vec_field_name(x)));
+            VecU8_append_span(res, cstr(", "));
+        VecU8_append_vec(res, VecU8_fmt(".%s = {", vec_field_name(x)));
         for (int y = 0; y < n; y++) {
             if (y)
-                VecU8_append_span(&res, cstr(", "));
-            VecU8_append(&res, '0' + (x == y));
+                VecU8_append_span(res, cstr(", "));
+            VecU8_append(res, '0' + (x == y));
         }
-        VecU8_append_span(&res, cstr(" }"));
+        VecU8_append_span(res, cstr(" }"));
     }
-    VecU8_append_span(&res, cstr(" };\n\n"));
+    VecU8_append_span(res, cstr(" };\n\n"));
 
+    VecU8_append_vec(res, VecU8_fmt(
+        "%s %s_new_for_proj(%s v){\n"     /* xmatn, xmatn, xvecn */
+        SPACE "return (%s){ ",      /* xmatn */
+        xmatn, xmatn, xvecn, xmatn));
+    for (int x = 0; x < n; x++) {
+        VecU8_append_span(res, cstr("\n" SPACE SPACE));
+        VecU8_append_vec(res, VecU8_fmt(".%s={", vec_field_name(x)));
+        for (int y = 0; y < n; y++) {
+            if (y)
+                VecU8_append_span(res, cstr(", "));
+            VecU8_append_vec(res, VecU8_fmt("v.%s * v.%s", vec_field_name(x), vec_field_name(y)));
+        }
+        VecU8_append_span(res, cstr("},"));
+    }
+    VecU8_append_span(res, cstr("};\n}\n\n"));
     VecU8_drop(g_xmatn);
-    return res;
+    VecU8_drop(g_xvecn);
 }
 
 NODISCARD VecU8 generate_xmat_inverse_methods(SpanU8 xmat, SpanU8 xvec, SpanU8 memb){
@@ -466,7 +533,7 @@ NODISCARD VecU8 generate_xvec234_structs_and_cool_methods(SpanU8 xvec, SpanU8 me
     return res;
 }
 
-NODISCARD VecU8 generate_xmat234x234_structs_methods(SpanU8 xmat, SpanU8 xvec, SpanU8 memb, int sizeof_member) {
+NODISCARD VecU8 generate_xmat234x234_structs_and_base_methods(SpanU8 xmat, SpanU8 xvec, SpanU8 memb, int sizeof_member){
     VecU8 res = VecU8_new();
     for (int cols = 2; cols <= 4; cols++) {
         for (int rows = 2; rows <= 4; rows++) {
@@ -479,7 +546,7 @@ NODISCARD VecU8 generate_xmat234x234_structs_methods(SpanU8 xmat, SpanU8 xvec, S
         }
     }
     for (int n = 2; n <= 4; n++) {
-        VecU8_append_vec(&res, generate_square_xmatn_methods(xmat, xvec, memb, n));
+        codegen_append_square_xmatn_methods(&res, xmat, xvec, memb, n);
     }
     for (int n = 2; n <= 4; n++) {
         for (int m = 2; m <= 4; m++) {
@@ -488,6 +555,13 @@ NODISCARD VecU8 generate_xmat234x234_structs_methods(SpanU8 xmat, SpanU8 xvec, S
             }
         }
     }
+    VecU8_append_vec(&res, generate_xmat234_det_method(xmat, xvec, memb));
+    return res;
+}
+
+NODISCARD VecU8 generate_xmat234x234_structs_and_cool_methods(SpanU8 xmat, SpanU8 xvec, SpanU8 memb, int sizeof_member){
+    VecU8 res = VecU8_new();
+    VecU8_append_vec(&res, generate_xmat234x234_structs_and_base_methods(xmat, xvec, memb, sizeof_member));
     VecU8_append_vec(&res, generate_xmat_inverse_methods(xmat, xvec, memb));
     return res;
 }
@@ -507,7 +581,10 @@ void generate_geom_header() {
 
     VecU8_append_vec(&res.result, generate_xvec234_structs_and_cool_methods(cstr("vec"), cstr("float"), cstr("sqrtf")));
     VecU8_append_vec(&res.result, generate_xvec234_structs_and_cool_methods(cstr("dvec"), cstr("double"), cstr("sqrt")));
-    VecU8_append_vec(&res.result, generate_xmat234x234_structs_methods(cstr("mat"), cstr("vec"), cstr("float"), sizeof(float)));
-    /* VecU8_append_vec(&res.result, generate_xmat234x234_structs_methods(cstr("dmat"), cstr("dvec"), cstr("double"), sizeof(double))); */
+
+    // todo: remove padding from matrix structure. VERY IMPORTANT!!! Add padding on the fly when transferring to vulkan
+    VecU8_append_vec(&res.result, generate_xmat234x234_structs_and_cool_methods(cstr("mat"), cstr("vec"), cstr("float"), sizeof(float)));
+    VecU8_append_vec(&res.result, generate_xmat234x234_structs_and_base_methods(cstr("s64mat"), cstr("s64vec"), cstr("S64"), sizeof(S64)));
+
     finish_header(res);
 }

src/l1/anne/precise_integers.h

@@ -0,0 +1,20 @@
+#pragma once
+
+// todo: write something here after I close КТ-1
+
+#include "../codegen/codegen.h"
+
+void generate_precise_addition_func(VecU8* res, U64 liA, U64 liB){
+    assert(liA < 1000 && liB < 1000);
+    assert(liA >= liB);
+    VecU8_append_vec(res, VecU8_fmt(
+        "void S%u_add_S%u(S%u* a, S%u b){\n", liA * 64, liB * 64, liA * 64, liB * 64));
+    // for (size_t i = 10; )
+    VecU8_append_span(res, cstr("}\n\n"));
+}
+
+void generate_l1_header_for_precise_integers(){
+    GeneratedHeader h = begin_header(cstr("l1/precise_integers.h"));
+    VecU8_append_span(&h.result, cstr("#include \"../../src/l1/core/int_primitives.h\"\n"));
+    finish_header(h);
+}

src/l1/core/VecU8_as_str.h

@@ -40,6 +40,7 @@ void VecU8_print(VecU8 self){
     VecU8_drop(self);
 }
 
+/* Wrapper around vsnprintf */
 NODISCARD VecU8 VecU8_format(const char *fmt, ...) {
     assert(fmt);
 
@@ -152,17 +153,10 @@ VecU8 U64_stringification(U64 x){
     return res;
 }
 
-/* %s - SpanU8
- * %v - VecU8
- * %u - U64
- * %c - int (one byte character)
- * %i - S64
- */
-NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
-    assert(fmt);
+NODISCARD VecU8 VecU8_fmt_va_list(const char* fmt, va_list args){
     size_t k = 0;
-    va_list args;
-    va_start(args, fmt);
+    va_list repeat;
+    va_copy(repeat, args);
     for (const char *ch = fmt; *ch; ) {
         if (*ch == '%') {
             ch++;
@@ -175,7 +169,7 @@ NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
                 va_arg(args, int);
                 k++;
             } else if (*ch == 'v') {
-                /* We had not taken ownership of YET (will take it the next iteration) */
+                /* We had not taken ownership of it YET (will take it the next iteration) */
                 VecU8 vs = va_arg(args, VecU8);
                 k += vs.len;
             } else if (*ch == 'i') {
@@ -184,6 +178,9 @@ NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
             } else if (*ch == 'u') {
                 U64 x = va_arg(args, U64);
                 k += U64_stringification_get_length(x);
+            } else if (*ch == 'r') {
+                VecU8 vs = va_arg(args, VecU8);
+                k += vs.len;
             } else
                 abortf("Format syntax error at pos %lu! Watch out, be careful", (size_t)(ch - fmt));
         } else {
@@ -191,29 +188,31 @@ NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
         }
         ch++;
     }
-    va_end(args);
     VecU8 res = VecU8_new_reserved(k);
-    va_start(args, fmt);
     for (const char *ch = fmt; *ch;) {
         if (*ch == '%') {
             ch++;
             if (*ch == '%') {
                 VecU8_append(&res, '%');
             } else if (*ch == 's') {
-                SpanU8 s = va_arg(args, SpanU8);
+                SpanU8 s = va_arg(repeat, SpanU8);
                 VecU8_append_span(&res, s);
             } else if (*ch == 'c') {
-                int byte = va_arg(args, int);
+                int byte = va_arg(repeat, int);
                 VecU8_append(&res, (U8)byte);
             } else if (*ch == 'v') {
-                VecU8 vs = va_arg(args, VecU8);
+                VecU8 vs = va_arg(repeat, VecU8);
                 VecU8_append_vec(&res, vs); /* Moved ownership of vs argument */
             } else if (*ch == 'i') {
-                S64 x = va_arg(args, S64);
+                S64 x = va_arg(repeat, S64);
                 S64_stringification_into_buf(x, &res);
             } else if (*ch == 'u') {
-                U64 x = va_arg(args, U64);
+                U64 x = va_arg(repeat, U64);
                 U64_stringification_into_buf(x, &res);
+            } else if (*ch == 'r') {
+                VecU8 vs = va_arg(repeat, VecU8);
+                VecU8_append_span(&res, VecU8_to_span(&vs));
+                /* But we won't deallocate `vs` because never took ownership of it */
             } else
                 assert(false);
         } else {
@@ -221,10 +220,40 @@ NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
         }
         ch++;
     }
+    va_end(repeat);
+    return res;
+}
+
+/* %s - SpanU8
+ * %v - VecU8
+ * %u - U64
+ * %c - int (one byte character)
+ * %i - S64
+ * %r - VecU8 (but taken by reference, it won't be deallocated)
+ * Yes, I am breaking the dumb ownership-syntax-crutch convention I set up myself
+ */
+NODISCARD VecU8 VecU8_fmt(const char* fmt, ...) {
+    assert(fmt);
+    va_list args;
+    va_start(args, fmt);
+    VecU8 res = VecU8_fmt_va_list(fmt, args);
     va_end(args);
     return res;
 }
 
+void VecU8_append_cstr(VecU8* res, const char* lit){
+    VecU8_append_span(res, SpanU8_from_cstr(lit));
+}
+
+// todo: add VecX_resize method. I really REALLY need to refactor this awful code
+void VecU8_append_fmt(VecU8* res, const char* fmt, ...){
+    va_list args;
+    va_start(args, fmt);
+    VecU8 t = VecU8_fmt_va_list(fmt, args);
+    VecU8_append_vec(res, t);
+    va_end(args);
+}
+
 // todo: generate a special span method to check equality of contents
 bool SpanU8_cont_equal(SpanU8 a, SpanU8 b) {
     if (a.len != b.len)

src/l1/marie/geom_alg_utils.h

@@ -36,60 +36,6 @@ bool marie_same_dir3(float A0, float A1, float A2, float B) {
     return A0 == 0 || A1 == 0 || A2 == 0 || (diff & 0x80000000u) == 0;
 }
 
-// This function uses floating point arithmetic, which is very bad
-/* cross product (vi - u) x (vj - u). If {vi, vj, vu} are points of CCW triangle, this will be it's positive surface */
-float marie_surface(vec2 vi, vec2 vj, vec2 u) {
-    return u.x * (vi.y - vj.y) + u.y * (vj.x - vi.x) + (vi.x * vj.y - vj.x * vi.y);
-}
-
-typedef vec4 MarieVertAttr;
-
-typedef struct {
-    vec2 pos;
-    MarieVertAttr attr;
-} MariePlaneVertAttr;
-
-typedef struct {
-    vec2 v0;
-    vec2 v1;
-    vec2 v2;
-} MarieTriangle;
-
-MarieTriangle MarieTriangle_mat3x2_mul_pos(MarieTriangle self, mat3x2 trop) {
-    return (MarieTriangle){
-        .v0 = mat3x2_mul_vec3(trop, vec2_and_one(self.v0)),
-        .v1 = mat3x2_mul_vec3(trop, vec2_and_one(self.v1)),
-        .v2 = mat3x2_mul_vec3(trop, vec2_and_one(self.v2))};
-}
-
-#include "../../../gen/l1/eve/marie/VecMarieTriangle.h"
-
-typedef struct {
-    MariePlaneVertAttr v0;
-    MariePlaneVertAttr v1;
-    MariePlaneVertAttr v2;
-} MarieTriangleAttr;
-
-MarieTriangleAttr MarieTriangle_goto_nat_cords_pres_par(MarieTriangle self, mat3x2 trop) {
-    return (MarieTriangleAttr){
-        {mat3x2_mul_vec3(trop, vec2_and_one(self.v0)), {self.v0.x, self.v0.y, 0, 0}},
-        {mat3x2_mul_vec3(trop, vec2_and_one(self.v1)), {self.v1.x, self.v1.y, 0, 0}},
-        {mat3x2_mul_vec3(trop, vec2_and_one(self.v2)), {self.v2.x, self.v2.y, 0, 0}},
-    };
-}
-
-#include "../../../gen/l1/eve/marie/VecMarieTriangleAttr.h"
-
-vec2 marie_intersect_lines(vec2 A1, vec2 B1, vec2 A2, vec2 B2) {
-    vec2 alpha = vec2_minus_vec2(B1, A1);
-    vec2 beta = vec2_minus_vec2(A2, B2);
-    vec2 gamma = vec2_minus_vec2(A2, A1);
-    float det_alpha_beta = alpha.x * beta.y - alpha.y * beta.x;
-    float det_gamma_beta = gamma.x * beta.y - gamma.y * beta.x;
-    float t1 = det_gamma_beta / det_alpha_beta;
-    return vec2_add_vec2(A1, vec2_mul_scal(alpha, t1));
-}
-
 float marie_clamp_float(float a, float l, float r) {
     return (a < l) ? l : ((a <= r) ? a : r);
 }

src/l1_5/core/quaternion.h

@@ -0,0 +1,59 @@
+#pragma once
+
+#include "../marie/geom_util.h"
+#include "math.h"
+
+typedef vec4 quaternion_t;
+
+vec4 quaternion_fl_mul_vec4(vec4 A, vec4 B){
+    return (vec4){
+        A.x * B.x - A.y * B.y - A.z * B.z - A.w * B.w,
+        A.x * B.y + A.y * B.x + A.z * B.w - A.w * B.z,
+        A.x * B.z - A.y * B.w + A.z * B.x + A.w * B.y,
+        A.x * B.w + A.y * B.z - A.z * B.y + A.w * B.x,
+    };
+}
+
+/* Assumes abs(a) is 1 */
+vec4 unit_quaternion_fl_inverse(vec4 a){
+    return (vec4){a.x, -a.y, -a.z, -a.w};
+}
+
+vec4 quaternion_fl_inverse(vec4 a){
+    float m = vec4_length(a);
+    m = 1 / (m * m);
+    return (vec4){a.x * m, -a.y * m, -a.z * m, -a.w * m};
+}
+
+mat3 quaternion_fl_to_rot3d_mat3(vec4 q){
+    float a = q.x;
+    float Rx = q.y;
+    float Ry = q.z;
+    float Rz = q.w;
+
+    return mat3_add_mat3(mat3_add_mat3(
+        mat3_mul_scal(mat3_E, 2*a*a - 1),
+        mat3_mul_scal(mat3_for_cross_product((vec3){Rx, Ry, Rz}), 2 * a)), mat3_E);
+}
+
+/* This thing requires some MATH */
+vec4 unit_quaternion_fl_power(vec4 q, float t){
+    float qva = vec3_length((vec3){q.y, q.z, q.w});
+    float fi = atan2f(qva, q.x);
+    const float rub = 0.001f;
+    float e_x = cosf(t * fi);
+    float e_vc;
+    // float t_to_3 = t*t*t;
+    // float t_to_5 = t_to_3*t*t;
+    // float fi_to_2 = fi * fi;
+    // float fi_to_4 = fi_to_2 * fi_to_2;
+    // Remember, qva = sin(fi)
+    if (qva < rub) {
+        /* I have a feeling this is an overkill... Well, only time __and 100 hours of testing__ will show ;) */
+        // e_vc = t + (-t_to_3+t) / 6 * fi_to_2 + (3 * t_to_5 - 10 * t_to_3 + 7 * t) / 360 * fi_to_4;
+        e_vc = t + (-t*t*t+t) / 6 * fi * fi ;
+    } else {
+        e_vc = sinf(t *fi) / qva;
+    }
+    return (vec4){e_x, e_vc * q.y, e_vc * q.z, e_vc * q.w};
+}

src/l1_5/marie/geom_util.h

@@ -0,0 +1,7 @@
+#pragma once
+
+#include "../../../gen/l1/geom.h"
+
+mat3 mat3_for_cross_product(vec3 r){
+    return (mat3){.x = {0, -r.z, r.y}, .y = {r.z, 0, -r.x}, .z = {-r.y, r.x, 0}};
+}

src/l1_5/marie/prim_shape_geom.h

@@ -0,0 +1,95 @@
+#pragma once
+
+#include "../../../gen/l1/geom.h"
+
+/* cross product (vi - u) x (vj - u). If {vi, vj, vu} are points of CCW triangle, this will be it's positive surface */
+float marie_surface(vec2 vi, vec2 vj, vec2 u) {
+    return u.x * (vi.y - vj.y) + u.y * (vj.x - vi.x) + (vi.x * vj.y - vj.x * vi.y);
+}
+
+/* 2 lines on a plane (of float type). Assumes lines are intersecting */
+vec2 marie_intersect_lines(vec2 A1, vec2 B1, vec2 A2, vec2 B2) {
+    vec2 alpha = vec2_minus_vec2(B1, A1);
+    vec2 beta = vec2_minus_vec2(A2, B2);
+    vec2 gamma = vec2_minus_vec2(A2, A1);
+    float det_alpha_beta = alpha.x * beta.y - alpha.y * beta.x;
+    float det_gamma_beta = gamma.x * beta.y - gamma.y * beta.x;
+    float t1 = det_gamma_beta / det_alpha_beta;
+    return vec2_add_vec2(A1, vec2_mul_scal(alpha, t1));
+}
+
+/* Triangle (ABC) and a ray (from origin P to point Q) in 3d space (of float type). Returns bool: do they intersect */
+bool marie_intersect_triangle_and_ray(vec3 A, vec3 B, vec3 C, vec3 P, vec3 Q, vec3* ret){
+    vec3 alpha = vec3_minus_vec3(A, C);
+
+    vec3 beta = vec3_minus_vec3(B, C);
+    vec3 gamma = vec3_minus_vec3(P, Q);
+    vec3 delta = vec3_minus_vec3(P, C);
+    float det_alpha_beta_gamma = mat3_det((mat3){.x = alpha, .y = beta, .z = gamma});
+    float det_delta_beta_gamma = mat3_det((mat3){.x = delta, .y = beta, .z = gamma});
+    float det_alpha_delta_gamma = mat3_det((mat3){.x = alpha, .y = delta, .z = gamma});
+    float det_alpha_beta_delta = mat3_det((mat3){.x = alpha, .y = beta, .z = delta});
+    if (det_alpha_beta_gamma < 0.0001)
+        return false;
+    float b_0 = det_delta_beta_gamma / det_alpha_beta_gamma;
+    float b_1 = det_delta_beta_gamma / det_alpha_beta_gamma;
+    float t = det_delta_beta_gamma / det_alpha_beta_gamma;
+    if (b_0 < 0 || b_0 > 1 || b_1 < 0 || b_1 > 1 || t < 0)
+        return false;
+    *ret = vec3_add_vec3(P, vec3_mul_scal(gamma, -t));
+    return true;
+}
+
+
+
+/* any inv exp q => any abs bitness a =>
+ * vec2#(a)#.(q) -> vec2#(a)#.(q) -> vec2#(a)#.(q) -> vec2#(2a+3)#.(2q) */
+S64 marie_precise_surface(s64vec2 vi, s64vec2 vj, s64vec2 u){
+    s64vec2 da = s64vec2_minus_s64vec2(vi, u);
+    s64vec2 db = s64vec2_minus_s64vec2(vj, u);
+    /* da, db are vec2#(a+1)#.(q) */
+    return da.x * db.y - da.y * db.x;
+}
+
+bool marie_precise_do_intersect_2d_interv_and_seg(s64vec2 a, s64vec2 b, s64vec2 c, s64vec2 d){
+    s64vec2 alpha = s64vec2_minus_s64vec2(b, a);
+    s64vec2 beta = s64vec2_minus_s64vec2(c, d);
+    s64vec2 gamma = s64vec2_minus_s64vec2(c, a);
+    S64 det_alpha_beta = alpha.x * beta.y - alpha.y * beta.x;
+    S64 det_gamma_beta = gamma.x * beta.y - gamma.y * beta.x;
+    if (det_alpha_beta == 0) {
+        if (det_gamma_beta != 0)
+            return false;
+        if (a.x < b.x) {
+            return !(c.x <= a.x && d.x <= a.x) && !(b.x <= c.x && b.x <= d.x);
+        } else if (b.x < a.x) {
+            return !(c.x <= b.x && d.x <= b.x) && !(a.x <= c.x && a.x <= d.x);
+        } else if  (a.y < b.y) {
+            return !(c.y <= a.y && d.y <= a.y) && !(b.y <= c.y && b.y <= d.y);
+        } else if (b.y < a.y) {
+            return !(c.y <= b.y && d.y <= b.y) && !(a.y <= c.y && a.y <= d.y);
+        } else
+            return false;  /* This code is unreachable, actually */
+    }
+    /* Return det_gamma_beta/det_alpha_beta in (0; 1) */
+    return 0 < det_gamma_beta && det_gamma_beta < det_alpha_beta;
+}
+
+/* Does not work for degenerate case where s(b_prev, b, b_next) = 0.
+ * Returns false if a is on the edge of angle(b_prev, b, b_next) */
+bool marie_precise_is_in_ccw_angle(s64vec2 b_prev, s64vec2 b, s64vec2 b_next, s64vec2 a){
+    S64 sx = marie_precise_surface(b_prev, b, b_next);
+    S64 sy = marie_precise_surface(b_prev, b, a);
+    S64 sz = marie_precise_surface(b, b_next, a);
+    if (sx < 0) {
+        return sy > 0 || sz > 0;
+    } else {
+        return sy > 0 && sz > 0;
+    }
+}
+
+/* Read source code to understand */
+bool marie_order_s64vec2_less(s64vec2 a, s64vec2 b){
+    return a.x < b.x || (a.x == b.x && a.y < b.y);
+}
+

src/l2/alice/assets.h

@@ -1,6 +1,6 @@
 #pragma once
 
-#include "../marie/graphics_geom.h"
+#include "../../l1_5/marie/graphics_geom.h"
 #include "../../../gen/l1/VecAndSpan_U32.h"
 #include "../../../gen/l1/VecAndSpan_U8.h"
 

src/l2/marie/rasterization.h

@@ -2,6 +2,46 @@
 
 #include "../../../gen/l1/geom.h"
 #include "../../l1/marie/geom_alg_utils.h"
+#include "../../l1_5/marie/prim_shape_geom.h"
+
+typedef vec4 MarieVertAttr;
+
+typedef struct {
+    vec2 pos;
+    MarieVertAttr attr;
+} MariePlaneVertAttr;
+
+typedef struct {
+    vec2 v0;
+    vec2 v1;
+    vec2 v2;
+} MarieTriangle;
+
+MarieTriangle MarieTriangle_mat3x2_mul_pos(MarieTriangle self, mat3x2 trop) {
+    return (MarieTriangle){
+        .v0 = mat3x2_mul_vec3(trop, vec2_and_one(self.v0)),
+        .v1 = mat3x2_mul_vec3(trop, vec2_and_one(self.v1)),
+        .v2 = mat3x2_mul_vec3(trop, vec2_and_one(self.v2))};
+}
+
+#include "../../../gen/l1/eve/marie/VecMarieTriangle.h"
+
+typedef struct {
+    MariePlaneVertAttr v0;
+    MariePlaneVertAttr v1;
+    MariePlaneVertAttr v2;
+} MarieTriangleAttr;
+
+MarieTriangleAttr MarieTriangle_goto_nat_cords_pres_par(MarieTriangle self, mat3x2 trop) {
+    return (MarieTriangleAttr){
+            {mat3x2_mul_vec3(trop, vec2_and_one(self.v0)), {self.v0.x, self.v0.y, 0, 0}},
+            {mat3x2_mul_vec3(trop, vec2_and_one(self.v1)), {self.v1.x, self.v1.y, 0, 0}},
+            {mat3x2_mul_vec3(trop, vec2_and_one(self.v2)), {self.v2.x, self.v2.y, 0, 0}},
+        };
+}
+
+#include "../../../gen/l1/eve/marie/VecMarieTriangleAttr.h"
+
 
 typedef struct {
     /* guest, x, y, attribute (custom) */

src/l2/marie/shape_geom.h

@@ -8,57 +8,8 @@
 #include "../../../gen/l1/VecAndSpan_VecU64.h"
 #include <stdbool.h>
 #include <string.h>
+#include "../../l1_5/marie/prim_shape_geom.h"
 
-/* any inv exp q => any abs bitness a =>
- * vec2#(a)#.(q) -> vec2#(a)#.(q) -> vec2#(a)#.(q) -> vec2#(2a+3)#.(2q) */
-S64 marie_precise_surface(s64vec2 vi, s64vec2 vj, s64vec2 u){
-    s64vec2 da = s64vec2_minus_s64vec2(vi, u);
-    s64vec2 db = s64vec2_minus_s64vec2(vj, u);
-    /* da, db are vec2#(a+1)#.(q) */
-    return da.x * db.y - da.y * db.x;
-}
-
-bool marie_precise_do_intersect_2d_interv_and_seg(s64vec2 a, s64vec2 b, s64vec2 c, s64vec2 d){
-    s64vec2 alpha = s64vec2_minus_s64vec2(b, a);
-    s64vec2 beta = s64vec2_minus_s64vec2(c, d);
-    s64vec2 gamma = s64vec2_minus_s64vec2(c, a);
-    S64 det_alpha_beta = alpha.x * beta.y - alpha.y * beta.x;
-    S64 det_gamma_beta = gamma.x * beta.y - gamma.y * beta.x;
-    if (det_alpha_beta == 0) {
-        if (det_gamma_beta != 0)
-            return false;
-        if (a.x < b.x) {
-            return !(c.x <= a.x && d.x <= a.x) && !(b.x <= c.x && b.x <= d.x);
-        } else if (b.x < a.x) {
-            return !(c.x <= b.x && d.x <= b.x) && !(a.x <= c.x && a.x <= d.x);
-        } else if  (a.y < b.y) {
-            return !(c.y <= a.y && d.y <= a.y) && !(b.y <= c.y && b.y <= d.y);
-        } else if (b.y < a.y) {
-            return !(c.y <= b.y && d.y <= b.y) && !(a.y <= c.y && a.y <= d.y);
-        } else
-            return false;  /* This code is unreachable, actually */
-    }
-    /* Return det_gamma_beta/det_alpha_beta in (0; 1) */
-    return 0 < det_gamma_beta && det_gamma_beta < det_alpha_beta;
-}
-
-/* Does not work for degenerate case where s(b_prev, b, b_next) = 0.
- * Returns false if a is on the edge of angle(b_prev, b, b_next) */
-bool marie_precise_is_in_ccw_angle(s64vec2 b_prev, s64vec2 b, s64vec2 b_next, s64vec2 a){
-    S64 sx = marie_precise_surface(b_prev, b, b_next);
-    S64 sy = marie_precise_surface(b_prev, b, a);
-    S64 sz = marie_precise_surface(b, b_next, a);
-    if (sx < 0) {
-        return sy > 0 || sz > 0;
-    } else {
-        return sy > 0 && sz > 0;
-    }
-}
-
-/* Read source code to understand */
-bool marie_order_s64vec2_less(s64vec2 a, s64vec2 b){
-    return a.x < b.x || (a.x == b.x && a.y < b.y);
-}
 
 typedef struct{
     U64 prev, next;

src/l2/tests/r2/r2a.c

@@ -5,7 +5,7 @@
 #include <stdio.h>
 #include "../../../l1/system/pthread.h"
 #include "../../margaret/time_utils.h"
-#include "../../marie/graphics_geom.h"
+#include "../../../l1_5/marie/graphics_geom.h"
 #include "../../marie/rasterization.h"
 
 #include <fcntl.h>