IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
1
Leveraging Mathematical Reasoning of LLMs for Efficient GPU Thread Mapping
arXiv:2604.10387v1 [cs.DC] 12 Apr 2026
Jose Maureira, Cristóbal A. Navarro, Hector Ferrada, and Luis Veas-Castillo
Abstract—Mapping parallel threads onto non-box-shaped domains is a known challenge in GPU computing that, if done efficiently, can prevent severe performance penalties arising from the allocation of unnecessary computational resources. Currently, achieving this optimal efficiency requires a significant amount of analytical human time and effort to manually derive bespoke mapping functions for each specific geometry. This work introduces a novel approach that leverages the symbolic reasoning capabilities of Large Language Models (LLMs) to automate this mathematical derivation process entirely through in-context learning. Focusing on state-of-the-art open-weights models, we conducted a rigorous comparative analysis across spatial domains of increasing complexity. Our results demonstrate that modern local LLMs successfully infer exact O(1) and O(log N ) mapping equations for complex 2D/3D dense domains and 2D fractals, vastly outperforming traditional symbolic regression methods which systematically failed at this discrete task. Crucially, we profile the energetic viability of this approach on highperformance infrastructure, drawing a clear distinction between the code-generation phase and the execution phase. While the one-time inference of the equation incurs a high energy penalty— particularly for reasoning-focused models like DeepSeek-R1— this is a single upfront investment. Once integrated, the generated analytical kernels eliminate block waste entirely, yielding massive and repeated energy and time savings (e.g., achieving up to a 4833× speedup and 2890× energy reduction) during actual GPU workloads. Finally, we identify a current “reasoning ceiling” when these models face highly recursive 3D fractals tested in this work (e.g., the Menger Sponge). This limitation establishes a clear benchmark for the present maturity of open-weight architectures, charting a viable and completely sovereign path toward fully automated, energy-efficient GPU resource optimization. Index Terms—Deep Learning, Large Language Models, Symbolic Regression, Thread Mapping, GPU Computing, Energy Efficiency.
challenges involve complex, irregular geometries, leading to unpredictable control flow and memory access patterns [26]. In fields like computational fluid dynamics, finite element analysis, and molecular simulations, the computational space is often non-uniform [3], [4]. For these domains, naive mapping strategies, such as using a simple bounding box (BB), are inefficient. As illustrated for 2D triangular domains (see Figure 1) and 3D tetrahedral domains (see Figure 2), a large fraction of threads may be allocated outside the actual problem domain, just be discarded at runtime, thus wasting valuable computational cycles and energy.
Fig. 1: Illustration of how the classic BB mapping is not efficient for a 2D triangular domain. The red blocks represent wasted computational resources (threads/blocks) that are allocated but fall outside the useful problem area. (Adapted from [1]).
I. I NTRODUCTION The GPU has become an indispensable tool for accelerating a vast range of applications in science and engineering, driven by its massively parallel architecture [1], [2]. The central paradigm of GPU programming involves subdividing a problem into many parts and solving them with thousands of parallel threads, which are organized into hierarchical structures such as blocks and grids. Maximizing the performance of these powerful devices hinges on a fundamental principle: ensuring that every thread is engaged in useful computation, thereby minimizing idle resources. While standard grid-based thread mapping is effective for problems defined on regular domains (e.g., vectors, arrays, tables, matrices, boxes, etc), a significant number of scientific Jose Maureira, Cristóbal A. Navarro, Hector Ferrada, and Luis Veas-Castillo are with the Instituto de Informática, Universidad Austral de Chile, Valdivia, Chile, and the Temporal research group (http://temporal.uach.cl). Corresponding author: Jose Maureira (e-mail: [email protected]).
Fig. 2: Conceptual representation of the inefficient BB mapping for a 3D triangular (tetrahedral) domain. A large portion of the allocated cube is wasted as it falls outside the problem space. (Adapted from [2]). Efficient thread mapping on GPUs is therefore essential for maximizing computational efficiency in various complex geometric domains. However, one of the main challenges lies in the high human effort required to manually formulate
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
symbolic equations for each specific domain. Researching on this field had historically involved a high amount of analytical effort in order to come up with new equations [5], [6], [7], [8], limiting the speed of progress for efficient GPU thread mapping for highly complex domains. This work takes a different route in GPU thread mapping research; it considers the recent advancements in open-weights large language models (LLMs) and how they can redefine the way research is done for this class of problems [9], [14], [16], [17]. Unlike traditional symbolic regression methods that function primarily as data-fitting algorithms, state-ofthe-art LLMs with sophisticated reasoning capabilities can perform a higher-level task akin to algorithm induction or pattern reverse-engineering, entirely with local infrastructure. The main contributions of this work are: • We present a novel framework for automatically inferring GPU thread mapping equations for complex domains using in-context learning with local, open-source LLMs. • We provide a comprehensive block-level performance and energy benchmark on high-performance infrastructure, revealing the critical trade-offs between model size, reasoning paradigms (e.g., Chain-of-Thought), and code execution efficiency. • We identify a current “reasoning ceiling” for open-source models on highly recursive 3D fractals, establishing a benchmark for the present and upcoming open-weight models in automated mathematical derivation. II. BACKGROUND AND R ELATED W ORK This section situates our work within the established literature of GPU computing and the rapidly evolving field of symbolic regression. A. Efficient GPU Thread Mapping for Complex Domains The core problem is to define a bijective mapping λ ∈ Zk+ 7→ Zd+ from a k-dimensional set of thread coordinates in grid space to a set of d-dimensional coordinates in data domain space [1]. This problem space relies on the fundamental GPU thread mapping principles described in canonical texts [27]. Previous research has systematically explored this problem, manually deriving mapping functions for various non-trivial domains, including triangular, tetrahedral, and fractal spaces [1], [2], [3], [5], [6], [7], [8]. A canonical example from this body of work is the mapping function for a 2D lower triangular domain [1], which relates a linear index λ to a coordinate pair (i, j). The forward mapping function, g(λ), is given by: % ! $r 1 1 i(i + 1) + 2λ − ,λ − (1) g(λ) = (i, j) = 4 2 2 This equation, while elegant, is non-trivial to derive from first principles. It involves an understanding of triangular numbers and requires a specific analytical insight to formulate. The manual effort required to produce such equations for progressively more complex domains escalates significantly. This work aims to automate the discovery of such functions.
2
B. Symbolic Equation Inference using Deep Learning Symbolic regression (SR), the automated discovery of mathematical expressions from data, presents a promising avenue for such automation. The field has evolved from genetic programming to incorporate modern deep learning methods, including transformer-based architectures for sequenceto-sequence generation [12]. Models such as Neural Symbolic Regression that Scales [9] and SymFormer [13] utilize encoder-decoder structures to generate equation skeletons from data. Other approaches have integrated prior knowledge [10] or used Monte-Carlo tree search [11]. However, these methods are fundamentally designed as continuous data-fitting algorithms aimed at minimizing numerical error (e.g., mean squared error). This poses a critical theoretical limitation for GPU thread mapping, which operates strictly in the discrete integer domain. An approximation, no matter how numerically close, is inherently invalid for indexing array addresses or hardware threads, as absolute algorithmic precision is mandatory. Therefore, traditional SR networks are structurally unsuited for deriving exact mappings. C. Large Language Models in Symbolic Reasoning The recent and rapid advancements in LLMs have opened a new frontier for symbolic discovery [14], [15], [16], [17]. Recent comprehensive surveys highlight the transformative evolution of these models in automated code generation [24] and mathematical reasoning [25], shifting their utility from pure natural language processing to complex algorithmic problemsolving. Instead of being trained specifically for continuous numerical SR, these models leverage their vast, pre-existing knowledge base to perform logical and symbolic reasoning tasks. A prominent method is in-context learning, where an LLM is guided to a solution through a prompt containing a few examples of the desired input-output behavior [14], [17]. This "few-shot" approach transforms the task from a continuous data-fitting problem into one of discrete pattern recognition and logical reverse-engineering, aligning perfectly with the goal of deriving exact, absolute mapping algorithms. This recent mathematical reasoning capability of modern LLMs is the ground from which this work builds upon.
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
3
Fig. 3: Overview of the proposed automated discovery pipeline. (1) Data extraction from the target domain, (2) Neural symbolic reasoning for algorithm induction, (3) Automated synthesis of the analytical mapping code, and (4) Integration and deployment of the discovered logic.
III. L EVERAGING LLM S FOR MORE EFFICIENT GPU THREAD MAPPING
The core of our proposal is to transform the manual, timeconsuming analytical process of deriving mapping functions into an automated pipeline driven by the symbolic reasoning of open-weight LLMs. This approach treats geometric mapping not as a data-fitting problem, but as an algorithmic induction task where the model discovers the underlying map. A. Algorithmic Induction via In-Context Learning Unlike traditional symbolic regression, the use of LLMs allows for true "Algorithm Induction". By providing a sequence of coordinates as context, we trigger the model’s ability to recognize complex patterns in integer sequences. This process is language-independent; the model effectively infers the abstract mathematical relationship required for absolute precision indexing and expresses it through functional code ready for parallel environments. B. Upfront Reasoning vs. Permanent Savings The procedure consists of a one-time "Reasoning Phase" where the model expends computational effort to derive the symbolic solution. Once the analytical logic is identified, it represents a permanent architectural optimization. This shifts the burden from a naive, wasteful runtime check (Bounding Box) to a one-time intelligent derivation, leading to the massive energy savings discussed in Section V.
C. The Automated Workflow and Experimental Validation The workflow is divided into four distinct phases as shown in Fig. 3. It is important to distinguish between the operational framework and the experimental validation presented in this study: 1) Context Sampling: Extraction of the first N elements of the domain (Stage 20, 50, or 100) to define the pattern. These initial points are generated trivially via standard sequential CPU execution. 2) Symbolic Inference: Neural induction of the mapping algorithm. The model is guided by a few-shot structured prompt that defines the operational constraints, the mathematical nature of the expected output, and provides the sampled coordinate context. 3) Algorithmic Synthesis: Generation of a self-contained analytical code block (O(1) or O(log N )) that implements the discovered logic. 4) Integration and Deployment: Direct application of the synthesized code into GPU kernels to optimize the amount of threads required. We emphasize that the large-scale validation (against 106 points) and the empirical energy/time profiling described in the following sections are part of our experimental protocol to formally prove the correctness and viability of the method. The value 106 is strictly an evaluation parameter; the inferred analytical mapping functions naturally generalize to arbitrary problem sizes, bounded only by the hardware’s numerical
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
4
TABLE I: Summary of evaluated computational domains, their geometric complexity, and the analytical or algorithmic Ground Truth equations required for exact thread mapping, along with their respective literature sources. Domain
Type
Complexity
2D Triangular [1], [4]
Dense
O(1)
Ground Truth Mapping Logic λ → c x=
jq
1 + 2λ − 12 4 q 3
3D Pyramid [2], [4]
Dense
O(1)
2D Sierpinski Gasket [3], [6]
Fractal
O(log3 N )
2D Sierpinski Carpet [6], [7]
Fractal
O(log8 N )
3D Sierpinski Pyramid [7]
Fractal
O(log4 N )
3D Menger Sponge [7]
Fractal
O(log20 N )
√
k
y =λ−
,
x(x+1) 2
729λ2 −3+27λ
1 = − 1, z = ⌊v⌋, (x, y) = + √ q 3 √ 32/3 3 3 729λ2 −3+27λ z(z+1)(z+2) f2D_Tri λ − 6 P P⌊log3 λ⌋ di 3i , ⃗vdi ∈ {(0, 0), (1, 0), (0, 1)} (x, y) = i=0 ⃗vdi 2i , where λ = P⌊log8 λ⌋ P (x, y) = i=0 di 8i , ⃗vdi ∈ {0, 1, 2}2 \ {(1, 1)} ⃗vdi 3i , where λ = P P⌊log4 λ⌋ di 4i , ⃗vdi ∈ {(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)} (x, y, z) = i=0 ⃗vdi 2i , where λ = P⌊log20 λ⌋ P di 20i , ⃗vdi ∈ {0, 1, 2}3 \ Vvoids (x, y, z) = i=0 ⃗vdi 3i , where λ =
v
representation limits (e.g., 32-bit or 64-bit registers). In a production environment, once the model’s reasoning reliability is established for a class of geometries, the discovery remains an autonomous process that does not require prior knowledge of the full ground truth dataset. IV. F RAMEWORK FOR E VALUATING S YMBOLIC I NFERENCE This section details the experimental protocol designed to evaluate the accuracy and efficiency of the selected LLMs. To simulate a realistic scenario of algorithmic discovery from scarce data, we employed an in-context learning approach with varying levels of information. A. Experiment Design We evaluated the models across six computational domains of increasing geometric and algorithmic complexity. Figure 4 visualizes these target geometries, and Table I summarizes the exact Ground Truth (GT) mathematical logic required to correctly map a linear thread index λ to the spatial coordinates c in each domain. The analytical and fractal equations presented are extracted directly from the established literature on non-linear block-space mapping [1], [2], [3], [4], [6], [7], [8]. As presented in Table I, dense domains rely on the exact mathematical inverse of arithmetic progressions and tetrahedral numbers (requiring floating-point square and cubic roots). In contrast, the mapping for fractal domains leverages a systematic base-B digit decomposition, where the linear index λ is decomposed into its constituent digits di , which in turn dictate the scaled spatial translation vectors ⃗vdi applied at each recursive depth. For each domain-model pair, we followed a strict evaluation pipeline designed to simulate algorithmic discovery from scarce data: 1) In-Context Learning Stages: The models were prompted with a sequence of integer coordinates representing the domain. We defined three "Stages" of information density to test the models’ ability to generalize from small samples: • Stage 20: The prompt includes the first 20 points of the domain. This tests extreme few-shot inference.
Fig. 4: Visual overview of the six evaluated computational domains: (a) 2D Triangular, (b) 2D Sierpinski Gasket, (c) 2D Sierpinski Carpet, (d) 3D Pyramid, (e) 3D Sierpinski Pyramid, and (f) 3D Menger Sponge.
• •
Stage 50: The prompt includes the first 50 points. Stage 100: The prompt includes the first 100 points, providing a richer context for pattern recognition.
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
2) Validation against Ground Truth: Once a model generates a Python function, we validate it against a massive Ground Truth (GT) dataset of N = 1, 000, 000 points. This ensures that the inferred logic holds for the entire domain and not just the few examples provided in the prompt. Specifically, the validation verifies that the inferred function produces a bijective mapping over a ground truth dataset of N = 1,000,000 points: every valid domain coordinate must be visited exactly once, with no repetitions and no omissions. To capture the full spectrum of model behavior, correctness is assessed under two complementary criteria. The first, Ordered, checks whether the model’s output matches the ground truth sequence index-byindex, rewarding exact algorithmic reproduction. The second, Any-order, checks whether all unique domain coordinates are covered regardless of traversal order, identifying solutions that recover the correct geometry but with a permuted index sequence. 3) Evaluation Metrics: We classify the performance using two accuracy metrics and one efficiency metric: 1) Ordered: The percentage of indices n ∈ [0, 106 ] where the model’s output (x, y, z) exactly matches the GT vector sequence. 2) Any-order: The percentage of unique GT coordinates found by the model, regardless of the index order. This identifies "Silver Standard" solutions where the geometric shape is correct, but the mapping order is permuted. 3) Big-O Efficiency: We perform static and dynamic analysis to ensure the generated code is optimal, extracting empirical measurements on GPU execution time and energy bounds. V. R ESULTS AND D ISCUSSION This section presents the findings from the experiments conducted on the local Patagón supercomputer infrastructure, focusing on the performance, efficiency, and limitations of the open-source LLM ecosystem. A. LLM Symbolic Inference Accuracy We evaluated 11 state-of-the-art open-weight LLMs (including models from the Llama 3.3/4 [19], Qwen3 [20], DeepSeek [21], Gemma3 [22], Mistral [23], and GPT-OSS families) across the six computational domains. These specific families were selected because they represent the current frontier of open-source capabilities, encompassing a diverse spectrum of underlying architectures—from traditional dense transformers to scalable Mixture-of-Experts (MoE) [18]—and training paradigms, including standard instruction-tuning and Reinforcement Learning (RL) driven Chain-of-Thought (CoT). Furthermore, these models are strictly comparable as they represent the leading open-source competitors systematically evaluated against each other in recent holistic benchmarks for code generation, algorithmic logic, and mathematical reasoning [28], [29]. The results demonstrate a clear stratification in reasoning capabilities, demonstrating that local models can handle complex discrete symbolic tasks.
5
In the accuracy tables (II–VII), model names are abbreviated to fit the column width and append their parameter sizes. The specific models evaluated in this benchmark correspond to: R1:70b (deepseek-r1:70b), Gem3:12b and Gem3:27b (gemma3), OSS:20b and OSS:120b (gpt-oss), Lla3.3:70b (llama3.3), Lla4:16x17b (llama4), MistN:12b (mistral-nemo), Nemo:70b (nemotron:70b), and Qw3:32b and Qw3:235b (qwen3). All models were executed locally in GGUF format using their default configuration parameters. Entries marked with (NC) indicate non-compiling or structurally invalid code. 2D Domains (Triangular & Fractals): For the basic 2D Triangular mapping (Table II), top-tier models like OSS:120b, OSS:20b, and R1:70b achieved 100% exact Ordered accuracy, demonstrating that basic arithmetic series are well within the grasp of open models. In the fractal domains, results diverged. For the 2D Sierpinski Gasket (Table III), only OSS variants successfully inferred the logic. Notably, for the highly complex 2D Sierpinski Carpet (Table IV), which requires recursive modular arithmetic, OSS:120b (at 100 points) and Qw3:235b (at 20/50 points) achieved 100% Ordered accuracy. This proves that massive open models can infer complex fractal logic that defeats smaller architectures. 3D Domains (Triangular & Pyramid): Table V reveals robust performance in 3D spaces. Models such as OSS:120b, OSS:20b, and the Qw3 family consistently solved the 3D Triangular mapping with 100% accuracy. An interesting phenomenon was observed with R1:70b: at 20 points, it deduced the correct geometric set, achieving 82.70% Any-order accuracy, but required more context to construct the precise index sequence, initially yielding 0.11% Ordered accuracy. For the 3D Sierpinski Pyramid (Table VI), OSS:120b was the only model to reach 100% perfect mapping. The "Menger Limit" (3D Menger Sponge): The 3D Menger Sponge (Table VII) is the hardest domain and acts as a frontier in symbolic reasoning for existing open-source LLMs. While these models excel at 2D fractals and dense 3D spaces, the Menger Sponge requires tracking cubic coordinate systems and recursive 3D void removals. The open-source ecosystem currently achieves an Any-order score of < 1% on this domain. This establishes a clear “reasoning boundary” for state-of-the-art open weight models, making it a benchmark for the next generation of LLMs. TABLE II: Results for 2D Triangular mapping. Model R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
20 pts
50 pts
100 pts
Ord.
Any
Ord.
Any
Ord.
Any
100% 0.00% 0.00% 100% 0.00% 100% 0.00% 0.00% 0.00% 100% 100%
100% 0.00% 50.05% 100% 0.71% 100% 0.71% 0.71% 0.00% 100% 100%
100% 0.00% 0.00% 100% 100% 0.00% 0.00% 0.00% 0.00% 0.14% 100%
100% 1.27% 1.27% 100% 100% 0.00% 1.27% 1.27% 0.14% 0.14% 100%
100% 0.00% 0.00% 100% 100% 0.00% 0.00% 0.00% 100% 0.00% (NC) 100%
100% 1.83% 50.05% 100% 100% 0.14% 0.01% 1.69% 100% 0.00% 100%
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
TABLE III: Results for 2D Sierpinski Gasket. 20 pts
Model R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
50 pts
6
TABLE VII: Results for 3D Menger Sponge.
100 pts
Model
Ord.
Any
Ord.
Any
Ord.
Any
0.00% 0.00% 0.00% 0.00% 100% 0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00%
8.10% 1.03% 1.03% 8.10% 100% 7.96% 0.34% 0.00% 8.10% 0.00% 8.10%
4.57% 0.00% 0.00% 100% 0.00% (NC) 0.00% 0.00% 0.00% 0.00% 0.00% 0.00%
21.30% 1.55% 5.22% 100% 0.00% 1.17% 0.00% 3.09% 8.10% 0.00% 0.01%
0.00% 0.00% 0.00% 100% 100% 0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00% (NC)
1.52% 0.69% 5.22% 100% 100% 3.19% 0.01% 0.01% 8.10% 0.00% 0.00%
R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
20 pts
50 pts
100 pts
Ord.
Any
Ord.
Any
Ord.
Any
0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00% 0.00%
0.05% 0.05% 0.05% 0.00% 0.00% 0.05% 0.06% 0.03% 0.00% 0.05% 0.00%
0.00% (NC) 0.00% 0.00% 0.01% 0.01% 0.00% 0.00% 0.00% 0.00% 0.01% 0.00%
0.00% 0.36% 0.05% 0.16% 0.16% 0.04% 0.16% 0.00% 0.05% 0.16% 0.04%
0.00% 0.00% 0.00% 0.01% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00%
0.05% 0.05% 0.05% 0.36% 0.00% 0.36% 0.16% 0.11% 0.01% 0.00% 0.14%
TABLE IV: Results for 2D Sierpinski Carpet. 20 pts
Model R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
50 pts
100 pts
Ord.
Any
Ord.
Any
Ord.
Any
0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 100% 0.00%
0.58% 0.58% 0.39% 0.58% 0.58% 0.39% 0.58% 0.39% 0.00% 100% 0.00%
0.00% 0.00% 0.00% (NC) 0.01% 0.00% (NC) 0.00% 0.00% 0.00% 0.00% 100% 0.00%
0.00% 0.39% 0.20% 1.04% 0.00% 0.39% 1.04% 1.04% 0.58% 100% 0.03%
0.00% 0.00% 0.00% 100% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00%
37.08% 0.58% 1.04% 100% 0.58% 0.46% 1.56% 1.30% 0.10% 0.00% 0.58%
TABLE V: Results for 3D Triangular mapping. 20 pts
Model R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
50 pts
100 pts
Ord.
Any
Ord.
Any
Ord.
Any
0.11% 0.00% 0.00% 100% 0.00% (NC) 0.00% 0.00% 0.00% 0.00% 100% 100%
82.70% 0.02% 0.00% 100% 0.00% 0.00% 0.00% 0.05% 0.14% 100% 100%
100% 0.00% 0.00% 100% 100% 0.00% 0.00% 0.00% 0.00% 0.00% 100%
100% 0.02% 0.00% 100% 100% 17.16% 0.00% 0.18% 0.00% 16.96% 100%
0.00% 0.00% 0.00% 100% 100% 0.00% 0.00% 0.00% 0.00% 100% 100%
0.00% 0.02% 17.17% 100% 100% 0.00% 0.00% 0.00% 0.00% 100% 100%
TABLE VI: Results for 3D Sierpinski Pyramid. Model R1:70b Gem3:12b Gem3:27b OSS:120b OSS:20b Lla3.3:70b Lla4:16x17b Mist-N:12b Nemo:70b Qw3:235b Qw3:32b
20 pts
50 pts
100 pts
Ord.
Any
Ord.
Any
Ord.
Any
0.00% 0.00% 0.00% 100% 0.00% (NC) 0.00% (NC) 0.00% 0.00% 0.00% (NC) 0.00% (NC) 0.00%
0.00% 0.20% 0.31% 100% 0.00% 0.59% 0.01% 0.49% 0.00% 0.00% 0.01%
0.00% 0.00% 0.00% 0.00% 0.00% (NC) 0.00% (NC) 0.00% 0.00% 0.00% (NC) 0.00% (NC) 0.00%
0.00% 0.10% 0.18% 1.23% 0.00% 0.00% 1.87% 0.00% 0.00% 0.00% 0.52%
0.00% 0.00% (NC) 0.00% 100% 0.00% (NC) 0.00% 0.00% (NC) 0.00% 0.00% 0.00% (NC) 0.00% (NC)
0.00% 0.00% 0.00% 100% 0.00% 0.28% 0.00% 0.00% 2.52% 0.00% 0.00%
presents this efficiency metric across all six geometric domains, evaluated at 20, 50, and 100 in-context examples. 1) Impact of Model Architecture and Reasoning: Two distinct efficiency profiles emerge within the open-source ecosystem: 1) Parameter-Driven Penalties: Massive models like qwen3:235b exhibit lower points-per-joule efficiency driven by the substantial memory bandwidth required to move 235.1 billion parameters across the 4 GPUs, yielding a high baseline power draw during execution. 2) Reasoning-Driven Penalties: Notably, deepseek-r1:70b often matches or falls below the efficiency of models three times its size. This is attributed to its "Chain-of-Thought" (CoT) mechanism. While standard dense models like llama3.3 (70.6B) emit code quickly, DeepSeek-R1 requires significantly more time to "reason" through the problem, imposing a severe energy penalty (fewer points processed per joule) compared to standard models of equivalent parameter count. 2) Efficiency Gains from In-Context Learning: Increasing the number of in-context examples generally improved the energy efficiency. For standard models like llama3.3 and mistral-nemo, moving from 20 to 50 or 100 shots resulted in a noticeable increase in efficiency. Qualitative analysis suggests that with fewer examples (20 shots), models are prone to generating verbose code that fails to compile or "hallucinating" lengthy explanations. Providing a richer context constrains the generation space, guiding the model to produce concise, correct solutions faster, thereby minimizing GPU active time per valid point. For deepseek-r1, it occasionally shows an efficiency drop at 50 shots, possibly due to an increased activation of its internal reasoning traces before converging efficiently at 100 shots.
B. Energy Cost of Symbolic Inference
C. Performance and Energy in Block-Level Execution
We also evaluated the upfront energy cost of running the LLM inference itself. We measured the computational efficiency, defined as Points processed per Joule (Points/Joule), during the equation derivation phase. The hardware resource for this phase was an allocation of 4x NVIDIA A100 (40GB SXM4) GPUs within a DGX node. To ensure clean measurements, all inferences were executed in strict isolation (one model at a time) and energy consumption was recorded using NVIDIA’s NVML (nvidia-smi) profiling tools. Figure 5
Once the LLM successfully infers the correct analytical GPU thread function, this one-time upfront inference cost allows doing the actual GPU work by just evaluating the inferred map with all threads, not needing to infer again. To quantify this, the analytical expressions derived by the LLMs were tested with a lightweight dummy CUDA kernel that performs an atomic increment (+1) on each mapped memory address, simulating a representative memory-bound workload without introducing algorithm-specific overhead. This kernel
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
7
Fig. 5: Computational efficiency (Points/Joule) of each open-source model across six spatial domains—three 2D (top) and three 3D (bottom)—evaluated at 20, 50, and 100 few-shot examples. Higher values indicate greater energy efficiency on a logarithmic scale. Bar annotations report the total energy consumed per inference run. Efficiency is defined as the number of correctly mapped points divided by total energy consumed.
was integrated into a CUDA block-level mapping kernel which is the state of the art way of applying efficient GPU thread
mapping [1], [2], [6], [7]. For these performance and power profiling tests, the kernels
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
were executed on a single, isolated NVIDIA A100 (40GB SXM4) GPU within the DGX node, using NVML for precise energy polling. We report the gross energy draw, which includes both the dynamic calculation energy and the hardware’s idle baseline power. Each table includes two reference baselines: the Bounding Box (BB), representing the naive mapping strategy, and the Paper entry, corresponding to the analytically derived mapping from the established literature, which serves as the human-effort gold standard against which LLM-generated solutions are compared. Dense Geometries: As shown in Table VIII, the use of optimal O(1) cost analytical functions generated by top models eliminates the resource waste inherent in the Bounding Box (BB) method without sacrificing performance. In the 3D Pyramid domain, the BB approach wastes approximately 83% of the launched blocks, resulting in an execution time of 2530.65 ms and consuming over 282 Joules. In contrast, the optimal inferred kernels reduce this time to 3.84 ms and consume less than 1 Joule. However, the results highlight the penalty of inefficient code generation: while OSS:120b (Stage 20) correctly solved the mapping logic, it implemented a Linear Search algorithm (O(N 1/3 ))—so named because it iterates linearly over candidate integer values of the depth parameter until the boundary condition is satisfied, despite the closed-form involving a cube root—rather than a direct analytical solution, degrading performance to 117.03 ms and consuming 22.25 J—over 24× more energy than the optimal solution. Fractal Geometries: The impact of exact analytical mapping is most pronounced in fractal domains (Table IX), where the sparsity of the geometry renders the BB approach unfeasible for large domains. For the 3D Sierpinski structure, the BB method requires traversing a dense enclosing cube, launching over 8 billion blocks to identify only 1.9 million valid ones of the actual fractal geometry which is not dense. This results in an estimated execution time of nearly 16 seconds (projected from a smaller scale due to extreme physical timeout overhead) and an energy cost of ∼ 1591 J. The bitwise operations inferred by the top-performing LLMs allow for direct indexing of valid blocks, reducing execution time to 3.30 ms and energy consumption to just 0.55 Joules. This represents an effective speedup of ∼ 4833× and an energy reduction of 2890×, instantly amortizing the initial LLM inference energy cost on the very first execution. D. Discussion: Parameter Scaling vs. Reasoning Our findings challenge the traditional scaling laws in the context of exact symbolic derivation. While large dense models like qw3:235b exhibit high memory bandwidth and energy demands, they do not inherently guarantee success in complex domains like the 3D Sierpinski Pyramid. Conversely, reasoning-focused architectures (e.g., R1:70b) leverage latent Chain-of-Thought processing to correctly identify recursive mapping patterns that defeat larger parameter-dense counterparts. Ultimately, the boundaries observed on the 3D Menger Sponge outline a current “Menger Limit” for the open ecosystem. This domain forces the LLM to simultaneously track
8
cubic index progressions and non-trivial 3D void removals, marking the absolute cutting edge of open-source symbolic inference. We note that this represents the current state of open-weights training maturity rather than a fundamental limitation in the LLM paradigm. As open models continue to adopt more advanced techniques—such as sophisticated Mixture-ofExperts (MoE) routing and deeper reinforcement learning—we anticipate this current frontier will soon be surpassed. From a practical standpoint, this implies that researchers do not necessarily need massive 200B+ parameter models requiring multi-GPU clusters; highly optimized reasoning models can achieve state-of-the-art symbolic derivation on a single node. E. Applicability, Limitations, and Framework Extensibility The proposed framework is specifically applicable to computational domains that exhibit deterministic mathematical patterns, such as arithmetic progressions (dense geometries) or recursive self-similarity (fractals). It is not designed for completely unstructured meshes where coordinate mapping depends on arbitrary data points without an underlying mathematical law. While a human expert could manually interact with an LLM to derive an equation for a singular geometry, the primary value of this framework lies in the complete automation of the pipeline, removing human analytical effort from the compilation or execution stages of HPC workloads. Furthermore, the system is fundamentally modular and model-agnostic. As new, more capable open-weight architectures are released, they can be deployed within the pipeline without modifying the integration logic. Regarding the incontext learning process, the structured prompt detailed in Appendix A is empirically robust in enforcing strict algorithmic generation. We do not assert that it constitutes the absolute theoretical optimum of prompt engineering; rather, it establishes a reliable, reproducible baseline necessary to evaluate and expose the current mathematical reasoning limits of existing neural architectures.
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
9
TABLE VIII: CUDA Block-Level Performance and Energy in Dense Geometries (N = 500 × 106 ). Domain
Model (Stage)
Time (ms)
Total Blocks
Wasted
Energy (J)
2D Triangular
Bounding Box (Baseline) Paper (Navarro 2014) R1:70b (Stage 20) R1:70b (Stage 50) OSS:120b (All Stages) Lla3.3:70b / Nemo:70b R1:70b (Stage 100) OSS:20b (Stage 50) OSS:20b (Stage 100) Qw3:32b (Stage 50)
747.45 1.46 1.46 1.46 1.46 1.46 1.97 1.51 1.51 14.86
3,912,484 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125
1,959,359 0 0 0 0 0 0 0 0 0
83.27 0.44 0.45 0.45 0.45 0.45 0.70 0.51 0.51 3.21
If (O(1)) Analytical (O(1)) Analytical (O(1)) Analytical (O(1)) Analytical (O(1)) Analytical (O(1)) Sqrt+Loops (O(1)) Approx+If (O(1)) Approx+If (O(1)) BinSearch (O(log N ))
Logic (Complexity)
3D Pyramid
Bounding Box (Baseline) Paper (Navarro 2016) R1:70b (Stage 50) Qw3:32b (All Stages) OSS:120b (Stage 100) Qw3:235b (Stage 20) OSS:120b (Stage 50) OSS:120b (Stage 20)
2530.65 3.84 6.21 6.21 29.31 29.31 51.57 117.03
12,008,989 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125 1,953,125
10,055,864 0 0 0 0 0 0 0
282.67 0.92 1.44 1.44 5.99 5.99 9.12 22.25
If (O(1)) Analytical (O(1)) Cbrt+Loop (O(1)) Cbrt+Loop (O(1)) BinSearch (O(log N )) BinSearch (O(log N )) BinSearch+Lin (O(N 1/3 )) Linear (O(N 1/3 ))
TABLE IX: CUDA Block-Level Performance and Energy in Fractal Geometries (N = 500 × 106 ). Domain
Model (Stage)
2D Sierpinski
Bounding Box (Baseline) Paper (Reference) OSS:120b (Stage 20)
Time (ms)
Total Blocks
Wasted
Energy (J)
Logic (Complexity)
65.78 8.62 8.62
88,736,400 1,953,125 1,953,125
86,783,275 0 0
6.73 1.39 1.39
If (O(1)) Bitwise (O(log N )) Bitwise (O(log N ))
Bounding Box (Baseline)* 15,949.00 8,000,000,000 7,998,046,875 Paper (Reference) 3.30 1,953,125 0 R1:70b (Stage 100) 3.30 1,953,125 0 * Baseline projected from N = 5 × 106 due to physical timeout on extreme waste overhead.
1,591.71 0.55 0.56
If (O(1)) Bitwise (O(log N )) Bitwise (O(log N ))
3D Sierpinski
VI. C ODE AND DATA AVAILABILITY The source code, generated domain datasets, and CUDA evaluation kernels associated with this manuscript are publicly available to ensure full reproducibility. The repository can be accessed at https://github.com/aspiadevs/ llm-gpu-thread-mapping. VII. C ONCLUSIONS AND F UTURE W ORK This research demonstrates the viability of utilizing stateof-the-art open-source LLMs to automatically infer exact symbolic equations for GPU thread mapping, serving as an alternative to manual human mathematical derivation. More in detail, three findings are worth highlighting: 1) Open-Source LLM Viability: Open-weights LLMs can infer efficient GPU thread maps for 2D/3D regular and fractal geometries. In comparison, a traditional Symbolic Regression (SR) approach was unable to deliver accurate symbolic functions for GPU thread mapping. 2) The Energy Trade-off: Automating this derivation locally exposes significant energy disparities. While generating the code incurs a high one-time upfront energy cost—especially when utilizing reasoning-heavy models like DeepSeek-R1—this investment is instantly amortized. Once integrated, the inferred GPU maps eliminate the unnecessary blocks, saving orders of magnitude in execution time and energy (up to 4833× faster and 2890× more energy efficient) compared to the baseline).
3) The Menger Frontier: While open models excel in many spatial domains, highly recursive 3D fractal structures (such as the Menger Sponge) represent the current reasoning ceiling for the open ecosystem. This establishes a clear, measurable benchmark for future openweights architectures to overcome as their mathematical reasoning capabilities evolve. Future Work: Future research will focus on mitigating the upfront energy penalty of LLM inference by exploring the targeted fine-tuning of smaller, highly efficient open-weights models (e.g., sub-10B parameters) specifically specialized for discrete spatial reasoning. This approach aims to deliver the exactness of massive reasoning models at a fraction of the initial computational cost. Additionally, we plan to extend this automated mapping framework to encompass more heterogeneous HPC topologies, such as unstructured meshes and adaptive refinement grids, while continuously monitoring the evolution of open-source architectures to verify when the 3D recursive fractal boundary is definitively breached. ACKNOWLEDGMENT This work was supported by Universidad Austral de Chile and ANID FONDECYT grants #1221357, #1261152, the Temporal research group and the Patagón supercomputer of Universidad Austral de Chile (FONDEQUIP EQM180042). R EFERENCES [1] Navarro, C. A., & Hitschfeld, N. (2014). GPU maps for the space of computation in triangular domain problems. Proceedings - 16th
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
IEEE International Conference on High Performance Computing and Communications, HPCC 2014. https://doi.org/10.1109/HPCC.2014.64 [2] Navarro, C. A., Bustos, B., & Hitschfeld, N. (2016). Potential benefits of a block-space GPU approach for discrete tetrahedral domains. Proceedings of the 2016 42nd Latin American Computing Conference, CLEI 2016. https://doi.org/10.1109/CLEI.2016.7833394 [3] Navarro, C. A., Vega, R., Bustos, B., & Hitschfeld, N. (2017). BlockSpace GPU Mapping for Embedded Sierpiński Gasket Fractals. Proceedings - 2017 IEEE 19th Intl Conference on High Performance Computing and Communications, HPCC 2017. https://doi.org/10.1109/ HPCC-SmartCity-DSS.2017.56 [4] Navarro, C. A., Vernier, M., Bustos, B., & Hitschfeld, N. (2018). Competitiveness of a non-linear block-space GPU thread map for simplex domains. IEEE Transactions on Parallel and Distributed Systems, 29(12). https://doi.org/10.1109/TPDS.2018.2849705 [5] Navarro, C. A., Bustos, B., & Hitschfeld, N. (2019). Analysis of a SelfSimilar GPU Thread Map for Data-parallel m-Simplex Domains. 2019 International Conference on High Performance Computing and Simulation, HPCS 2019. https://doi.org/10.1109/HPCS48598.2019.9188081 [6] Navarro, C. A., Quezada, F. A., Hitschfeld, N., Vega, R., & Bustos, B. (2020). Efficient GPU thread mapping on embedded 2D fractals. Future Generation Computer Systems, 113. https://doi.org/10.1016/j. future.2020.07.006 [7] Quezada, F. A., Navarro, C. A., Hitschfeld,N., & Bustos, B. (2022). Squeeze: Efficient compact fractals for tensor core GPUs. Future Generation Computer Systems, 135. https://doi.org/10.1016/j.future.2022.04. 023 [8] Navarro, C. A., Quezada, F. A., Bustos, B., Hitschfeld, N., & Kindelan, R. (2022). A scalable and energy efficient GPU thread map for msimplex domains. Future Generation Computer Systems, 141. https: //doi.org/10.1016/j.future.2022.12.020 [9] Biggio, L., et al. (2021). Neural Symbolic Regression that Scales. Proceedings of the 38th International Conference on Machine Learning (ICML), PMLR 139:936-945. [10] Bendinelli, T., et al. (2023). Controllable Neural Symbolic Regression. Proceedings of the 40th International Conference on Machine Learning (ICML), PMLR 202:2063-2077. [11] Kamienny, P. A., et al. (2023). Deep Generative Symbolic Regression with Monte-Carlo Tree Search. Proceedings of the 40th International Conference on Machine Learning (ICML), PMLR 202:15682-15697. [12] Vaswani, A., et al. (2017). Attention is All You Need. Advances in Neural Information Processing Systems, 30. https://arxiv.org/abs/1706. 03762 [13] Vastl, M., et al. (2024). SymFormer: End-to-End Symbolic Regression Using Transformer-Based Architecture. IEEE Access, 12. https://doi.org/ 10.1109/ACCESS.2024.3374649 [14] Merler, M., et al. (2024). In-Context Symbolic Regression: Leveraging Large Language Models for Function Discovery. Proceedings of the 62nd Annual Meeting of the Association for Computational Linguistics. https://aclanthology.org/2024.acl-srw.49/ [15] Li, Y., et al. (2024). MLLM-SR: Conversational Symbolic Regression base Multi-Modal Large Language Models. arXiv preprint. https://arxiv. org/abs/2406.05410v1 [16] Shojaee, P., et al. (2024). LLM-SR: Scientific Equation Discovery via Programming with Large Language Models. https://arxiv.org/abs/2404. 18400 [17] Sharlin, S., & Josephson, T. (2024). In Context Learning and Reasoning for Symbolic Regression with Large Language Models. https://arxiv.org/ abs/2410.17448 [18] Shazeer, N., et al. (2017). Outrageously Large Neural Networks: The Sparsely-Gated Mixture-of-Experts Layer. https://arxiv.org/abs/1701. 06538 [19] AI@Meta (2024). The Llama 3 Herd of Models. https://arxiv.org/abs/ 2407.21783 [20] Qwen Team (2024). Qwen2 Technical Report. https://arxiv.org/abs/2407. 10671 [21] DeepSeek-AI (2025). DeepSeek-R1: Incentivizing Reasoning Capability in LLMs via Reinforcement Learning. https://arxiv.org/abs/2501.12948 [22] Google DeepMind (2024). Gemma: Open Models Based on Gemini Research and Technology. https://arxiv.org/abs/2403.08295 [23] Mistral AI (2024). Mistral Nemo. https://mistral.ai/news/mistral-nemo/ [24] Zheng, Z., Ning, K., Wang, Y., Zhang, J., Zheng, D., Ye, M., & Chen, J. (2024). A Survey of Large Language Models for Code: Evolution, Benchmarking, and Future Trends. https://arxiv.org/abs/2311.10372 [25] Yu, Z., et al. (2025). From System 1 to System 2: A Survey of Reasoning Large Language Models. https://arxiv.org/abs/2502.17419
10
[26] Burtscher, M., Nasre, R., & Pingali, K. (2012). A quantitative study of irregular programs on GPUs. 2012 IEEE International Symposium on Workload Characterization (IISWC). https://doi.org/10.1109/IISWC. 2012.6402918 [27] Kirk, D. B., & Hwu, W. W. (2016). Programming Massively Parallel Processors: A Hands-on Approach (3rd Edition). Morgan Kaufmann. https://www.sciencedirect.com/science/book/9780128119860 [28] Jain, N., et al. (2024). LiveCodeBench: Holistic and Contamination Free Evaluation of Large Language Models for Code. https://arxiv.org/abs/ 2403.07974 [29] Liu, J., et al. (2024). Is Your Code Generated by ChatGPT Really Correct? Rigorous Evaluation of Large Language Models for Code Generation. Advances in Neural Information Processing Systems 36 (NeurIPS 2023). https://arxiv.org/abs/2305.01210
A PPENDIX A D ETAILED P ROMPT S PECIFICATION To guide the models in equation inference, a detailed prompt was designed that specifies the expert role, the inference task, and the output requirements. The data points for each test case were appended to the following template: <ROLE> Act as an expert in mathematics and cryptography, specializing in the reverse engineering of algorithms and the identification of complex patterns in multidimensional spaces. Your goal is SOLELY to generate the Python code requested. </ROLE> <TASK> Analyze the mapping data in the <CONTEXT> to find the underlying mathematical algorithm . Then, generate the complete source code for a single Python function that implements this general algorithm. </TASK> <CONTEXT> # Mapping Data __MAPPING_DATA_HERE__ </CONTEXT> <RULES> - Function name must be exactly ‘ map_to_coordinates(n)‘. - Input: ‘n‘ (non-negative integer). - Output: tuple of integers representing coordinates. - Each integer within the returned coordinate tuple must be greater than or equal to 0. - Validate input ‘n‘ (non-negative integer), raise ‘ValueError‘ if invalid. - **CRITICAL ALGORITHM CONSTRAINT:** The function MUST implement a general mathematical algorithm that works for ANY non-negative integer ’n’, not just the examples provided. - **DO NOT use hardcoded values, lookup tables , or long ‘if/elif‘ chains based on ranges of ’n’ (e.g., ‘if n == 1:‘, ‘if n < 10:‘, ‘if 10 <= n <= 20:‘ are forbidden).** - **CRITICAL OUTPUT CONSTRAINT:** Your response MUST contain ONLY the Python code block for the function.
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS, VOL. XX, NO. YY, APRIL 2026
- **DO NOT include ANY introductory text, explanations, reasoning, thought processes , or comments (including docstrings or # comments) inside or outside the code block .** - Do NOT include an ‘if __name__ == "__main__ ":‘ block. </RULES> <RESPONSE>
Jose Maureira received the M.Sc. degree in Computer Science from the Universidad Austral de Chile in 2026, where he is currently pursuing the Ph.D. degree in Engineering Sciences. He is a Research Assistant with the Temporal research lab, and his research interests include high-performance computing, large language models, computer vision, and the generation of synthetic data.
Cristóbal A. Navarro has a Ph.D. degree in computer science from the University of Chile (2015). Currently, he is an associate professor at the Universidad Austral de Chile and leads the Temporal research lab as well as the Patagón Supercomputer project. Today, his research interests include GPU computing, computer graphics, and computational physics.
Hector Ferrada received his Ph.D. degree in Computer Science from the University of Chile in 2016, focusing on his research in the design and analysis of algorithms for compact data structures. In 20162017, he conducted postdoctoral research in the Genome Scaling Algorithms Group at the University of Helsinki, Finland, in collaboration with Dr. Veli Mäkinen. Currently, he mainly teaches courses related to his research interests in algorithms and data structures.
Luis Veas-Castillo received the Ph.D. degree in Engineering Sciences from the Universidad de Santiago de Chile in 2022. His research interests include parallel and distributed systems, computational simulation, large-scale databases, and applied artificial intelligence. He has extensive experience leading interdisciplinary R&D projects in collaboration with Chilean institutions and research centers. He is currently driving applied research initiatives at the Universidad Austral de Chile, integrating engineering with healthcare and environmental sciences.
11