[Paper Review] An O(M(n) log n) algorithm for the Jacobi symbol

Richard P. Brent, Paul Zimmermann|arXiv (Cornell University)|Apr 13, 2010

Algorithms and Data Compression2 references4 citations

TL;DR

This paper presents an $O(M(n)"log n)$ algorithm for computing the Jacobi symbol of two $n$-bit integers using a binary recursive GCD approach, avoiding reliance on continued fractions or Gauss's identity. The method achieves subquadratic time complexity by efficiently handling quotient sequences and is shown to outperform GMP's quadratic implementation for inputs larger than ~10,000 decimal digits.

ABSTRACT

The best known algorithm to compute the Jacobi symbol of two n-bit integers runs in time O(M(n) log n), using Schönhage's fast continued fraction algorithm combined with an identity due to Gauss. We give a different O(M(n) log n) algorithm based on the binary recursive gcd algorithm of Stehlé and Zimmermann. Our implementation - which to our knowledge is the first to run in time O(M(n) log n) - is faster than GMP's quadratic implementation for inputs larger than about 10000 decimal digits.

Motivation & Objective

• To develop a subquadratic algorithm for computing the Jacobi symbol of two $n$-bit integers.
• To avoid dependence on Schönhage's continued fraction algorithm or Gauss's identity, offering a simpler alternative.
• To achieve an $O(M(n)\log n)$ time complexity using only binary GCD techniques.
• To implement and benchmark the algorithm, showing performance gains over existing quadratic implementations for large inputs.
• To provide a provably efficient and practical method for high-performance number-theoretic computation.

Proposed method

• Adapts the binary recursive GCD algorithm of Stehlé and Zimmermann to compute the Jacobi symbol directly.
• Uses binary division with positive quotient to maintain positivity and ensure termination.
• Introduces a modified stopping criterion to allow the algorithm to terminate despite using positive quotients.
• Employs a symbolic transformation of the Jacobi symbol using identities: $(bc|a)=(b|a)(c|a)$, $(2|a)=(-1)^{(a^2-1)/8}$, and $(b|a)=(-1)^{(a-1)(b-1)/4}(a|b)$ for odd $a,b$.
• Classifies iterations as good, bad, or ugly and combines sequences of 'ugly' iterations into a single 'harmless' step to reduce worst-case complexity to $O(M(n)\log n)$.

Experimental results

Research questions

• RQ1 Can the Jacobi symbol be computed in $O(M(n)\log n)$ time using only binary GCD techniques, without relying on continued fractions or Gauss's identity?
• RQ2 What is the worst-case time complexity of a binary GCD-based algorithm for the Jacobi symbol, and can it be reduced to $O(M(n)\log n)$?
• RQ3 How does the performance of the new algorithm compare to existing implementations like GMP’s quadratic algorithm for large inputs?
• RQ4 What is the average-case bit reduction per iteration in the proposed algorithm, and how does it compare to standard binary GCD?
• RQ5 Can the algorithm be implemented efficiently in practice, and at what input size does it outperform existing methods?

Key findings

• The proposed algorithm achieves an $O(M(n)\log n)$ time complexity, matching the best-known theoretical bound.
• The algorithm is the first to achieve this bound without using Schönhage’s continued fraction algorithm or Gauss’s identity.
• The implementation outperforms GMP’s quadratic implementation for inputs larger than approximately 10,300 decimal digits (34,240 bits).
• For inputs of one million 64-bit words, the algorithm takes about 83.1 seconds and 47.5 million iterations, with a theoretical bit reduction rate of 0.651993 bits per iteration.
• Empirical results show that about 53.33% of iterations are good, 26.67% bad, and 20% harmless, converging toward the expected heuristic distribution.
• The break-even point between the $O(n^2)$ GMP implementation and the new $O(M(n)\log n)$ implementation occurs at around 535 words (~34,240 bits or 10,300 decimal digits).