1+ // Cooley–Tukey FFT (radix-2, iterative) and polynomial/big-integer multiplication
2+ // Exports: fft, ifft, multiplyPolynomials, convolveReal, multiplyBigIntegers
3+
4+ function isPowerOfTwo ( n ) {
5+ return n && ( n & ( n - 1 ) ) === 0
6+ }
7+
8+ function nextPowerOfTwo ( n ) {
9+ let p = 1
10+ while ( p < n ) p <<= 1
11+ return p
12+ }
13+
14+ function bitReverse ( n , bits ) {
15+ let rev = 0
16+ for ( let i = 0 ; i < bits ; i ++ ) {
17+ rev = ( rev << 1 ) | ( n & 1 )
18+ n >>= 1
19+ }
20+ return rev
21+ }
22+
23+ function fft ( re , im , invert = false ) {
24+ const n = re . length
25+ if ( ! isPowerOfTwo ( n ) ) {
26+ throw new Error ( 'fft input length must be a power of two' )
27+ }
28+ if ( im . length !== n ) throw new Error ( 're and im lengths must match' )
29+
30+ // Bit-reverse permutation
31+ const bits = Math . floor ( Math . log2 ( n ) )
32+ for ( let i = 0 ; i < n ; i ++ ) {
33+ const j = bitReverse ( i , bits )
34+ if ( i < j ) {
35+ ; [ re [ i ] , re [ j ] ] = [ re [ j ] , re [ i ] ]
36+ ; [ im [ i ] , im [ j ] ] = [ im [ j ] , im [ i ] ]
37+ }
38+ }
39+
40+ // Iterative FFT
41+ for ( let len = 2 ; len <= n ; len <<= 1 ) {
42+ const ang = 2 * Math . PI / len * ( invert ? 1 : - 1 )
43+ const wLenRe = Math . cos ( ang )
44+ const wLenIm = Math . sin ( ang )
45+
46+ for ( let i = 0 ; i < n ; i += len ) {
47+ let wRe = 1
48+ let wIm = 0
49+ for ( let j = 0 ; j < len / 2 ; j ++ ) {
50+ const uRe = re [ i + j ]
51+ const uIm = im [ i + j ]
52+ const vRe = re [ i + j + len / 2 ]
53+ const vIm = im [ i + j + len / 2 ]
54+
55+ // v * w
56+ const tRe = vRe * wRe - vIm * wIm
57+ const tIm = vRe * wIm + vIm * wRe
58+
59+ // butterfly
60+ re [ i + j ] = uRe + tRe
61+ im [ i + j ] = uIm + tIm
62+ re [ i + j + len / 2 ] = uRe - tRe
63+ im [ i + j + len / 2 ] = uIm - tIm
64+
65+ // w *= wLen
66+ const nwRe = wRe * wLenRe - wIm * wLenIm
67+ const nwIm = wRe * wLenIm + wIm * wLenRe
68+ wRe = nwRe
69+ wIm = nwIm
70+ }
71+ }
72+ }
73+
74+ if ( invert ) {
75+ for ( let i = 0 ; i < n ; i ++ ) {
76+ re [ i ] /= n
77+ im [ i ] /= n
78+ }
79+ }
80+ }
81+
82+ function ifft ( re , im ) {
83+ fft ( re , im , true )
84+ }
85+
86+ function convolveReal ( a , b ) {
87+ const need = a . length + b . length - 1
88+ const n = nextPowerOfTwo ( need )
89+
90+ const reA = new Array ( n ) . fill ( 0 )
91+ const imA = new Array ( n ) . fill ( 0 )
92+ const reB = new Array ( n ) . fill ( 0 )
93+ const imB = new Array ( n ) . fill ( 0 )
94+
95+ for ( let i = 0 ; i < a . length ; i ++ ) reA [ i ] = a [ i ]
96+ for ( let i = 0 ; i < b . length ; i ++ ) reB [ i ] = b [ i ]
97+
98+ fft ( reA , imA )
99+ fft ( reB , imB )
100+
101+ const re = new Array ( n )
102+ const im = new Array ( n )
103+ for ( let i = 0 ; i < n ; i ++ ) {
104+ // (reA + i imA) * (reB + i imB)
105+ re [ i ] = reA [ i ] * reB [ i ] - imA [ i ] * imB [ i ]
106+ im [ i ] = reA [ i ] * imB [ i ] + imA [ i ] * reB [ i ]
107+ }
108+
109+ ifft ( re , im )
110+
111+ const res = new Array ( need )
112+ for ( let i = 0 ; i < need ; i ++ ) {
113+ res [ i ] = Math . round ( re [ i ] ) // round to nearest integer to counter FP errors
114+ }
115+ return res
116+ }
117+
118+ function multiplyPolynomials ( a , b ) {
119+ return convolveReal ( a , b )
120+ }
121+
122+ function trimLSD ( arr ) {
123+ // Remove trailing zeros in LSD-first arrays
124+ let i = arr . length - 1
125+ while ( i > 0 && arr [ i ] === 0 ) i --
126+ return arr . slice ( 0 , i + 1 )
127+ }
128+
129+ function multiplyBigIntegers ( A , B , base = 10 ) {
130+ // A, B are LSD-first arrays of digits in given base
131+ if ( ! Array . isArray ( A ) || ! Array . isArray ( B ) ) {
132+ throw new Error ( 'Inputs must be digit arrays' )
133+ }
134+ const conv = convolveReal ( A , B )
135+
136+ // Carry handling
137+ const res = conv . slice ( )
138+ let carry = 0
139+ for ( let i = 0 ; i < res . length ; i ++ ) {
140+ const total = res [ i ] + carry
141+ res [ i ] = total % base
142+ carry = Math . floor ( total / base )
143+ }
144+ while ( carry > 0 ) {
145+ res . push ( carry % base )
146+ carry = Math . floor ( carry / base )
147+ }
148+ const trimmed = trimLSD ( res )
149+ return trimmed . length === 0 ? [ 0 ] : trimmed
150+ }
151+
152+ export { fft , ifft , convolveReal , multiplyPolynomials , multiplyBigIntegers }
0 commit comments