Cycle-Stride Evaluation of Finite Fractional Fourier Orbits over Dyadic-Symmetry-Complete Prime Fields

1. Introduction

The finite fractional Fourier transform (FrFT) over a prime field ${\mathbb{F}}_p$ generalizes the discrete Fourier transform by interpreting it as one step of a cyclic rotation in representation space: the ordinary Fourier transform occupies a quarter-turn position in a four-cycle whose other cardinal points are identity, parity, and the inverse Fourier transform. Between these cardinal positions sit intermediate fractional transforms indexed by a parameter s in a cyclic group. The continuum FrFT has well-established applications in optics, time-frequency signal analysis, and quantum mechanics [1,2]; finite-field constructions have been developed for discrete signal processing [3,4] and, more recently, for qudit-based quantum circuit implementations [5].

For primes $p\equiv 1(mod4)$, the field ${\mathbb{F}}_p$ contains a primitive fourth root of unity — a square root of $-1$ — and the finite Fourier matrix admits an exact normalization with the four-cycle structure of the classical FrFT. This setting is convenient because the four-cycle, the projector decomposition onto eigenspaces, and the polynomial form of the fractional powers all transfer cleanly from the continuum theory. The price is that the cycle length $n=p-1$ is constrained to $4t$ for some integer t, and the standard radix-2 Cooley–Tukey FFT does not in general apply.

This paper makes three contributions.

• Output-optimal full-orbit streaming.

After computing the four spectral projections $u_{\ell }={\Pi }_{\ell }v$ of an input signal $v\in {\mathbb{F}}_p^n$ ($n=p-1$), the entire orbit $\{F^{\left[s\right]}v:s\in {\mathbb{Z}}_n\}$ is recovered by a streaming procedure performing $O\left(n\right)$ field operations per output. Total streaming cost is $O\left(n^2\right)$, which matches the output-size lower bound up to constant factors: the orbit consists of n vectors each of length n, so any algorithm that materializes the full orbit must perform $\Omega \left(n^2\right)$ output operations. This result holds for every prime $p=4t+1$ and every choice of primitive generator $e_t\in {\mathbb{F}}_p^{\times }$. We refer to such primes as symmetry-complete throughout, the condition simply being that ${\mathbb{F}}_p$ contains a primitive fourth root of unity and the normalized finite Fourier operator admits the four-cycle $F_t^4=I$.

• Dyadic-symmetry-complete primes and shift-based streaming.

We isolate the class of primes for which the streaming phase admits a particularly hardware-friendly implementation. Call p dyadic-symmetry-complete (DSC) if $p\equiv 1(mod4)$ and 2 is a primitive root of ${\mathbb{F}}_p^{\times }$. For these primes, the streaming phase can be reformulated using elementwise contribution arrays whose updates between meridian steps consist of multiplications by $2,4,8$ modulo p — which are bit-shifts followed by conditional reductions, requiring no general modular multiplier. We establish two structural facts: every DSC prime satisfies $p\equiv 5(mod8)$ (because 2 being primitive forces 2 to be a quadratic non-residue), and among Fermat primes only $p=5$ is DSC (because no other Fermat prime has 2 as a primitive root).

• Output optimality interpretation.

The full-orbit streaming cost of $O\left(n^2\right)$ matches the output-size lower bound. It also matches, up to constants, the cost of the most natural NTT-based alternative: precompute v, $F_{t}v$, $Jv$, $F_t^{-1}v$ (four transforms in $O(nlogn)$ each under any applicable NTT), assemble the projector components, and stream the orbit by recombination in $O\left(n^2\right)$. Our contribution is therefore not an asymptotic speedup over all alternatives but rather a unified streaming construction valid for every symmetry-complete prime, with a shift-only inner loop in the DSC subclass that we believe is novel.

1.1. Relation to Existing Work

The eigenstructure-based discretization of the FrFT was introduced by Pei, Yeh, and Tseng [4] and by Candan, Kutay, and Ozaktas [3], who used the projector decomposition to define canonical discrete fractional powers of the Fourier matrix. Our use of the same projectors as a precomputed decomposition enabling streaming evaluation across the full meridian cycle does not appear in those references, which focus on defining a single fractional transform rather than enumerating the orbit. The arithmetic FrFT of Floratos and Pavlidis [5] constructs discrete fractional Fourier families adapted to qudit quantum circuits.

The number-theoretic transform literature, beginning with Pollard [6] and now applied widely in lattice-based cryptography [7,8], exploits the fast modular arithmetic of specific prime fields. Fermat primes $p=2^{2^k}+1$ admit radix-2 NTTs of length a power of two, and the property “multiplication by 2 is a bit-shift modulo p” is used as an implementation optimization for the twiddle-factor updates. Crucially, in ${\mathbb{F}}_{2^{2^k}+1}$ the element 2 has order $2^{k+1}$ rather than $p-1$, so it is not a primitive root of ${\mathbb{F}}_p^{\times }$ for $k\ge 2$. The class of primes we identify — where 2 is primitive — is therefore distinct from the Fermat primes (with the exception of $p=5$, the only common member).

For prime-length DFTs more generally, Rader’s algorithm [9] and Bluestein’s chirp-z transform [10] provide $O(nlogn)$ alternatives when n is prime or has an inconvenient factorization. The cycle-stride streaming we develop here applies after a single setup phase regardless of how that phase is implemented; in particular, the post-projection streaming cost does not depend on whether the projections were obtained by direct matrix multiplication, by an NTT-based method, or by any other route.

1.2. Organization

Section 2 establishes the algebraic setting and proves the four-cycle and projector identities from first principles. Section 3 states and proves the general streaming theorem (Theorem 1) valid for every symmetry-complete prime. Section 4 introduces DSC primes, proves the congruence and Fermat-prime separation results (Theorem 2), and gives the contribution-array reformulation (Theorem 3). Section 5 analyzes complexity in operation counts. Section 6 reports empirical validation. Section 7 discusses applications and limitations.

2. Preliminaries

2.1. Setting and Notation

Throughout, p is a prime with $p\equiv 1(mod4)$. Write $p=4t+1$ and $n=p-1=4t$. The multiplicative group ${\mathbb{F}}_p^{\times }$ is cyclic of order n, and we fix a primitive generator $e_t\in {\mathbb{F}}_p^{\times }$. Set

$$i_t:=e_t^t\in {\mathbb{F}}_p^{\times }.$$

Let ${\Phi }_t:={\mathbb{Z}}_n$ and $V_t:={\mathbb{F}}_p^{{\Phi }_t}$, the space of functions $v:{\Phi }_t\to {\mathbb{F}}_p$ (equivalently, n-vectors over ${\mathbb{F}}_p$ indexed by ${\Phi }_t$).

2.2. The Finite Fourier Matrix

Definition 1

(Shell Fourier matrix). The shell Fourier matrix $W_t:V_t\to V_t$ has entries

$${\left(W_t\right)}_{k,j}:=e_t^{jk},j,k\in {\Phi }_t.$$

Definition 2

(Reversal operator). The reversal operator $J:V_t\to V_t$ acts by ${\left(Jv\right)}_k:=v_{-k}$ (indices mod n). Equivalently $J_{k,\ell }={\delta }_{k,-\ell }$.

It is immediate that $J^2=I$ and that J has trace equal to the number of indices $k\in {\Phi }_t$ with $k\equiv -k(modn)$, namely $k\in \{0,2t\}$.

Lemma 1

(Squared Fourier identity). $i_t^2=-1$ in ${\mathbb{F}}_p$, and

$$W_t^2=-J.$$

Proof. 

First claim. Since $e_t$ is a primitive generator of ${\mathbb{F}}_p^{\times }$, we have $e_t^n=1$ and $e_t^{n/2}\ne 1$. The element $e_t^{n/2}$ squares to 1, so it equals either 1 or $-1$; not equal to 1, hence equal to $-1$. Therefore $i_t^2=e_t^{2t}=e_t^{n/2}=-1$.

Second claim. For $k,\ell \in {\Phi }_t$, the $(k,\ell )$ entry of $W_t^2$ is

$${\left(W_t^2\right)}_{k,\ell }=\sum _{j\in {\Phi }_t}{e}_t^{jk}·e_t^{j\ell }={\displaystyle \sum _{j=0}^{n-1}}{e}_t^{j(k+\ell )}.$$

If $k+\ell \equiv 0(modn)$, every summand equals 1 and the sum is $n=p-1=-1$ in ${\mathbb{F}}_p$. Otherwise $e_t^{k+\ell }\ne 1$ and the finite geometric sum is $(e_t^{n(k+\ell )}-1)/(e_t^{k+\ell }-1)=0$. Hence ${\left(W_t^2\right)}_{k,\ell }=-{\delta }_{\ell ,-k}=-J_{k,\ell }$.    □

2.3. The Normalized Fourier Operator and Its Four-Cycle

Definition 3

(Normalized shell Fourier operator). Thenormalized shell Fourier operator is

$$F_t:=i_t{W}_t.$$

Proposition 1

(Four-cycle). $F_t^2=J$ and $F_t^4=I$.

Proof. 

By Lemma 1, $F_t^2=i_t^2{W}_t^2=(-1)(-J)=J$. Squaring, $F_t^4=J^2=I$.    □

The four-cycle structure is the algebraic skeleton on which the finite FrFT is built. The four cardinal points are

$$F_t^0=I,F_t^1=F_t\left(\mathrm{Fourier}\right),F_t^2=J\left(\mathrm{parity}\right),F_t^3=F_t^{-1}(\mathrm{inverse}\mathrm{Fourier}).$$

2.4. Eigenspace Projectors

The eigenvalues of $F_t$ are the four roots of $x^4=1$ in ${\mathbb{F}}_p$, namely $\{1,i_t,-1,-i_t\}=\{i_t^{\ell }:\ell =0,1,2,3\}$. Each eigenvalue’s algebraic multiplicity is the dimension of the corresponding eigenspace, and the projector onto eigenspace ℓ is given by the standard polynomial formula in $F_t$.

Definition 4

(Spectral projectors). For $\ell \in \{0,1,2,3\}$, define

$${\Pi }_{\ell }:={\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{-\ell r}{F}_t^r.$$

The scalar $1/4\in {\mathbb{F}}_p$ denotes the multiplicative inverse of 4; since p is odd prime with $p>2$, this is well-defined.

Proposition 2

(Projector identities). The operators ${\Pi }_{\ell }$ satisfy:

(i) 

${\Pi }_{\ell }{\Pi }_m={\delta }_{\ell m}{\Pi }_{\ell }$;

(ii) 

${\sum }_{\ell =0}^3{\Pi }_{\ell }=I$;

(iii) 

$F_t{\Pi }_{\ell }=i_t^{\ell }{\Pi }_{\ell }$.

Proof. 

The proof uses the discrete orthogonality identity over ${\mathbb{Z}}_4$: for $a,b\in \{0,1,2,3\}$,

$${\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{r(a-b)}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}a\equiv b(mod4),\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

which follows from $i_t^4=1$ and the same geometric-sum argument used in Lemma 1.

(iii): Multiplying ${\Pi }_{\ell }$ by $F_t$ on the left and using $F_t^4=I$:

$$F_t{\Pi }_{\ell }={\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{-\ell r}{F}_t^{r+1}={\displaystyle \frac{1}{4}}{\displaystyle \sum _{r^{\prime }=1}^4}{i}_t^{-\ell (r^{\prime }-1)}{F}_t^{r^{\prime }}=i_t^{\ell }·{\displaystyle \frac{1}{4}}{\displaystyle \sum _{r^{\prime }=0}^3}{i}_t^{-\ell r^{\prime }}{F}_t^{r^{\prime }}=i_t^{\ell }{\Pi }_{\ell },$$

where we used $i_t^{-\ell ·4}=1$ to reindex.

(ii):${\sum }_{\ell }{\Pi }_{\ell }={\displaystyle \frac{1}{4}}{\sum }_{r=0}^3\left({\sum }_{\ell =0}^3{i}_t^{-\ell r}\right)F_t^r$. By the orthogonality identity with $a=0,b=r$, the inner sum is 4 for $r=0$ and 0 otherwise. So ${\sum }_{\ell }{\Pi }_{\ell }=F_t^0=I$.

(i): Using (iii) and the formula for ${\Pi }_m$,

$${\Pi }_{\ell }{\Pi }_m={\Pi }_{\ell }·{\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{-mr}{F}_t^r={\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{-mr}{\left(i_t^{\ell }\right)}^r{\Pi }_{\ell }=\left({\displaystyle \frac{1}{4}}{\displaystyle \sum _{r=0}^3}{i}_t^{r(\ell -m)}\right){\Pi }_{\ell }={\delta }_{\ell m}{\Pi }_{\ell }.$$

   □

2.5. The Finite FrFT Family

Definition 5

(Finite fractional Fourier transform). For $s\in {\Phi }_t$,

$$F^{\left[s\right]}:={\displaystyle \sum _{\ell =0}^3}{e}_t^{\ell s}{\Pi }_{\ell }.$$

Proposition 3

(FrFT properties). The family $\{F^{\left[s\right]}:s\in {\Phi }_t\}$ satisfies:

(i) 

$F^{[s+r]}=F^{\left[s\right]}{F}^{\left[r\right]}$ for all $s,r\in {\Phi }_t$.

(ii) 

$F^{\left[0\right]}=I$, $F^{\left[t\right]}=F_t$, $F^{\left[2t\right]}=J$, $F^{\left[3t\right]}=F_t^{-1}$, $F^{\left[4t\right]}=F^{\left[0\right]}=I$.

Proof. 

(i): Using the orthogonality of the projectors,

$$F^{\left[s\right]}{F}^{\left[r\right]}=\left(\sum _{\ell }{e}_t^{\ell s}{\Pi }_{\ell }\right)\left(\sum _m{e}_t^{mr}{\Pi }_m\right)=\sum _{\ell }{e}_t^{\ell (s+r)}{\Pi }_{\ell }=F^{[s+r]}.$$

(ii): For $s=0$: $F^{\left[0\right]}={\sum }_{\ell }{\Pi }_{\ell }=I$ by Proposition 2(ii). For $s=t$: $e_t^{\ell t}=i_t^{\ell }$, so $F^{\left[t\right]}={\sum }_{\ell }{i}_t^{\ell }{\Pi }_{\ell }$. Using $F_t{\Pi }_{\ell }=i_t^{\ell }{\Pi }_{\ell }$ and ${\sum }_{\ell }{\Pi }_{\ell }=I$, this equals $F_t$. For $s=2t$: $e_t^{\ell ·2t}=i_t^{2\ell }={(-1)}^{\ell }$, and $F^{\left[2t\right]}={\sum }_{\ell }{(-1)}^{\ell }{\Pi }_{\ell }=F_t^2=J$. For $s=3t$: $e_t^{3\ell t}=i_t^{3\ell }=i_t^{-\ell }$, giving $F_t^{-1}$. For $s=4t$: $e_t^{4\ell t}=1$, recovering ${\sum }_{\ell }{\Pi }_{\ell }=I$.    □

This is the algebraic setup we use throughout. We will return to the question of when the map $s\mapsto F^{\left[s\right]}$ is injective (faithfulness) in Section 7; for the algorithmic results below, what matters is the projector identity (2) below.

3. The Cycle-Stride Streaming Algorithm

3.1. The Projector-Orbit Identity

The central observation of this paper is the following identity, which follows directly from Definition 5.

Theorem 1

(Projector-orbit identity; Theorem A). Let $p=4t+1$ be prime and $e_t\in {\mathbb{F}}_p^{\times }$ a primitive generator. For any signal $v\in V_t$, set

$$u_{\ell }:={\Pi }_{\ell }v,\ell \in \{0,1,2,3\}.$$

(1)

Then for every $s\in {\Phi }_t$,

$$F^{\left[s\right]}v=u_0+e_t^s{u}_1+e_t^{2s}{u}_2+e_t^{3s}{u}_{3}in{V}_t.$$

(2)

Consequently, after the projections ${\left\{u_{\ell }\right\}}_{\ell =0}^3$ have been computed, the entire orbit $\{F^{\left[s\right]}v:s\in {\Phi }_t\}$ can be evaluated in $O\left(n^2\right)$ field operations and $O\left(n\right)$ working memory beyond the output stream.

Proof. 

By Definition 5,

$$F^{\left[s\right]}v={\displaystyle \sum _{\ell =0}^3}{e}_t^{\ell s}{\Pi }_{\ell }v={\displaystyle \sum _{\ell =0}^3}{e}_t^{\ell s}{u}_{\ell },$$

which is (2). The cost analysis: there are n values of s, and for each, computing $u_0+e_t^s{u}_1+e_t^{2s}{u}_2+e_t^{3s}{u}_3$ requires $3n$ scalar multiplications and $3n$ additions in ${\mathbb{F}}_p$, plus a constant number of operations to update the powers $e_t^s,e_t^{2s},e_t^{3s}$ between successive values of s. The total is $O\left(n^2\right)$ field operations. Working memory beyond the output stream is the four arrays $u_0,u_1,u_2,u_3$, totalling $4n$ field elements, plus $O\left(1\right)$ coefficient scalars.    □

Remark 1

(Output-size optimality). The full orbit ${\left\{F^{\left[s\right]}v\right\}}_{s\in {\Phi }_t}$ consists of n vectors in $V_t$, each of length n, totalling $n^2$ field elements. Any algorithm that explicitly materializes or streams every element of the full orbit must therefore perform $\Omega \left(n^2\right)$ output operations. The cycle-stride streaming phase matches this output-size lower bound up to constant factors: it performs $O\left(n^2\right)$ field operations and produces $n^2$ output entries.

3.2. Algorithm A: Scalar-Coefficient Streaming

The most direct realization of Theorem 1 maintains three scalar coefficients $c_1,c_2,c_3$ representing $e_t^s,e_t^{2s},e_t^{3s}$, and updates them multiplicatively between meridian steps.

This algorithm performs the projector setup in $O\left(n^2\right)$ field operations and the streaming phase in $O\left(n^2\right)$ field operations (with $3n$ multiplications and $3n$ additions per meridian step, plus three scalar updates). The cost is dominated by general modular multiplications in the inner loop — the per-output term $c_{\ell }{u}_{\ell }\left[j\right]$ is a full modular multiply regardless of the value of $e_t$.

4. Dyadic-Symmetry-Complete Primes

We now isolate a class of primes for which the streaming phase admits a structurally cleaner implementation.

4.1. Definition and Structural Characterization

Definition 6

(Dyadic-symmetry-complete prime). A prime p is dyadic-symmetry-complete (DSC) if $p\equiv 1(mod4)$ and 2 is a primitive root of ${\mathbb{F}}_p^{\times }$.

Theorem 2

(Necessary congruence and Fermat-prime separation; Theorem C). Let p be a DSC prime.

(i) 

$p\equiv 5(mod8)$.

(ii) 

If $p-1$ is a power of 2, then $p=5$.

(iii) 

Among the known Fermat primes $\{3,5,17,257,65537\}$, only $p=5$ is DSC.

Proof. 

(i): If 2 is a primitive root of ${\mathbb{F}}_p^{\times }$, then 2 has order $p-1$, in particular order strictly greater than $(p-1)/2$. Hence $2^{(p-1)/2}\ne 1(modp)$. By Euler’s criterion, $2^{(p-1)/2}\equiv \left(\frac{2}{p}\right)(modp)$, where $\left(\frac{2}{p}\right)$ is the Legendre symbol. Therefore $\left(\frac{2}{p}\right)=-1$, i.e., 2 is a quadratic non-residue mod p.

The second supplementary law of quadratic reciprocity states $\left(\frac{2}{p}\right)={(-1)}^{(p^2-1)/8}$, equivalently $\left(\frac{2}{p}\right)=1$ iff $p\equiv \pm 1(mod8)$ and $\left(\frac{2}{p}\right)=-1$ iff $p\equiv \pm 3(mod8)$ [11]. Combining with $p\equiv 1(mod4)$, the only admissible residue class is $p\equiv 5(mod8)$.

(ii): Suppose $p>5$ is DSC and $p-1=2^N$ for some $N\ge 2$. Then $p=2^N+1$ is a prime of the form $2^N+1$, which forces N itself to be a power of 2: indeed, if N has an odd factor $d>1$, writing $N=d{N}^{\prime }$, then $2^N+1={\left(2^{N^{\prime }}\right)}^d+1$ is divisible by $2^{N^{\prime }}+1>1$, contradicting primality. Hence $N=2^k$ for some $k\ge 1$, and p is a Fermat prime $F_k=2^{2^k}+1$.

In ${\mathbb{F}}_p$ we have $2^{2^k}=2^N=p-1\equiv -1(modp)$, so $2^{2·2^k}=1$. Hence the order of 2 in ${\mathbb{F}}_p^{\times }$ divides $2·2^k=2^{k+1}$; and the order does not divide $2^k$ (since $2^{2^k}=-1\ne 1$). Therefore ${ord}_p\left(2\right)=2^{k+1}$. For 2 to be primitive, we need ${ord}_p\left(2\right)=p-1=2^{2^k}$, i.e., $k+1=2^k$. This equation holds at $k=1$ (giving $p=5$) but at no $k\ge 2$ (since $2^k$ grows much faster than $k+1$). Hence $p\le 5$, contradicting $p>5$.

(iii): For each Fermat prime, one verifies directly: in ${\mathbb{F}}_3$, $2=-1$ has order $2=p-1$, but $p=3$ does not satisfy $p\equiv 1(mod4)$. In ${\mathbb{F}}_5$, 2 has order $4=p-1$ and $p\equiv 1(mod4)$, so $p=5$ is DSC. In ${\mathbb{F}}_{17}$, ${\mathbb{F}}_{257}$, ${\mathbb{F}}_{65537}$, the order of 2 is $2^{k+1}$ for $k=2,3,4$, namely $8,16,32$, far below $p-1$. None of these is DSC.    □

Corollary 1

(Structural separation from Fermat-prime NTT). The DSC primes and Fermat primes (other than $p=5$) are disjoint classes. In particular, the dyadic shift property used in Fermat-prime NTT implementations [6,8] — where “multiplication by 2 is bit-shift” refers to twiddle-factor update for non-primitive 2 — is structurally different from the property exploited in this paper, where 2 is a primitive generator and powers of 2 traverse all of ${\mathbb{F}}_p^{\times }$.

Remark 2

(Effective search class). By Theorem 2(i), every DSC prime satisfies $p\equiv 5(mod8)$. The effective search space for DSC primes is therefore the arithmetic progression $\{8q+5:q\in {\mathbb{Z}}_{\ge 0}\}$ intersected with the primes; within this class one then tests whether 2 is a primitive root. The empirical density results below (Proposition 4) are reported relative to this effective class.

4.2. Density: Empirical and Conjectural

The density of DSC primes within primes $p\equiv 5(mod8)$ is governed heuristically by Artin’s conjecture on primitive roots, conditional on the Generalized Riemann Hypothesis for certain Dedekind zeta functions [12]. Artin’s constant

$$C_{\mathrm{Artin}}=\prod _{q\mathrm{prime}}\left(1-{\displaystyle \frac{1}{q(q-1)}}\right)\approx 0.3739558$$

governs the global density of primes for which 2 is a primitive root. Restricted to the congruence class $p\equiv 5(mod8)$, the conjectural density may differ; we are not aware of a closed form, but the heuristic and the empirical data agree that it is bounded away from zero.

Unconditionally, the deepest known result is Heath-Brown’s theorem [13]: at least one of the integers $2,3,5$ is a primitive root for infinitely many primes. This does not single out 2, so it does not unconditionally establish the infinitude of DSC primes. We therefore record the existence statement as conjectural and the algorithm’s applicability as empirically demonstrated for the explicit examples enumerated below.

Proposition 4

(Empirical density at small p). Among the 43 primes $p\le 1000$ with $p\equiv 5(mod8)$, exactly 36 are DSC, an empirical density of $0.837$. The first thirty DSC primes are

$$\begin{array}{c}\hfill 5,13,29,37,53,61,101,149,173,181,197,269,293,317,349,\\ \hfill 373,389,421,461,509,541,557,613,653,661,677,701,709,757,773.\end{array}$$

The empirical density within $p\equiv 5(mod8)$ is notably higher than the global Artin constant. The discrepancy is consistent with the heuristic that restricting to the congruence class removes the primes where 2 must be a quadratic residue, and is in line with conditional results on primitive roots in arithmetic progressions surveyed in [14].

4.3. Bit-shift Multiplication Modulo P

Lemma 2

(Bit-shift multiplication). For any prime p and any $a\in \{0,1,\dots ,p-1\}$, the operation $a\mapsto 2a(modp)$ is realized by:

(1) 

left-shift a by one bit to obtain $2a$ as an integer in $\{0,\dots ,2p-2\}$;

(2) 

if $2a\ge p$, subtract p.

On a register of $\lceil {log}_{2}p\rceil +1$ bits, the operation requires one shift and one conditional subtraction. Iterated k times, the operation realizes $a\mapsto 2^{k}a(modp)$ using k shifts and at most k conditional subtractions.

This is elementary but worth recording. The DSC condition ensures that iterating this operation generates all of ${\mathbb{F}}_p^{\times }$ from the starting value 1.

4.4. Algorithm B: Contribution-Array Streaming for DSC Primes

The scalar-coefficient streaming of Algorithm 1 performs three general modular multiplications per output element ($c_{\ell }{u}_{\ell }\left[j\right]$). For DSC primes with $e_t=2$, this can be restructured to eliminate general modular multiplications from the streaming phase entirely.

Algorithm 1 Scalar-coefficient streaming (Algorithm A)

Require: Prime $p=4t+1$; primitive $e_t\in {\mathbb{F}}_p^{\times }$; signal $v\in {\mathbb{F}}_p^n$; precomputed projectors ${\Pi }_0,{\Pi }_1,{\Pi }_2,{\Pi }_3$ Ensure: The orbit $\{F^{\left[s\right]}v:s=0,1,\dots ,n-1\}$   1: for $\ell =0,1,2,3$ do   2:     $u_{\ell }\leftarrow {\Pi }_{\ell }v$   3: end for   4: $c_1,c_2,c_3\leftarrow 1,1,1$   5: ${\mu }_1,{\mu }_2,{\mu }_3\leftarrow e_t,e_t^2,e_t^3(modp)$                   ▹ Per-step multipliers   6: for $s=0,1,\dots ,n-1$ do   7:     for $j=0,1,\dots ,n-1$ do   8:         $w\left[j\right]\leftarrow u_0\left[j\right]+c_1{u}_1\left[j\right]+c_2{u}_2\left[j\right]+c_3{u}_3\left[j\right](modp)$   9:     end for 10:     output w as $F^{\left[s\right]}v$ 11:     $c_1\leftarrow c_1{\mu }_{1}modp$ 12:     $c_2\leftarrow c_2{\mu }_{2}modp$ 13:     $c_3\leftarrow c_3{\mu }_{3}modp$ 14: end for

The reformulation maintains contribution arrays

$$a_{\ell }^{\left(s\right)}:=e_t^{\ell s}{u}_{\ell }\in V_t,\ell =1,2,3,$$

which absorb the scalar coefficient into the vector. The recombination at each meridian step is then an additive expression with no scalar-vector multiplications, and the update between meridian steps is an elementwise multiplication of each $a_{\ell }$ by $e_t^{\ell }\in \{2,4,8\}$ — which by Lemma 2 is a sequence of shifts.

Theorem 3

(Shift-based streaming over DSC primes; Theorem B). Algorithm 2 is correct: it outputs $F^{\left[s\right]}v$ at step s for every $s\in {\Phi }_t$. The streaming phase after the projection setup performs no general modular multiplications; the inner-loop operations consist entirely of:

• modular additions (in the output recombination);
• left bit-shifts and conditional subtractions (in the contribution-array updates).

Per meridian step, Algorithm 2 performs $3n$ modular additions to form the output $w=u_0+a_1+a_2+a_3$, and updates the three contribution arrays by multiplication by $2,4,8$ modulo p. Implemented by repeated doubling, the three updates require $6n$ elementary shift-reduce operations per meridian step. Equivalently, if multiplication by $2,4,8$ is implemented as fixed dyadic shift-reduce primitives, the updates consist of $3n$ fixed dyadic modular updates per step. In either interpretation, the streaming inner loop contains no general modular multiplication.

Algorithm 2 Contribution-array streaming for DSC primes (Algorithm B)

Require: DSC prime p; signal $v\in {\mathbb{F}}_p^n$; precomputed projectors ${\Pi }_0,{\Pi }_1,{\Pi }_2,{\Pi }_3$ Ensure: The orbit $\{F^{\left[s\right]}v:s=0,1,\dots ,n-1\}$   1: for $\ell =0,1,2,3$ do   2:     $u_{\ell }\leftarrow {\Pi }_{\ell }v$   3: end for   4: $a_1\leftarrow u_1$; $a_2\leftarrow u_2$; $a_3\leftarrow u_3$                                         ▹ Contribution arrays   5: for $s=0,1,\dots ,n-1$ do   6:     for $j=0,1,\dots ,n-1$ do   7:         $w\left[j\right]\leftarrow u_0\left[j\right]+a_1\left[j\right]+a_2\left[j\right]+a_3\left[j\right](modp)$   8:     end for   9:     output w as $F^{\left[s\right]}v$ 10:     for $j=0,1,\dots ,n-1$ do 11:           $a_1\left[j\right]\leftarrow 2a_1\left[j\right]modp$                                                       ▹ One shift 12:           $a_2\left[j\right]\leftarrow 4a_2\left[j\right]modp$                                                     ▹ Two shifts 13:           $a_3\left[j\right]\leftarrow 8a_3\left[j\right]modp$                                                  ▹ Three shifts 14:     end for 15: end for

Proof. 

We prove correctness by induction on s. The invariant maintained at the beginning of meridian step s is $a_{\ell }^{\left(s\right)}=2^{\ell s}{u}_{\ell }$ for $\ell =1,2,3$. Base case ($s=0$):$a_{\ell }^{\left(0\right)}=u_{\ell }=2^0{u}_{\ell }$. Inductive step: After the output at step s, line 11 updates $a_{\ell }$ to $2^{\ell }{a}_{\ell }^{\left(s\right)}=2^{\ell }·2^{\ell s}{u}_{\ell }=2^{\ell (s+1)}{u}_{\ell }=a_{\ell }^{(s+1)}$.

The output at step s is

$$w=u_0+a_1^{\left(s\right)}+a_2^{\left(s\right)}+a_3^{\left(s\right)}=u_0+2^s{u}_1+2^{2s}{u}_2+2^{3s}{u}_3,$$

which equals $F^{\left[s\right]}v$ by Theorem 1 applied with $e_t=2$.

Operation count: the recombination $w\left[j\right]=u_0\left[j\right]+a_1\left[j\right]+a_2\left[j\right]+a_3\left[j\right]$ contains three modular additions per index, hence $3n$ modular additions per meridian step. The contribution-array updates $a_1\leftarrow 2a_1$, $a_2\leftarrow 4a_2$, $a_3\leftarrow 8a_3$, each applied elementwise to a length-n array, require $1+2+3=6$ elementary shift-reduce operations per index under the repeated-doubling implementation, for $6n$ shift-reduce operations per meridian step. Under the alternative interpretation in which multiplication by $2,4,8$ modulo p is treated as a primitive dyadic operation, each contribution update is one dyadic modular update per index, for $3n$ such updates per step. Neither interpretation involves a general modular multiplication. □

Remark 3

(Comparison of Algorithms A and B). Algorithm 1 performs $3n$ general modular multiplications per meridian step in the inner loop. Algorithm 2 replaces these with elementwise multiplications by the fixed small constants $2,4,8$ modulo p, which decompose into $6n$ single-bit shift-reduce operations under repeated doubling or, equivalently, $3n$ fixed dyadic modular updates if these constants are implemented as dedicated shift-reduce primitives. On hardware where general modular multiplication is more expensive than shift-reduce — typically FPGA, ASIC, and microcontroller targets without hardware multipliers — Algorithm 2 is preferable. On a modern x86 or ARM CPU with single-cycle 64-bit integer multiplication, multiplication by a small constant is no cheaper than general multiplication, and the two algorithms have similar wall-clock cost; we report this explicitly in Section 6.

5. Complexity Analysis

We compare four methods for computing the full FrFT orbit ${\left\{F^{\left[s\right]}v\right\}}_{s\in {\Phi }_t}$ over a symmetry-complete prime $p=4t+1$ with $n=p-1$. Let $\mathsf{M},\mathsf{A},\mathsf{S}$ denote the cost of a modular multiplication, modular addition, and shift-with-reduction in ${\mathbb{F}}_p$ respectively.

• Method 1: Direct.

Construct $F^{\left[s\right]}$ as a matrix for each s via Definition 5 (cost $4n^2$ multiplications), then apply to v ($n^2$ multiplications). Total over n values of s: $O\left(n^3\mathsf{M}\right)$.

• Method 2: Iterated NTT.

If an NTT of length n is available at cost $O(nlogn)$, then $F^{\left[s\right]}v$ for a single s can be evaluated in $O(nlogn)$ field operations by leveraging the matrix form of $F^{\left[s\right]}$. Iterating across n values of s: $O(n^{2}logn)$. This baseline applies only when an NTT of length $n=4t$ is available; by Theorem 2(ii), n is not a power of 2 for any DSC prime $p>5$, so radix-2 NTT does not apply directly.

• Method 3: NTT-based projector setup, then streaming.

Compute v, $F_{t}v$, $Jv$, $F_t^{-1}v$ once via NTT (cost $O(nlogn)$ each), then form the projector components $u_{\ell }$ as linear combinations. Stream the orbit via Algorithm 1 or 2. Total: $O(nlogn)$ setup + $O\left(n^2\right)$ streaming = $O\left(n^2\right)$.

• Method 4: Cycle-stride (this paper).

Compute the projector components via direct matrix-vector multiplication (cost $4n^2\mathsf{M}$) or via Method 3’s setup (cost $O(nlogn)$ if NTT applies). Stream via Algorithm 1 (general primitive $e_t$) or Algorithm 2 (DSC, $e_t=2$). Total: $O\left(n^2\right)$.

The principal observations from Table 1:

(a)

The post-projection streaming phase is $O\left(n^2\right)$, which is output-optimal: the orbit contains $n^2$ field elements, so no algorithm producing the full orbit can have lower asymptotic cost.

(b)

Cycle-stride is not an asymptotic improvement over Method 3 (NTT setup + streaming); both achieve $O\left(n^2\right)$ total. The advantages of cycle-stride are: (i) it does not require NTT to be available at length n, which by Theorem 2(ii) is significant for DSC primes; (ii) its streaming phase has a particularly simple structure; and (iii) in the DSC subclass, Algorithm 2 replaces modular multiplications with shifts in the inner loop.

(c)

For single fractional transforms (s fixed), an NTT-based method computes one $F^{\left[s\right]}v$ in $O(nlogn)$ field operations when applicable, beating both cycle-stride variants. Cycle-stride is preferable only when the full orbit (or a substantial fraction of it) is required.

6. Empirical Validation

We implemented both Algorithm 1 (scalar-coefficient streaming) and Algorithm 2 (contribution-array streaming for DSC primes) in Python with NumPy, using 64-bit integer arithmetic throughout. The reference baseline is direct construction of $F^{\left[s\right]}$ as a matrix from the precomputed projectors, followed by matrix-vector multiplication. All implementations share the same projector precomputation; only the per-s work differs.

We emphasize the scope of these experiments at the outset. The measurements below validate correctness of Algorithms 1 and 2 against an independent dense reference, and validate the predicted $O\left(n^2\right)$ scaling of the streaming phase against the $O\left(n^3\right)$ scaling of the dense baseline. They do not benchmark Algorithm 2’s intended shift-reduce hardware advantage — which is a property of FPGA, ASIC, and microcontroller implementations and is not exposed on a general-purpose CPU running Python with NumPy — nor do they compare against optimized NTT pipelines, where a fair head-to-head requires careful implementation in a common framework. The shift-reduce hardware advantage and the NTT comparison are both left as future work.

Table 2 reports per-orbit wall-clock timings on a random signal $v\in {\mathbb{F}}_p^n$, averaged over 30 runs, on a single CPU core. The rightmost three columns show speedups relative to direct evaluation.

• Correctness validation.

Both streaming algorithms produce outputs identical to the direct baseline at all n meridian indices for every prime in Table 2. We additionally verified correctness of Algorithm 1 on a non-DSC symmetry-complete prime ($p=41$, with primitive root $e_t=6$); the general streaming structure applies regardless of whether $e_t=2$.

• Scaling validation.

The ratio between direct and streaming costs grows approximately linearly in n, consistent with the predicted $\Theta \left(n\right)$ ratio between $O\left(n^3\right)$ direct evaluation and $O\left(n^2\right)$ streaming. The empirical speedup of Algorithm 1 over direct at $p=197$ is approximately $97\times$, with $n=196$.

• Algorithm A vs Algorithm B in Python.

A notable observation is that the contribution-array streaming (Algorithm B) runs approximately $30\%$ slower than scalar-coefficient streaming (Algorithm A) in our Python/NumPy implementation. The reason is implementational rather than asymptotic: NumPy’s elementwise scalar multiplication (used in Algorithm B to update the contribution arrays) carries per-call overhead comparable to a small matrix operation, while scalar coefficient updates in Algorithm A use native Python integers with low overhead. The crossover in Algorithm B’s favor occurs on hardware where general modular multiplication is substantially more expensive than shift-and-reduce, which is not the case for 64-bit integer arithmetic on modern x86. The contribution-array variant’s principal advantage is therefore on FPGA, microcontroller, or custom-ASIC targets where multiplier area or cycle cost is significant.

7. Discussion

7.1. Faithfulness and Signal-Orbit Degeneracy

The streaming algorithms (Algorithms 1 and 2) compute $F^{\left[s\right]}v$ correctly for every $s\in {\Phi }_t$ regardless of whether the operator family $s\mapsto F^{\left[s\right]}$ is faithful or whether the realized signal orbit is degenerate. We separate two distinct questions about the orbit and clarify the role each plays in the present work.

Operator-level faithfulness concerns whether the map $s\mapsto F^{\left[s\right]}$ is injective on ${\Phi }_t$ as a map into $\mathrm{End}\left(V_t\right)$. From Definition 5, $F^{\left[s\right]}=F^{\left[r\right]}$ iff ${\sum }_{\ell }(e_t^{\ell s}-e_t^{\ell r}){\Pi }_{\ell }=0$, which (using ${\Pi }_{\ell }{\Pi }_m={\delta }_{\ell m}{\Pi }_{\ell }$) holds iff $(e_t^{\ell s}-e_t^{\ell r}){\Pi }_{\ell }=0$ for every ℓ. Hence operator-level faithfulness reduces to whether the spectral projectors ${\Pi }_1$ and ${\Pi }_3$ are nonzero: if ${\Pi }_1\ne 0$, then $e_t^s=e_t^r$ in ${\mathbb{F}}_p$ forces $s=r$ modulo n by primitivity of $e_t$.

Signal-orbit distinctness concerns whether, for a fixed vector v, the realized outputs $\{F^{\left[s\right]}v:s\in {\Phi }_t\}$ are pairwise distinct as vectors. By the projector-orbit identity (Theorem 1), $F^{\left[s\right]}v=F^{\left[r\right]}v$ iff $(e_t^{\ell s}-e_t^{\ell r})u_{\ell }=0$ for every ℓ, where $u_{\ell }={\Pi }_{\ell }v$ is the realized projected component. A particular vector with $u_1=0$ has a signal orbit of period at most $2t$ in s, regardless of whether ${\Pi }_1$ is nonzero in operator terms.

These are different questions, and the algorithms do not require either: they produce $F^{\left[s\right]}v$ correctly for every s, with repeated outputs when the signal orbit is degenerate. The full operator-level faithfulness result for $t\ge 2$ — which can be obtained by computing the projector ranks $\mathrm{rank}{\Pi }_{\ell }$ via $\mathrm{tr}{\Pi }_{\ell }$ and Gauss-sum evaluations of $\mathrm{tr}{F}_t^r$ — is not required by Theorems 1 and 3.

7.2. Applications

The principal use cases for streaming over the FrFT orbit are settings that require exact arithmetic over a finite field combined with sweeping over the fractional-rotation parameter. Three concrete examples:

Time-frequency signal analysis with exact arithmetic.

The continuum FrFT is widely used in chirp detection, time-frequency filtering, and time-frequency localization [2]. Implementations over a finite field via NTT-style arithmetic avoid floating-point roundoff in applications that demand reproducibility (e.g., regulatory signal processing or symbolic computation). When the application requires sweeping the fractional-rotation parameter — as in chirp parameter estimation or time-frequency contour search — the cycle-stride algorithm computes the entire sweep at output-optimal cost.

Qudit quantum circuit verification.

Floratos and Pavlidis [5] use the finite FrFT family as a one-parameter family of qudit gates. Classical simulation and verification of such circuits requires computing the action of each gate on a state vector; the cycle-stride streaming computes the action of all n fractional Fourier gates on a fixed state vector in $O\left(n^2\right)$ field operations, producing the n-vector orbit of total output size $n^2$. Explicitly materializing the full family of n gate matrices would have output size $\Omega \left(n^3\right)$ and is outside the vector-orbit complexity claimed here; for circuit-simulation applications that need only gate actions on specific input states, the vector-orbit complexity is the relevant figure.

Hardware implementations.

On targets without hardware modular multipliers (low-power microcontrollers, area-constrained FPGAs), Algorithm 2 is structurally attractive: the streaming inner loop requires no multiplication unit beyond what is needed to compute the initial projections. For applications that perform many full-orbit evaluations on the same prime field — e.g., batch processing of different signals through the same fractional Fourier sweep — the projector precomputation amortizes and the streaming cost dominates, making the multiplier-free inner loop the relevant figure of merit.

7.3. Limitations and Future Work

We note several limitations of the present work.

First, the empirical comparison is against a dense matrix-vector baseline rather than against an optimized NTT pipeline. A fair comparison against, for example, hand-tuned Bluestein chirp-z implementations or production NTT libraries (NFLlib, Microsoft SEAL) is necessary to establish the practical regime in which cycle-stride dominates. We expect cycle-stride to remain competitive for full-orbit workloads but to lose to NTT for single-transform use cases when NTT applies.

Second, the cycle-stride algorithm produces the full orbit. Some applications may want a contiguous sub-arc of the orbit (e.g., a small window of fractional parameters around the standard Fourier transform). For such applications, a tailored algorithm restricting the streaming to the desired sub-arc is natural, with proportionally reduced cost; we have not developed this variant in detail.

Third, the DSC condition that 2 is primitive could be relaxed to $2^k$ for small k, broadening the prime class at the cost of slightly more expensive contribution updates (each step requires k shifts rather than ℓ shifts for the ℓ-th contribution array). The full characterization of which prime classes admit a similar shift-only streaming phase is an open question.

Fourth, the structural separation of the DSC primes from the Fermat-prime NTT class (Theorem 2(ii)) means that the two methods address different prime classes and are not in direct competition for the same applications. A hybrid system that uses standard radix-2 NTT for Fermat primes and cycle-stride for DSC primes would cover a broader range of cycle lengths than either method alone, though the cycle lengths involved (n a power of 2 versus $n=4·\left(\mathrm{odd}\mathrm{composite}\right)$) typically correspond to different signal sizes.

8. Conclusion

We have identified a streaming structure for the finite fractional Fourier transform over symmetry-complete prime fields. The full orbit of $p-1$ fractional powers on a signal of length $p-1$ can be generated in $O\left({(p-1)}^2\right)$ post-projection field operations, matching the $\Omega \left({(p-1)}^2\right)$ output-size lower bound up to constant factors. The streaming procedure separates cleanly into a projection setup that depends only on the field and the input signal, and a streaming phase whose cost is independent of the implementation of the setup. For the dyadic-symmetry-complete primes — those with $p\equiv 5(mod8)$ for which 2 is a primitive root — the streaming phase admits a contribution-array reformulation that eliminates general modular multiplications from the inner loop, leaving only modular additions and shift updates.

The class of dyadic-symmetry-complete primes is structurally distinct from the Fermat primes used in classical radix-2 NTT settings: among the Fermat primes only $p=5$ is DSC, and for every DSC prime with $p>5$, the cycle length $n=p-1$ is not a power of two. The two prime classes are therefore complementary rather than competing, and a finite-field signal-processing toolkit that includes both broadens the range of admissible cycle lengths beyond what either method covers individually.

The result is a clean algorithmic note rather than a dramatic asymptotic improvement: cycle-stride achieves the same $O\left(n^2\right)$ full-orbit total cost as the natural NTT-setup-and-stream alternative, and its principal contribution is in the structure of the streaming phase rather than in the asymptotic complexity. We expect the contribution-array variant to be most useful in hardware implementations where modular multiplication carries significant area or energy cost relative to shift-and-reduce; a head-to-head comparison against tuned NTT pipelines is left to future work.