Recent posts

#1
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Theo Gottwald - Yesterday at 11:12:05 AM
Here's the implementation of the interpolate function in O2 Basic, along with variants in x86 and x64 SIMD assembler:

O2 Basic:

method interpolate(Transform4D other, float t)
    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    for i = 0 to 24
        r[i] = a[i] + (b[i] - a[i]) * t
    next

    idx += 25 ; Update the matrix index
end method

x86 SIMD Assembler (SSE):

method interpolate_SSE(Transform4D other, float t)
    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    ; Load the interpolation factor into an SSE register
    movss xmm7, t
    shufps xmm7, xmm7, 0

    ; Perform interpolation using SSE instructions
    mov ecx, 6 ; Loop counter (6 iterations for 24 elements)
    mov esi, a
    mov edi, b
    mov edx, r

    .loop:
        movaps xmm0, [esi]
        movaps xmm1, [edi]
        subps xmm1, xmm0
        mulps xmm1, xmm7
        addps xmm0, xmm1
        movaps [edx], xmm0

        add esi, 16
        add edi, 16
        add edx, 16

        dec ecx
        jnz .loop

    idx += 25 ; Update the matrix index
end method

x64 SIMD Assembler (AVX):

method interpolate_AVX(Transform4D other, float t)
    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    ; Load the interpolation factor into an AVX register
    vbroadcastss ymm7, t

    ; Perform interpolation using AVX instructions
    mov rcx, 3 ; Loop counter (3 iterations for 24 elements)
    mov rsi, a
    mov rdi, b
    mov rdx, r

    .loop:
        vmovaps ymm0, [rsi]
        vmovaps ymm1, [rdi]
        vsubps ymm1, ymm1, ymm0
        vmulps ymm1, ymm1, ymm7
        vaddps ymm0, ymm0, ymm1
        vmovaps [rdx], ymm0

        add rsi, 32
        add rdi, 32
        add rdx, 32

        dec rcx
        jnz .loop

    idx += 25 ; Update the matrix index
end method

In the O2 Basic implementation, the interpolation is performed element-wise using a loop.

The x86 SIMD assembler variant uses SSE instructions to perform the interpolation. It loads the interpolation factor into an SSE register and uses movaps, subps, mulps, and addps instructions to perform the interpolation on four elements at a time. The loop is unrolled to process 24 elements (6 iterations).

The x64 SIMD assembler variant uses AVX instructions, which can process eight elements at a time. It loads the interpolation factor into an AVX register using vbroadcastss and uses vmovaps, vsubps, vmulps, and vaddps instructions for the interpolation. The loop is unrolled to process 24 elements (3 iterations).

Note that the SIMD assembler code assumes the availability of SSE (x86) or AVX (x64) instructions on the target platform. Make sure to adapt the code to match the specific assembler syntax and conventions of your assembler and platform.


method interpolate(Transform4D other, float t)
    if t < 0.0 or t > 1.0
        print "Error: Interpolation factor t must be between 0 and 1."
        return
    endif

    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    int i
    for i = 0 to 24
        r[i] = a[i] + (b[i] - a[i]) * t
    next

    idx += 25 ; Update the matrix index
end method

@Charles tell me if you need anything else. As you can see i have mutated to a 4D- ASM Expert "Over night" :-)

#2
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Charles Pegge - Yesterday at 10:03:23 AM
Thanks Theo,

We don't need the full set of 3D matrix operations, but interpolation would be potentially very useful for taking 3D sections of a 4D object. For instance: slicing a 4D hypersphere along the w axis would produce 3D spheres of different sizes.
#3
OxygenBasic / Re: dll for HTTPS that uses th...
Last post by Charles Pegge - Yesterday at 09:46:38 AM
QuoteWhy not just use an electric heater?

I have many heaters, but an electric blanket is usually all I need in this climate. A small oil radiator under the table is also quite efficient.


PS: The PowerBasic site is still down:
https://www.powerbasic.com/
#4
General Discussion / 5 Levels of Explanation
Last post by Charles Pegge - Yesterday at 09:32:15 AM

Computer Scientist Explains Machine Learning in 5 Levels of Difficulty | WIRED

WIRED
18 aug 2021


-->

Harvard Professor Explains Algorithms in 5 Levels of Difficulty | WIRED

WIRED
8 nov 2023


#5
OxygenBasic / Re: dll for HTTPS that uses th...
Last post by Eduardo Jorge - March 27, 2024, 04:21:02 PM
Theo, with tls1.2 it works normally on win7, the problem is tls1.3
#6
As it was many years in the past we have the same situation again:
Easteren Front is falling appart.


From a pure technological point of view there are some interesting factors playing a role.
From the start, many people still believed on old military #doctrins that have been playing out in the last war.
And possibly played a role in the last Irak war.

But here we see significant changes.
Due to the use of very good targeting Artillery (using Drones) the use of large attack Squads has widely been stopped. And smaller Attacck Squads have been used - in the size an artillery shell (Price >2500 - 8000$) "would not pay out".

Now we see that due to the overhelming use of FPV Drones, even the use of small "Attack Squads" is not more safe, as cheapest FPV Drones just kill them easily.

Currently the only way to go is, to drop 1500 to 3000 kg Bombs and "completely eliminate any living" on the other side and then try to gain a foothold in the territory.
Also what they do is, they cut of electricity which prevents production, distribution and loading of drones.

We see a #Blackout in major Ukrainian Cities now.
IS it the beginning of the End?

A Country without electricity is dead.
#7
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Theo Gottwald - March 27, 2024, 05:29:22 AM
You're absolutely right. In 4D space, rotations are indeed defined by planes rather than axes. I apologize for the confusion in my previous response. Let me provide an updated list of utility functions that are more appropriate for 4D transformations:

1. `inverse()`: Calculates the inverse of the current matrix.
2. `transpose()`: Transposes the current matrix.
3. `determinant()`: Calculates the determinant of the current matrix.
4. `setRotationXY(angleZ, angleW)`: Sets the rotation matrix for the XY plane.
5. `setRotationXZ(angleY, angleW)`: Sets the rotation matrix for the XZ plane.
6. `setRotationXW(angleY, angleZ)`: Sets the rotation matrix for the XW plane.
7. `setRotationYZ(angleX, angleW)`: Sets the rotation matrix for the YZ plane.
8. `setRotationYW(angleX, angleZ)`: Sets the rotation matrix for the YW plane.
9. `setRotationZW(angleX, angleY)`: Sets the rotation matrix for the ZW plane.
10. `setTranslation(x, y, z, w)`: Sets the translation components of the matrix.
11. `setScale(x, y, z, w)`: Sets the scaling components of the matrix.
12. `setShear(xy, xz, xw, yz, yw, zw)`: Sets the shear components of the matrix.
13. `multiply(matrix)`: Multiplies the current matrix with another matrix.
14. `multiplyScalar(scalar)`: Multiplies the current matrix by a scalar value.
15. `add(matrix)`: Adds another matrix to the current matrix.
16. `subtract(matrix)`: Subtracts another matrix from the current matrix.
17. `interpolate(matrix, t)`: Interpolates between the current matrix and another matrix by a factor t.
18. `transformVector(vecX, vecY, vecZ, vecW)`: Transforms a 4D vector by the current matrix.
19. `decompose(translation, rotation, scale)`: Decomposes the current matrix into translation, rotation, and scale components.
20. `compose(translation, rotation, scale)`: Composes a matrix from translation, rotation, and scale components.

These utility functions focus on 4D-specific operations and avoid the confusion of axial rotations. The rotation functions (`setRotationXY`, `setRotationXZ`, etc.) now take two angles corresponding to the two rotational axes of each plane.

PS: #AnthropicAI


method inverse()
    float* m = @f + idx * 4
    float* inv = @f + (idx + 25) * 4

    ; Calculate the adjugate matrix
    inv[0] = m[6] * m[12] * m[18] - m[6] * m[13] * m[17] - m[11] * m[7] * m[18] + m[11] * m[8] * m[17] + m[16] * m[7] * m[13] - m[16] * m[8] * m[12]
    inv[1] = -m[1] * m[12] * m[18] + m[1] * m[13] * m[17] + m[11] * m[2] * m[18] - m[11] * m[3] * m[17] - m[16] * m[2] * m[13] + m[16] * m[3] * m[12]
    inv[2] = m[1] * m[7] * m[18] - m[1] * m[8] * m[17] - m[6] * m[2] * m[18] + m[6] * m[3] * m[17] + m[16] * m[2] * m[8] - m[16] * m[3] * m[7]
    inv[3] = -m[1] * m[7] * m[13] + m[1] * m[8] * m[12] + m[6] * m[2] * m[13] - m[6] * m[3] * m[12] - m[11] * m[2] * m[8] + m[11] * m[3] * m[7]
    inv[4] = -m[5] * m[12] * m[18] + m[5] * m[13] * m[17] + m[10] * m[7] * m[18] - m[10] * m[8] * m[17] - m[15] * m[7] * m[13] + m[15] * m[8] * m[12]
    inv[5] = m[0] * m[12] * m[18] - m[0] * m[13] * m[17] - m[10] * m[2] * m[18] + m[10] * m[3] * m[17] + m[15] * m[2] * m[13] - m[15] * m[3] * m[12]
    inv[6] = -m[0] * m[7] * m[18] + m[0] * m[8] * m[17] + m[5] * m[2] * m[18] - m[5] * m[3] * m[17] - m[15] * m[2] * m[8] + m[15] * m[3] * m[7]
    inv[7] = m[0] * m[7] * m[13] - m[0] * m[8] * m[12] - m[5] * m[2] * m[13] + m[5] * m[3] * m[12] + m[10] * m[2] * m[8] - m[10] * m[3] * m[7]
    inv[8] = m[5] * m[11] * m[18] - m[5] * m[13] * m[16] - m[10] * m[6] * m[18] + m[10] * m[8] * m[16] + m[15] * m[6] * m[13] - m[15] * m[8] * m[11]
    inv[9] = -m[0] * m[11] * m[18] + m[0] * m[13] * m[16] + m[10] * m[1] * m[18] - m[10] * m[3] * m[16] - m[15] * m[1] * m[13] + m[15] * m[3] * m[11]
    inv[10] = m[0] * m[6] * m[18] - m[0] * m[8] * m[16] - m[5] * m[1] * m[18] + m[5] * m[3] * m[16] + m[15] * m[1] * m[8] - m[15] * m[3] * m[6]
    inv[11] = -m[0] * m[6] * m[13] + m[0] * m[8] * m[11] + m[5] * m[1] * m[13] - m[5] * m[3] * m[11] - m[10] * m[1] * m[8] + m[10] * m[3] * m[6]
    inv[12] = -m[5] * m[11] * m[17] + m[5] * m[12] * m[16] + m[10] * m[6] * m[17] - m[10] * m[7] * m[16] - m[15] * m[6] * m[12] + m[15] * m[7] * m[11]
    inv[13] = m[0] * m[11] * m[17] - m[0] * m[12] * m[16] - m[10] * m[1] * m[17] + m[10] * m[2] * m[16] + m[15] * m[1] * m[12] - m[15] * m[2] * m[11]
    inv[14] = -m[0] * m[6] * m[17] + m[0] * m[7] * m[16] + m[5] * m[1] * m[17] - m[5] * m[2] * m[16] - m[15] * m[1] * m[7] + m[15] * m[2] * m[6]
    inv[15] = m[0] * m[6] * m[12] - m[0] * m[7] * m[11] - m[5] * m[1] * m[12] + m[5] * m[2] * m[11] + m[10] * m[1] * m[7] - m[10] * m[2] * m[6]
    inv[16] = m[4] * m[12] * m[18] - m[4] * m[13] * m[17] - m[9] * m[7] * m[18] + m[9] * m[8] * m[17] + m[14] * m[7] * m[13] - m[14] * m[8] * m[12]
    inv[17] = -m[0] * m[12] * m[18] + m[0] * m[13] * m[17] + m[9] * m[2] * m[18] - m[9] * m[3] * m[17] - m[14] * m[2] * m[13] + m[14] * m[3] * m[12]
    inv[18] = m[0] * m[7] * m[18] - m[0] * m[8] * m[17] - m[4] * m[2] * m[18] + m[4] * m[3] * m[17] + m[14] * m[2] * m[8] - m[14] * m[3] * m[7]
    inv[19] = -m[0] * m[7] * m[13] + m[0] * m[8] * m[12] + m[4] * m[2] * m[13] - m[4] * m[3] * m[12] - m[9] * m[2] * m[8] + m[9] * m[3] * m[7]
    inv[20] = -m[4] * m[11] * m[18] + m[4] * m[13] * m[16] + m[9] * m[6] * m[18] - m[9] * m[8] * m[16] - m[14] * m[6] * m[13] + m[14] * m[8] * m[11]
    inv[21] = m[0] * m[11] * m[18] - m[0] * m[13] * m[16] - m[9] * m[1] * m[18] + m[9] * m[3] * m[16] + m[14] * m[1] * m[13] - m[14] * m[3] * m[11]
    inv[22] = -m[0] * m[6] * m[18] + m[0] * m[8] * m[16] + m[4] * m[1] * m[18] - m[4] * m[3] * m[16] - m[14] * m[1] * m[8] + m[14] * m[3] * m[6]
    inv[23] = m[0] * m[6] * m[13] - m[0] * m[8] * m[11] - m[4] * m[1] * m[13] + m[4] * m[3] * m[11] + m[9] * m[1] * m[8] - m[9] * m[3] * m[6]
    inv[24] = m[4] * m[11] * m[17] - m[4] * m[12] * m[16] - m[9] * m[6] * m[17] + m[9] * m[7] * m[16] + m[14] * m[6] * m[12] - m[14] * m[7] * m[11]
    inv[25] = -m[0] * m[11] * m[17] + m[0] * m[12] * m[16] + m[9] * m[1] * m[17] - m[9] * m[2] * m[16] - m[14] * m[1] * m[12] + m[14] * m[2] * m[11]
    inv[26] = m[0] * m[6] * m[17] - m[0] * m[7] * m[16] - m[4] * m[1] * m[17] + m[4] * m[2] * m[16] + m[14] * m[1] * m[7] - m[14] * m[2] * m[6]
    inv[27] = -m[0] * m[6] * m[12] + m[0] * m[7] * m[11] + m[4] * m[1] * m[12] - m[4] * m[2] * m[11] - m[9] * m[1] * m[7] + m[9] * m[2] * m[6]

    ; Calculate the determinant
    float det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12] + m[4] * inv[16]

    ; Check if the determinant is zero (matrix is not invertible)
    if det == 0
        print "Error: Matrix is not invertible"
        return
    endif

    ; Multiply the adjugate matrix by the reciprocal of the determinant
    float invDet = 1.0 / det
    for i = 0 to 24
        inv[i] *= invDet
    next

    ; Update the matrix index
    idx += 25
end method


method inverse()
    float* m = @f + idx * 4
    float* inv = @f + (idx + 25) * 4

    ; Create an identity matrix
    for i = 0 to 24
        inv[i] = 0.0
    inv[0] = 1.0
    inv[6] = 1.0
    inv[12] = 1.0
    inv[18] = 1.0
    inv[24] = 1.0

    ; Perform Gaussian elimination
    for i = 0 to 4
        ; Find pivot element
        float pivot = m[i * 5 + i]
        int pivotRow = i

        for j = i + 1 to 4
            if abs(m[j * 5 + i]) > abs(pivot)
                pivot = m[j * 5 + i]
                pivotRow = j

        ; Swap rows if necessary
        if pivotRow != i
            for k = 0 to 4
                float temp = m[i * 5 + k]
                m[i * 5 + k] = m[pivotRow * 5 + k]
                m[pivotRow * 5 + k] = temp

                temp = inv[i * 5 + k]
                inv[i * 5 + k] = inv[pivotRow * 5 + k]
                inv[pivotRow * 5 + k] = temp

        ; Divide row by pivot element
        float invPivot = 1.0 / pivot
        for k = 0 to 4
            m[i * 5 + k] *= invPivot
            inv[i * 5 + k] *= invPivot

        ; Subtract row from other rows to eliminate pivot element
        for j = 0 to 4
            if j != i
                float factor = m[j * 5 + i]
                for k = 0 to 4
                    m[j * 5 + k] -= factor * m[i * 5 + k]
                    inv[j * 5 + k] -= factor * inv[i * 5 + k]

    ; Copy the inverted matrix back to the original matrix
    for i = 0 to 24
        m[i] = inv[i]

    idx += 25 ; Update the matrix index
end method

[quote]This method uses Gaussian elimination to calculate the inverse of the current matrix. Here's how it works:

It starts by creating an identity matrix of the same size as the current matrix. This identity matrix will be transformed into the inverse matrix during the elimination process.
It performs Gaussian elimination on the current matrix and the identity matrix simultaneously. The steps are as follows:
For each row, it finds the pivot element (the element with the largest absolute value in the current column).
If necessary, it swaps rows to bring the pivot element to the diagonal position.
It divides the current row by the pivot element to make the pivot element equal to 1.
It subtracts multiples of the current row from the other rows to eliminate the pivot element from those rows.
After the elimination process, the identity matrix has been transformed into the inverse matrix.
Finally, it copies the inverted matrix back to the original matrix and updates the matrix index.
Note that this method assumes that the matrix is invertible. If the matrix is singular or non-invertible, the method may produce incorrect results or encounter division by zero.

Add to Conversation[/quote]
#8
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Charles Pegge - March 26, 2024, 02:16:31 PM
Once you go beyond 3D, axial rotations don't really exist, so you cant define a w axial rotation. Rotation must be treated as a transform around a fixed point on a plane. Whereas in 3D there are 3 basic planes of rotation: xy xz yz, with corresponding axes z y z, in 4D there are 6 planes:

xy xz xw yz yw zw

Each plane would have 2 rotational axes, for instance the xy plane needs 2 axes z and w, which does not make sense
#9
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Theo Gottwald - March 26, 2024, 01:35:09 PM
These methods set the corresponding transformation matrices for rotation around specific axes, translation, scaling, shearing, perspective projection, and orthographic projection.

Each method follows a similar pattern:

It pushes the current matrix onto the stack using pushmatrix().
It increments the matrix index (idx) by 25 to allocate space for the new matrix.
It sets the matrix to an identity matrix using identity().
It retrieves a pointer to the current matrix using float g at @f + (idx * 4).
It sets the appropriate matrix elements based on the provided parameters.
Finally, it multiplies the current matrix with the newly created matrix using mulmatrix().
The setPerspective() and setOrthographic() methods set up the projection matrices based on the provided field of view, aspect ratio, and clipping plane distances.

Make sure to integrate these methods into the Transform4D class and update the method signatures and variable names to match your existing code style.

method setRotationX(float angle)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    float c = cos(rad(angle))
    float s = sin(rad(angle))
    g[6] = c
    g[7] = s
    g[11] = -s
    g[12] = c
    mulmatrix()
end method

method setRotationY(float angle)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    float c = cos(rad(angle))
    float s = sin(rad(angle))
    g[0] = c
    g[2] = -s
    g[10] = s
    g[12] = c
    mulmatrix()
end method

method setRotationZ(float angle)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    float c = cos(rad(angle))
    float s = sin(rad(angle))
    g[0] = c
    g[1] = s
    g[5] = -s
    g[6] = c
    mulmatrix()
end method

method setRotationW(float angle)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    float c = cos(rad(angle))
    float s = sin(rad(angle))
    g[0] = c
    g[3] = s
    g[15] = -s
    g[18] = c
    mulmatrix()
end method

method setTranslation(float x, y, z, w)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    g[20] = x
    g[21] = y
    g[22] = z
    g[23] = w
    mulmatrix()
end method

method setScale(float x, y, z, w)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    g[0] = x
    g[6] = y
    g[12] = z
    g[18] = w
    mulmatrix()
end method

method setShear(float xy, xz, yz, xw, yw, zw)
    pushmatrix()
    idx += 25
    identity()
    float g at @f + (idx * 4)
    g[1] = xy
    g[2] = xz
    g[7] = yz
    g[3] = xw
    g[8] = yw
    g[13] = zw
    mulmatrix()
end method

method setPerspective(float fov, aspect, nearPlane, farPlane)
    pushmatrix()
    idx += 25
    clear()
    float g at @f + (idx * 4)
    float f = 1.0 / tan(rad(fov) / 2.0)
    g[0] = f / aspect
    g[6] = f
    g[12] = (farPlane + nearPlane) / (nearPlane - farPlane)
    g[13] = (2.0 * farPlane * nearPlane) / (nearPlane - farPlane)
    g[14] = -1.0
    mulmatrix()
end method

method setOrthographic(float left, right, bottom, top, nearPlane, farPlane)
    pushmatrix()
    idx += 25
    clear()
    float g at @f + (idx * 4)
    g[0] = 2.0 / (right - left)
    g[6] = 2.0 / (top - bottom)
    g[12] = -2.0 / (farPlane - nearPlane)
    g[20] = -(right + left) / (right - left)
    g[21] = -(top + bottom) / (top - bottom)
    g[22] = -(farPlane + nearPlane) / (farPlane - nearPlane)
    g[24] = 1.0
    mulmatrix()
end method
#10
OxygenBasic / Re: Tools for 4 Dimensional Ge...
Last post by Theo Gottwald - March 26, 2024, 01:28:52 PM
How about this?

method mulMatrixSIMD(Transform4D other)
    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    ; Assuming SSE instructions are available
    ; Perform matrix multiplication using SIMD instructions
    ; Unroll the loops for better performance

    ; Example using SSE intrinsics (pseudocode)
    for i = 0 to 4
        for j = 0 to 4
            __m128 sum = _mm_setzero_ps()
            for k = 0 to 4
                __m128 a_vec = _mm_load_ps(a + i * 5 + k)
                __m128 b_vec = _mm_load_ps(b + k * 5 + j)
                __m128 prod = _mm_mul_ps(a_vec, b_vec)
                sum = _mm_add_ps(sum, prod)
            _mm_store_ps(r + i * 5 + j, sum)

    idx += 25 ; Update the matrix index
end method

These utility functions provide the following functionality:

transformNormal: Transforms a normal vector by the current matrix. It multiplies the input normal vector by the upper-left 3x3 portion of the matrix, which represents the rotation and scale components. The translation component is ignored since normals are only affected by rotation and scale.

transformPoint: Transforms a point by the current matrix. It multiplies the input point by the entire 4x4 matrix, including the translation component. If the w-component of the transformed point is not zero, it performs perspective division to normalize the point.

decompose: Decomposes the current matrix into its translation, rotation, and scale components. It extracts the translation values from the fourth column of the matrix. It calculates the scale factors by computing the magnitude of each row vector in the upper-left 3x3 portion of the matrix. Finally, it extracts the rotation values by dividing the upper-left 3x3 portion of the matrix by the corresponding scale factors.

compose: Composes a matrix from the given translation, rotation, and scale components. It sets the translation values in the fourth column of the matrix. It multiplies the rotation values by the corresponding scale factors and sets them in the upper-left 3x3 portion of the matrix. Finally, it scales the matrix by the given scale factors.

Note: The code assumes that the input arrays for translation, rotation, and scale have the appropriate sizes (4 elements for translation, 12 elements for rotation, and 4 elements for scale).

Remember to integrate these methods into the Transform4D class and adapt them to match your existing code style and variable names.

method inverse()
    float* m = @f + idx * 4
    float* inv = @f + (idx + 25) * 4

    ; Calculate the adjugate matrix
    inv[0] = m[6] * m[12] * m[18] - m[6] * m[13] * m[17] - m[11] * m[7] * m[18] + m[11] * m[8] * m[17] + m[16] * m[7] * m[13] - m[16] * m[8] * m[12]
    inv[1] = -m[1] * m[12] * m[18] + m[1] * m[13] * m[17] + m[11] * m[2] * m[18] - m[11] * m[3] * m[17] - m[16] * m[2] * m[13] + m[16] * m[3] * m[12]
    inv[2] = m[1] * m[7] * m[18] - m[1] * m[8] * m[17] - m[6] * m[2] * m[18] + m[6] * m[3] * m[17] + m[16] * m[2] * m[8] - m[16] * m[3] * m[7]
    inv[3] = -m[1] * m[7] * m[13] + m[1] * m[8] * m[12] + m[6] * m[2] * m[13] - m[6] * m[3] * m[12] - m[11] * m[2] * m[8] + m[11] * m[3] * m[7]
    inv[4] = m[21] * m[12] * m[18] - m[21] * m[13] * m[17] - m[11] * m[22] * m[18] + m[11] * m[23] * m[17] + m[16] * m[22] * m[13] - m[16] * m[23] * m[12]
    inv[5] = -m[5] * m[12] * m[18] + m[5] * m[13] * m[17] + m[10] * m[7] * m[18] - m[10] * m[8] * m[17] - m[15] * m[7] * m[13] + m[15] * m[8] * m[12]
    inv[6] = m[5] * m[11] * m[18] - m[5] * m[13] * m[16] - m[10] * m[6] * m[18] + m[10] * m[8] * m[16] + m[15] * m[6] * m[13] - m[15] * m[8] * m[11]
    inv[7] = -m[5] * m[11] * m[17] + m[5] * m[12] * m[16] + m[10] * m[6] * m[17] - m[10] * m[7] * m[16] - m[15] * m[6] * m[12] + m[15] * m[7] * m[11]
    inv[8] = -m[21] * m[7] * m[18] + m[21] * m[8] * m[17] + m[6] * m[22] * m[18] - m[6] * m[23] * m[17] - m[16] * m[22] * m[8] + m[16] * m[23] * m[7]
    inv[9] = m[20] * m[7] * m[18] - m[20] * m[8] * m[17] - m[6] * m[21] * m[18] + m[6] * m[23] * m[17] + m[16] * m[21] * m[8] - m[16] * m[23] * m[7]
    inv[10] = -m[20] * m[7] * m[13] + m[20] * m[8] * m[12] + m[1] * m[22] * m[18] - m[1] * m[23] * m[17] - m[16] * m[22] * m[3] + m[16] * m[23] * m[2]
    inv[11] = m[20] * m[7] * m[13] - m[20] * m[8] * m[12] - m[1] * m[21] * m[18] + m[1] * m[23] * m[17] + m[11] * m[21] * m[3] - m[11] * m[23] * m[2]
    inv[12] = m[20] * m[2] * m[13] - m[20] * m[3] * m[12] - m[1] * m[22] * m[13] + m[1] * m[23] * m[12] + m[11] * m[22] * m[3] - m[11] * m[23] * m[2]
    inv[13] = -m[20] * m[2] * m[8] + m[20] * m[3] * m[7] + m[1] * m[22] * m[8] - m[1] * m[23] * m[7] - m[6] * m[22] * m[3] + m[6] * m[23] * m[2]
    inv[14] = m[5] * m[2] * m[18] - m[5] * m[3] * m[17] - m[1] * m[7] * m[18] + m[1] * m[8] * m[17] + m[16] * m[7] * m[3] - m[16] * m[8] * m[2]
    inv[15] = -m[5] * m[2] * m[13] + m[5] * m[3] * m[12] + m[1] * m[7] * m[13] - m[1] * m[8] * m[12] - m[11] * m[7] * m[3] + m[11] * m[8] * m[2]
    inv[16] = -m[20] * m[12] * m[18] + m[20] * m[13] * m[17] + m[10] * m[22] * m[18] - m[10] * m[23] * m[17] - m[15] * m[22] * m[13] + m[15] * m[23] * m[12]
    inv[17] = m[5] * m[12] * m[18] - m[5] * m[13] * m[17] - m[10] * m[2] * m[18] + m[10] * m[3] * m[17] + m[15] * m[2] * m[13] - m[15] * m[3] * m[12]
    inv[18] = -m[5] * m[11] * m[18] + m[5] * m[13] * m[16] + m[10] * m[1] * m[18] - m[10] * m[3] * m[16] - m[15] * m[1] * m[13] + m[15] * m[3] * m[11]
    inv[19] = m[5] * m[11] * m[17] - m[5] * m[12] * m[16] - m[10] * m[1] * m[17] + m[10] * m[2] * m[16] + m[15] * m[1] * m[12] - m[15] * m[2] * m[11]
    inv[20] = m[20] * m[2] * m[18] - m[20] * m[3] * m[17] - m[1] * m[22] * m[18] + m[1] * m[23] * m[17] + m[16] * m[22] * m[3] - m[16] * m[23] * m[2]
    inv[21] = -m[20] * m[2] * m[13] + m[20] * m[3] * m[12] + m[1] * m[22] * m[13] - m[1] * m[23] * m[12] - m[11] * m[22] * m[3] + m[11] * m[23] * m[2]
    inv[22] = m[20] * m[2] * m[8] - m[20] * m[3] * m[7] - m[1] * m[22] * m[8] + m[1] * m[23] * m[7] + m[6] * m[22] * m[3] - m[6] * m[23] * m[2]
    inv[23] = -m[20] * m[2] * m[8] + m[20] * m[3] * m[7] + m[1] * m[21] * m[8] - m[1] * m[23] * m[7] - m[6] * m[21] * m[3] + m[6] * m[23] * m[2]
    inv[24] = -m[5] * m[2] * m[8] + m[5] * m[3] * m[7] + m[1] * m[7] * m[8] - m[1] * m[8] * m[7] - m[6] * m[7] * m[3] + m[6] * m[8] * m[2]

    ; Calculate the determinant
    float det = m[0] * inv[0] + m[1] * inv[5] + m[2] * inv[10] + m[3] * inv[15] + m[4] * inv[20]

    ; Check if the determinant is zero (matrix is not invertible)
    if det == 0
        print "Matrix is not invertible"
        return
    endif

    ; Multiply the adjugate matrix by the reciprocal of the determinant
    float invDet = 1.0 / det
    int i
    for i = 0 to 24
        inv[i] *= invDet
    next

    ; Update the matrix index
    idx += 25
end method

method interpolate(Transform4D other, float t)
    float* a = @f + idx * 4
    float* b = @other.f + other.idx * 4
    float* r = @f + (idx + 25) * 4

    int i
    for i = 0 to 24
        r[i] = a[i] + (b[i] - a[i]) * t
    next

    idx += 25 ; Update the matrix index
end method

method transformDirection(float dirX, float dirY, float dirZ)
    float* m = @f + idx * 4
    float rx, ry, rz

    rx = m[0] * dirX + m[1] * dirY + m[2] * dirZ
    ry = m[5] * dirX + m[6] * dirY + m[7] * dirZ
    rz = m[10] * dirX + m[11] * dirY + m[12] * dirZ

    ; Return the transformed direction vector
    float result[3]
    result[0] = rx
    result[1] = ry
    result[2] = rz

    return result
end method


method transformNormal(float normalX, normalY, normalZ)
    float* m = @f + idx * 4
    float nx = normalX
    float ny = normalY
    float nz = normalZ

    normalX = m[0] * nx + m[1] * ny + m[2] * nz
    normalY = m[5] * nx + m[6] * ny + m[7] * nz
    normalZ = m[10] * nx + m[11] * ny + m[12] * nz

    return normalX, normalY, normalZ
end method

method transformPoint(float pointX, pointY, pointZ)
    float* m = @f + idx * 4
    float px = pointX
    float py = pointY
    float pz = pointZ

    pointX = m[0] * px + m[1] * py + m[2] * pz + m[3]
    pointY = m[5] * px + m[6] * py + m[7] * pz + m[8]
    pointZ = m[10] * px + m[11] * py + m[12] * pz + m[13]
    float pw = m[15] * px + m[16] * py + m[17] * pz + m[18]

    if pw != 0
        pointX /= pw
        pointY /= pw
        pointZ /= pw
    endif

    return pointX, pointY, pointZ
end method

method decompose(float* translation, float* rotation, float* scale)
    float* m = @f + idx * 4

    ; Extract translation
    translation[0] = m[3]
    translation[1] = m[8]
    translation[2] = m[13]
    translation[3] = m[18]

    ; Extract scale
    float sx = sqrt(m[0] * m[0] + m[1] * m[1] + m[2] * m[2])
    float sy = sqrt(m[5] * m[5] + m[6] * m[6] + m[7] * m[7])
    float sz = sqrt(m[10] * m[10] + m[11] * m[11] + m[12] * m[12])
    float sw = sqrt(m[15] * m[15] + m[16] * m[16] + m[17] * m[17])

    scale[0] = sx
    scale[1] = sy
    scale[2] = sz
    scale[3] = sw

    ; Extract rotation
    if sx != 0
        m[0] /= sx
        m[1] /= sx
        m[2] /= sx
    endif
    if sy != 0
        m[5] /= sy
        m[6] /= sy
        m[7] /= sy
    endif
    if sz != 0
        m[10] /= sz
        m[11] /= sz
        m[12] /= sz
    endif
    if sw != 0
        m[15] /= sw
        m[16] /= sw
        m[17] /= sw
    endif

    ; Store rotation values
    rotation[0] = m[0]
    rotation[1] = m[1]
    rotation[2] = m[2]
    rotation[3] = m[5]
    rotation[4] = m[6]
    rotation[5] = m[7]
    rotation[6] = m[10]
    rotation[7] = m[11]
    rotation[8] = m[12]
    rotation[9] = m[15]
    rotation[10] = m[16]
    rotation[11] = m[17]
end method

method compose(float* translation, float* rotation, float* scale)
    float* m = @f + idx * 4

    ; Set translation
    m[3] = translation[0]
    m[8] = translation[1]
    m[13] = translation[2]
    m[18] = translation[3]

    ; Set rotation
    m[0] = rotation[0] * scale[0]
    m[1] = rotation[1] * scale[0]
    m[2] = rotation[2] * scale[0]
    m[5] = rotation[3] * scale[1]
    m[6] = rotation[4] * scale[1]
    m[7] = rotation[5] * scale[1]
    m[10] = rotation[6] * scale[2]
    m[11] = rotation[7] * scale[2]
    m[12] = rotation[8] * scale[2]
    m[15] = rotation[9] * scale[3]
    m[16] = rotation[10] * scale[3]
    m[17] = rotation[11] * scale[3]

    ; Set scale
    m[0] *= scale[0]
    m[6] *= scale[1]
    m[12] *= scale[2]
    m[18] *= scale[3]
end method