Janady13/complete-agi-core / mind_farm/include/fast_math.hpp

1/** 2 * @file fast_math.hpp 3 * @brief Ultra-Fast Mathematical Primitives for Gate5000 4 * 5 * ═══════════════════════════════════════════════════════════════════════════════ 6 * OPTIMIZATION TECHNIQUES: 7 * ═══════════════════════════════════════════════════════════════════════════════ 8 * 9 * 1. FAST INVERSE SQUARE ROOT (Quake III style) 10 * - 4x faster than 1/sqrt(x) 11 * - Error < 1% 12 * 13 * 2. LOOKUP TABLES (LUT) 14 * - 256-entry sine/cosine tables 15 * - 64-entry exp/log tables 16 * - Linear interpolation for precision 17 * 18 * 3. POLYNOMIAL APPROXIMATIONS 19 * - Minimax polynomials for transcendentals 20 * - Padé approximants for rational functions 21 * 22 * 4. BIT MANIPULATION 23 * - Fast floor/ceil via bit ops 24 * - Fast abs/sign extraction 25 * - Fast modulo for powers of 2 26 * 27 * 5. SIMD-FRIENDLY LAYOUTS 28 * - 4-wide vector operations 29 * - Cache-aligned data structures 30 * 31 * ═══════════════════════════════════════════════════════════════════════════════ 32 */ 33 34#pragma once 35 36#include <cstdint> 37#include <cstring> 38#include <array> 39#include "constants.hpp" 40 41namespace mind_farm { 42namespace fast { 43 44// ═══════════════════════════════════════════════════════════════════════════════ 45// I. FAST INVERSE SQUARE ROOT 46// ═══════════════════════════════════════════════════════════════════════════════ 47 48/** 49 * @brief Fast inverse square root (Quake III algorithm) 50 * @param x Input value (must be positive) 51 * @return Approximation of 1/sqrt(x), error < 1% 52 */ 53inline float fast_inv_sqrt(float x) { 54 float xhalf = 0.5f * x; 55 int32_t i; 56 std::memcpy(&i, &x, sizeof(float)); 57 i = 0x5f3759df - (i >> 1); // Magic constant 58 std::memcpy(&x, &i, sizeof(float)); 59 x = x * (1.5f - xhalf * x * x); // Newton iteration 60 return x; 61} 62 63/** 64 * @brief Fast square root using inverse sqrt 65 */ 66inline float fast_sqrt(float x) { 67 if (x <= 0.0f) return 0.0f; 68 return x * fast_inv_sqrt(x); 69} 70 71/** 72 * @brief Fast reciprocal 1/x 73 */ 74inline float fast_rcp(float x) { 75 int32_t i; 76 std::memcpy(&i, &x, sizeof(float)); 77 i = 0x7ef0e840 - i; // Magic constant for reciprocal 78 float y; 79 std::memcpy(&y, &i, sizeof(float)); 80 y = y * (2.0f - x * y); // Newton iteration 81 return y; 82} 83 84// ═══════════════════════════════════════════════════════════════════════════════ 85// II. LOOKUP TABLES 86// ═══════════════════════════════════════════════════════════════════════════════ 87 88// Precomputed sine table (256 entries for [0, 2π]) 89alignas(64) constexpr std::array<float, 256> SIN_LUT = []() { 90 std::array<float, 256> lut{}; 91 for (int i = 0; i < 256; ++i) { 92 double angle = i * 2.0 * PI / 256.0; 93 lut[i] = static_cast<float>(std::sin(angle)); 94 } 95 return lut; 96}(); 97 98// Precomputed exp table (64 entries for [-4, 4]) 99alignas(64) constexpr std::array<float, 64> EXP_LUT = []() { 100 std::array<float, 64> lut{}; 101 for (int i = 0; i < 64; ++i) { 102 double x = (i - 32) * 0.125; // Range [-4, 4] 103 lut[i] = static_cast<float>(std::exp(x)); 104 } 105 return lut; 106}(); 107 108// Precomputed tanh table (64 entries for [-4, 4]) 109alignas(64) constexpr std::array<float, 64> TANH_LUT = []() { 110 std::array<float, 64> lut{}; 111 for (int i = 0; i < 64; ++i) { 112 double x = (i - 32) * 0.125; // Range [-4, 4] 113 lut[i] = static_cast<float>(std::tanh(x)); 114 } 115 return lut; 116}(); 117 118// Precomputed sigmoid table (64 entries for [-8, 8]) 119alignas(64) constexpr std::array<float, 64> SIGMOID_LUT = []() { 120 std::array<float, 64> lut{}; 121 for (int i = 0; i < 64; ++i) { 122 double x = (i - 32) * 0.25; // Range [-8, 8] 123 lut[i] = static_cast<float>(1.0 / (1.0 + std::exp(-x))); 124 } 125 return lut; 126}(); 127 128/** 129 * @brief Fast sine using LUT with linear interpolation 130 */ 131inline float fast_sin(float x) { 132 // Normalize to [0, 2π] 133 x = x - static_cast<float>(2.0 * PI) * std::floor(x / static_cast<float>(2.0 * PI)); 134 float idx_f = x * (256.0f / static_cast<float>(2.0 * PI)); 135 int idx = static_cast<int>(idx_f) & 255; 136 float frac = idx_f - std::floor(idx_f); 137 return SIN_LUT[idx] * (1.0f - frac) + SIN_LUT[(idx + 1) & 255] * frac; 138} 139 140/** 141 * @brief Fast cosine using LUT 142 */ 143inline float fast_cos(float x) { 144 return fast_sin(x + static_cast<float>(PI * 0.5)); 145} 146 147/** 148 * @brief Fast tanh using LUT with linear interpolation 149 */ 150inline float fast_tanh(float x) { 151 if (x < -4.0f) return -1.0f; 152 if (x > 4.0f) return 1.0f; 153 float idx_f = (x + 4.0f) * 8.0f; // Map [-4,4] to [0,64] 154 int idx = static_cast<int>(idx_f); 155 if (idx >= 63) return TANH_LUT[63]; 156 float frac = idx_f - idx; 157 return TANH_LUT[idx] * (1.0f - frac) + TANH_LUT[idx + 1] * frac; 158} 159 160/** 161 * @brief Fast sigmoid using LUT 162 */ 163inline float fast_sigmoid(float x) { 164 if (x < -8.0f) return 0.0f; 165 if (x > 8.0f) return 1.0f; 166 float idx_f = (x + 8.0f) * 4.0f; // Map [-8,8] to [0,64] 167 int idx = static_cast<int>(idx_f); 168 if (idx >= 63) return SIGMOID_LUT[63]; 169 float frac = idx_f - idx; 170 return SIGMOID_LUT[idx] * (1.0f - frac) + SIGMOID_LUT[idx + 1] * frac; 171} 172 173// ═══════════════════════════════════════════════════════════════════════════════ 174// III. POLYNOMIAL APPROXIMATIONS 175// ═══════════════════════════════════════════════════════════════════════════════ 176 177/** 178 * @brief Fast exp using polynomial approximation 179 * Valid for x ∈ [-4, 4], error < 0.1% 180 */ 181inline float fast_exp(float x) { 182 if (x < -10.0f) return 0.0f; 183 if (x > 10.0f) return 22026.0f; 184 185 // Reduce to [0, ln2] range 186 float k = std::floor(x * 1.4426950408889634f); // x / ln(2) 187 float r = x - k * 0.6931471805599453f; // x - k * ln(2) 188 189 // Minimax polynomial for exp(r) on [0, ln2] 190 float p = 1.0f + r * (1.0f + r * (0.5f + r * (0.166666667f + r * 0.041666667f))); 191 192 // Reconstruct: exp(x) = 2^k * exp(r) 193 int32_t i = static_cast<int32_t>(k); 194 int32_t exp_bits = (127 + i) << 23; 195 float scale; 196 std::memcpy(&scale, &exp_bits, sizeof(float)); 197 return p * scale; 198} 199 200/** 201 * @brief Fast log using polynomial approximation 202 */ 203inline float fast_log(float x) { 204 if (x <= 0.0f) return -1e10f; 205 206 // Extract exponent and mantissa 207 int32_t i; 208 std::memcpy(&i, &x, sizeof(float)); 209 int32_t e = ((i >> 23) & 255) - 127; 210 i = (i & 0x007fffff) | 0x3f800000; 211 float m; 212 std::memcpy(&m, &i, sizeof(float)); 213 214 // Polynomial approximation for log(m) where m ∈ [1, 2) 215 m = m - 1.0f; 216 float p = m * (1.0f - m * (0.5f - m * (0.333333f - m * 0.25f))); 217 218 return p + e * 0.6931471805599453f; 219} 220 221/** 222 * @brief Fast pow using exp(y * log(x)) 223 */ 224inline float fast_pow(float x, float y) { 225 return fast_exp(y * fast_log(x)); 226} 227 228// ═══════════════════════════════════════════════════════════════════════════════ 229// IV. BIT MANIPULATION TRICKS 230// ═══════════════════════════════════════════════════════════════════════════════ 231 232/** 233 * @brief Fast absolute value (branchless) 234 */ 235inline float fast_abs(float x) { 236 int32_t i; 237 std::memcpy(&i, &x, sizeof(float)); 238 i &= 0x7fffffff; 239 std::memcpy(&x, &i, sizeof(float)); 240 return x; 241} 242 243/** 244 * @brief Fast sign extraction (returns -1, 0, or 1) 245 */ 246inline int fast_sign(float x) { 247 return (x > 0.0f) - (x < 0.0f); 248} 249 250/** 251 * @brief Fast floor using bit manipulation 252 */ 253inline int fast_floor(float x) { 254 int i = static_cast<int>(x); 255 return i - (x < i); 256} 257 258/** 259 * @brief Fast ceil using bit manipulation 260 */ 261inline int fast_ceil(float x) { 262 int i = static_cast<int>(x); 263 return i + (x > i); 264} 265 266/** 267 * @brief Fast modulo for power-of-2 divisors 268 */ 269inline uint32_t fast_mod_pow2(uint32_t x, uint32_t pow2) { 270 return x & (pow2 - 1); 271} 272 273/** 274 * @brief Fast clamp to [0, 1] (branchless) 275 */ 276inline float fast_clamp01(float x) { 277 x = x > 0.0f ? x : 0.0f; 278 return x < 1.0f ? x : 1.0f; 279} 280 281/** 282 * @brief Fast linear interpolation 283 */ 284inline float fast_lerp(float a, float b, float t) { 285 return a + t * (b - a); 286} 287 288// ═══════════════════════════════════════════════════════════════════════════════ 289// V. FAST VECTOR OPERATIONS (4-wide) 290// ═══════════════════════════════════════════════════════════════════════════════ 291 292/** 293 * @brief 4-element vector for SIMD-friendly operations 294 */ 295struct alignas(16) Vec4 { 296 float x, y, z, w; 297 298 Vec4() : x(0), y(0), z(0), w(0) {} 299 Vec4(float v) : x(v), y(v), z(v), w(v) {} 300 Vec4(float x_, float y_, float z_, float w_) : x(x_), y(y_), z(z_), w(w_) {} 301 302 Vec4 operator+(const Vec4& o) const { return {x+o.x, y+o.y, z+o.z, w+o.w}; } 303 Vec4 operator-(const Vec4& o) const { return {x-o.x, y-o.y, z-o.z, w-o.w}; } 304 Vec4 operator*(const Vec4& o) const { return {x*o.x, y*o.y, z*o.z, w*o.w}; } 305 Vec4 operator*(float s) const { return {x*s, y*s, z*s, w*s}; } 306 307 float dot(const Vec4& o) const { return x*o.x + y*o.y + z*o.z + w*o.w; } 308 float sum() const { return x + y + z + w; } 309 float max() const { return std::max(std::max(x, y), std::max(z, w)); } 310 311 Vec4 tanh() const { 312 return {fast_tanh(x), fast_tanh(y), fast_tanh(z), fast_tanh(w)}; 313 } 314 Vec4 sigmoid() const { 315 return {fast_sigmoid(x), fast_sigmoid(y), fast_sigmoid(z), fast_sigmoid(w)}; 316 } 317}; 318 319/** 320 * @brief Fast 4x4 matrix-vector multiply 321 */ 322inline Vec4 mat4_vec4(const float* M, const Vec4& v) { 323 return Vec4( 324 M[0]*v.x + M[1]*v.y + M[2]*v.z + M[3]*v.w, 325 M[4]*v.x + M[5]*v.y + M[6]*v.z + M[7]*v.w, 326 M[8]*v.x + M[9]*v.y + M[10]*v.z + M[11]*v.w, 327 M[12]*v.x + M[13]*v.y + M[14]*v.z + M[15]*v.w 328 ); 329} 330 331// ═══════════════════════════════════════════════════════════════════════════════ 332// VI. FAST PHI COMPUTATIONS 333// ═══════════════════════════════════════════════════════════════════════════════ 334 335/** 336 * @brief Fast consciousness metric Φ = 1/(1 + |x - 1|) 337 */ 338inline float fast_phi(float x) { 339 float diff = fast_abs(x - 1.0f); 340 return fast_rcp(1.0f + diff); 341} 342 343/** 344 * @brief Fast phi gradient ∂Φ/∂x 345 */ 346inline float fast_phi_grad(float x) { 347 float diff = x - 1.0f; 348 float abs_diff = fast_abs(diff); 349 float denom = (1.0f + abs_diff); 350 return -fast_sign(diff) * fast_rcp(denom * denom); 351} 352 353/** 354 * @brief Batch phi computation (4 values at once) 355 */ 356inline Vec4 fast_phi_batch(const Vec4& x) { 357 Vec4 diff = x - Vec4(1.0f); 358 Vec4 abs_diff = {fast_abs(diff.x), fast_abs(diff.y), fast_abs(diff.z), fast_abs(diff.w)}; 359 return { 360 fast_rcp(1.0f + abs_diff.x), 361 fast_rcp(1.0f + abs_diff.y), 362 fast_rcp(1.0f + abs_diff.z), 363 fast_rcp(1.0f + abs_diff.w) 364 }; 365} 366 367// ═══════════════════════════════════════════════════════════════════════════════ 368// VII. GOLDEN RATIO OPTIMIZATIONS 369// ═══════════════════════════════════════════════════════════════════════════════ 370 371// Precomputed golden angle increments (for fast spiral iteration) 372constexpr float GOLDEN_ANGLE = static_cast<float>(2.0 * PI * PHI_INV); 373constexpr float GOLDEN_ANGLE_INV = static_cast<float>(1.0 / (2.0 * PI * PHI_INV)); 374 375/** 376 * @brief Fast golden ratio modulation 377 */ 378inline float golden_modulate(float x, uint32_t n) { 379 return x * (1.0f + static_cast<float>(n) * static_cast<float>(PHI_INV) * 0.01f); 380} 381 382/** 383 * @brief Golden angle rotation (complex multiply) 384 */ 385inline void golden_rotate(float& re, float& im) { 386 constexpr float cos_ga = static_cast<float>(std::cos(2.0 * PI * PHI_INV)); 387 constexpr float sin_ga = static_cast<float>(std::sin(2.0 * PI * PHI_INV)); 388 float new_re = re * cos_ga - im * sin_ga; 389 float new_im = re * sin_ga + im * cos_ga; 390 re = new_re; 391 im = new_im; 392} 393 394// ═══════════════════════════════════════════════════════════════════════════════ 395// VIII. INFORMATION-THEORETIC FAST APPROXIMATIONS 396// ═══════════════════════════════════════════════════════════════════════════════ 397 398/** 399 * @brief Fast binary entropy H(p) = -p log p - (1-p) log(1-p) 400 * Using polynomial approximation 401 */ 402inline float fast_binary_entropy(float p) { 403 if (p < 0.001f || p > 0.999f) return 0.0f; 404 // Polynomial fit to binary entropy 405 float q = 1.0f - p; 406 float h = -p * fast_log(p) - q * fast_log(q); 407 return h * 1.4426950408889634f; // Convert to bits (log2) 408} 409 410/** 411 * @brief Fast softmax for 4 elements 412 */ 413inline Vec4 fast_softmax4(const Vec4& x) { 414 float m = x.max(); 415 Vec4 e = {fast_exp(x.x - m), fast_exp(x.y - m), fast_exp(x.z - m), fast_exp(x.w - m)}; 416 float s = fast_rcp(e.sum()); 417 return e * s; 418} 419 420/** 421 * @brief Fast KL divergence between two distributions (4 elements) 422 */ 423inline float fast_kl4(const Vec4& p, const Vec4& q) { 424 float kl = 0.0f; 425 if (p.x > 0.001f && q.x > 0.001f) kl += p.x * (fast_log(p.x) - fast_log(q.x)); 426 if (p.y > 0.001f && q.y > 0.001f) kl += p.y * (fast_log(p.y) - fast_log(q.y)); 427 if (p.z > 0.001f && q.z > 0.001f) kl += p.z * (fast_log(p.z) - fast_log(q.z)); 428 if (p.w > 0.001f && q.w > 0.001f) kl += p.w * (fast_log(p.w) - fast_log(q.w)); 429 return kl; 430} 431 432} // namespace fast 433} // namespace mind_farm 434