Repository: unum-cloud/usearch (language: rust)

Issue:
Bug: [Rust] all distances are zero
### Describe the bug

The computed `L2sq` distances seem to be always 0.

### Steps to reproduce

```rust
use usearch::ffi::{IndexOptions, MetricKind};
use usearch::new_index;

fn main() {
    let options = IndexOptions {
        dimensions: 3,
        metric: MetricKind::L2sq,
        ..Default::default()
    };

    let index = new_index(&options).unwrap();
    index.reserve(10).unwrap();
    index.add(0, &[0.4, 0.1, 0.1]).unwrap();
    index.add(1, &[0.5, 0.1, 0.1]).unwrap();
    index.add(2, &[0.6, 0.1, 0.1]).unwrap();
    println!("{:?}", index.search(&[0.05, 0.1, 0.1], 2).unwrap());
}
```
prints
```
Matches { keys: [2, 1], distances: [0.0, 0.0] }
```

### Expected behavior

Expected (and this is what `usearch` 2.9.2, which seems to be the last working version, outputs):
```
Matches { keys: [0, 1], distances: [0.122499995, 0.20249999] }
```

### USearch version

2.12.0

### Operating System

CentOS Stream 9

### Hardware architecture

x86

### Which interface are you using?

Other bindings

### Contact Details

_No response_

### Are you open to being tagged as a contributor?

- [X] I am open to being mentioned in the project `.git` history as a contributor

### Is there an existing issue for this?

- [X] I have searched the existing issues

### Code of Conduct

- [X] I agree to follow this project's Code of Conduct


Discussion / hints:
Interesting! Is it true for both `f32` and `f16` representations? Can you please share the exact model of the CPU printed from `lscpu` and the `index.hardware_acceleration()`. Thanks, @embe-pw 🤗 
`index.hardware_acceleration()` returns `"serial"`.
CPU model is `AMD EPYC 7413 24-Core Processor`.

`quantization: ScalarKind::F32` and `quantization: ScalarKind::F64` seem to work correctly.

I can also reproduce this on an M1 Mac (with `index.hardware_acceleration() == "neon"`

Fix patch:
diff --git a/.vscode/settings.json b/.vscode/settings.json
index 6b8d56f6..3b5d8520 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -128,6 +128,7 @@
     "arange",
     "ashvardanian",
     "astype",
+    "autovec",
     "Availible",
     "bidict",
     "BLAS",
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 6948ec98..adbc0719 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -174,8 +174,8 @@ Following options are enabled:
 - The `-p no:warnings` option will suppress and allow warnings.
 
 ```sh
-pip install pytest pytest-repeat # for repeated fuzzy tests
-pytest # if you trust the default settings
+pip install pytest pytest-repeat            # for repeated fuzzy tests
+pytest                                      # if you trust the default settings
 pytest python/scripts/ -s -x -p no:warnings # to overwrite the default settings
 ```
 
diff --git a/README.md b/README.md
index ab501244..3aed89d2 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 <h1 align="center">USearch</h1>
 <h3 align="center">
 Smaller & <a href="https://www.unum.cloud/blog/2023-11-07-scaling-vector-search-with-intel">Faster</a> Single-File<br/>
-Similarity Search Engine for <a href="https://github.com/ashvardanian/simsimd">Vectors</a> & 🔜 <a href="https://github.com/ashvardanian/stringzilla">Texts</a>
+Similarity Search & Clustering Engine for <a href="https://github.com/ashvardanian/simsimd">Vectors</a> & 🔜 <a href="https://github.com/ashvardanian/stringzilla">Texts</a>
 </h3>
 <br/>
 
@@ -71,13 +71,10 @@ Linux • MacOS • Windows • iOS • WebAssembly •
 
 __Technical Insights__ and related articles:
 
-- [Uses Horner's method for polynomial approximations, beating GCC 12 by 119x](https://ashvardanian.com/posts/gcc-12-vs-avx512fp16/).
 - [Uses Arm SVE and x86 AVX-512's masked loads to eliminate tail `for`-loops](https://ashvardanian.com/posts/simsimd-faster-scipy/#tails-of-the-past-the-significance-of-masked-loads).
-- [Uses AVX-512 FP16 for half-precision operations, that few compilers vectorize](https://ashvardanian.com/posts/simsimd-faster-scipy/#the-challenge-of-f16).
-- [Substitutes LibC's `sqrt` calls with bithacks using Jan Kadlec's constant](https://ashvardanian.com/posts/simsimd-faster-scipy/#bonus-section-bypassing-sqrt-and-libc-dependencies).
+- [Uses Horner's method for polynomial approximations, beating GCC 12 by 119x](https://ashvardanian.com/posts/gcc-12-vs-avx512fp16/).
 - [For every language implements a custom separate binding](https://ashvardanian.com/posts/porting-cpp-library-to-ten-languages/).
-- [For Python avoids slow PyBind11, and even `PyArg_ParseTuple` for speed](https://ashvardanian.com/posts/pybind11-cpython-tutorial/).
-- [For JavaScript uses typed arrays and NAPI for zero-copy calls](https://ashvardanian.com/posts/javascript-ai-vector-search/).
+
 
 ## Comparison with FAISS
 
@@ -119,17 +116,15 @@ USearch is compact and broadly compatible without sacrificing performance, prima
 Base functionality is identical to FAISS, and the interface must be familiar if you have ever investigated Approximate Nearest Neighbors search:
 
 ```py
-# pip install numpy usearch
+# pip install usearch
 
 import numpy as np
 from usearch.index import Index
 
-index = Index(ndim=3)
-
-vector = np.array([0.2, 0.6, 0.4])
-index.add(42, vector)
-
-matches = index.search(vector, 10)
+index = Index(ndim=3)               # Default settings for 3D vectors
+vector = np.array([0.2, 0.6, 0.4])  # Can be a matrix for batch operations
+index.add(42, vector)               # Add one or many vectors in parallel
+matches = index.search(vector, 10)  # Find 10 nearest neighbors
 
 assert matches[0].key == 42
 assert matches[0].distance <= 0.001
@@ -137,13 +132,13 @@ assert np.allclose(index[42], vector, atol=0.1) # Ensure high tolerance in mixed
 ```
 
 More settings are always available, and the API is designed to be as flexible as possible.
-The default storage/quantization level is hardware-dependant for efficiency, but `f16` is recommended for most modern CPUs.
+The default storage/quantization level is hardware-dependant for efficiency, but `bf16` is recommended for most modern CPUs.
 
 ```py
 index = Index(
     ndim=3, # Define the number of dimensions in input vectors
     metric='cos', # Choose 'l2sq', 'ip', 'haversine' or other metric, default = 'cos'
-    dtype='f16', # Store as 'f64', 'f32', 'f16', 'i8', 'b1'..., default = None
+    dtype='bf16', # Store as 'f64', 'f32', 'f16', 'i8', 'b1'..., default = None
     connectivity=16, # Optional: Limit number of neighbors per graph node
     expansion_add=128, # Optional: Control the recall of indexing
     expansion_search=64, # Optional: Control the quality of the search
diff --git a/golang/lib.go b/golang/lib.go
index d54f6502..6b61b1f1 100644
--- a/golang/lib.go
+++ b/golang/lib.go
@@ -63,6 +63,7 @@ type Quantization uint8
 // Different quantization kinds supported by the USearch library.
 const (
 	F32 Quantization = iota
+	BF16
 	F16
 	F64
 	I8
@@ -72,6 +73,8 @@ const (
 // String returns the string representation of the Quantization.
 func (a Quantization) String() string {
 	switch a {
+	case BF16:
+		return "BF16"
 	case F16:
 		return "F16"
 	case F32:
diff --git a/include/usearch/index.hpp b/include/usearch/index.hpp
index e9e940e1..dcf8b52f 100644
--- a/include/usearch/index.hpp
+++ b/include/usearch/index.hpp
@@ -2870,16 +2870,16 @@ class index_gt {
         result.computed_distances = context.computed_distances;
         result.visited_members = context.iteration_cycles;
 
-        // If we are updating the entry node itself, it won't contain any neighbors,
-        // so we should traverse a level down to find the closest match.
-        if (updated_node_level == max_level_copy)
-            updated_node_level--;
-
         // Go down the level, tracking only the closest match;
         // It may even be equal to the `updated_slot`
-        compressed_slot_t closest_slot = search_for_one_( //
-            value, metric, prefetch,                      //
-            entry_slot_copy, max_level_copy, updated_node_level, context);
+        compressed_slot_t closest_slot =
+            // If we are updating the entry node itself, it won't contain any neighbors,
+            // so we should traverse a level down to find the closest match.
+            updated_node_level == max_level_copy //
+                ? entry_slot_copy
+                : search_for_one_(             //
+                      value, metric, prefetch, //
+                      entry_slot_copy, max_level_copy, updated_node_level, context);
 
         // From `updated_node_level` down - perform proper extensive search
         for (level_t level = (std::min)(updated_node_level, max_level_copy); level >= 0; --level) {
@@ -2943,7 +2943,7 @@ class index_gt {
             config.expansion = default_expansion_search();
 
         // Using references is cleaner, but would result in UBSan false positives
-        context_t* context_ptr = contexts_.data() + config.thread;
+        context_t* context_ptr = contexts_.data() ? contexts_.data() + config.thread : nullptr;
         top_candidates_t* top_ptr = context_ptr ? &context_ptr->top_candidates : nullptr;
         search_result_t result{*this, top_ptr};
         if (!nodes_count_.load(std::memory_order_relaxed))
@@ -3697,6 +3697,7 @@ class index_gt {
     inline neighbors_ref_t neighbors_base_(node_t node) const noexcept { return {node.neighbors_tape()}; }
 
     inline neighbors_ref_t neighbors_non_base_(node_t node, level_t level) const noexcept {
+        usearch_assert_m(level > 0 && level <= node.level(), "Linking to missing level");
         return {node.neighbors_tape() + pre_.neighbors_base_bytes + (level - 1) * pre_.neighbors_bytes};
     }
 
diff --git a/include/usearch/index_dense.hpp b/include/usearch/index_dense.hpp
index 9fcd8234..62bef01a 100644
--- a/include/usearch/index_dense.hpp
+++ b/include/usearch/index_dense.hpp
@@ -2240,6 +2240,7 @@ class index_dense_gt {
         case scalar_kind_t::f64_k: return make_casts_<f64_t>();
         case scalar_kind_t::f32_k: return make_casts_<f32_t>();
         case scalar_kind_t::f16_k: return make_casts_<f16_t>();
+        case scalar_kind_t::bf16_k: return make_casts_<bf16_t>();
         case scalar_kind_t::i8_k: return make_casts_<i8_t>();
         case scalar_kind_t::b1x8_k: return make_casts_<b1x8_t>();
         default: return {};
diff --git a/include/usearch/index_plugins.hpp b/include/usearch/index_plugins.hpp
index 68904a8d..84f628bd 100644
--- a/include/usearch/index_plugins.hpp
+++ b/include/usearch/index_plugins.hpp
@@ -70,19 +70,26 @@ struct uuid_t {
 };
 
 class f16_bits_t;
+class bf16_bits_t;
 
 #if !USEARCH_USE_FP16LIB
 #if USEARCH_USE_SIMSIMD
 using f16_native_t = simsimd_f16_t;
+using bf16_native_t = simsimd_bf16_t;
 #elif defined(USEARCH_DEFINED_ARM)
 using f16_native_t = __fp16;
+using bf16_native_t = __fp16; // No better choice on most compilers!
 #else
 using f16_native_t = _Float16;
+using bf16_native_t = _Float16; // No better choice on most compilers!
 #endif
 using f16_t = f16_native_t;
+using bf16_t = bf16_native_t;
 #else
 using f16_native_t = void;
+using bf16_native_t = void;
 using f16_t = f16_bits_t;
+using bf16_t = bf16_bits_t;
 #endif
 
 using f64_t = double;
@@ -132,6 +139,7 @@ enum class scalar_kind_t : std::uint8_t {
     b1x8_k = 1,
     u40_k = 2,
     uuid_k = 3,
+    bf16_k = 4,
     // Common:
     f64_k = 10,
     f32_k = 11,
@@ -167,6 +175,8 @@ template <typename scalar_at> scalar_kind_t scalar_kind() noexcept {
         return scalar_kind_t::f32_k;
     if (std::is_same<scalar_at, f16_t>())
         return scalar_kind_t::f16_k;
+    if (std::is_same<scalar_at, bf16_t>())
+        return scalar_kind_t::bf16_k;
     if (std::is_same<scalar_at, i8_t>())
         return scalar_kind_t::i8_k;
     if (std::is_same<scalar_at, u64_t>())
@@ -208,6 +218,7 @@ inline std::size_t bits_per_scalar(scalar_kind_t scalar_kind) noexcept {
     switch (scalar_kind) {
     case scalar_kind_t::uuid_k: return 128;
     case scalar_kind_t::u40_k: return 40;
+    case scalar_kind_t::bf16_k: return 16;
     case scalar_kind_t::b1x8_k: return 1;
     case scalar_kind_t::u64_k: return 64;
     case scalar_kind_t::i64_k: return 64;
@@ -229,6 +240,7 @@ inline std::size_t bits_per_scalar_word(scalar_kind_t scalar_kind) noexcept {
     switch (scalar_kind) {
     case scalar_kind_t::uuid_k: return 128;
     case scalar_kind_t::u40_k: return 40;
+    case scalar_kind_t::bf16_k: return 16;
     case scalar_kind_t::b1x8_k: return 8;
     case scalar_kind_t::u64_k: return 64;
     case scalar_kind_t::i64_k: return 64;
@@ -250,6 +262,7 @@ inline char const* scalar_kind_name(scalar_kind_t scalar_kind) noexcept {
     switch (scalar_kind) {
     case scalar_kind_t::uuid_k: return "uuid";
     case scalar_kind_t::u40_k: return "u40";
+    case scalar_kind_t::bf16_k: return "bf16";
     case scalar_kind_t::b1x8_k: return "b1x8";
     case scalar_kind_t::u64_k: return "u64";
     case scalar_kind_t::i64_k: return "i64";
@@ -291,10 +304,14 @@ inline expected_gt<scalar_kind_t> scalar_kind_from_name(char const* name, std::s
         parsed.result = scalar_kind_t::f64_k;
     else if (str_equals(name, len, "f16"))
         parsed.result = scalar_kind_t::f16_k;
+    else if (str_equals(name, len, "bf16"))
+        parsed.result = scalar_kind_t::bf16_k;
     else if (str_equals(name, len, "i8"))
         parsed.result = scalar_kind_t::i8_k;
+    else if (str_equals(name, len, "b1"))
+        parsed.result = scalar_kind_t::b1x8_k;
     else
-        parsed.failed("Unknown type, choose: f32, f16, f64, i8");
+        parsed.failed("Unknown type, choose: f64, f32, f16, bf16, i8, b1");
     return parsed;
 }
 
@@ -327,11 +344,13 @@ inline expected_gt<metric_kind_t> metric_from_name(char const* name, std::size_t
                       "tanimoto, sorensen");
     return parsed;
 }
-
 inline expected_gt<metric_kind_t> metric_from_name(char const* name) {
     return metric_from_name(name, std::strlen(name));
 }
 
+/**
+ *  @brief Convenience function to upcast a half-precision floating point number to a single-precision one.
+ */
 inline float f16_to_f32(std::uint16_t u16) noexcept {
 #if !USEARCH_USE_FP16LIB
     f16_native_t f16;
@@ -342,6 +361,9 @@ inline float f16_to_f32(std::uint16_t u16) noexcept {
 #endif
 }
 
+/**
+ *  @brief Convenience function to downcast a single-precision floating point number to a half-precision one.
+ */
 inline std::uint16_t f32_to_f16(float f32) noexcept {
 #if !USEARCH_USE_FP16LIB
     f16_native_t f16 = f16_native_t(f32);
@@ -353,6 +375,36 @@ inline std::uint16_t f32_to_f16(float f32) noexcept {
 #endif
 }
 
+/**
+ *  @brief Convenience function to upcast a brain-floating point number to a single-precision one.
+ *  https://github.com/ashvardanian/SimSIMD/blob/ff51434d90c66f916e94ff05b24530b127aa4cff/include/simsimd/types.h#L395-L410
+ */
+inline float bf16_to_f32(std::uint16_t u16) noexcept {
+    union float_or_unsigned_int_t {
+        float f;
+        unsigned int i;
+    };
+    union float_or_unsigned_int_t result_union;
+    result_union.i = u16 << 16; // Zero extends the mantissa
+    return result_union.f;
+}
+
+/**
+ *  @brief Convenience function to downcast a single-precision floating point number to a brain-floating point one.
+ *  https://github.com/ashvardanian/SimSIMD/blob/ff51434d90c66f916e94ff05b24530b127aa4cff/include/simsimd/types.h#L412-L425
+ */
+inline std::uint16_t f32_to_bf16(float f32) noexcept {
+    union float_or_unsigned_int_t {
+        float f;
+        unsigned int i;
+    };
+    union float_or_unsigned_int_t value;
+    value.f = f32;
+    value.i >>= 16;
+    value.i &= 0xFFFF;
+    return (unsigned short)value.i;
+}
+
 /**
  *  @brief  Numeric type for the IEEE 754 half-precision floating point.
  *          If hardware support isn't available, falls back to a hardware
@@ -412,6 +464,65 @@ class f16_bits_t {
     }
 };
 
+/**
+ *  @brief  Numeric type for brain-floating point half-precision floating point.
+ *          If hardware support isn't available, falls back to a hardware
+ *          agnostic in-software implementation.
+ */
+class bf16_bits_t {
+    std::uint16_t uint16_{};
+
+  public:
+    inline bf16_bits_t() noexcept : uint16_(0) {}
+    inline bf16_bits_t(bf16_bits_t&&) = default;
+    inline bf16_bits_t& operator=(bf16_bits_t&&) = default;
+    inline bf16_bits_t(bf16_bits_t const&) = default;
+    inline bf16_bits_t& operator=(bf16_bits_t const&) = default;
+
+    inline operator float() const noexcept { return bf16_to_f32(uint16_); }
+    inline explicit operator bool() const noexcept { return bf16_to_f32(uint16_) > 0.5f; }
+
+    inline bf16_bits_t(std::int8_t v) noexcept : uint16_(f32_to_bf16(v)) {}
+    inline bf16_bits_t(bool v) noexcept : uint16_(f32_to_bf16(v)) {}
+    inline bf16_bits_t(float v) noexcept : uint16_(f32_to_bf16(v)) {}
+    inline bf16_bits_t(double v) noexcept : uint16_(f32_to_bf16(static_cast<float>(v))) {}
+
+    inline bool operator<(bf16_bits_t const& other) const noexcept { return float(*this) < float(other); }
+
+    inline bf16_bits_t operator+(bf16_bits_t other) const noexcept { return {float(*this) + float(other)}; }
+    inline bf16_bits_t operator-(bf16_bits_t other) const noexcept { return {float(*this) - float(other)}; }
+    inline bf16_bits_t operator*(bf16_bits_t other) const noexcept { return {float(*this) * float(other)}; }
+    inline bf16_bits_t operator/(bf16_bits_t other) const noexcept { return {float(*this) / float(other)}; }
+    inline float operator+(float other) const noexcept { return float(*this) + other; }
+    inline float operator-(float other) const noexcept { return float(*this) - other; }
+    inline float operator*(float other) const noexcept { return float(*this) * other; }
+    inline float operator/(float other) const noexcept { return float(*this) / other; }
+    inline double operator+(double other) const noexcept { return float(*this) + other; }
+    inline double operator-(double other) const noexcept { return float(*this) - other; }
+    inline double operator*(double other) const noexcept { return float(*this) * other; }
+    inline double operator/(double other) const noexcept { return float(*this) / other; }
+
+    inline bf16_bits_t& operator+=(float v) noexcept {
+        uint16_ = f32_to_bf16(v + bf16_to_f32(uint16_));
+        return *this;
+    }
+
+    inline bf16_bits_t& operator-=(float v) noexcept {
+        uint16_ = f32_to_bf16(v - bf16_to_f32(uint16_));
+        return *this;
+    }
+
+    inline bf16_bits_t& operator*=(float v) noexcept {
+        uint16_ = f32_to_bf16(v * bf16_to_f32(uint16_));
+        return *this;
+    }
+
+    inline bf16_bits_t& operator/=(float v) noexcept {
+        uint16_ = f32_to_bf16(v / bf16_to_f32(uint16_));
+        return *this;
+    }
+};
+
 /**
  *  @brief  An STL-based executor or a "thread-pool" for parallel execution.
  *          Isn't efficient for small batches, as it recreates the threads on every call.
@@ -1009,18 +1120,22 @@ template <typename to_scalar_at> struct cast_from_i8_gt {
 };
 
 template <> struct cast_gt<i8_t, f16_t> : public cast_from_i8_gt<f16_t> {};
+template <> struct cast_gt<i8_t, bf16_t> : public cast_from_i8_gt<bf16_t> {};
 template <> struct cast_gt<i8_t, f32_t> : public cast_from_i8_gt<f32_t> {};
 template <> struct cast_gt<i8_t, f64_t> : public cast_from_i8_gt<f64_t> {};
 
 template <> struct cast_gt<f16_t, i8_t> : public cast_to_i8_gt<f16_t> {};
+template <> struct cast_gt<bf16_t, i8_t> : public cast_to_i8_gt<bf16_t> {};
 template <> struct cast_gt<f32_t, i8_t> : public cast_to_i8_gt<f32_t> {};
 template <> struct cast_gt<f64_t, i8_t> : public cast_to_i8_gt<f64_t> {};
 
 template <> struct cast_gt<b1x8_t, f16_t> : public cast_from_b1x8_gt<f16_t> {};
+template <> struct cast_gt<b1x8_t, bf16_t> : public cast_from_b1x8_gt<bf16_t> {};
 template <> struct cast_gt<b1x8_t, f32_t> : public cast_from_b1x8_gt<f32_t> {};
 template <> struct cast_gt<b1x8_t, f64_t> : public cast_from_b1x8_gt<f64_t> {};
 
 template <> struct cast_gt<f16_t, b1x8_t> : public cast_to_b1x8_gt<f16_t> {};
+template <> struct cast_gt<bf16_t, b1x8_t> : public cast_to_b1x8_gt<bf16_t> {};
 template <> struct cast_gt<f32_t, b1x8_t> : public cast_to_b1x8_gt<f32_t> {};
 template <> struct cast_gt<f64_t, b1x8_t> : public cast_to_b1x8_gt<f64_t> {};
 
@@ -1602,6 +1717,7 @@ class metric_punned_t {
         case scalar_kind_t::f32_k: datatype = simsimd_datatype_f32_k; break;
         case scalar_kind_t::f64_k: datatype = simsimd_datatype_f64_k; break;
         case scalar_kind_t::f16_k: datatype = simsimd_datatype_f16_k; break;
+        case scalar_kind_t::bf16_k: datatype = simsimd_datatype_bf16_k; break;
         case scalar_kind_t::i8_k: datatype = simsimd_datatype_i8_k; break;
         case scalar_kind_t::b1x8_k: datatype = simsimd_datatype_b8_k; break;
         default: break;
@@ -1654,9 +1770,10 @@ class metric_punned_t {
         switch (metric_kind_) {
         case metric_kind_t::ip_k: {
             switch (scalar_kind_) {
-            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<f32_t>>; break;
-            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::bf16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<bf16_t, f32_t>>; break;
             case scalar_kind_t::i8_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<i8_t, f32_t>>; break;
+            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<f32_t>>; break;
             case scalar_kind_t::f64_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_ip_gt<f64_t>>; break;
             default: metric_ptr_ = 0; break;
             }
@@ -1664,9 +1781,10 @@ class metric_punned_t {
         }
         case metric_kind_t::cos_k: {
             switch (scalar_kind_) {
-            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<f32_t>>; break;
-            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::bf16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<bf16_t, f32_t>>; break;
             case scalar_kind_t::i8_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<i8_t, f32_t>>; break;
+            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<f32_t>>; break;
             case scalar_kind_t::f64_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_cos_gt<f64_t>>; break;
             default: metric_ptr_ = 0; break;
             }
@@ -1674,9 +1792,10 @@ class metric_punned_t {
         }
         case metric_kind_t::l2sq_k: {
             switch (scalar_kind_) {
-            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<f32_t>>; break;
-            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::bf16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<bf16_t, f32_t>>; break;
             case scalar_kind_t::i8_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<i8_t, f32_t>>; break;
+            case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<f16_t, f32_t>>; break;
+            case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<f32_t>>; break;
             case scalar_kind_t::f64_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_l2sq_gt<f64_t>>; break;
             default: metric_ptr_ = 0; break;
             }
@@ -1684,6 +1803,9 @@ class metric_punned_t {
         }
         case metric_kind_t::pearson_k: {
             switch (scalar_kind_) {
+            case scalar_kind_t::bf16_k:
+                metric_ptr_ = (uptr_t)&equidimensional_<metric_pearson_gt<bf16_t, f32_t>>;
+                break;
             case scalar_kind_t::i8_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_pearson_gt<i8_t, f32_t>>; break;
             case scalar_kind_t::f16_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_pearson_gt<f16_t, f32_t>>; break;
             case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_pearson_gt<f32_t>>; break;
@@ -1694,7 +1816,8 @@ class metric_punned_t {
         }
         case metric_kind_t::haversine_k: {
             switch (scalar_kind_) {
-            case scalar_kind_t::f16_k: metric_ptr_ = 0; break; //< Having half-precision 2D coordinates is a bit silly.
+            case scalar_kind_t::bf16_k: metric_ptr_ = 0; break; //< Half-precision 2D vectors are silly.
+            case scalar_kind_t::f16_k: metric_ptr_ = 0; break;  //< Half-precision 2D vectors are silly.
             case scalar_kind_t::f32_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_haversine_gt<f32_t>>; break;
             case scalar_kind_t::f64_k: metric_ptr_ = (uptr_t)&equidimensional_<metric_haversine_gt<f64_t>>; break;
             default: metric_ptr_ = 0; break;
@@ -1703,6 +1826,9 @@ class metric_punned_t {
         }
         case metric_kind_t::divergence_k: {
             switch (scalar_kind_) {
+            case scalar_kind_t::bf16_k:
+                metric_ptr_ = (uptr_t)&equidimensional_<metric_divergence_gt<bf16_t, f32_t>>;
+                break;
             case scalar_kind_t::f16_k:
                 metric_ptr_ = (uptr_t)&equidimensional_<metric_divergence_gt<f16_t, f32_t>>;
                 break;
diff --git a/javascript/usearch.ts b/javascript/usearch.ts
index 89e8f073..a04ca1aa 100644
--- a/javascript/usearch.ts
+++ b/javascript/usearch.ts
@@ -70,6 +70,7 @@ export enum ScalarKind {
   F32 = "f32",
   F64 = "f64",
   F16 = "f16",
+  BF16 = "bf16",
   I8 = "i8",
   B1 = "b1",
 }
@@ -610,12 +611,12 @@ function exactSearch(
 }
 
 const usearch = {
-    Index,
-    MetricKind,
-    ScalarKind,
-    Matches,
-    BatchMatches,
-    exactSearch,
+  Index,
+  MetricKind,
+  ScalarKind,
+  Matches,
+  BatchMatches,
+  exactSearch,
 };
 export default usearch;
 
diff --git a/python/lib.cpp b/python/lib.cpp
index 112a1a97..fa1665da 100644
--- a/python/lib.cpp
+++ b/python/lib.cpp
@@ -905,6 +905,7 @@ PYBIND11_MODULE(compiled, m) {
         .value("B1", scalar_kind_t::b1x8_k)
         .value("U40", scalar_kind_t::u40_k)
         .value("UUID", scalar_kind_t::uuid_k)
+        .value("BF16", scalar_kind_t::bf16_k)
         .value("F64", scalar_kind_t::f64_k)
         .value("F32", scalar_kind_t::f32_k)
         .value("F16", scalar_kind_t::f16_k)
diff --git a/python/scripts/bench.py b/python/scripts/bench.py
index cac7a8e7..976c518f 100644
--- a/python/scripts/bench.py
+++ b/python/scripts/bench.py
@@ -26,7 +26,7 @@ def bench_speed(
     # Build various indexes:
     indexes = []
     jit_options = [False, True] if jit else [False]
-    dtype_options = [ScalarKind.F32, ScalarKind.F16, ScalarKind.I8]
+    dtype_options = [ScalarKind.F32, ScalarKind.F16, ScalarKind.BF16, ScalarKind.I8]
     for jit, dtype in itertools.product(jit_options, dtype_options):
         metric = MetricKind.IP
         if jit:
diff --git a/python/scripts/bench_exact.py b/python/scripts/bench_exact.py
index c1062d11..87d7d2fa 100644
--- a/python/scripts/bench_exact.py
+++ b/python/scripts/bench_exact.py
@@ -34,7 +34,7 @@ def run(
     q: int = 10,
     k: int = 100,
     ndim: int = 256,
-    dtype: Literal["b1", "i8", "f16", "f32", "f64"] = "f32",
+    dtype: Literal["b1", "i8", "f16", "bf16", "f32", "f64"] = "f32",
     metric: Literal["cos", "ip", "l2sq"] = "ip",
 ):
 
@@ -62,7 +62,7 @@ def run(
 
     if metric not in [MetricKind.L2sq, MetricKind.IP]:
         return
-    if dtype not in [ScalarKind.I8, ScalarKind.F16, ScalarKind.F32, ScalarKind.F64]:
+    if dtype not in [ScalarKind.I8, ScalarKind.F16, ScalarKind.BF16, ScalarKind.F32, ScalarKind.F64]:
         return
 
     start = time()
diff --git a/python/usearch/eval.py b/python/usearch/eval.py
index 09b803d9..4baee713 100644
--- a/python/usearch/eval.py
+++ b/python/usearch/eval.py
@@ -124,7 +124,7 @@ def self_recall(index: Index, sample: Union[float, int] = 1.0, **kwargs) -> Sear
     if "vectors" in kwargs:
         vectors = kwargs.pop("vectors")
     else:
-        vectors = index.get(keys, index.dtype)
+        vectors = index.get(keys)
 
     matches = index.search(vectors, **kwargs)
     count_matches: int = (
diff --git a/python/usearch/index.py b/python/usearch/index.py
index 2323501b..aeda1771 100644
--- a/python/usearch/index.py
+++ b/python/usearch/index.py
@@ -93,6 +93,8 @@ def _normalize_dtype(
     if dtype is None or dtype == "":
         if metric in MetricKindBitwise:
             return ScalarKind.B1
+        if _hardware_acceleration(dtype=ScalarKind.BF16, ndim=ndim, metric_kind=metric):
+            return ScalarKind.BF16
         if _hardware_acceleration(dtype=ScalarKind.F16, ndim=ndim, metric_kind=metric):
             return ScalarKind.F16
         return ScalarKind.F32
@@ -106,12 +108,14 @@ def _normalize_dtype(
     _normalize = {
         "f64": ScalarKind.F64,
         "f32": ScalarKind.F32,
+        "bf16": ScalarKind.BF16,
         "f16": ScalarKind.F16,
         "i8": ScalarKind.I8,
         "b1": ScalarKind.B1,
         "b1x8": ScalarKind.B1,
         "float64": ScalarKind.F64,
         "float32": ScalarKind.F32,
+        "bfloat16": ScalarKind.BF16,
         "float16": ScalarKind.F16,
         "int8": ScalarKind.I8,
         np.float64: ScalarKind.F64,
@@ -124,6 +128,8 @@ def _normalize_dtype(
 
 
 def _to_numpy_dtype(dtype: ScalarKind):
+    if dtype == ScalarKind.BF16:
+        return None
     _normalize = {
         ScalarKind.F64: np.float64,
         ScalarKind.F32: np.float32,
@@ -734,10 +740,15 @@ def get(
         """
         if not dtype:
             dtype = self.dtype
+            view_dtype = _to_numpy_dtype(dtype)
+            if view_dtype is None:
+                dtype = ScalarKind.F32
+                view_dtype = np.float32
         else:
             dtype = _normalize_dtype(dtype)
-
-        view_dtype = _to_numpy_dtype(dtype)
+            view_dtype = _to_numpy_dtype(dtype)
+            if view_dtype is None:
+                raise NotImplementedError("The requested representation type is not supported by NumPy")
 
         def cast(result):
             if result is not None:
diff --git a/rust/README.md b/rust/README.md
index 3ae58294..56690e2a 100644
--- a/rust/README.md
+++ b/rust/README.md
@@ -42,8 +42,8 @@ use usearch::{Index, IndexOptions, MetricKind, ScalarKind, new_index};
 
 let options = IndexOptions {
     dimensions: 3, // necessary for most metric kinds
-    metric: MetricKind::IP, // or MetricKind::L2sq, MetricKind::Cos ...
-    quantization: ScalarKind::F16, // or ScalarKind::F32, ScalarKind::I8, ScalarKind::B1x8 ...
+    metric: MetricKind::IP, // or ::L2sq, ::Cos ...
+    quantization: ScalarKind::BF16, // or ::F32, ::F16, ::I8, ::B1x8 ...
     connectivity: 0, // zero for auto
     expansion_add: 0, // zero for auto
     expansion_search: 0, // zero for auto
diff --git a/rust/lib.cpp b/rust/lib.cpp
index b480a596..87ee1d2c 100644
--- a/rust/lib.cpp
+++ b/rust/lib.cpp
@@ -28,6 +28,7 @@ metric_kind_t rust_to_cpp_metric(MetricKind value) {
 scalar_kind_t rust_to_cpp_scalar(ScalarKind value) {
     switch (value) {
     case ScalarKind::I8: return scalar_kind_t::i8_k;
+    case ScalarKind::BF16: return scalar_kind_t::bf16_k;
     case ScalarKind::F16: return scalar_kind_t::f16_k;
     case ScalarKind::F32: return scalar_kind_t::f32_k;
     case ScalarKind::F64: return scalar_kind_t::f64_k;
diff --git a/rust/lib.rs b/rust/lib.rs
index 9637e710..603d3bc7 100644
--- a/rust/lib.rs
+++ b/rust/lib.rs
@@ -280,6 +280,8 @@ pub mod ffi {
         F32,
         /// 16-bit half-precision IEEE 754 floating-point number (different from `bf16`).
         F16,
+        /// 16-bit brain floating-point number.
+        BF16,
         /// 8-bit signed integer.
         I8,
         /// 1-bit binary value, packed 8 per byte.
@@ -515,7 +517,7 @@ impl Default for ffi::IndexOptions {
         Self {
             dimensions: 256,
             metric: MetricKind::Cos,
-            quantization: ScalarKind::F16,
+            quantization: ScalarKind::BF16,
             connectivity: 0,
             expansion_add: 0,
             expansion_search: 0,


Test patch:
diff --git a/python/scripts/test_distances.py b/python/scripts/test_distances.py
index 5569739d..ec637189 100644
--- a/python/scripts/test_distances.py
+++ b/python/scripts/test_distances.py
@@ -25,6 +25,7 @@
     [
         ScalarKind.F32,
         ScalarKind.F16,
+        ScalarKind.BF16,
         ScalarKind.I8,
     ],
 )
diff --git a/python/scripts/test_index.py b/python/scripts/test_index.py
index 43987481..e3a99592 100644
--- a/python/scripts/test_index.py
+++ b/python/scripts/test_index.py
@@ -27,6 +27,7 @@
     ScalarKind.F32,
     ScalarKind.F64,
     ScalarKind.F16,
+    ScalarKind.BF16,
     ScalarKind.I8,
 ]
 dtypes = [np.float32, np.float64, np.float16]