Comput. Math. Applic. Vol. 15, No. 9, pp. 757-768, 1988 0097-4943/88 $3.00+0.

00
Printed in Great Britain. All rights reserved Copyright © 1988 Pergamon Press pie

EXACT DECONVOLUTION USING

NUMBER-THEORETIC TRANSFORMS
G. DRAUSCHKE and M. TASCHE
Wilhelm-Pieck-Universit/it Rostock, Sektion Mathematik, Universit~itsplatz l,
DDR-2500 Rostock, G.D.R.

(Received 5 October 1987)

Communicated by E. Y. Rodin

Abstract--Using residue arithmetic and number-theoretic transforms, new algorithms of the one- and
two-dimensional deconvolution are described. These algorithms work without round-off errors and yield
the exact solution under natural conditions. By the same method, the exact value of the determinant of
a circulant integer matrix is computed, too.

1. D E C O N V O L U T I O N PROBLEMS

Let us denote by Z the ring o f integers, by [~ the set o f positive integers and by Q the field o f
rational numbers. Let n ~ [~, n > 2. For given n-dimensional integer vectors
a = [a(0) . . . . . a ( n - 1)]TeT/" and b = [ b ( 0 ) . . . . . b(n - 1)]Te Z ", b # o.'=[0 . . . . . 0] T, we want to
solve numerically without r o u n d - o f f errors the following special system o f n linear equations:
(circ a) x = b, (1)
where
F a(0) a(n -- 1) ... a(1)]
circa.-= / all) a!0) a<:2)/ E Z ~x~

La(n--1) a(n -- 2)
... a(O)J
is a nonsingular circulant matrix and x = [x(0) . . . . . x ( n - 1)]T e Q" is an u n k n o w n vector. We
r e m a r k that the one-dimensional deconvolution p r o b l e m (1) can be written as
a * x = b, (2)

where • denotes the one-dimensional cyclic convolution, i.e.

n--I
b(k)'.= ~ a ( k - i ) x(i), k=O ..... n-1.
i=0

N o t e that the difference k - i must be calculated m o d u l o n.

T h e two-dimensional deconvolution p r o b l e m considered in Section 6 reads as follows. Let two
n x n integer matrices A = [a ( i ,J)]ia=0
" " - i ~ 7/" ×" and B = [b(t,j)],.j=
' ' " - Io ~ 7/" ×", B ~ O , = [0]7.fl0
. - be given.
Find the exact solution X - [x(t,J)]~j=oe
' " " - 1 Q "×" of the relation
A * X = B, (3)
where • denotes the two-dimensional cyclic convolution, i.e.
n-In-1
b(k,l),= ~ ~ a(k-i,l-j) x(i,j), k,l=O ..... n-1.
~=0j=0
Here, both k - i and 1 - j must be calculated m o d u l o n. Of course, we assume that (3) is uniquely
solvable.
C.A.M.W.A.,5/9--D 757
758 G. DRAUSCHKE and M. TASCHE

Sampling an integer matrix X • 7/" ×" by rows, we obtain the n2-dimensional vector
row X.'=[x(0, 0) . . . . . x(0, n - 1) . . . . . x ( n - 1,0) . . . . . x ( n - l , n - l)]X • 7/"2.
N o t e that x = [ x ( 0 ) , . . . , x ( n 2 - 1)]T = row X if and only if x ( n i + j ) = x ( i , j ) , i , j = 0 . . . . . n - 1.
F o r A • 7/n ×", let

circ A.'= I AA1

I A0
A._l •
...
•.• A2
At 1 77": × n2

-I A._2 -.- A0

denote a circulant block matrix consisting of circulant blocks A i : : circ[a(i, 0) . . . . .

a(i, n - 1)]Te ~7" ×". Then (3) is equivalent to the special system o f n 2 linear equations
(circ A ) r o w X = row B. (4)
We see that (3) is uniquely solvable if and only if circ A is a nonsingular matrix•
Using residue arithmetic and number-theoretic transforms, we treat the above deconvolution
problems. F o r the m o r e general problem, the c o m p u t a t i o n of the exact solution o f an integer system
o f linear equations using residue arithmetic, see [1] and [2].

2. RESIDUE ARITHMETIC

Let m • N, m > 1 be an odd integer• Given any x • Z, if x = r ( m o d m) and if Iris< (m - 1)/2,

then we write r = / x / m and we say t h a t / x / m is the s y m m e t r i c r e s i d u e o f x ~ 2~ m o d u l o m [3, p. 1 13]•
F o r arbitrary x, y • 7/, we have
I x +_ y l m = I l x l ~ + _ / Y / m / m ,

I x "Ylm = I l x l , ~ ' l Y l m l m .

Further, if gcd(x, m ) = 1, then there exists a unique integer x ' w i t h / x x ' / m = 1, Ix'l-< (m - 1)/2.
This integer x ' is called m u l t i p l i c a t i v e i n v e r s e o f x ~ 7_ m o d u l o m and x ' is denoted by
/1/
The multiplicative inverse m o d u l o m can be c o m p u t e d via the extended Euclidean algorithm [1].
The symmetric residue o f an integer vector x = [x(0) . . . . . x ( n - 1)]T e Z" m o d u l o m is defined
by
/x/re'.= [Ix(O)~ . . . . . . Ix(n - 1)/m]r • 7/".
Analogously, we explain the symmetric residue o f an integer matrix m o d u l o m. F o r arbitrary
x = Ix(0) . . . . . x ( n - 1)]T• Z" and y = [y(0) . . . . . y ( n - 1)IX• Z ", we have

/x + y/m =/lx/~ +/Y/m/m,

/X * y /~ = / / X / . */y/m/~,
Ix o y/~ =//x/m o/y/m/m,
where o signifies the H a d a m a r d p r o d u c t
x o y..= Ix(0) y(0) ..... x(n -- 1)y(n - 1)]T.
Further, if gcd(x(j), m ) = 1 for a l l j = 0 . . . . . n - 1, then there exists a unique vector x ' • 7/n with
/x ° x'/m = [1 . . . . . 1]r, x' =/X'/m. We have
1 T
. . . . .

This vector x' is called H a d a m a r d inverse of x • Z" modulo m.

 Exact deeonvolution using number-theoretic transforms 759

• ' n-I ~.7/nxn, i " n-I

If X - [x(t,J)],j=0 Y = [ Y ( , J ) ] , 4 = 0 • Z ~ ~" and x • Z ~, then it holds

ix x/. =//x/./x/~/.,
/x _+vim = //X/m "~-/Y/m/m,
/x.v/~ m //X/ ./V/m/m ,

/x • v/~ = //X/m */Y/m/m,

/x v / . =//x/~ o/y/~/~,
o

/X (~ Y/m=//X/m ®/Y/m/m,
where o signifies the Hadamard product
X o Y,= [x(i,j) y(i,j)]~i~o,

and where ® denotes the K r o n e c k e r product. Further, if g c d ( x ( i , j ) , m) = 1 for all i , j = 0 . . . . .

n - 1, then there exists a unique matrix X ' e Z "×" with /X o X'/m = [1]i,:= " - 0, X ' = / X ' / , , . We have

Me ~
i,j=0

This matrix X' is called H a d a m a r d inverse o f X • 7/"×" modulo m.

N o w we consider a finite n u m b e r o f pairwise relatively prime, odd moduli mi • N, mi > 1,
i = 1. . . . . s. Very often, these moduli are different odd primes. Let m = m ~ . . . m s. Using the
Chinese Remainder T h e o r e m (see [4]), we can reconstruct uniquely an integer x with
Ix l~< (m - 1)/2 f r o m its symmetric residues x : = / X / m : i = 1 . . . . , S.

T h e o r e m 2.1
U n d e r above assumptions, any x • Z with Ixl ~ (m - 1)/2 and with given symmetric residues x~,
i = 1. . . . . s can be expressed uniquely in the f o r m

X = [Xl] + [ X I X2] ml + [xl x2x3] m l m 2 + " " + [ x l . . . x~] m l . . . ms_ i (5)

•- .
with moaeular divided differences [xi],= xi, i = 1. . . . . s,

[XiXj]:=/[XJ]~ii[Xi!/" mj. l. <~i. < j . <~s,[x~, xj,_lxj,=/[x~l 1 *

" " "xi'-2XJm~,--[Xil °Xir-I]/mir' (6)

where 1-%<ij< . . . < i , ~ < s .

Proof. It is k n o w n (see [4]) that every x • Z with Ix l~< (m - 1)/2 can be represented uniquely
in the form
X = X 0) + x ( 2 ) m l + x O ) m t m 2 + • • • + x ( ~ ) m ~ . . , m~_ ~ (7)

with x t/) • Z, Ix°~l ~ ( m ~ - 1)/2, i = 1. . . . . s. By the Chinese Remainder Theorem, the residue-class
ring 7//m7/ and the direct p r o d u c t Z / m ~ Z x . . . x Z / m s Z are isomorphic. Hence, x • 7 / with
Ix l~< (m - 1)/2 is uniquely determined by its symmetric residues xi, i = 1. . . . . s.
N o w let x • 7/with Ix l~< (m - 1)/2 be given by its symmetric residues x;, i = 1. . . . . s. We have
to show that x ~° = [x~ x2. • .xi], i = 1 , . . . , s. We employ induction on s. The case s = 1 is trivial.
Assume that the assertion is true for s - 1. Setting

y "= [xl] + [xl x2] m~ + . . " + [xl ... Xs_2]ml... m~-3 + [xl . . . xs_ I]ml . . . m s - 2 ,
Z'=[Xl]+[XlX2]mI +'''+[Xl...X~_z]ml...m,_3+[xl...x,_2x,]ml...m,_2, (8)
we have by the induction hypothesis that

/ y /,,~ = xi, i = l ..... s-l,

/Z/m, = Xi, i = 1 . . . . , S -- 2, S. (9)
760 G. DRAUSCHKEand M. TASCHE

Then we conclude from (8) and (9) that

w:=y + (z - y ) m , _ l
m s

= Y Jl-/_ l / ,,,~(Ix[ " . X s. 2 X s ]. - - [ X. l . . . X s .l ] ) m l .ms_

satisfies /W/m ~ = X~, i = 1. . . . . S -- 1. Hence, w -- x ( m o d m) by the Chinese Remainder Theorem.

F r o m Ix I ~< (m - 1)/2 and (7), we conclude that

x = y + /[XI'''Xs-2Xs]--[XI'''Xs-I]/---- m , . .. m~_,.
m s - i m s

By (6)-(8), it follows that x ~° = [ x t . . . x~], i = 1. . . . . s. This completes the proof. [ ]

It is convenient to exhibit the m o d u l a r divided differences in the form o f a triangular array called
a m o d u l a r divided difference table. This will look for s = 4 as follows:
ml m2 m3 m4
[xt] [x2] [x3] [x41
[x, xd [x, xA [xlx,]
[Xl X 2 2 3 ] [Xl X 2 X 4 ]

[Xl X2X3X4]
In this way, we can reconstruct also integer vectors or integer matrices f r o m the corresponding
symmetric residues.

Corollary 2.2
Let m i ~ N, m i > 1, i = 1. . . . . s, be pairwise relatively prime, odd moduli. Further let
m = m l . . . m s-
(i) Then any x = [x(0) . . . . . x ( n - 1)]xe Z" with Ix(i)l ~< ( m - 1)/2, i = 0 . . . . . n - 1 , and with
given x~:=/X/m,, i = 1. . . . . S, can be represented uniquely in the form

x = [x,] + [x~ x2]m~ + " " + [Xl.. • X s ] m l . . . m,_ i

with [x/].-= xi, i = 1. . . . . s,

[x~xj],= / [ x j ] - - [ x i ] / , l ~ < i < j ~ < S,

/ mi /,,j
x, x j , = / [ x " " " ' x~'-2 x j - - [ x ~ ' " " " x i r - ' ] /
Xil • . .
r I mi, -- I mi r ,

where l ~ < i ~ < ' ' ' < i r < ~ S .

(ii) A n y X = [x (,J)]t.j=0~
i " n-I 7/n×, with I x ( i , j ) l <<,(m- 1)/2, i , j = 0 , . . . , n - 1, and with given
X;.'=/X/m,, i = 1. . . . . s, can be expressed uniquely in the form

x = pi,] + [x,x2lm, + . . . + I X , . . . X j m , . . . m._,

with [X~].-= X~, i = 1. . . . . s,
[X~Xj],=/[XJ]~[Xi]// / , 1 ~i <j ~s,
/ m~ /mj
.X~, Xi,].'=/[xq'''xt'-2x~r]-[x~''''x~r ']/
Xil • .
- I mi, I mir ~

where 1 ~< i~ < • • • < ir <~S.

3. N U M B E R - T H E O R E T I C TRANSFORMS

Let n e N, n > 2, and let m ~ N, m > 1, be an o d d m o d u l u s in the residue arithmetic. T h e n e ~ 7/,

2~< [el ~ ( m - 1)/2, is called a primitive nth root o f unity modulo m, if e " = 1 ( m o d m) and
 Exact deconvolution using number-theoretic transforms 761

gcd(e k - 1, m ) = 1, k = 1. . . . . n - 1. T h e n one can prove the following result:

Theorem 3.1 [5]

Let n e N, n > 2, and let m E N, m > 1, be an odd integer. A n u m b e r e e 7/, 2 ~< l e I < (m - 1)/2,
is a primitive n t h r o o t o f unity m o d u l o m if and only if ~ . ( e ) =- 0 (rood m ) and gcA (n, m ) = 1,
where 4,. is the n t h cyclotomic polynomial.
The concept o f a primitive n th root o f unity m o d u l o m is f u n d a m e n t a l in the following context.
N o t e that the equality o f two integer vectors x, y e 7/" is defined in the residue arithmetic by
/X/m =/Y/re. Analogously, the equality o f integer matrices is explained in the residue arithmetic.
The one-dimensional number-theoretic transform ( N T T ) ~- o f length n with e as primitive n t h root
o f unity m o d u l o m, and its inverse ~ ' - * , are defined for x e 7/" by
.~x = i = [:~(0) . . . . . :~(n - 1)] T,

O~--Ix = y = [y(0) . . . . . y(n -- 1)]T

with

Yc(k),=
iffiO
x(i) e a'/ m~
k = O. . . . . n - 1,

y(i),= n
= x ( k ) e-U' / m, i = O . . . . . n --1.

N o t e that there exists such an integer

by gcd (n, m ) = 1 (see T h e o r e m 3.1). F o r arbitrary x e Z " , we have #'-~3r-x = ~ - # - - I x =/X/m.

Introducing the integer matrices F = [/e • /,.]i, , -kl= 0 and

F_ 1 = //n/ra
1 [/e --ik/m]~,k=0
n-lira ,
we get ~ - x = / F x / m and ~ - t x =/F-*X/m for all x ~ Z". In the case n = 2 '+ *, t ~ N, and e = 2, we
obtain the so-called Fermat number transform. By T h e o r e m 3.1, we can choose as possible moduli
m ~ • only divisors > 1 o f ~2, +*(2) = 22` + 1.
The N T T was introduced in [6] and [7] as a natural generalization o f the discrete Fourier
t r a n s f o r m ( D F T ) in order to p e r f o r m fast cyclic convolutions without round-off errors. The N T T
possesses properties resembling those o f the D F T [8, pp. 211-216], particularly the cyclic
convolution property
~'(x • y) = / ~ o ~/,. (10)
for all x, y ~ 7/".

Lemma 3.2
Let x e 7/". Then we have
/F(circ x) F-t/m = diag i , 01)
Met (circ x)/m = 2 (k)/m. 02)
Ik=O
Proof. Let x, y E Z". By (10) and x , y = ( c i r c x ) y , it follows that
3r(circ x) y = 3r(x * y) = / i o Y/re.
Setting z , = ~ , i.e. 3r-~z = / y / = , this yields
oar(circ x) ~ r - tz = / i o z/~ = / ( d i a g i ) z/m
762 G. DRAUSCHKEand M. TASCHE

for all z • Z". Then we obtain (11). Hence we have for the determinant of a circulant matrix
/ - I

/det(circ x)/m = / d e t ( d i a g i)/,, =/"IZI ' ~(k)/m.

/k=O
This completes the proof. []
From numerical point of view, the following three essential conditions on NTTs are required:
(1) the length n has to be large enough and highly factorizable in order to implement fast algorithms
(see [8, pp. 85-94, 116-120, and 125-144]); (2) the primitive nth root e of unity modulo m should
have a simple binary representation such that the arithmetic modulo m is easy to perform; (3) the
modulus m has to be large enough to avoid overflow, but on the other hand small enough such
that the machine word length is not exceeded; furthermore, m should have a simple binary
representation.
For instance, the pseudo-Fermat number transform with n = T + l, t e ~, e = 2 q, q E ~, and m e ~,
m > 1, m l2qV+ 1, is a compromise between these various conditions. By Theorem 3.1, we can
determine all possible moduli m e [~, m > 1 for given length n e [~, n > 2, and given e ~ 7/, l e[ >~ 2,
such that e is a primitive n th root of unity modulo m (see [5]).
The two-dimensional number-theoretic transform ~-2 of size n x n with e as primitive nth root
of unity modulo m, and its inverse , ~ 1 are defined for all X = t rx~t i,J)li./=0
:,1,-~ E Z" ×" by
#'2X = ,'K,=/FXF/m,
~-IX = y..=/F-1XF-I/~.
More precisely, we have X = [~(k, l)]k,J=0e
"-t 77"×" and Y = [y(,J)]~.j=0~
i " "-~ 7/"×" with
/,-1,-I /
.~(k, l).'= // X-'z., ~X-"xti, ,J):"eik+itl/m,
/ i=O j = O /
/ / 1 / n--ln-I /
y(i,j),=//-~/m ~ 0 t =~0 x(k, l)e-'ik+/t),,,./
For arbitrary X e 77"×", we get ~-{l .~2X = . ~ 2 . ~ -l X =/X/m. Introducing the Kronecker product
F ® F, we have the known relation between the two-dimensional N T T of size n × n and a
one-dimensional transform of length n 2
row R = / ( F ® F) row X/m (1 3)
for all X • Z" ×". The two-dimensional N T T has properties resembling those of the two-dimensional
D F T , particularly the cyclic convolution property
/F(X • Y) rim = / ( F X F) o (F Y e)/m = / f ( o '~/m (14)
for all X, Y e 77"×".

Lemma 3.3
Let X • Z" ×". Then we have
/(F ® F) (circ X) (F- ' ® F - l)/,, = diag (row X), (15)
/n-I n 1
/det(circ X)/m =/,~I=o too :c(k, l)/m . (16)

Proof Let X, Y E 7/" ×". By (4) we see that row(X • Y) = (circ X) row Y. Using (13) and (14), we
have then
/(F ® F)(tire X) row Y/m = / ( F ® F) row (X * Y)/m = / r o w F(X * Y)F/m = / ( r o w X) o (row Y)/m.
Setting Z,=¢¢, this yields row Z = / ( F ® F ) row Y/,, by (13). Hence, we obtain by
/row Y/,, = / ( F -1 ® F -l) row Z/m the relation
/(F ® F)(circ X)(F -I ® F -I) row Z/m = / ( r o w X) o (row Z)/,. =/(diag(row X)) row Z/m
for all Z • Z" ×". Thus it follows (15). Calculating the determinants on both sides of (15), we get
the equation (16). This completes the proof. []
 Exact deconvolution using number-theoretic transforms 763

4. A L G O R I T H M FOR ONE-DIMENSIONAL DECONVOLUTION

Let n e N, n > 2. F o r given integer vectors a, b e Z n, we consider the one-dimensional de-

convolution problem (1) and (2), respectively. Assume that A , = c i r c a e Z n × n is nonsingular. Let
d , = det A # 0 denote the determinant o f A. By the well-known H a d a m a r d ' s inequality, it holds
]tn- 1 ~n/2
[di <...cii=tk~=oa(k)2 ) (17)

The formal solution x e Qn o f (1) reads

x.'=A-Ib = 1 y (18)

with y = [y(0), . . . , y ( n - 1)IT.'= A~dJb. F o r A ~¢.-.- tfaadJti~,J/~,S=o;~ln-I~ 7in ×n, it follows again by Hada-
mard's inequality that
(~£, )(~-,)/2
max la~dJ(i,j)[ <~ a(k) 2
i,j \k =0
Consequently, we obtain that
n- I '~ i n - I \(n- I)/2
m aI x l,y ( i ) ~ c2,= C ~°'= . ~ ( j ) , ) t ~ o a ( k ) ') (19)

Applying T h e o r e m 3.1, we choose a finite n u m b e r o f N T T s ~,., i = 1. . . . . s, with e, e L l e , I >t 2,

as primitive n th r o o t o f unity m o d u l o m~ e N, mi > 1. Denote the corresponding transformed vector
o f arbitrary w e Z n by ~ = [~,~(0). . . . . ~,i(n - 1)]r , = ~ w, i = 1. . . . . s. We assume that the moduli
m~, i = 1. . . . . s, are pairwise relatively prime and that m = m o . . . ms fulfils the condition
max{c~, c2} ~< (m - 1)/2. (20)
Further we suppose that
gcd(d, m) = 1. (21)
By (12) we obtain that
/~_=J
4'=idA.,-- I l l ,=l .... (22)
/k=O
Then one can see by (22) that (21) is equivalent to
gcd(~i(j), mi) = I, j = 0 ..... n - I; i = I..... s. (23)
F r o m (17) and (20) it follows that Idl ~< (m - I)/2. Using Theorem 2.1, we get the exact value
d =[dl] + [dldz]m~ + ... + [dl... ds]m~... ms_~ (24)
of the determinant of the circulant matrix A.
Applying the N T T s :;, i = I,..., s, to (2), we obtain by the cyclic convolution property (10)
that
/'~i ° i,/,., = 6. i = 1. . . . . s. (25)
By (23) there exist the H a d a m a r d inverses (ii)' o f i~ m o d u l o m~, i = 1. . . . . s, which can be
computed via the extended Euclidean algorithm. Then from (25) it follows that
~; =/($~)' o ~/,~, i = 1. . . . . s. (26)
Using the inverse N T T s ,~'7 I, i = 1. . . . . s, we obtain
/x/~, = / F T ' i , / m , , i = 1. . . . . s. (27)
Then by (18) we know the symmetric residues y:=/y/m~ o f y = d x ~ ~,n m o d u l o m~, i = 1. . . . . s. By
(19) and (20) it follows from Corollary 2.2 that

Y = [Yl] + [YlYz]ml + " " + [Yl-.-Ys] m l . . . m , _ l . (28)

764 G. DRAUSCHKEand M. TASCHE

Finally, by (21) the unique solution of (1) and (2), respectively reads
1
x=2y~Q".

We s u m m a r i z e this deconvolution algorithm to c o m p u t e the exact solution x e Q" of (1) and (2),
respectively. U n d e r a b o v e mentioned assumptions, the algorithm o f one-dimensional deconvolution
can be expressed as follows:
1. Calculate ai, i = 1 , . . . , s .
2. Calculate di=/d/m,, i = 1. . . . . s by (22).
3. Determine d by (24). Prove that (21) is fulfilled.
4. Determine the H a d a m a r d inverses (~)' o f hi m o d u l o mi, i = 1. . . . . s, via the extended
Euclidean algorithm.
. C o m p u t e G~, i = 1. . . . . s.
6. C o m p u t e xi, i = 1. . . . . s by (26).
7. F o r m / x / , , , , i = 1. . . . . s by (27).
8. Calculate Yi = / d x / m , , i = l , . . . , s .
9. F o r m y by (28).
10. Set
1
x=~y.

N o t e that all operations o f steps 1, 2 and 4-8 can be implemented in parallel manner.

5. D E C O N V O L U T I O N ALGORITHM WITH KNOWN DETERMINANT

Recently, in [9] a precise deconvolution algorithm using the F e r m a t n u m b e r t r a n s f o r m was

described. N o w we generalize this m e t h o d . We use the same notations as before.
Let n e N, n > 2, not necessarily a 2-radix. Let us consider the convolution equation (2) under
the a s s u m p t i o n that d = det A :~ 0 is known. Further, let ~- be a N T T of length n with e e 77,
2 ~< lel <~ (m - 1)/2, as primitive n t h root o f unity m o d u l o an odd integer m > 1. Suppose that
gcd(d, m ) = 1.
F r o m (2) and (18) we conclude that a • y = db. Using the symmetric m - r a d i x representation of
y ~ 77",
y = y~O)+ y~t)m + . . . + y~t)mt

with some t ~ I~ and y~J~=/YO~/m ~ Z", y~0 =/=o, j = 0 . . . . . t, we obtain

(a * y~0)) + (a * y~l)) m + .. • + (a * y~O) m t = d b . (29)

Introducing the remainder vectors r <°).'=db ~ 77", r ~°) =/=o,

r (j+ I).-= 1 (r ('~ -- (a * y(J~)) ~ 7/", j = 0 . . . . . t, (30)

m
we observe that r~t÷ t) = o. This indicates the end o f the calculation.
By (29) and (30), the u n k n o w n vectors y(J~, j = 0 . . . . . t, can be calculated on line from
/a * y~J~/,. =/r~J~/,., j = 0 ..... t. (31)

Then we solve (31) step by step by the k n o w n procedure

Y~ =/Y(J)/m = / F - l ( ( ~ ) , o r(J))/m, j = 0..... t. (32)
Summarizing this m e t h o d , we obtain following deconvolution algorithm:
1. Calculate i.
2. Determine the H a d a m a r d inverse (fi)' o f i m o d u l o m via the extended Euclidean algorithm.
3. Set j = 0.
 Exact deconvolution using number-theoretictransforms 765

4. Calculate r o~ by (30). If r o~ = o, then set t = j - 1 and go to the step 9. If r o~ # o, then go

to the next step.
5. Calculate r ~ , = / F r°~/m.
6. C o m p u t e / ( i ) ' o l ' ~ .
7. Form yO~ by (32).
8. Increment j and repeat from step 4 on.
9. Set
1
x = at (y(0) + y(l) m + . . . + y(,)mr).

An essential assumption of this deconvolution algorithm is that the value of d is known. We have
seen that d can be calculated via the first three steps of our algorithm described in Section 4. By
the following Lemma, we can easily determine the sign of d.

Lemma 5.1 [10, p. 76]

Let n E N, n > 2 and let a = [a(0) . . . . . a(n - 1)]+ • Z". Then we have:
(i) If n is odd, then d/> 0 if and only if
.--1
E a(i) >10.
i=O

(ii) If n = 2r, r > l, then d t> 0 if and only if

,=~a(2i-1) >~ ; ~ i a ( 2 i ) .

If we want to calculate only the value of d, where the sign of d is known, via the first three steps
of the algorithm in Section 4, then it is sufficient to choose a finite number of NTTs with pairwise
relatively prime moduli rn~• I~, m~ > 1, i = 1. . . . . s, such that m~... m, I> el + 1 is satisfied instead
of (20).

6. A L G O R I T H M FOR T W O - D I M E N S I O N A L D E C O N V O L U T I O N

Now we extend the first algorithm for the one-dimensional deconvolution problem to the
two-dimensional case.
Let n • N, n > 2. For given integer matrices A, B E Z "×", we consider the two-dimensional
deconvolution problem (3) and (4), respectively. Assume that C..= circ A • 77"2×,2 is nonsingular, i.e.
d.'= det C ~ 0. By Hadamard's inequality, it holds

Idl <~c3'=(~] "~'

t=0 a(k' l)2) (33)

The formal solution x = row X • 7/'~ of (4) reads as follows:

1
x = C-~b = ~ y (34)

with b r o w B • Z "2 and y Iv(0), .,y(n 2

= = " " I
1)]T:=cadJb. For cadj:=rcadJQ;~ln2-1fZ
t I ,J]Ji, j=O "2×n2, it
follows by Hadamard's inequality that
"~(.2 -- I)/2
m ax'c'dj(i,j)l <~ (.~, .~1 a(k. 1)2) .
t,j \k=0 1=0

Consequently, we obtain that

{n - I n-- I "~/n-- I n- I ~(n 2 - I)/2

maxyl()i<
l ~c4=' (k~=' ~ Ib(k'l)l)~k~= o ~=oa(k,l) =)
0,=0 (35)

Now we choose a finite number of NTTs a~,, i = 1 . . . . ,s, of the length n with e ; • Z ,
2 ~< le, I ~< (m~- 1)/2, as primitive nth root of unity modulo an odd integer m~ > 1. Denote the
corresponding transformed matrix of W • Z"*" by "~/~= [~,~(k, I)]~5~=o,=/F~WF~/=,, i = 1. . . . . s.
766 G. DRAUSCHKE and M. TASCHE

Assume that mr, i = 1. . . . . s, are pairwise relatively prime and that m = m j . . . mr fulfils the
condition
max{c3, c4} ~< (m - 1)/2. (36)
Further, we suppose that

gcd(d, m ) = 1. (37)
By (16) we obtain that

drl=ldlm,=lHH,~,(k,l)/, i = 1 . . . . . s. (38)
l k = O /=0 ] mi

One can see by (38) that (37) is equivalent to

gcd(~ir(k, l), mr) = 1, k, l = 0 . . . . . n - 1; i = 1. . . . . s. (39)
F r o m (33) and (36) it follows that Idl <~ ( m - 1 ) / 2 . Using (24), we get the exact value o f the
d e t e r m i n a n t o f C = circ A.
A p p l y i n g two-dimensional N T T s to (3), we obtain by the cyclic convolution p r o p e r t y (14) that
^ ^

/Fr(A * X) Fi/ra ' = / A i o Xi/m ' =/~r/mi , i = 1..... S. (40)

By (39) there exist the H a d a m a r d inverses (~,~)' o f -g-r m o d u l o mi, i = 1. . . . . s. Then from (40) we
conclude that
Xr = / ( A , ) " ° B,/m,, i = l ..... s. (41)
Using the two-dimensional inverse N T T s , we obtain

/X/m i =/FT' $,FT'/m,, i = 1..... s. (42)

Then by (34) we k n o w the symmetric residues Yr"=/Y/m, o f Y = d X 6 7/n ×" m o d u l o mr, i = 1. . . . . s.
N o t e that y = row Y ~ 7/"2. By (35) and (36) it follows f r o m the Corollary 2.2 that

Y = [YI] + [YI Y2] m l + " " + [YI • • • YJ m l - . . m,_l. (43)

Finally, by (37) the unique solution o f (3) reads

1 Q.
X=~YE ×".

W e s u m m a r i z e the a l g o r i t h m o f t w o - d i m e n s i o n a l deconvolution. U n d e r a b o v e mentioned assump-

tions, this algorithm can be expressed as follows:
1. Calculate ,i, r, i = 1. . . . , s.
2. Calculate d,. = ~d/m,, i = 1 . . . . . s by (38).
3. D e t e r m i n e d by (24). Prove that (37) is fulfilled.
4. Determine the H a d a m a r d inverses (J-r)' o f ~,t m o d u l o mi, i = 1. . . . . s, via the extended
Euclidean algorithm.
. C o m p u t e ]~i, i = 1. . . . . s.
6. C o m p u t e $~, i = 1. . . . . s by (41).
7. F o r m / X / m , , i = 1 , . . . , s by (42).
8. Calculate Y~ =/dX/m,, i = 1 . . . . . s.
9. F o r m Y by (43).
10. Set
1
X =~Y.

Finally, we note that the second deconvolution algorithm described in Section 5 can be extended
in obvious m a n n e r to the two-dimensional case, too.
 Exact deconvolution using number-theoretic transforms 767

7. E X A M P L E

We illustrate the a b o v e algorithm by a simple example. Let n = 3,

100 00] [0 i]
-1 0 --1 -1 --1
A= 1 , B= 0 •
-1 0 -
We want to solve (3) without r o u n d - o f f errors. By (33) and (35) we estimate that Idl ~< 512 and

m a x l y ( i ) l ~< 2048.
i

Thus we have to choose m >/4097 by (36). Let ~'~ be the N T T o f the length 3 with el = 2 as
primitive third root o f unity m o d u l o m~ = 7. Then the corresponding integer matrices o f this N T T
and its inverse read as follows:

Fl= 2 -- , F? 1=/-2 -3 /7.

-3 2 -
Let ~a~"2 be the N T T o f the length 3 with e 2 = 3 as primitive third root o f unity m o d u l o m 2 = | 3.
Thus we have

F2 -- 3 - , F21 = / - 4 -4 /13"
--4 3 --
Further let ~3 be the N T T o f the length 3 with e3 = 8 as primitive third root o f unity m o d u l o
m3 = 73. Hence it holds

F3 = 8 -- , F31 = / _ 2 4 --9 /73"

-9 8 -

N o t e that m~, m2, m 3 are pairwise relatively prime and that m = ml m2m3 ---- 6643 > 4097. Then we
obtain the exact solution X ~ ~ a × 3 of (3) by following steps:

. ~,l = 2 - , A2 = 3 , A,3 = 8 - 3 •
- 2 - 2 - 21

2. dl = / d / 7 = - 2 , dE =/d/13 = 6, d3 =/d/73 = 19.

3. [d,]= - 2 , [dtd~=[d~d3]=3, [dl d2 d3] = 0, d = 19.

. (.~,)'= -3
2
-3
-3
- , (-~2)' = [il !1
-- -4
6
(-i-3)' = -9
8
1
-9
7
2 •

. 1~1=

[°°!1
-3

-3
-2

3
- , B2 = [i s :il
--
--
--4
1
B3 --
L°
-3
-3
--
13
14
34
-41
-34
20
•

[ 0 13 -~).
. 51 = -1

-2
, 52 = -
:i]
-
-4
1
53 =
1-24
27
--2019 -- 14
768 G. DRAUSCHKEand M. TASCHE

. /X/7= ]o 1 !1
3
0
0
0
, /XI,3=
E i4 i] -2
5
, /X/73=
F 5
28
-30
2 - 32
26 .

.
YI = Ei2ij 0
0
-- , Y2 = L,2!I
--4
--5
1
4
-- , Y3 = [,423j
22
21
14
--35
-- 17 •
--24

. [Y1] = Yl,

1F
[Y, Y2] = [YI Y3] =

: 31
E! °:!1-5
2 , [Yl Y2Y3] = O, Y=Y3.

10. x=~122 14 -17

1_21 --35 -24 J
8. C O N C L U S I O N

T h e c o n s i d e r e d d e c o n v o l u t i o n p r o b l e m s are typical inverse p r o b l e m s . It is k n o w n t h a t using the

D F T , the s o l u t i o n o f a d e c o n v o l u t i o n p r o b l e m is very sensitive to r o u n d - o f f errors.
T h e d e s c r i b e d d e c o n v o l u t i o n a l g o r i t h m s using residue a r i t h m e t i c a n d N T T s eliminate the
influence o f r o u n d - o f f e r r o r s a n d yield the exact s o l u t i o n o f the one- a n d t w o - d i m e n s i o n a l
d e c o n v o l u t i o n p r o b l e m , respectively. These a l g o r i t h m s a r e new. T h e first a l g o r i t h m w o r k s m a i n l y
parallel a n d affords p a r a l l e l p r o c e s s o r architectures for real time a p p l i c a t i o n s . F u r t h e r , the exact
value o f the d e t e r m i n a n t o f a c i r c u l a n t integer m a t r i x can be c o m p u t e d by the first steps o f this
algorithm.
By g e n e r a l i z a t i o n o f a recent result [9], we o b t a i n the s e c o n d d e c o n v o l u t i o n a l g o r i t h m , which
w o r k s on line. U n l i k e [9], we s h o w t h a t it is n o t necessary to c o n s i d e r o n l y 2-radix lengths a n d
F e r m a t n u m b e r t r a n s f o r m s . F u r t h e r m o r e , we i m p r o v e k n o w n results c o n c e r n i n g the Chinese
R e m a i n d e r T h e o r e m b y i n t r o d u c t i o n o f new m o d u l a r d i v i d e d differences.

REFERENCES

I. J. A. Howell, Exact solution of linear equations using residue arithmetic. Algorithm 406. Commun. A C M 14, 180-184
(1971).
2. S. Cabay and T. P. L. Lam, Algorithm 522 ESOLVE, Congruence techniques for the exact solution of integer systems
of linear equations. A C M Trans. math. Software 3, 404-410 (1977).
3. N. S. Szab6 and R. I. Tanaka, Residue Arithmetic and Its Applications to Computer Technology. McGraw-Hill, New
York 0967).
4. J. A. Howell and R. T. Gregory, An algorithm for solving linear algebraic equations using residue arithmetic. I, II.
BIT 9, 200-224, 324-337 (1969).
5. R. Creutzburg and M. Tasche, Number-theoretic transforms of prescribed length. Math. Comput. 47, 693-701 0986).
6. P. J. Nicholson, Algebraic theory of finite Fourier transforms. J. Comput. System Sci. 5, 524-547 (1971).
7. J. M. Pollard, The fast Fourier transform in a finite field. Math. Comput. 25, 365-374 (1971).
8. H. J. Nussbaumer, Fast Fourier Transform and Convolution Algorithms. Springer, Berlin (1981).
9. M. Mother6, Precise deconvolution using the Fermat number transform. Comput. Math. Applic. 12A, 319-329 (1986).
10. P. J. Davis, Cireulant Matrices. Wiley, New York (1979).