Update documentation

docs/algorithm/add.md

@@ -0,0 +1,11 @@
+# Addition
+The addition algorithm operates in Base $2^{64}$, using 64-bit unsigned integers (`uint64_t`) to represent each limb of the arbitrary-precision integer. Overflows beyond the 64-bit boundary are detected using standard unsigned arithmetic comparisons.
+
+1. **Iteration:** The algorithm traverses from the least significant limb (LSB, `index 0`) to the most significant limb (MSB, `index N`).
+2. **Carry Propagation:** For each limb, represented by `i`, it calculates:
+   `sum = aVal + bVal`
+   An overflow is detected if the intermediate `sum` is smaller than `aVal`.
+3. **Carry Handling:**
+   * The carry for the next iteration (`nextCarry`) is initially flagged as `1` if an overflow occurred in the initial sum.
+   * The carry from the previous iteration is then added to the `sum`. If this addition causes a secondary overflow (detected by checking `sum < carry`), `nextCarry` is incremented.
+4. **Finalization:** If a carry remains after the final limb is processed, an additional limb containing the value `1` is appended to the vector. Finally, the vector is normalized to prune trailing zeros.

docs/algorithm/barret-reducer.md

@@ -0,0 +1,62 @@
+# Barrett Reduction
+
+This module implements **Barrett reduction**, a mathematical technique used to speed up modulo operations on large integers.
+
+---
+
+## 1. What is Barrett Reduction?
+
+In computer science, when working with very large integers (arbitrary-precision arithmetic), traditional division is one of the slowest operations. Calculating a simple remainder $x \pmod m$ normally requires a slow long-division step.
+
+Barrett reduction solves this by replacing the slow division with a combination of fast **multiplication** and **digit shifting** (limb-shifting).
+
+### The Core Concept:
+To find the remainder of $x$ divided by $m$, we want to find the quotient $q = \lfloor x / m \rfloor$. Instead of dividing $x$ by $m$, we can multiply $x$ by the "inverse" of $m$ (i.e., $1/m$).
+
+Since we work with integers, we scale this inverse up to make it an integer. Let $b = 2^{64}$ represent our base (since each "limb" or array element in our code is a 64-bit integer), and let $k$ be the number of limbs in the modulus $m$.
+
+1. **Precomputation:**
+   We precalculate a scaling constant, $\mu$ (mu):
+   $$\mu = \lfloor \frac{b^{2k}}{m} \rfloor$$
+   This is the only actual division we need to perform, and it is done once when initializing the reducer.
+
+2. **Estimating the Quotient:**
+   To calculate the modulo, we estimate the quotient $q$. Instead of full multiplication, we use a slightly adjusted formula to keep the numbers from growing too large:
+   $$q \approx \lfloor \frac{\lfloor x / b^{k-1} \rfloor \cdot \mu}{b^{k+1}} \rfloor$$
+   In code, dividing by $b^{k-1}$ and $b^{k+1}$ is done using simple right-shifts of the array elements (limbs).
+
+3. **Calculating the Remainder and Correcting:**
+   Once we have our estimated quotient $q$, we calculate the remainder:
+   $$r = x - (q \cdot m)$$
+   Because we used an approximation, our estimated quotient might be slightly smaller than the actual quotient (by at most 2). Consequently, the temporary remainder $r$ might still be larger than our modulus $m$. We correct this by subtracting $m$ from $r$ until $r < m$ (which is guaranteed to happen in at most two steps).
+
+---
+
+## 2. Estimation and Final Correction
+
+To understand how the mathematical approximation interacts with the final adjustment loop, we can look at how the estimate behaves.
+
+### Visualizing the Process:
+The following diagram illustrates how the algorithm behaves:
+
+![Quotient Estimation and Remainder Progress](img/barrett_reduction.png)
+
+* **Top Graph (Quotient Estimation):** Shows that the estimated quotient (due to rounding down during limb-shifts) occasionally falls slightly below the actual quotient.
+* **Bottom Graph (Remainder Progress):** Shows the preliminary remainder $r$. Whenever the estimated quotient is slightly too small, the temporary remainder $r$ crosses the red line (the modulus $m$). In these instances, the correction loop subtracts $m$ until the remainder is brought within the correct bounds.
+
+---
+
+## 3. Algorithm Flowchart
+
+The step-by-step control flow is represented in the diagram below:
+
+![Barrett Reduction Flowchart](img/barrett_flowchart.png)
+
+### Summary of Steps in Code:
+1. **Shortcut check:** If the input $x$ is already smaller than the modulus $m$, return $x$ immediately.
+2. **First shift:** Shift $x$ right by $k-1$ limbs to get $q_1$.
+3. **Multiply:** Multiply $q_1$ by our precomputed value $\mu$ to get $q_2$.
+4. **Second shift:** Shift $q_2$ right by $k+1$ limbs to get our final estimated quotient, $q_3$.
+5. **Find remainder:** Calculate $r = x - (q_3 \cdot m)$.
+6. **Adjust:** While $r \geq m$, subtract $m$ from $r$ (runs at most twice).
+7. **Return:** Return the final reduced value $r$.

docs/algorithm/div.md

@@ -0,0 +1,166 @@
+# Division
+This documentation explains how the **Bitwise Long Division** (Binary Long Division) is implemented.
+
+## Definitions
+* **Dividend (`data`):** The number being divided.
+* **Divisor (`divisor`):** The number you are dividing by.
+* **Quotient (`quotient`):** This is the result number. How many times the divisor fits entirely into the dividend.
+* **Remainder / Modulo (`remaining`):** The leftover amount that is strictly less than the divisor.
+* **Dividend Mask (`dividendMask`):** The active "working remainder" used dynamically within the long division loop to evaluate the current subset of bits.
+
+## The Basic Concept
+Binary Long Division works exactly like the long division taught in elementary school, but it is simpler because the quotient at any step can only ever be `1` or `0` (it either fits, or it doesn't).
+1. Isolate the highest piece of the dividend.
+2. Check if the divisor fits into it.
+3. If it fits, write `1` to the quotient, subtract the divisor, and "bring down" the next bit.
+4. If it doesn't fit, write `0` to the quotient, don't subtract, and "bring down" the next bit.
+
+**Quick Binary Example `7 / 2`:**
+* `7` (`111`)
+* `2` (`10`)
+
+**The Setup:**
+```text
+             [Quotient]
+           _____________
+Divisor: 10 | 1   1   1  <-- Dividend
+              ^       ^
+             MSB     LSB
+         (Highest) (Lowest)
+```
+
+**Step 1: Isolate the highest bit (MSB)**
+```text
+           0          <-- 10 does not fit into 1. Write 0 to quotient.
+         _______
+     10 |  1 1 1
+         - 0          <-- No subtraction occurs.
+         ---
+           1          <-- Working remainder is 1.
+```
+
+**Step 2: Bring down the next bit**
+```text
+           0 1        <-- 10 fits into 11! Write 1 to quotient.
+         _______
+     10 |  1 1 1
+           . |
+           1 1        <-- Bring down the 1. Working remainder is now 11.
+         - 1 0        <-- Subtract divisor (10).
+         -----
+             1        <-- New working remainder is 1.
+```
+
+**Step 3: Bring down the last bit (LSB)**
+```text
+           0 1 1      <-- 10 fits into 11! Write 1 to quotient.
+         _______
+     10 |  1 1 1
+             . |
+             1 1      <-- Bring down the final 1. Working remainder is 11.
+           - 1 0      <-- Subtract divisor (10).
+           -----
+               1      <-- Final Remainder!
+```
+
+**Result:**
+The quotient is `001` (`3` in decimal). The leftover value at the bottom is our remainder: **`1`**.
+*(7 / 2 = 3 R 1)*.
+
+
+## Implementation
+The division algorithm is implemented within 5 steps.
+
+### 1. Initialization
+1. **Division by Zero:** If the divisor is `0` or the dividend is `0` (yielding an invalid MSB index), we abort and default the quotient and remainder to `0`.
+2. **Finding the MSB:** Using the `getStartBitIndex` function, we get the absolute highest active bit index across the 64-bit limb array and assign it to `initialDividendIndex`. This ensures that we skip processing leading zero-limbs or "ghost" zeros.
+3. **Preallocation:** Because the exact bit-length of the quotient correlates directly to the `initialDividendIndex`, the `quotient` vector is pre-allocated to the maximum required size: `(initialDividendIndex / 64) + 1` 64-bit limbs.
+
+### 2. Working Remainder (`dividendMask`)
+`dividendMask` is initialized by taking a value of `0` and using the `addBitFromNumber` helper to copy down the absolute highest bit from the **dividend** at `initialDividendIndex`.
+
+### 3. Bitwise Evaluation Loop
+In a loop, bit by bit is processed from `initialDividendIndex` down to `-1`.
+For every bit position, the mathematical power index (`currentQBitIndex`) is evaluated.
+
+At each step, we check if the `dividendMask` (working remainder) is $\ge$ the `divisor`.
+
+**Case A: The Mask is $\ge$ the Divisor**
+* The divisor "fits" inside the working remainder.
+* A `1` is written to the exact corresponding bit inside the pre-allocated `quotient` vector using bitwise operations: `quotient[currentQBitIndex / 64] |= (1ULL << (currentQBitIndex % 64))`.
+* If we have reached the end of the dividend (`dividendIndex < 0`), the final modulo subtraction occurs, and the loop terminates.
+* Otherwise, the `divisor` is subtracted from the `dividendMask` in-place, the next bit from the dividend is appended to the mask using `addBitFromNumberInPlace`, and `dividendIndex` is decremented.
+
+**Case B: The Mask is $<$ the Divisor**
+* The divisor does not fit. The quotient bit remains `0` (its default pre-allocated state).
+* No subtraction occurs.
+* If we have reached the end of the dividend (`dividendIndex < 0`), the `dividendMask` is preserved as-is as the remainder, and the loop terminates.
+* Otherwise, the next bit from the dividend is appended to the mask using `addBitFromNumberInPlace` to increase its value for the next loop iteration, and `dividendIndex` is decremented.
+
+### 4. Modulo (The Remainder)
+Because this is integer division, there is often a fractional remainder. The `div` function accepts an optional pointer to a `remaining` vector (`ByteArray* remaining`).
+When the loop finishes processing the final bit (`dividendIndex < 0`), whatever mathematical value is left inside the `dividendMask` represents the modulo. If the pointer is provided, the mask (or subtracted mask) is copied into it.
+
+*(Note: Because of this architecture, evaluating `A % B` requires the same computational effort as `A / B`. Therefore, if both the quotient and remainder are needed, they are extracted simultaneously).*
+
+### 5. Final Normalization
+Even though pre-allocation is tightly bound to the `initialDividendIndex`, the final quotient might have leading zeros depending on the magnitude of the divisor. The `div` function concludes by stripping any trailing zero-limbs from the little-endian vector using `normalizeVector` to maintain strict Base $2^{64}$ normalization guarantees.
+
+---
+
+## Visualization: `13 / 3`
+* **Dividend:** `13` (`1101`)
+* **Divisor:** `3` (`0011`)
+
+**Initialization:**
+* `getStartBitIndex(1101)` returns `3` (the 0-based index of the highest `1` bit).
+* `dividendMask` (the working remainder) is initialized to `0000`.
+* `quotient` is pre-allocated and initialized to `0000`.
+
+### Step 1: Processing Bit Index 3
+* **Action:** Shift `dividendMask` left by 1, and bring down Bit 3 of the dividend (`1`).
+* **Mask Status:** `dividendMask` becomes `0001` (Decimal: 1).
+* **Comparison:** Is `0001` $\ge$ `0011` (Divisor)? $\rightarrow$ **FALSE**
+* **Result:**
+  * Divisor does not fit. No subtraction.
+  * Bit 3 of `quotient` remains `0`.
+  * **Current Quotient:** `0000`
+  * **Current Mask:** `0001`
+
+### Step 2: Processing Bit Index 2
+* **Action:** Shift `dividendMask` left by 1 (`0001` $\rightarrow$ `0010`), and bring down Bit 2 of the dividend (`1`).
+* **Mask Status:** `dividendMask` becomes `0011` (Decimal: 3).
+* **Comparison:** Is `0011` $\ge$ `0011` (Divisor)? $\rightarrow$ **TRUE**
+* **Result:**
+  * Divisor fits!
+  * Set Bit 2 of `quotient` to `1` using bitwise OR.
+  * Subtract divisor from mask: `0011` - `0011` = `0000`.
+  * **Current Quotient:** `0100`
+  * **Current Mask:** `0000`
+
+### Step 3: Processing Bit Index 1
+* **Action:** Shift `dividendMask` left by 1 (`0000` $\rightarrow$ `0000`), and bring down Bit 1 of the dividend (`0`).
+* **Mask Status:** `dividendMask` becomes `0000` (Decimal: 0).
+* **Comparison:** Is `0000` $\ge$ `0011` (Divisor)? $\rightarrow$ **FALSE**
+* **Result:**
+  * Divisor does not fit. No subtraction.
+  * Bit 1 of `quotient` remains `0`.
+  * **Current Quotient:** `0100`
+  * **Current Mask:** `0000`
+
+### Step 4: Processing Bit Index 0 (LSB)
+* **Action:** Shift `dividendMask` left by 1 (`0000` $\rightarrow$ `0000`), and bring down Bit 0 of the dividend (`1`).
+* **Mask Status:** `dividendMask` becomes `0001` (Decimal: 1).
+* **Comparison:** Is `0001` $\ge$ `0011` (Divisor)? $\rightarrow$ **FALSE**
+* **Result:**
+  * Divisor does not fit. No subtraction.
+  * Bit 0 of `quotient` remains `0`.
+  * **Current Quotient:** `0100`
+  * **Current Mask:** `0001`
+
+### Final Output Evaluation:
+The loop terminates because we have processed bit index `0`.
+1. **The Quotient:** The pre-allocated quotient vector holds `0100`, which mathematically evaluates to **`4`**.
+2. **The Remainder:** The `dividendMask` is left holding `0001`, which mathematically evaluates to **`1`**. If a pointer for the modulo was provided, this value is safely copied over.
+
+**Conclusion:** `1101 / 0011` = `0100` with a remainder of `0001` ($13 / 3 = 4\text{ R }1$).

docs/algorithm/mul.md

@@ -0,0 +1,31 @@
+# Multiplication Algorithm
+
+This document outlines the sequential multiplication logic implemented in the `BigInt` class.
+
+---
+
+### 1. Decimal Analogy
+The multiplication process starts with a standard positional approach. Using a decimal analogy (e.g., $23 \times 45$), the operands are broken down into their positional parts (tens and ones), cross-multiplied, scaled by their place value, and summed to produce the final product.
+
+![Dezimales Multiplikationsmodell](img/dezimales_modell.png)
+
+---
+
+### 2. 128-Bit Alignment Layout
+In the binary domain, multiplying two 64-bit inputs (each split into 32-bit halves) results in four 64-bit partial products ($r1$, $r2$, $r3$, $r4$). This layout diagram illustrates how these partial products are shifted and placed within the target 128-bit result space.
+
+![128-Bit-Alignment der Teilprodukte](img/teilprodukte_alignment_128bit.png)
+
+---
+
+### 3. C++ Register Model
+The sequential C++ register model processes these partial products to form a 128-bit product. It splits the calculation into a lower 64-bit register (`fl`) and a higher 64-bit register (`fh`), sequentially accumulating the terms and tracking any intermediate overflows using a running carry.
+
+![C++ Registermodell](img/cpp_registermodell_sequentiell.png)
+
+---
+
+### 4. Result Vector Integration
+Finally, the calculated 128-bit product (`fl` and `fh`) is integrated back into the target array elements (`result`). This step handles addition with existing cell values, tracks two potential overflows, and generates a carry-over value to propagate into the subsequent step.
+
+![Sequentielle Überführung ins ByteArray](img/cpp_ueberfuehrung_bytearray_sequentiell.png)

docs/algorithm/sub.md

@@ -0,0 +1,7 @@
+# Subtraction
+Subtraction implements a "borrow" propagation system to ensure mathematical correctness within an unsigned 64-bit integer space:
+
+1. **Underflow Protection:** Before processing, `isBigger` is checked. If the subtrahend `b` is larger than the minuend `a`, the operation clamps and returns `0` (or assigns `{0}` in-place) because negative integers are not supported.
+2. **Borrow Propagation:** During limb-by-limb iteration, the algorithm calculates `diff = aVal - bVal`. An underflow is flagged (`nextBorrow = true`) if `aVal < bVal`.
+3. **Regrouping:** If an underflow occurs, the algorithm effectively borrows $2^{64}$ (the base value) from the next higher limb (`index i + 1`). When subtracting the previous borrow from the current difference, a secondary underflow check (`diff < borrowVal`) ensures the borrow state is correctly propagated to the next iteration.
+4. **Normalization:** Post-calculation, `normalizeVector` is called to prune trailing zero-limbs, keeping the internal vector size proportional to the actual magnitude of the result.