Pith. sign in

REVIEW 4 major objections 6 minor 28 references

GPU Implementation of Zippel Method for Feynman Integral Reconstruction

T0 review · 4 major / 6 minor · reviewed 2026-08-07 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗

Pith's one-line read The paper reports the first GPU implementation of the Zippel algorithm for sparse rational-function reconstruction, with measured speedups up to 14.57x over an Intel Xeon CPU on a 4.8-million-monomial Feynman integral reduction.

desk verdict Solid GPU-engineering paper with plausible speedups, but the missing GPU-vs-CPU correctness check is the kind of thing a referee must ask for before this speedup claim is trusted. read the letter →

arxiv 2505.24256 v1 pith:RIWJWQOL submitted 2025-05-30 hep-ph hep-th

Alexander V. Smirnov , Boris I. Rozhnov , Vadim V. Voevodin This is my paper
classification hep-phhep-th MSC 68W3065Y05
keywords ZippelalgorithmrationalreconstructionsparsepolynomialinterpolationFeynmanintegralsGPUcomputingCUDAmodulararithmeticparallelalgorithms
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper aims to show that the Zippel algorithm, a standard method for reconstructing sparse multivariate polynomials and rational functions from evaluations, can be ported to graphics processing units and run much faster than on a CPU. This matters because Zippel reconstruction of intermediate expressions can become the bottleneck in automating Feynman integral reductions in particle physics, where millions of monomial terms must be handled. The authors implement both the classical Zippel algorithm and its balanced variant for rational functions in a GPU programming model, porting low-level modular arithmetic routines to GPU kernels. On the largest test case, a reconstruction with 4.8 million skeleton monomials, the fastest GPU is reported 14.57 times faster than a 12-core Intel Xeon CPU, with GPU utilization above 90 percent.

What carries the argument

The central object is the Zippel reconstruction step itself: given a reconstructed polynomial in $k-1$ variables, the algorithm assumes the same monomial support for the next slice $x_k = c_i$, builds a Vandermonde system from evaluations at specially chosen points, and solves it in $O(t^2)$ time, then runs univariate Newton reconstruction in $x_k$. The GPU port carries this out with ported modular arithmetic kernels for multiplication, inversion, and reduction over machine-word primes, evaluating polynomials with Horner's scheme, generating power vectors in parallel, and using warp-level shuffle reductions for summation. Key engineering devices are aligned memory layouts, a shared-memory accumulator that defers modular reduction until the end, emulating a 192-bit integer, and fast exponentiation arrays for balanced coefficient preparation.

What would settle it

Run the CPU and GPU implementations on the same four input datasets from the paper's benchmarks and compare every reconstructed coefficient (or the resulting rational functions) symbolically; any mismatch would falsify the implied correctness, while a full match would confirm the speedups are genuine.

Watch

Extended reading notes

Core claim

On its own terms, the paper establishes that a GPU implementation of the balanced Zippel reconstruction is feasible and scales with problem size. The implementation reconstructs in modular arithmetic over large primes, using a skeleton polynomial's support to set up Vandermonde systems whose solution reveals the coefficients for the next variable, and the balanced variant handles rational functions by reconstructing numerator and denominator separately. Measured speedups over the CPU version grow monotonically with problem size, reaching 3.84x on a P100, 8.85x on a V100, and 14.57x on an A100 for the 4.8-million-monomial input, while relative GPU-vs-GPU speedups rise from under 1.5x for small inputs to 3.8x between A100 and P100. The authors report GPU load above 90 percent for the largest case and note that the gain comes from parallelizing the costliest reconstruction steps and from replacing repeated modular reductions with deferred reductions over a wider accumulator.

Load-bearing premise

The load-bearing premise is that the GPU kernels compute exactly the same modular arithmetic as the CPU version; the paper reports wall-clock benchmarks but never compares reconstructed coefficients between the two, so a port that silently changed results would still appear as a big speedup.

Editorial extensions

If this is right

  • Feynman integral reductions whose reconstruction step takes hours on a CPU can be completed in a fraction of the time on GPU-equipped systems; the paper's largest CPU run took about 20.7 hours and the A100 run about 1.4 hours.
  • The speedup grows with skeleton size, so the method is most useful precisely for the largest, most expensive reconstructions that previously strained supercomputer time limits.
  • The GPU implementation frees CPU cores for other reduction tasks and can be scaled horizontally by distributing independent coefficients across GPUs or nodes.
  • The ported modular arithmetic layer is reusable beyond this algorithm; the authors plan to release it separately, which could accelerate other modular polynomial computations.
  • Both the plain Zippel algorithm and its balanced rational-function variant are covered, so the approach applies to full Feynman integral reduction output, not just polynomial interpolation.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the correctness gap is closed, the same CUDA modular-arithmetic kernels could accelerate other reconstruction schemes that are not Zippel-based, such as Newton interpolation or p-adic approaches, since they share the same low-level operations.
  • The observed pattern, with speedup increasing with problem size and small-input differences between GPUs being statistically insignificant, suggests the implementation is limited by memory latency or launch overhead on small problems and by arithmetic throughput on large ones; profiling kernel occupancy would test this.
  • Because independent coefficients can be reconstructed in parallel, the algorithm could be extended to multi-GPU or distributed-memory settings where the skeleton is partitioned by coefficient, potentially scaling beyond the 40 GB memory of a single A100.
  • Releasing the code as a standalone library would let other Feynman integral reduction tools adopt the GPU reconstruction without adopting the full reduction framework, possibly broadening the impact beyond the authors' own package.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 6 minor

Summary. The paper describes a GPU implementation of the Zippel algorithm for sparse multivariate rational function reconstruction, aimed at Feynman integral reduction. The authors port FLINT modular arithmetic routines to CUDA, implement the Zippel and balanced Zippel reconstruction on GPUs, and benchmark four problem sizes (125k, 300k, 750k, 4.8m monomials) on P100, V100, and A100 GPUs against an Intel Xeon Gold 6126 CPU baseline. They report speedups up to 3.84x, 8.85x, and 14.57x on the largest case for P100, V100, and A100, respectively, and claim this is the first GPU implementation of the Zippel algorithm.

Significance. If the reported speedups are real and the GPU implementation produces mathematically identical results to the CPU version, this is a practically valuable contribution: Zippel reconstruction is indeed a bottleneck in some Feynman integral reductions, and a 14-fold speedup on a single GPU could significantly shorten production runs. The paper is clearly written and gives a useful description of the porting choices, including the use of shared memory, aligned memory, and warp-level reductions. The explicit comparison of four modular-multiplication implementations is also informative. However, the paper is a benchmark-driven implementation report, and its central quantitative claims currently rest on unverified correctness and a very narrow set of benchmark problems; the absence of code or raw data further limits reproducibility.

major comments (4)
  1. [Section 3, Tables 3-4] The paper never verifies that the GPU reconstruction produces the same output as the CPU version. All reported metrics are timings; a faster but incorrect port would produce identical tables. Please add a correctness check: for each of the four benchmark sizes, compare the full set of reconstructed coefficients (or the final rational functions) between the CPU and GPU implementations, and also validate the ported modular primitives against known-answer tests.
  2. [Section 2.1, Variant 2] The code computing pow = 2^64 mod n is incorrect for primes n > 2^63: (1ULL << 63) % n equals 2^63, and the subsequent pow << 1 overflows to 0 before the modulo operation, giving pow = 0 rather than 2^64 mod n. This concrete arithmetic error in a variant the paper benchmarks illustrates that timing-only fuzzing is insufficient validation of modular arithmetic; the authors should either correct the code or clearly label it as illustrative, and should validate all arithmetic used in the final implementation against the CPU reference.
  3. [Section 3.1] The statement that 'the computational complexity (and therefore the calculation time) does not depend much on the structure of the datasets, only on their size' is an unsupported assumption. Zippel reconstruction time depends on the sparsity pattern and the number of probes, not merely the monomial count. The benchmarks use only one physical problem, so the reported speedups may not generalize. Please provide evidence for the complexity claim, e.g., benchmarks on several problems with the same term counts but different monomial or skeleton structures.
  4. [Section 3.1] The CPU baseline is not fully specified. The text gives the Xeon Gold 6126 with 12 cores and INT64 peak, but does not state whether the original CPU version of FIRE uses a single thread or all 12 cores. If the CPU baseline is serial, the speedup numbers misrepresent the comparison against a typical multi-core CPU run. Please state the threading configuration and, if possible, compare against a multi-threaded CPU baseline.
minor comments (6)
  1. [Section 2.1, Table 1] There are typos: 'straightforwad' should be 'straightforward', and the variant name is inconsistently written as 'CUDMODP' in the table and 'CUMODP' in the text.
  2. [Section 2.2, Eq. (4)] The symbols C, hi, and n' in the balanced-coefficient formula are not defined in the text; please define them or refer to the exact definitions in [26].
  3. [Section 3.2] The paper states that 'confidence intervals are overlapping' but does not report any confidence intervals, standard deviations, or per-run times. Please include these statistics or raw timing data to support the claim of statistical significance.
  4. [Section 2.2] The word 'discributed' should be 'distributed'.
  5. [Section 3.2] The GPU load metric is described as 'preliminary information obtained not for all types of experiments'; this caveat should be stated in the abstract or conclusion as well, since the text later presents the >90% figure prominently.
  6. [Introduction] The claim of being the first GPU implementation of the Zippel algorithm would be more persuasive with a brief survey of existing GPU polynomial-arithmetic or sparse-interpolation libraries, to rule out prior work explicitly.

Circularity Check

0 steps flagged · score 0.0 of 10

Benchmark/implementation paper with no circular derivation; speedups are measured against an external CPU baseline.

full rationale

The paper's central claim is an engineering result: a GPU port of the Zippel/balanced-Zippel reconstruction with measured wall-clock speedups compared to a CPU version of FIRE. No result is fitted to data or defined in terms of its own output. The speedups in Tables 3 and 4 are comparisons against an external CPU implementation running the same algorithm, so they do not presuppose the GPU version's correctness or the paper's conclusions. The reliance on the authors' prior work [26] for the balanced Zippel formulation is a citation to previously published material, not a derivation that reduces to the present claim; the present paper contributes an implementation and benchmarks, not a new mathematical reconstruction theorem. The absence of a direct check that GPU and CPU produce identical reconstructed coefficients is a genuine validation gap, and the variant-2 modular multiplication code shown in Section 2.1 does overflow for primes larger than 2^63, but these are correctness risks, not circular reasoning: a wrong-but-fast port would be incorrect rather than circularly justified. No circular step can be exhibited where an equation or fitted parameter is renamed as a prediction. Therefore the appropriate circularity score is 0.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

No fitted parameters appear in this engineering paper. The central speedup claim is a measurement, not a fit, but its validity depends on the Zippel support-stability assumption, the prior balanced reconstruction from the authors' own [26], correctness of the FLINT port without explicit output validation, and standard modular reconstruction.

assumptions (4)
  • domain assumption The set of nonzero monomials in f(x1,...,x_{k-1}, c_i) is assumed to be the same as in the skeleton at c_0 for all probe values c_i.
    Core to Zippel; if a coefficient vanishes at a probe, the assumed support is wrong and reconstruction can fail or require extra probes. This paper inherits the assumption without new justification (Section 1, Eq. 3).
  • domain assumption The balanced numerator/denominator reconstruction from Smirnov and Zeng [26] produces correct rational functions.
    The GPU code implements this prior algorithm as part of the balanced Zippel method. The paper relies on this earlier published work without re-deriving it.
  • domain assumption The ported FLINT modular arithmetic functions (nmod_mul, NMOD_RED2, NMOD_RED3, etc.) are correct on GPU across operand ranges.
    Only fuzzing timing results are shown (Table 1), no correctness validation of the port is reported, and full reconstruction outputs are not checked against CPU.
  • standard math Large-prime modular arithmetic followed by rational reconstruction yields the true rational function.
    Standard technique in the modular Feynman integral reduction literature; the paper applies it as background methodology.

how reviews work

0 comments
Cite this review

Pith. "Pith review of GPU Implementation of Zippel Method for Feynman Integral Reconstruction." pith.science (2026). https://pith.science/paper/RIWJWQOL

@misc{pith2026250524256,
  author       = {Pith},
  title        = {Pith review of: GPU Implementation of Zippel Method for Feynman Integral Reconstruction},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/RIWJWQOL}},
  note         = {Machine review of arXiv:2505.24256}
}
read the original abstract

The Zippel algorithm performs a rational reconstruction of multivariate polynomials and aims specifically at the sparse case. It is applied in different fields of science, lately becoming an important step in Feynman integral reduction in elementary particle physics. In some cases with multiple variables it might become a bottleneck for the whole evaluation so that different optimizations are required. In this paper we describe how we ported the classical Zippel algorithm together with its balanced version for rational functions to graphical processor units and perform its evaluation on several GPUs.

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

28 extracted references · 14 canonical work pages

  1. [1]

    JHEP 07, 046 (2004).https://doi.org/10.1088/1126-6708/2004/07/046

    Anastasiou, C., Lazopoulos, A.: Automatic integral reduction for higher order perturbative calculations. JHEP 07, 046 (2004).https://doi.org/10.1088/1126-6708/2004/07/046

    doi:10.1088/1126-6708/2004/07/046 2004

  2. [2]

    Belitsky, A.V., Smirnov, A.V., Yakovlev, R.V.: Balancing act: Multivariate rational re- construction for IBP. Nucl. Phys. B 993, 116253 (2023). https://doi.org/10.1016/j. nuclphysb.2023.116253

    arXiv 2023

  3. [3]

    In: Proceedings of the Twentieth Annual ACM Symposium on Theory of Computing

    Ben-Or, M., Tiwari, P.: Adeterministicalgorithmforsparsemultivariatepolynomialinterpo- lation. In: Proceedings of the Twentieth Annual ACM Symposium on Theory of Computing. p. 301–309. STOC ’88, Association for Computing Machinery, New York, NY, USA (1988). https://doi.org/10.1145/62212.62241

    arXiv 1988

  4. [4]

    Chetyrkin, K.G., Tkachov, F.V.: Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops. Nucl. Phys. B 192, 159–204 (1981). https://doi.org/10.1016/ 0550-3213(81)90199-1

    work page 1981

  5. [5]

    JHEP 12, 140 (2022).https://doi.org/10.1007/JHEP12(2022)140

    De Laurentis, G., Page, B.: Ans¨ atze for scattering amplitudes from p-adic numbers and algebraic geometry. JHEP 12, 140 (2022).https://doi.org/10.1007/JHEP12(2022)140

    work page doi:10.1007/jhep12(2022)140 2022

  6. [6]

    In: Fukuda, K., van der Hoeven, J., Joswig, M., Takayama, N

    Hart, W.: Fast Library for Number Theory: An introduction. In: Fukuda, K., van der Hoeven, J., Joswig, M., Takayama, N. (eds.) Mathematical Software – ICMS 2010. Lecture Notes in Computer Science, vol. 6327, pp. 88–91. Springer (2010).https://doi.org/10. 1007/978-3-642-15582-6_15

    work page 2010

  7. [7]

    Hart, W., et al.: FLINT: Fast Library for Number Theory.http://www.flintlib.org/, version 2.9 (or your version), accessed May 2025

    work page 2025

  8. [8]

    In: Proceedings of the 2000 International Symposium on Symbolic and Algebraic Computation

    Kaltofen, E., Lee, W.s., Lobo, A.A.: Early termination in ben-or/tiwari sparse interpolation and a hybrid of zippel’s algorithm. In: Proceedings of the 2000 International Symposium on Symbolic and Algebraic Computation. p. 192–201. ISSAC ’00, Association for Computing Machinery, New York, NY, USA (2000).https://doi.org/10.1145/345542.345629

    work page arXiv 2000

Show all 28 references
  1. [9]

    Klappert, J., Lange, F.: Reconstructing rational functions with FireFly. Comput. Phys. Commun. 247, 106951 (2020).https://doi.org/10.1016/j.cpc.2019.106951

    2020

  2. [10]

    Klappert, J., Lange, F., Maierh¨ ofer, P., Usovitsch, J.: Integral reduction with Kira 2.0 and finite field methods. Comput. Phys. Commun. 266, 108024 (2021).https://doi.org/10. 1016/j.cpc.2021.108024

    2021

  3. [11]

    Lange, F., Usovitsch, J., Wu, Z.: Kira 3: integral reduction with efficient seeding and optimized equation selection (2025),https://arxiv.org/abs/2505.20197

    2025 arXiv

  4. [12]

    Laporta, S.: High precision calculation of multiloop Feynman integrals by difference equations. Int. J. Mod. Phys. A 15, 5087–5159 (2000). https://doi.org/10.1142/ S0217751X00002159

    2000

  5. [13]

    JHEP 07, 123 (2019).https://doi.org/10.1007/JHEP07(2019)123

    Laurentis, G., Maˆ ıtre, D.: Extracting analytical one-loop amplitudes from numerical evalu- ations. JHEP 07, 123 (2019).https://doi.org/10.1007/JHEP07(2019)123

    2019 doi

  6. [14]

    Lee, R.N.: Presenting LiteRed: a tool for the Loop InTEgrals REDuction (12 2012)

    2012

  7. [15]

    Lee, R.N.: LiteRed 1.4: a powerful tool for reduction of multiloop integrals. J. Phys. Conf. Ser. 523, 012059 (2014).https://doi.org/10.1088/1742-6596/523/1/012059

    2014 doi

  8. [16]

    Magerya, V.: Rational Tracer: a Tool for Faster Rational Function Reconstruction (11 2022)

    2022

  9. [17]

    Maierh¨ ofer, P., Usovitsch, J.: Kira 1.2 Release Notes (12 2018)

    2018

  10. [18]

    Com- put

    Maierh¨ ofer, P., Usovitsch, J., Uwer, P.: Kira—A Feynman integral reduction program. Com- put. Phys. Commun. 230, 99–112 (2018).https://doi.org/10.1016/j.cpc.2018.04.012

    2018 doi

  11. [19]

    von Manteuffel, A., Studerus, C.: Reduze 2 - Distributed Feynman Integral Reduction (1 2012)

    2012

  12. [20]

    von Manteuffel, A., Schabinger, R.M.: A novel approach to integration by parts reduction. Phys. Lett. B 744, 101–104 (2015).https://doi.org/10.1016/j.physletb.2015.03.029

    2015 doi

  13. [21]

    JHEP 12, 030 (2016).https://doi.org/10.1007/JHEP12(2016)030

    Peraro, T.: Scattering amplitudes over finite fields and multivariate functional reconstruc- tion. JHEP 12, 030 (2016).https://doi.org/10.1007/JHEP12(2016)030

    2016 doi

  14. [22]

    JHEP 07, 031 (2019).https://doi.org/10.1007/JHEP07(2019)031

    Peraro, T.: FiniteFlow: multivariatefunctionalreconstructionusingfinitefieldsanddataflow graphs. JHEP 07, 031 (2019).https://doi.org/10.1007/JHEP07(2019)031

    2019 doi

  15. [23]

    Smirnov, A.V., Chuharev, F.S.: FIRE6: Feynman Integral REduction with Modular Arith- metic. Comput. Phys. Commun. 247, 106877 (2020). https://doi.org/10.1016/j.cpc. 2019.106877

    2020

  16. [24]

    Smirnov, A.V., Smirnov, V.A.: FIRE4, LiteRed and accompanying tools to solve integration by parts relations. Comput. Phys. Commun. 184, 2820–2827 (2013).https://doi.org/10. 1016/j.cpc.2013.06.016

    2013

  17. [25]

    Smirnov, A.V.: FIRE5: a C++ implementation of Feynman Integral REduction. Comput. Phys. Commun. 189, 182–191 (2015).https://doi.org/10.1016/j.cpc.2014.11.024

    2015 doi

  18. [26]

    Smirnov, A.V., Zeng, M.: Feynman integral reduction: balanced reconstruction of sparse rationalfunctionsandimplementationonsupercomputersinaco-designapproach.Numerical Methods and Programming 25(Special issue), 30–45 (2024).https://doi.org/10.26089/ NumMet.2024s03

    2024

  19. [27]

    large scale, deep monitoring and fine analytics for the user community

    Voevodin, V., Antonov, A., Nikitenko, D., Shvets, P., Sobolev, S., Sidorov, I., Stefanov, K., Voevodin, V., Zhumatiy, S.: SupercomputerLomonosov-2, “large scale, deep monitoring and fine analytics for the user community”. Supercomputing Frontiers and Innovations 6(2), 4–11 (20...

    2019 doi

  20. [28]

    In: Ng, E.W

    Zippel, R.: Probabilistic algorithms for sparse polynomials. In: Ng, E.W. (ed.) Symbolic and Algebraic Computation. pp. 216–226. Springer Berlin Heidelberg, Berlin, Heidelberg (1979)

    1979

Pith tools

Reviewed August 7, 2026 · model on record in the stance chip above.