Building Blocks for Symmetric Cryptography Statistics and Skipjack TU-Decomposition and Kuznyechik The Butterfly Permutations and Functions Conclusion

⊕ x x ⊕ a a E κ E κ E κ (x) E κ (x ⊕ a) b ⊕ 4 / 46
⊕ x x ⊕ a a E κ E κ E κ (x) E κ (x ⊕ a) b ⊕

Differential Attack

If there are many x such that E κ (x) ⊕ E κ (x ⊕ a) = b, then the cipher is not secure. Motivation (1/3)

A malicious designer can easily hide a structure in an S-Box.

To keep an advantage in implementation (WB crypto)...

Motivation (1/3)

A malicious designer can easily hide a structure in an S-Box.

To keep an advantage in implementation (WB crypto)... ... or an advantage in cryptanalysis (backdoor). Summary

We can recover parts of the design process of an S-Box using some statistics.

1 The two tables (basics of Boolean functions for cryptography)

2 A satistical tool based on the two tables The Two Tables Let S : F n 2 → F n 2 be an S-Box.

Definition (DDT)

The

Difference Distribution Table of S is a matrix of size 2 n × 2 n such that DDT[a, b] = #{x ∈ F n 2 | S (x ⊕ a) ⊕ S(x) = b}.
Building Blocks for Symmetric Cryptography Statistics and Skipjack TU-Decomposition and Kuznyechik The Butterfly Permutations and Functions Conclusion

The Two Tables Statistical Analysis of the Two Tables Application to Skipjack

The Two Tables Let S : F n 2 → F n 2 be an S-Box.

Definition (DDT)

The

Difference Distribution Table of S is a matrix of size 2 n × 2 n such that DDT[a, b] = #{x ∈ F n 2 | S (x ⊕ a) ⊕ S(x) = b}.

Definition (LAT)

The Example [4,2,1,6,0,5,7,3] The DDT of S. Coefficient Distribution in the DDT If an n-bit S-Box is bijective, then its DDT coefficients behave like independent and identically distributed random variables following a Poisson distribution:

Linear Approximations Table of S is a matrix of size 2 n × 2 n such that LAT[a, b] = #{x ∈ F n 2 | x • a = S(x) • b} -2 n-
S =
        8 0 0 0 0 0 0 0 0 0 0 0 2 2 2 2 0 0 0 0 2 2 2 2 0 0 4 4 0 0 0 0 0 0 0 0 2 2 2 2 0 4 4 0 0 0 0 0 0 4 0 4 0 0 0 0 0 0 0 0 2 2 2 2         The LAT of S.         4 0 0 0 0 0 0 0 0 0 2 2 0 0 2 -2 0 2 2 0 0 2 -2 0 0 2 0 2 0 -2 0 2 0 2 0 -2 0 -2 0 -2 0 -2 2 0 0 -2 -2 0 0 0 -2 2 0 0 -2 -2 0 0 0 0 -4 0 0 0        
Pr [DDT[a, b] = 2z] = e -1/2 2 z z
.

Always even, ≥ 0 Typically between 0 and 16.

Lower is better. Coefficient Distribution in the LAT If an n-bit S-Box is bijective, then its LAT coefficients behave like independent and identically distributed random variables following this distribution:

Pr [LAT[a, b] = 2z] = ( 2 n-1 2 n-2+z ) ( 2 n 2 n-1
) .

The Two Tables Statistical Analysis of the Two Tables Application to Skipjack Coefficient Distribution in the LAT If an n-bit S-Box is bijective, then its LAT coefficients behave like independent and identically distributed random variables following this distribution:

Pr [LAT[a, b] = 2z] = ( 2 n-1 2 n-2+z ) ( 2 n 2 n-1
) .

Always even, signed.

Typically between -40 and 40.

Lower absolute value is better. Reverse-Engineering F For Skipjack's F, max(LAT) = 28 and #28 = 3. What Can We Deduce?

F has not been picked uniformly at random.

F has not been picked among a feasibly large set of random S-Boxes.

Its linear properties were optimized (though poorly). What Can We Deduce?

F has not been picked uniformly at random.

F has not been picked among a feasibly large set of random S-Boxes.

Its linear properties were optimized (though poorly).

The S-Box of Skipjack was built using a dedicated algorithm.

Common ground

Both are standard symmetric primitives in Russia.

Both were designed by the FSB (TC26).

Both use the same 8 × 8 S-Box, π. The TU-Decomposition

Definition

The TU-decomposition is a decomposition algorithm working against S-Boxes with vector spaces of zeroes in their LAT.

ω σ ϕ ⊙ ν 1 ν 0 I ⊙ α ⊙ Multiplication in F 2 4 α Linear permutation I Inversion in F 2 4 ν 0 , ν 1 , σ 4 × 4 permutations ϕ 4 × 4

Belarussian inspiration

The last standard of Belarus (BelT) uses an 8-bit S-box, somewhat similar to π... The Big APN Problem Definition (APN function)

ω ′ ⊗ -1 ⊞ q ′ log w,16
A function S :

F n 2 → F n 2 is Almost Perfect Non-linear (APN) if S(x ⊕ a) ⊕ S(x) = b
has 0 or 2 solutions for all a ̸ = 0 and for all b.

The Big APN Problem Definition (APN function)

A function S :

F n 2 → F n 2 is Almost Perfect Non-linear (APN) if S(x ⊕ a) ⊕ S(x) = b
has 0 or 2 solutions for all a ̸ = 0 and for all b.

Big APN Problem

Are there APN permutations operating on F n 2 where n is even?

Dillon et al.'s Permutation

Only One Known Solution!

For n = 6, Dillon et al. found an APN permutation.

Dillon et al.'s Permutation

Only One Known Solution!

For n = 6, Dillon et al. found an APN permutation.

Dillon et al.'s Permutation

Only One Known Solution!

For n = 6, Dillon et al. found an APN permutation.

Dillon et al.'s Permutation

Only One Known Solution!

For n = 6, Dillon et al. found an APN permutation.

It is possible to make a TU-decomposition!

On the Butterfly Structure

βx 3 x 1/3 ⊙ α ⊕ ⊕ βx 3 x 3 ⊙ α ⊕ ⊕ Definition (Open Butterfly H 3 α,β )
This permutation is an open butterfly.

On the Butterfly Structure

βx 3 x 1/3 ⊙ α ⊕ ⊕ βx 3 x 3 ⊙ α ⊕ ⊕ T U Definition (Open Butterfly H 3 α,β )
This permutation is an open butterfly.

Lemma

Dillon's permutation is affine-equivalent to H 3 w,1 , where Tr (w) = 0.

Closed Butterflies

⊙ α ⊕ x 3 βx 3 ⊕ ⊙ α ⊕ x 3 βx 3 ⊕ Definition (Closed butterfly V 3 α,β )
This quadratic function is a closed butterfly.

Closed Butterflies

⊙ α ⊕ x 3 βx 3 ⊕ ⊙ α ⊕ x 3 βx 3 ⊕ Definition (Closed butterfly V 3 α,β )
This quadratic function is a closed butterfly.

Lemma (Equivalence)

Open and closed butterflies with the same parameters are CCZ-equivalent.

Some Properties of Butterflies

Theorem (Properties of butterflies)

Let V 3 α,β and H 3 α,β be butterflies operating on 2n bits, n odd. Then:

deg ( V 3 α,β ) = 2, if n = 3, Tr (α) = 0 and β + α 3 ∈ {α, 1/α}, then max(DDT) = 2, max(W) = 2 n+1 and deg ( H 3 α,β ) = n + 1 , if β = (1 + α) 3 , then max(DDT) = 2 n+1 , max(W) = 2 (3n+1

  s F, max(LAT) = 28 and #28 = 3. = 28, #28 ≤ N28] Pr [max(LAT) = 28 and #28 ≤ 3] ≈ 2 -55

T

  and U are mini-block ciphers ; µ and η are linear permutations. Decomposing the Mysterious S-Box Final Decomposition Number 1

  b c d e T f 0 f 1 2 3 4 5 6 7 8 9 a b c d e Decomposing the Mysterious S-Box Conclusion on Kuznyechik/Stribog π Decomposing the Mysterious S-Box Conclusion on Kuznyechik/Stribog π Decomposing the Mysterious S-Box Conclusion on Kuznyechik/Stribog π Decomposing the Mysterious S-Box

  1 .

	Building Blocks for Symmetric Cryptography
	Statistics and Skipjack	The Two Tables
	TU-Decomposition and Kuznyechik	Statistical Analysis of the Two Tables
	The Butterfly Permutations and Functions	Application to Skipjack
	Conclusion
		14 / 46

  function ω Linear permutation Decomposing the Mysterious S-Box

	Building Blocks for Symmetric Cryptography Building Blocks for Symmetric Cryptography
	Statistics and Skipjack TU-Decomposition and Kuznyechik The Butterfly Permutations and Functions Statistics and Skipjack TU-Decomposition and Kuznyechik The Butterfly Permutations and Functions	Streebog and Kuznyechik Decomposing the Mysterious S-Box Streebog and Kuznyechik Decomposing the Mysterious S-Box Decomposing the Mysterious S-Box
	Conclusion Conclusion
	Hardware Performance Conclusion for Kuznyechik/Stribog? Conclusion for Kuznyechik/Stribog? Conclusion for Kuznyechik/Stribog?
	The Russian S-Box was built like a The Russian S-Box was built like a The Russian S-Box was built like a
	strange Feistel... strange Feistel... strange Feistel...
	Structure	Area (µm 2 ) ... or was it? ... or was it?	Delay (ns)
	Naive implementation	3889.6	362.52
	Feistel-like	1534.7	61.53
	Multiplications-first	1530.3	54.01
	Feistel-like (with tweaked MUX)	1530.1	46.11

  )/2 and deg If you choose a cipher, demand a full design explanation.

	Building Blocks for Symmetric Cryptography Building Blocks for Symmetric Cryptography Building Blocks for Symmetric Cryptography Building Blocks for Symmetric Cryptography Building Blocks for Symmetric Cryptography
	Statistics and Skipjack Statistics and Skipjack Statistics and Skipjack Statistics and Skipjack Statistics and Skipjack
	TU-Decomposition and Kuznyechik TU-Decomposition and Kuznyechik TU-Decomposition and Kuznyechik TU-Decomposition and Kuznyechik TU-Decomposition and Kuznyechik	Conclusion Conclusion Conclusion Conclusion Conclusion
	The Butterfly Permutations and Functions The Butterfly Permutations and Functions The Butterfly Permutations and Functions The Butterfly Permutations and Functions The Butterfly Permutations and Functions
	Conclusion Conclusion Conclusion Conclusion Conclusion
	Outline Open Problem Conclusion Conclusion Conclusion The Last S-Box
	14 11 60 6d e9 10 e3 2 b 90 d 17 c5 b0 9f c5
	d8 da be 22 8 f3 4 a9 fe f3 f5 fc bc 30 be 26
	bb 88 85 46 f4 2e e fd 76 fe b0 11 4e de 35 bb Building Blocks for Symmetric Cryptography Statistics and Skipjack 1 Cryptographers use mathematics but mathematicians could also use crypto! 1 Cryptographers use mathematics but mathematicians could also use crypto! 30 4b 30 d6 dd df df d4 90 7a d8 8c 6a 89 30 39 1 Cryptographers use mathematics but mathematicians could e9 1 da d2 85 87 d3 d4 ba 2b d4 9f 9c 38 8c 55 also use crypto! d3 86 bb db ec e0 46 48 bf 46 1b 1c d7 d9 1b e0
	23 d4 d7 7f 16 3f 3 3 44 c3 59 10 2a da ed e9
	TU-Decomposition and Kuznyechik 8e d8 d1 db cb cb c3 c7 38 22 34 3d db 85 23 7c
	( 24 d1 d8 2e fc 44 8 38 c8 c7 39 4c 5f 56 2a cf H 3 α,β ) = n , The Butterfly Permutations and Functions d0 e9 d2 68 e4 e3 e9 13 e2 c 97 e4 60 29 d7 9b
	otherwise, d9 16 24 94 b3 e3 4c 4c 4f 39 e0 4b bc 2c d3 94
	max(DDT) = 4, max(W) = 2 n+1 and deg ) 81 96 93 84 91 d0 2e d6 d2 2b 78 ef d6 9e 7b 72 ( H 3 α,β ) ∈ {n, n + 1} and deg ( Conclusion A hidden structure! ad c4 68 92 7a d2 5 2b 1e d0 dc b1 22 3f c3 c3 H 3 α,β = n if and only if 88 b1 8d b5 e3 4e d7 81 3 15 17 25 4e 65 88 4e
	e4 3b 81 81 fa 1 1d 4 22 0 6 1 27 68 27 2e
	3b 83 c7 cc 25 9b d8 d5 1c 1f e5 59 7f 3f 3f ef
		44 / 46 44 / 46

1 + αβ + α 4 = (β + α + α 3 ) 2 .

CMEA uses an 8-bit (non-bijective) S-Box... With a TU-decomposition! What is its actual structure?

2 If you design a cipher, justify every step of your design.

2 If you design a cipher, justify every step of your design. 3

S-Box Design

be an S-Box.


Versatile Block Cipher A Lightweight

T Suzaki

K Minematsu

S Morioka

E Kobayashi

Observations on GFNs Multi-Rotating Feistel Network (MRFN) Possible Applications Conclusion

Conclusion

Di usion in Generalized Feistel networks How long does it take for each input word to influence each output word?

The state consists of 2b branches.

Observations on GFNs

Multi-Rotating Feistel Network (MRFN)

Possible Applications Conclusion

Di usion in Generalized Feistel networks How long does it take for each input word to influence each output word?

The state consists of 2b branches. Nyberg/Type-II GFN: ≈ 2b rounds Observations on GFNs Multi-Rotating Feistel Network (MRFN) Possible Applications Conclusion

Di usion in Generalized Feistel networks How long does it take for each input word to influence each output word?

The state consists of 2b branches.

General Vue

π X i 0 f ⊕ X i 4 X i 1 f ⊕ X i 5 X i 2 f ⊕ X i 6 X i 3 f ⊕ X i 7
Optimal Di usion

The best we can achieve is for X 0 0 to influence ϕ i+2 branches at round i, where ϕ 0 = 0, ϕ 1 = 1, ϕ i+2 = ϕ i+1 + ϕ i .

Di usion in GFNs

(i < r , j < b) Sequence s i of rotations of b words.
The round i of a MRFN with b = 4 and s i = 1 is: Fibonacci Case

f i 0 f i 1 f i 2 f i 3 Introduction Observations
At round 0, X 0 0 has touched the first ϕ 1 = 1 branches of one side.

ϕ i+1 ϕ i X i X i-1 ϕ i ϕ i + ϕ i+1 ϕ i+2 ϕ i+1 X i+1 X i F i ⊕ ϕ 0 = ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ϕ 1 = ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ϕ 2 = ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ϕ 3 = ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ϕ 4 = ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ w b b
VRound function operating on 2bw bit internal state.

w b b f i 1 2. parallel layer of f i f i 2 2. parallel layer of f i f i 3 2. parallel layer of f i f i 4 2. parallel layer of f i f i 5 2. parallel layer of f i f i 6 2. parallel layer of f i f i 7 2. parallel layer of f i f i 8 2. parallel layer of f i f i 9 2. parallel layer of f i f i 10 2. parallel layer of f i w b b ≪ s i 3. rotations ≪ s i 3. rotations ≪ s i 3. rotations w b b ⊕ 4. XOR ⊕ 4. XOR ⊕ 4. XOR s i and s i + (-) i mod b are equivalent if gcd(s i , b) 1 for all i, no full di usion! Importance of the choice of {s i } i ≥0 Introduction Observations on GFNs Multi-Rotating Feistel Network (MRFN) Possible Applications Conclusion Security If s i = ϕ i , then full di usion in ≈ Λ(n) rounds, where Λ(x ) = i if ϕ i-1 < x ≤ ϕ i (optimal).
If s 2i = 0 and i 2i+1 = 2 i , then full di usion in ≈ 2 log 2 (n) rounds (like TWINE).

Both are quickly safe from miss-in-the-middle based impossible di erential a acks and MitM!

If s i = ϕ i , then full di usion in ≈ Λ(n) rounds, where Λ(x ) = i if ϕ i-1 < x ≤ ϕ i (optimal).
If s 2i = 0 and i 2i+1 = 2 i , then full di usion in ≈ 2 log 2 (n) rounds (like TWINE).

Both are quickly safe from miss-in-the-middle based impossible di erential a acks and MitM! When s i = ϕ i , bad truncated di erential with 2 active S-Boxes/round.

If s i = ϕ i , then full di usion in ≈ Λ(n) rounds, where Λ(x ) = i if ϕ i-1 < x ≤ ϕ i (optimal).
If s 2i = 0 and i 2i+1 = 2 i , then full di usion in ≈ 2 log 2 (n) rounds (like TWINE).

Both are quickly safe from miss-in-the-middle based impossible di erential a acks and MitM! When s i = ϕ i , bad truncated di erential with 2 active S-Boxes/round.

Open

GFN-based Linear Layers

Use linear { f i } i ≥0 ; s i = ϕ i n-bit block divided into 2b branches of w bits uses:

w 2 2 f i j ×b f layer × 2 log 2 (b) r XORs .

GFN-based Linear Layers

Use linear { f i } i ≥0 ; s i = ϕ i n-bit block divided into 2b branches of w bits uses:

w 2 2 f i j ×b f layer × 2 log 2 (b) r
XORs .

If we fix w to a small value, then the number of XORs scales with n log 2 (n) rather than n 2 .

GFN-based Linear Layers

Use linear { f i } i ≥0 ; s i = ϕ i n-bit block divided into 2b branches of w bits uses:

w 2 2 f i j ×b f layer × 2 log 2 (b) r
XORs .

If we fix w to a small value, then the number of XORs scales with n log 2 (n) rather than n 2 .

Practical gains even for n = 256:

f i j (x ) = χ 3 (x ⊕ c i j )
s 2i = 0, s 2i+1 = 2 i for 0 ≤ i < 2 log 2 (b) = 12, then repeat (4? times): s = {0, 1, 0, 2, 0, 4, 0, 8, 0, 16, 0, 32} n = 384, with b = 64 and w = 3

f i j (x ) = χ 3 (x ⊕ c i j )
s 2i = 0, s 2i+1 = 2 i for 0 ≤ i < 2 log 2 (b) = 12, then repeat (4? times): s = {0, 1, 0, 2, 0, 4, 0, 8, 0, 16, 0, 32}

E iciency estimates

On 64-bit processors, for each round:

3

Conclusion

Fun stu happens when we allow the use of di erent permutations in each round!

Open problems

What are good sequences of rotations?

How to bound number of active f functions?

What can we use it for?

What happens in other structures (SPN? ARX?) when the linear layers are round-dependent?

Thank you!

n

  Improvements to the Linear Layer of LowMC: A Faster Picnic, with Angela Promitzer, Sebastian Ramacher and Christian Rechberger (2017/448) = 384, with b = 64 and w = 3


Michael Burr
email: burr2@clemson.edu

Shuhong Gao
email: sgao@clemson.edu

Elias Tsigaridas
email: elias.tsigaridas@inria.fr

The Complexity of Subdivision for Diameter-Distance Tests

Keywords: Algorithmic complexity, worst case bit-complexity, subdivision algorithms, continuous amortization, adaptive complexity, separation bounds

We present a general framework for analyzing the complexity of subdivision-based algorithms whose tests are based on the sizes of regions and their distance to certain sets (often varieties) intrinsic to the problem under study. We call such tests diameter-distance tests. We illustrate that diameter-distance tests are common in the literature by proving that many interval arithmeticbased tests are, in fact, diameter-distance tests. For this class of algorithms, we provide both non-adaptive bounds for the complexity, based on separation bounds, as well as adaptive bounds, by applying the framework of continuous amortization.

Using this structure, we provide the first complexity analysis for the algorithm by Plantinga and Vegeter for approximating real implicit curves and surfaces. We present both adaptive and non-adaptive a priori worst-case bounds on the complexity of this algorithm both in terms of the number of subregions constructed and in terms of the bit complexity for the construction. Finally, we construct families of hypersurfaces to prove that our bounds are tight.

Introduction

Subdivision-based algorithms are adaptive methods that start with a domain of interest (often an axis-aligned box) and recursively split it into sub-domains until each sub-domain either isolates or does not contain an interesting feature of the problem at hand. The output is a partition of (a) (b)

Figure 1: (a) The real curve given by the zeros of f = 3y 3 +3xy 2 -2x 3 -3y 2 + xy+3x 2 -3y+3x+2. (b) The approximation produced by the PV algorithm [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] as well as the regions constructed by the algorithm.

Related Work

For univariate problems, the analysis of subdivision-based algorithms is well-understood, and there are several results, especially for the case of approximating the roots of polynomials, see, e.g., [START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF]; [START_REF] Sagraloff | A simple but exact and efficient algorithm for complex root isolation and its complexity analysis[END_REF]; [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF]; [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]; [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF]; [START_REF] Du | Amortized bound for root isolation via Sturm sequences[END_REF], and the references therein. Moreover, it also possible to modify the subdivision process by applying the Newton operator, see [START_REF] Sagraloff | When newton meets descartes: A simple and fast algorithm to isolate the real roots of a polynomial[END_REF]; [START_REF] Pan | Approximating complex polynomial zeros: modified weyl's quadtree construction and improved newton's iteration[END_REF], and considerably improve both the complexity and the actual running time of the corresponding algorithms. However, in higher dimensions, very little is known. For example, there are no explicit complexity results for pure subdivision-based algorithms for approximating curves and surfaces.

The design of efficient subdivision-based algorithms that are output-sensitive, precisionsensitive, certified, and exploit the underlying structure of the problem is an important challenge and an active area of research. An important step in this direction was the introduction of soft tests, see [START_REF] Wang | On soft predicates in subdivision motion planning[END_REF]; [START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF], that, roughly speaking, replace harder exact tests (usually comparisons with zero) with approximate computations which are exact in the limit. They introduce a new notion of correctness called resolution-exactness. In this context, it is exactly the continuous amortization tool [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF][START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] that captures the complexity of the soft predicates. Therefore, continuous amortization is a key tool for the analysis of such algorithms.

The previous work on subdivision-based methods and inclusion-exclusion predicates is quite extensive, so, we can only scratch its surface. For work that focuses on classical inclusionexclusion algorithms for the isolation of roots of algebraic and analytic functions, but without bit-complexity bounds, we refer the interested reader to [START_REF] Giusti | On location and approximation of clusters of zeros of analytic functions[END_REF]; [START_REF] Yakoubsohn | Approximating the zeros of analytic functions by the exclusion algorithm[END_REF]; [START_REF] Dedieu | Computing the real roots of a polynomial by the exclusion algorithm[END_REF]; [START_REF] Yap | Analytic root clustering: A complete algorithm using soft zero tests[END_REF], and the references therein. For other approaches for approximating curves and surfaces, we refer the interested reader to [START_REF] Cheng | Sampling and meshing a surface with guaranteed topology and geometry[END_REF]; [START_REF] Boissonnat | Isotopic implicit surface meshing[END_REF][START_REF] Boissonnat | Meshing of surfaces[END_REF]; [START_REF] Cheng | Certified rational parametric approximation of real algebraic space curves with local generic position method[END_REF], and the references therein. For the problem of isolating the roots of polynomials with subdivision-based methods, we refer the inter-ested reader to [START_REF] Mantzaflaris | On continued fraction expansion of real roots of polynomial systems, complexity and condition numbers[END_REF]; [START_REF] Krawczyk | Newton-Algorithmen zur Bestimmung von Nullstellen mit Fehlerschranken[END_REF]; [START_REF] Henrici | Methods of search for solving polynomial equations[END_REF]; [START_REF] Collins | Polynomial real root isolation using Descarte's rule of signs[END_REF]; [START_REF] Yakoubsohn | Numerical analysis of a bisection-exclusion method to find zeros of univeriate analytic functions[END_REF]; [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF]; [START_REF] Burr | Sqfreeeval: An (almost) optimal real-root isolation algorithm[END_REF]; [START_REF] Sagraloff | A simple but exact and efficient algorithm for complex root isolation and its complexity analysis[END_REF]; [START_REF] Cheng | Root isolation of zero-dimensional polynomial systems with linear univariate representation[END_REF], and the references therein. There are also approaches, see, e.g., [START_REF] Mourrain | Subdivision methods for solving polynomial equations[END_REF][START_REF] Mantzaflaris | On continued fraction expansion of real roots of polynomial systems, complexity and condition numbers[END_REF], that achieve locally quadratic convergence towards the simple roots of polynomial systems, and they are very efficient in practice. Another interesting direction for the application of subdivision-based algorithms, of a more geometric nature, concerns the approximation of algebraic varieties [START_REF] Snyder | Interval analysis for computer graphics[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]Burr et al., 2012;[START_REF] Sharma | Isotopic arrangement of simple curves: an exact numerical approach based on subdivision[END_REF][START_REF] Yap | Adaptive isotopic approximation of nonsingular curves: the parameterizability and nonlocal isotopy approach[END_REF][START_REF] Lin | Non-local isotopic approximation of nonsingular surfaces[END_REF] and the computation of the approximate Voronoi diagrams [START_REF] Yap | Towards exact numerical voronoi diagrams[END_REF]. There are also important applications of these algorithms to the problem of robotic motion planning (Wang et al., 2015a).

For the related problem of computing the topology of an implicitly defined curve in the plane, we refer the reader to [START_REF] Bouzidi | Solving bivariate systems using rational univariate representations[END_REF] for state-of-the-art results. Nevertheless, we emphasize, that even though analyzing the topology of an implicitly defined curve is related to the problem we consider in this paper, the problems and approaches are different and the complexity estimates are not directly comparable. Our approach is a general one that we can use for the analysis of any subdivision-based algorithm that uses diameter-distance tests and it is not a dedicated one for computing the topology of curves. This argument also holds for methods and algorithms based on cylindrical algebraic decomposition, which can be used as a black-box tool to solve similar problems with curves, see [START_REF] Hong | Variant quantifier elimination[END_REF] and references therein.

Main Results

We introduce diameter-distance tests, which formalize a type of test that is frequently used in subdivision-based algorithms, see Section 2.2. We then present both straightforward nonadaptive complexity bounds for such tests based on separation bounds, see Proposition 3, and adaptive bounds, based on continuous amortization, which exploit the local features of the problem at hand in Proposition 7. The diameter-distance tests are quite generic in nature and we illustrate this by formulating classical exclusion predicates, in any dimension, in conjunction with interval arithmetic as diameter-distance tests, see Section 3.2.

We provide the first complexity analysis for (a slightly modified version of) the PV algorithm for approximating curves and surfaces from Plantinga andVegter (2007, 2004). We extend the predicates of the PV algorithm to all dimensions and bound the number of regions and bit-complexity of these algorithms in two-and higher-dimensions, see Theorems 24 and 27. Moreover, using continuous amortization, first developed by [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF], we provide adaptive bounds on the number of regions and the bit-complexity of the PV algorithm in arbitrary dimensions, see Theorems 25 and 29. These results consist of the first application of continuous amortization to a pure high-dimensional problem. We provide examples that show that our bounds are tight in Lemma 30.

We anticipate that diameter-distance tests and the tools for the complexity analysis of the underlying subdivision-based algorithms that we develop in this paper will be applicable to many other algorithms and in related contexts.

Overview of Paper

The rest of this paper is organized as follows: In the next section, we present a general description of subdivision-based algorithms, we introduce diameter-distance tests, and we derive adaptive and non-adaptive complexity bounds for subdivision-based algorithms that use these tests. In Section 3, we show that exclusion tests based on interval arithmetic are diameterdistance tests. This illustrates that many algorithms in the literature can be analyzed with the techniques of this paper. In Section 4, we present a slight modification of the PV algorithm for curve and surface approximation. We then exhibit the tests in the PV algorithm as diameterdistance tests. In Section 5, we present both adaptive and non-adaptive bounds on the number of subdivisions that the PV algorithm performs and the bit-complexity for the overall algorithm. Finally, in Section 6 we present examples to demonstrate the tightness of our bounds.

Subdivision-based Methods and Diameter-Distance Tests

In Section 2.1, we present the general form of a subdivision-based method which is studied in this paper. In Section 2.2, we define the diameter-distance tests, which form the class of predicates studied in this paper. Even though our motivating example is the Plantinga and Vegter algorithm (Plantinga andVegter (2004, 2007)), we present this material in a general setting. Additional, related, background on this approach for the study of subdivision-based methods in this section can be found in [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] and [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF].

Throughout this section, we assume that X is both a measure space with measure µ and a metric space with distance function d. We note that we do not require any compatibility between µ and d. Additionally, we assume the technical condition that X is proper6 , i.e., closed balls are compact. Moreover, we let S be a collection of subsets of X which have finite measure (with respect to µ) and are compact (with respect to d). We call predicates (boolean functions) on S stopping criteria. In the motivating case, X = R n , µ is the Lebesgue measure, d is the Euclidean distance, and S is the collection of n-dimensional cubes in X.

Subdivision-based Methods

In this section, we provide the general form of a subdivision-based method considered in this paper. Let C 1 , . . . , C be stopping criteria, i.e., each C i is a function from subsets in S to {True, False}. When C i (J) is True for some i, then we do not split J, and, when C i (J) is False for all i, then we must subdivide J.

We use the following simple abstract algorithm to describe subdivision-based tests. Fix C 1 , . . . , C to be stopping criteria, and consider an input region I ∈ S. The output of the algorithm is a partition P of I such that for each element J in P, there is some i so that C i (J) = True. Initially, P = {I}.

Algorithm 1. Abstract Subdivision-based Algorithm

For each J ∈ P,

If there exists 1 ≤ i ≤ such that C i (J) = True, report J.
Otherwise, subdivide J and replace J with its children.

To subdivide a region J means to replace J in P with regions J 1 , . . . , J k , where k ≥ 2, each J j ∈ S, J = ∪ j J j , and the pairwise intersections of the J j are measure zero subsets. In this paper, we add two mild additional assumptions to the subdivisions under consideration: Let 0 < ε 1 , ε 2 < 1. Then, we add the following assumptions: Assumption 1: µ(J j ) ≥ ε 1 µ(J) and Assumption 2: Diam(J j ) ≤ ε 2 Diam(J).

The first condition prevents J from splitting into two many regions at any step, while the second condition generalizes the idea that the aspect ratio of J i should not be too large. For additional details on the first assumption, see Burr (2016, Lemma 3.5 and Remark 3.6). The subdivision tree is the tree whose root is I, whose internal nodes represent sub-domains J that are processed during the subdivision, whose leaves are the terminal regions, and where the parent-child relationship is given by subdivision. We observe that, in the motivating example for this paper,

ε 1 = 2 -n and ε 2 = 2 -1 .

Diameter-distance Tests

In this section, we define distance-diameter tests. These tests form the class of predicates that we consider in this paper. Many tests which have been developed, such as the one-circle condition in Descartes' rule of signs [START_REF] Krandick | New bounds for the Descartes method[END_REF][START_REF] Alesina | Vincent's theorem from a modern point of view[END_REF], are diameter-distance tests. At first glance, it might appear that the definition of diameter-distance tests is very specialized; this is not the case. In fact, in Section 3, we provide a nontrivial example of these tests which appear frequently in applications.

Definition 2. Let X and S be defined as above. Let C be a stopping criterion on S. C is a diameter-distance test if there exists a closed set V ⊆ X and a positive constant K such that for any J ∈ S,

If Diam(J) < K max x∈J d(x, V) , then C(J) = True, where d(x, V) = min v∈V d(x, v).
The extra conditions, such as compactness and properness allow one to use the minimum and maximum in the definition above instead of infimum and supremum. Loosely speaking, this definition states that C(J) must be True when J includes a point sufficiently far away from V and J isn't too large.

We note that this definition does not state that the stopping criterion C must be a distancebased test or even that V is known to C. Instead, the only assumption is that the criterion is less conservative than the conditional in the definition above. In particular, stopping criteria whose theory is based upon condition numbers are frequently diameter-distance tests because the condition number can be rewritten in terms of the inverse of the distance to the set of illconditioned inputs [START_REF] Bürgisser | Condition: The geometry of numerical algorithms[END_REF].

Throughout the remainder of this section, we assume that all stopping criteria are diameterdistance tests.

Non-adaptive Bounds

In this section, we provide a lower bound on the number of regions produced by Algorithm 1. This analysis is not adaptive, so it assumes the worst-case behavior everywhere. We include this approach for comparison because the adaptive bounds are based on the ideas of the non-adaptive bounds and, in some cases, the adaptive bounds may be too complicated to compute. In the next section, we provide an adaptive bound based on continuous amortization [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF].

Proposition 3. Suppose that the stopping criteria C 1 , . . . , C in Algorithm 1 are all diameterdistance tests, with associated positive constants K i and closed subsets V i . Furthermore, assume that the intersection V i is empty. Let K = min K i , and let I be the initial input region. Define the separation bound δ as 0 < δ ≤ min

x∈I max i d(x, V i ).
Then, the number of regions constructed by the algorithm is at most

max 1, ε -1+ ln Diam(I)-ln(Kδ) ln(ε 2 ) 1 .
Before presenting the proof, we note that δ is a lower bound on the smallest distance from any x ∈ I to the furthest V i . We call δ a separation bound because if the V i 's are pairwise disjoint and ∆ is the minimum distance between them, then ∆ 2 satisfies the conditions for δ. Proof. If no subdivisions occur, then the only region is the initial region. In this case, 1 region is constructed and the bound holds. We, therefore, assume that subdivisions occur. Let J j be a terminal region and J its parent; moreover, let x ∈ J j . Since x ∈ J and J was subdivided, we know that Diam(J) ≥ K max i d(x, V i ), since, otherwise, by the definition of a diameter-distance test, for some i, C i (J) would be True, contradicting the assumption that J was subdivided. Since this maximum is larger than δ, it follows that Diam(J) ≥ Kδ.

Suppose that the depth of J j in the tree is k, then the depth of J is k-1, and, by the assumption on diameters under subdivisions, we know that Diam(J) ≤ ε k-1 2 Diam(I). Therefore, Kδ ≤ ε k-1 2 Diam(I). Taking the logarithm of both sides (and recalling that ln(ε 2 ) < 0), it follows that k ≤ 1 + ln(Kδ)ln(Diam(I)) ln(ε 2 ) .

By the assumption on volumes for subdivisions, it follows that µ(J j ) ≥ ε k 1 µ(I). Substituting in our expression for k, we can conclude that

µ(J j ) ≥ ε 1+ ln(Kδ)-ln(Diam(I)) ln(ε 2 ) 1 µ(I).
This lower bound applies to the measure of every terminal region. Moreover, since the pairwise intersection of terminal regions has zero measure, we know that, in the worse-case, I is subdivided into regions of this size, which results in the desired bound.

We observe that in the motivating case, i.e., where

ε 1 = 2 -n and ε 2 = 2 -1 , the maximum above simplifies to max 1, 2 √ n Kδ n µ(I) (1)
since Diam(I) n = n n/2 µ(I). In many cases, however, this bound is much larger than necessary as the analysis assumes that the worst-case situation occurs everywhere. An adaptive bound is necessary to account for this non-uniformity.

Adaptive Bounds

In this section, we present adaptive bounds for the number of regions produced by Algorithm 1. This adaptive bound is based on the continuous amortization technique [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF], which we briefly review here.

Continuous amortization was introduced in [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF] as a way to adaptively analyze the complexity of subdivision-based algorithms. In [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF], the theory of continuous amortization was extended to apply to measure spaces and to evaluate functions on the regions of the partition, and we recall this technique here. The key to continuous amortization is a function on X, called a local size bound, which is a point estimate, locally describing the worst-case amount of work that is required at each point.

Definition 4. Let X and S be defined as above and C a stopping criterion on S. A local size bound for C is a function F : X → R ≥0 with the property that

F(x) ≤ inf J∈S J x C(J)=False µ(J).
In other words, F(x) is a lower bound on the measure of a region which contains x, but fails the stopping criterion. The local size bound provides the link between the algorithm and a quantity that we can compute.

Theorem 5 [START_REF] Burr | Continuous amortization: A non-probabilistic adaptive analysis technique[END_REF]; [START_REF] Burr | Continuous amortization and extensions: With applications to bisection-based root isolation[END_REF]). Let X and S be defined as above, C a stopping criterion on S, and F a local size bound for C. Let h : R ≥0 → R be a non-increasing function, and let P be the final partition formed by Algorithm 1, which recursively subdivides the input region I, subject to Assumption 1. The sum of h applied to the regions in P is bounded as follows:

J∈P h(µ(J)) ≤ max h(µ(I)), I h(ε 1 F(x)) ε 1 F(x) dµ .
If h ≡ 1, i.e., h is the constant function, then this integral counts the number of regions formed by Algorithm 1. In addition, if the algorithm does not terminate, then the integral is infinite.

We observe that we can use the continuous amortization integral to express the complexity of Algorithm 1 in the particular case where each stopping criterion is a diameter-distance test. Lemma 6. Let X and S be as above and C a stopping criterion on S with constant K and closed set V. Suppose that the subdivisions by Algorithm 1 are subject to the two additional conditions following Algorithm 1. Then,

F(x) = ε 1+ ln(Kd(x,V))-ln(Diam(I)) ln(ε 2 ) 1 µ(I) is a local size bound for C.
Proof Sketch. The proof is identical to that of Proposition 3 except that we begin with the condition that Diam(J) ≥ Kd(x, V) from the definition of a diameter-distance test.

In Algorithm 1, we have multiple stopping criteria. Therefore, for each C i , we can define a local size bound F i : X → R ≥0 . Moreover, since at least one of the stopping criteria must be true, we can take the maximum of all of them for the local size bound for Algorithm 1. In particular, we have the following result: Proposition 7. Suppose that the stopping criteria C 1 , . . . , C in Algorithm 1 are all diameterdistance tests with associated positive constants K i and closed subsets V i . Furthermore, assume that the intersection V i is empty. Let I be the initial input region. Then, the number of regions constructed by the algorithm is at most max 1, µ(I) -1

I min i ε -1+ ln(Diam(I))-ln(K i d(x,V i )) ln(ε 2 ) 1 dµ .
We observe that in the motivating case, i.e., where ε 1 = 2 -n and ε 2 = 2 -1 , the continuous amortization integral simplifies to max 1,

I min i 2 √ n K i d(x, V i ) n dµ .
We apply both the adaptive and non-adaptive bounds to the PV algorithm as a specific example in Section 5.

Interval Methods and Diameter-Distance Tests

In this section, we show that a common exclusion test which is based on the standard centered form is a diameter-distance test. We begin with a brief review of the standard centered form, for more details, see, for example, [START_REF] Ratschek | Computer methods for the range of functions[END_REF]; [START_REF] Moore | Introduction to interval analysis[END_REF].

Let Y be any set, S a collection of subsets of Y, and consider the function f : Y → R. An interval method for f is an algorithm f such that for any subset J ∈ S, f (J) ⊇ f (J), where f (J) is the image of J under f . In other words, f (J) is an over-approximation for the image f on J. In most applications, Y is a metric space, and we add the convergence condition for f , i.e., that for a sequence of domains {J k } which converge to a point p, then { f (J k )} converges to f (p).

In our applications, we consider the case where Y is R n , and the regions in S are axis aligned n-dimensional boxes, i.e., for J ∈ S, J = i [a i , b i ]. In this case, most interval methods use interval arithmetic, i.e., arithmetic operations on intervals that produce the set-theoretic image as an interval. In this section, we focus on the standard centered form for multivariate polynomials [START_REF] Ratschek | Computer methods for the range of functions[END_REF][START_REF] Moore | Introduction to interval analysis[END_REF]. Let f ∈ R[x 1 , . . . , x n ] be a multivariate polynomial of total degree d and J an axis aligned box. Let m = m(J) be the midpoint of J, then the standard centered form for f applied to J is

f (J) = f (m) + d |α|=1 ∂ α f (m) α! (J -m) α ,
where α ∈ N n and the notation is multi-index notation, i.e., |α|

= α i , ∂ α f (m) = ∂ α 1 1 . . . ∂ α n n f (m), α! = (α i )!, and (J -m) α = a i -b i 2 , b i -a i 2 α i . Since (J -m) α is a product of intervals centered at zero, using interval arithmetic, this product simplifies to b i -a i 2 α i [-1, 1]. In the special case
where J is an axis-aligned, n-dimensional cube, i.e., J is a product of n intervals all of the same width w, then all of the factors in (Jm) α are identical, and the standard centered form can be rewritten as:

f (J) = f (m) +         d |α|=1 |∂ α f (m)| α! w 2 |α|         [-1, 1], (2)
The standard centered form is an interval version of a multivariate Taylor expansion centered at m, and the standard centered form has several nice properties including a very structured expression and fast convergence.

In the remainder of this section, we consider the following predicate:

C(J) = True if and only if 0 f (J).
If 0 f (J), then we can directly conclude that 0 f (J); we observe that the converse does not hold in general, but converges in the limit, i.e., if {J k } is a sequence of n-dimensional boxes whose limit is p, then either C(J k ) = False for some k or f (p) = 0. It is often more efficient, in practice, to test C(J) and subdivide J, if necessary, rather than to compute f (J) directly. In the remainder of this section, we prove that C(J) is a diameter-distance test.

Bounds on Coefficients of Powers of Sines and Cosines

In this section, we prove a technical lemma on the magnitudes of the coefficients of sines and cosines. The main result in this section is used in the following section to prove that C(I) is a diameter-distance test.

Lemma 8. Suppose that for all θ, k j=0 a j cos j (θ) sin k-j (θ) ≤ C.

(3)

Then, |a j | ≤ 2 k+1 C.
Proof. Let f (θ) = k j=0 a j cos j (θ) sin k-j (θ). We observe that this is a square-integrable and 2πperiodic function. Moreover, its Fourier coefficients are bounded by 2C since

| f (θ) cos(nθ)| ≤ C and | f (θ) sin(nθ)| ≤ C.
We now observe that cos(x) = 1 2 (e ix + e -ix ) and sin(x) = 1 2i (e ixe -ix ). Therefore,

cos j (θ) sin k-j (θ) = 1 2 k i k-j (e ix + e -ix ) j (e ix -e -ix ) k-j = 1 2 k i k-j k l=0 b l e i(k-2l)x ,
where the b l 's are sums and products of binomial coefficients. Since e i(k-2l)x = cos((k -2l)x) + i sin((k -2l)x), it follows that cos j (θ) sin k-j (θ) has a finite Fourier series whose nonzero terms are of the form cos(nθ) and sin(nθ) where 0 ≤ n ≤ k and kn is even. Therefore, the Fourier series of f (θ) is can be written as follows:

f (θ) = c 0 2 + (k-1)/2 l=0 (c k-2l cos((k -2l)θ) + d k-2l sin((k -2l)θ)) (4)
where all the constants are bounded by 2C.

Suppose that k ≥ n and kn is even. Then, since cos(nx) = ((cos(x) + i sin(x)) n ) and sin(nx) = ((cos(x) + i sin(x)) n ), we have the following7 :

cos(nx) = (sin 2 (x) + cos 2 (x)) k-n 2         n/2 m=0 (-1) m n 2m cos n-2m (x) sin 2m (x)         sin(nx) = (sin 2 (x) + cos 2 (x)) k-n 2         (n-1)/2 m=0 (-1) m n 2m + 1 cos n-1-2m (x) sin 2m+1 (x)         .
We observe that in these expansions, cos(nx) and sin(nx) are written as a linear combination of products of sines and cosines of degree k. Reorganizing these sums, we find that

cos(nx) = k 2 m=0           min{ n 2 ,m} p=max{0,m-k-n 2 } (-1) p n 2p k-n 2 m -p           cos k-2m (x) sin 2m (x) (5)
and

sin(nx) = k-1 2 m=0           min{ n-1 2 ,m} p=max{0,m-k-n 2 } (-1) p n 2p + 1 k-n 2 m -p           cos k-2m-1 (x) sin 2m+1 (x). (6)
We observe that in the formula for cos(nx), the coefficient of cos k-2l (x) sin 2l (x) can be bounded as follows:

min{ n/2 ,m} p=max{0,m-(k-n)/2} (-1) p n 2p k-n 2 m -p ≤ n/2 p=0 n 2p (k-n)/2 q=0 k-n 2 q = n/2 p=0 n -1 2p + n -1 2p -1 (k-n)/2 q=0 k-n 2 q ≤ 2 n-1 2 (k-n)/2 = 2 n+k 2 -1 (7)
Similarly, the coefficient of cos k-2l-1 (x) sin 2l+1 (x) in the formula for sin(nx) is bounded by

min{ (n-1)/2 ,m} p=max{0,m-(k-n)/2} (-1) p n 2p + 1 k-n 2 m -p ≤ (n-1)/2 p=0 n 2p + 1 (k-n)/2 q=0 k-n 2 q = (n-1)/2 p=0 n -1 2p + 1 + n -1 2p (k-n)/2 q=0 k-n 2 q ≤ 2 n-1 2 (k-n)/2 = 2 n+k 2 -1 . (8)
Moreover, we observe that these bounds are independent of m and p, depending only on n and k In order to bound the coefficients a j , we substitute the formulas above for cos(nx) and sin(nx) into the Fourier series for f . In particular, we substitute n = k -2l into 2 n+k 2 -1 to get 2 k-l-1 . Moreover, by considering the powers in Equations ( 5) and ( 6), we conclude that if kj is even, then cos j (θ) sin k-j (θ) only appears in the expansion of the cosine terms (perhaps including the constant term) in the Fourier series for f (θ), while if kj is odd, then cos j (θ) sin k-j (θ) only appears in the expansion of the sine terms in the Fourier series for f (θ). We can then isolate the occurrences of cos j (θ) sin k-j (θ) (there are four cases, depending on the parity of k and j). Then, using the triangle inequality and the upper bounds in Inequalities ( 7) and ( 8), we find that

|a j | ≤ 2C k 2 l=0 2 k-l-1 < 2 k+1 C,
which completes the proof.

Corollary 9. Fix k 0 ∈ N, and suppose that for all θ 1 , . . . , θ m ,

k 0 k 1 =0 k 1 k 2 =0 • • • k m-1 k m =0 a (k 1 ,...,k m ) m j=1 sin k j-1 -k j (θ j ) cos k j (θ j ) ≤ C. Then, a (k 1 ,...,k m ) ≤ 2 m(k 0 +1) C.
Proof. Proof by induction on m; the base case is Lemma 8. For the inductive case, we fix θ 2 , . . . , θ m . For each k 1 , we define

a k 0 -k 1 = k 1 k 2 =0 • • • k m-1 k m =0 a (k 1 ,...,k m ) m j=2 sin k j-1 -k j (θ j ) cos k j (θ j ) .
Then, the given inequality simplifies to

k 0 k 1 =0 a k 0 -k 1 sin k 0 -k 1 (θ 1 ) cos k 1 (θ 1 ) < C. By Lemma 8, |a k 0 -k 1 | ≤ 2 k 0 +1 C. Since θ 2 , .
. . , θ m are fixed, but arbitrary, and the bound does not depend on θ 2 , . . . , θ m , we can apply the inductive hypothesis to

a k 0 -k 1 to give that |a (k 1 ,...,k m ) | ≤ 2 (m-1)(k 1 +1) |a k 0 -k 1 | ≤ 2 (m-1)(k 1 +1)+(k 0 +1) C. Since k 1 ≤ k 0 , the claim follows.

Exclusion Interval Arithmetic Tests are Distance-Diameter Tests

In this section, we use the results of Section 3.1 to prove that the predicate on n-dimensional cubes8 J where C(J) = True if and only if 0 f (J) is a distance-diameter test. In this case, the set in the definition of a distance-diameter test is the complex variety V C ( f ). As our first step, we reduce a higher-dimensional problem to a collection of one-dimensional problems as follows:

Definition 10. Let f ∈ R[x 1 , . . . , x n ], p ∈ R n ,
and v ∈ S n-1 . We define f v (t) to be the univariate polynomial passing through p and in the direction v, i.e., f v (t) = f (p + tv), see Figure 2. Next, we define Σ f v to be the sum of the reciprocals of the complex roots of f v , i.e.,

Σ f v (p) = s∈V C ( f v ) 1 |s| . p v α 3 α 1 α 2 Figure 2: For a polynomial f ∈ R[x,
y] and a point p ∈ R 2 , we consider the roots of f , α 1 , α 2 , α 3 , in the direction of a unit vector v.

In Burr and Krahmer (2012, Lemma 2.1), it was shown that for a univariate polynomial g ∈ R[x] with complex roots V C (g),

g (k) (x) g(x) ≤         α∈V C (g) 1 |x -α|         k .
This link between the Taylor coefficients of g to the geometry of the zero set of g can be extended to the current setting since f v is a univariate polynomial. We introduce the notation dist C (p, f ) to represent the complex distance between the point p and the variety V C ( f ). Explicitly, we have the following lemma:

Lemma 11. Let f ∈ R[x 1 , . . . , x n ], p ∈ R n , and v ∈ S n-1 . Then 1 f (p) d k f (p + tv) dt k t=0 ≤ Σ f v (p) k ≤ deg( f ) dist C (p, f ) k .
Proof. The claim is trivial when k = 0. Since f v is a univariate polynomial, the first inequality follows directly from Burr and Krahmer (2012, Lemma 2.1). The second inequality follows from the fact that in the sum for Σ f v (p), there are at most deg( f ) terms, and each element of the sum is the inverse of the distance between p and a point on V C ( f ), each of which is, in turn, bounded above by the inverse of the distance to the closest point on V C ( f ).

We now use this upper bound along with the results from Section 3.1 to bound individual Taylor coefficients in the multivariate Taylor expansion.

Proposition 12. Let f ∈ R[x 1 , . . . , x n ] and p ∈ R n . Then, for all multi-indices α ∈ N n with k = |α|, 1 f (p) k α ∂ k f ∂x α (p) ≤ 2 (n-1)(|α|+1) deg( f ) dist C (p, f ) |α| .
where k α is the multinomial coefficient.

Proof. Let (θ 1 , . . . , θ n-1 ) ∈ S 1 n-1 . Consider the surjective map S 1 n-1 → S n-1 given by (θ 1 , . . . , θ n-1 ) → x = (x 0 , . . . , x n-1 ) where

x i =         i j=1 cos θ j         sin θ i+1 for 0 ≤ i < n -1 and x n-1 = n-1 j=1 cos θ j . Then, for k 0 = k ≥ k 1 ≥ k 2 ≥ • • • ≥ k n-1 ≥ 0 = k n , x k 0 -k 1 0 x k 1 -k 2 1 . . . x k n-1 -k n n-1 = n-1 j=1 (sin k j-1 -k j θ j cos k j θ j ). ( 9
)
Observe that, by the chain rule for v ∈ S n-1 , 1 f (p)

d k f (p + tv) dt k = 1 f (p) |α|=k k α ∂ k f ∂x α (p)v α .
By Lemma 11, we know that the magnitude of these expressions are bounded above by deg( f )

dist C (p, f ) k . Moreover, since v is a unit vector, there exist k 0 = k ≥ k 1 ≥ k 2 ≥ • • • ≥ k n-1 ≥ 0 = k n
so that v α can be written in the form of Equation ( 9),

v α = n-1 j=1 (sin k j-1 -k j θ j cos k j θ j ). Therefore, 1 f (p) |α|=k k α ∂ k f ∂x α (p)v α
is of the form for Corollary 9 where m = n -1 and k 0 = k. Therefore, the individual terms are bounded by 2 (n-1)(k+1) deg( f )

dist C (p, f ) k , as desired.
With Proposition 12 in hand, we now prove a lower bound on an n-dimensional cube J of width w = w(J) that fails the predicate C(J) in the following corollary:

Corollary 13. Let f ∈ R[x 1 , . . . , x n ] and p ∈ R n . Suppose that 0 < w ≤ dist C (p, f ) ln(1+2 2-2n ) 2 n-1 deg( f ) . Then deg( f ) k=1 |α|=k 1 k! k α 1 f (p) ∂ |α| f ∂x α (p) w 2 k ≤ 1.
Proof. Observe that by the triangle inequality,

deg( f ) k=1 |α|=k 1 k! k α 1 f (p) ∂ |α| f ∂x α w 2 k ≤ deg( f ) k=1 |α|=k 1 k! k α 1 f (p) ∂ |α| f ∂x α w 2 k .
Now, we can substitute the bound on the derivatives in Proposition 12 as well as the assumed bound on w, resulting in an upper bound of

deg( f ) k=1 |α|=k 1 k! 2 (n-1)(k+1) ln(1 + 2 2-2n ) 2 n k .
Since there are n+k-1 k possibilities for α when |α| = k, which can be trivially bounded from above by 2 n+k-1 , we can bound the expression above by

2 2n-2 deg( f ) k=1 1 k! (ln(1 + 2 2-2n )) k .
Since this sum is a truncated version of the Taylor series expansion of e x -1 centered at 0 with x = ln(1 + 22-2n ), this sum is bounded above by 2 2-2n , and, hence, the entire expression is bounded above by 1.

Using Corollary 13, we can develop bounds on the size of a region which guarantees the success of the given predicate. We make this explicit in the following corollary:

Corollary 14. Let f ∈ R[x 1 , . . . ,
x n ] and J ⊆ R n an n-dimensional cube with midpoint m = m(J) and width w = w(J) such that w

≤ dist C (m, f ) ln(1+2 2-2n ) 2 n-1 deg( f )
. Then, C(J) is true.

Proof. From Equation (2), we see that 0 f (J) is equivalent to

        d |α|=1 |∂ α f (m)| | f (m)|α! w 2 |α|         < 1.
This inequality arises because, in Equation ( 2), f (J) is an interval centered at the origin shifted by f (m). In order for 0 to be excluded from this interval, the shift by f (m) must be larger than the half-width of the interval. Dividing both sides by | f (m)|, we get exactly the expression in Corollary 13.

While Corollary 14 gives a test for C(J) = True for an n-dimensional cube, this test is not enough to prove that C is a distance-diameter test because both sides of the inequality involve region J. In particular, the midpoint of the region J appears on the right-hand-side of the inequality. The following lemma changes the right-hand-side of the inequality to depend on any point within J, instead of the midpoint.

Lemma 15. Let f ∈ R[x 1 , . . . , x n ] and J ⊆ R n an n-dimensional cube with midpoint m = m(J) and width w = w(J). Suppose that x ∈ J and C and k are positive constants. If

w ≤ 2Cdist C (x, h) 2k + C √ n , then w ≤ Cdist C (m, h) k .
Proof. We follow the ideas of the argument in [START_REF] Burr | Sqfreeeval: An (almost) optimal real-root isolation algorithm[END_REF]. We observe that

w = 1 + C √ n 2k w - C √ n 2k w ≤ Cdist C (x, h) k - C √ n 2k w = C k dist C (x, h) - √ n 2 w , ( 10
)
where the inequality follows from the assumed upper bound on w. Suppose that α is the closest point of V(h) to m, then, by the triangle inequality,

dist C (m, h) ≥ dist C (x, α) -dist C (x, m). The distance dist C (x, m) is at most the radius of J, which is √ n 2 w. Moreover, the closest point on V(h) to x is distance at most the distance to α, so dist C (x, α) ≥ dist C (x, h). Hence, dist C (m, h) ≥ dist C (x, h) - √ n Corollary 16. Let f ∈ R[x 1 , . . . ,
x n ] and J ⊆ R n an n-dimensional cube with width w = w(J). If there is a point x ∈ J such that

w ≤ 2 ln 1 + 2 2-2n dist C (x, f ) 2 n deg( f ) + √ n ln 1 + 2 2-2n , then, C(J) is true.
Proof. This result follows from Corollary 14 with Lemma 15 by letting C = ln 1 + 2 2-2n and k = 2 n-1 deg( f ).

Since the diameter of an n-dimensional cube whose side is of length w is scaled by √ n, we have the following corollary:

Corollary 17. Let f ∈ R[x 1 , . . . , x n ] and J ⊆ R n an n-dimensional cube. If there is a point x ∈ J such that Diam(J) ≤ 2 √ n ln 1 + 2 2-2n dist C (x, f ) 2 n deg( f ) + √ n ln 1 + 2 2-2n .
Then, C(J) is true. Therefore, C is a diameter-distance test.

The Modified Plantinga-Vegter Algorithm

In this section, we provide a modified form of the PV algorithm [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]) for curve and surface approximation. Our version of the algorithm uses a slightly stronger variant of their, but makes both tests of the appropriate form for the application of Corollary 17. We begin by reviewing the original PV algorithm and then discuss our generalization and adaptation.

The Original PV Algorithm

Let f ∈ R[x, y] or R[x, y, z] be a square-free polynomial such that its real zero set V R ( f ) is smooth (see Footnote 5 for a brief discussion of the requirement for the curve to be bounded for correctness of the approximation). The PV algorithm recursively subdivides an initial input square or cube I with a quad-tree or oct-tree data structure until at least one of the following tests holds on each subregion J (in the literature, these tests are often referred to as C 0 and C 1 ):

C 0 (J) = True if and only if 0 f (J) C 1 (J) = True if and only if 0 ∇ f (J), ∇ f (J) .
For the purposes of curve approximation, when C 0 (J) is True, the variety does not enter the region J, and so J can be discarded. On the other hand, when C 1 (J) holds, the curve or surface does not bend much within the region J.

The PV algorithm is an instance of an Abstract Subdivision-based Algorithm that uses bisection and the two tests C 0 and C 1 , see Algorithm 1. We explicitly include the algorithm here for completeness and to illustrate the simplicity of the approach. Given an input polynomial f and region I, the PV algorithm constructs a partition P of I so that for all regions J of the partition, either C 0 (J) = True or C 1 (J) = True. Initially, P = {I}.

Algorithm 18. Main subdivision of PV algorithm

For each J ∈ P, If there C 0 (J) = True or C 1 (J) = True, report J. Otherwise, subdivide J and replace J with its children.

After every sub-region J satisfies C 0 (J) or C 1 (J), the authors of [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] perform post-processing steps, which include balancing the tree, evaluating the sign of f on the corners of each J in P, and using sign changes along the sides of regions J to detect and approximate the curve or surface. Their approximation is topologically correct for bounded curves as there is an ambient isotopy between the approximation and the variety V R ( f ). Additionally, by further subdivision, the isotopy can be made sufficiently small so that the Hausdorff distance between the approximation and the variety is as small as desired. We note that it is possible to extend the PV algorithm in the plane to provide a topologically correct approximation even when V R ( f ) is unbounded, when V R ( f ) is singular, and when I is not a bounding box, see Burr et al. (2012). In this paper, however, we focus on the original PV algorithm, but without the restriction of a bounded curve.

Our main target is to compute the number of regions that the PV algorithms construct, and not to approximate the curve or surface, per se. Therefore, we focus, exclusively, on the C 0 and C 1 tests and apply them in arbitrary dimensions. More precisely, let f ∈ R[x 1 , . . . , x n ] be such that its real zero set V R ( f ) is smooth. Let I ⊆ R n be an n-dimensional real cube. Then, we can generalize the tests C 0 and C 1 , along with Algorithm 18, to n dimensions, where the subdivision splits an n-dimensional cube into 2 n children. We mention that, in this case, we no longer use the output of the algorithm to construct an approximation to V R ( f ).

Modifying the C 1 test

As presented above, the C 0 test is of the form considered in Corollary 17, so it is a diameterdistance test. On the other hand, the C 1 test is not of this form; therefore, it is not clear if the C 1 test is a diameter-distance test. The difficulty in applying the corollary in this case is that arithmetic operations are performed on intervals after an application of interval methods. In this section, we describe an alternate C 1 test that satisfies the assumptions of Corollary 17.

The predicate C 1 (J) has the following two consequences that are fundamental in the proof of the correctness of the PV algorithm in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF]:

1. If a region J satisfies the C 1 condition, then, in J, there cannot be any pair of gradient vectors of f which are orthogonal to each other.

2. The variety V R ( f ) is parametrizable in the direction of at least one of the coordinate axes.

Fact 2 is a direct consequence of Fact 1, but it is used so frequently in the proofs in [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], that it is worthwhile to mention it explicitly. We now modify the C 1 test in arbitrary dimensions so that it has the form in the assumptions of Corollary 17. Consider the function g : R n × R n → R, defined as

g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) .
It follows that, if, for a region J, 0 g(J × J), then there is no pair of gradient vectors of f in J which are orthogonal to each other. Thus the modified C 1 test, briefly denoted C 1 , is as follows:

C 1 (J) = True if and only if 0 g(J × J).

Therefore, when C 1 (J) is true, we can conclude the truth of Facts 1 and 2, and the application of an interval method appears as the last step as opposed to in an intermediate step. Therefore, C 1 satisfies the assumptions in Corollary 17, and, therefore, is a diameter-distance test. For the rest of the paper, all references to the C 1 test refer to this new C 1 test. In particular, all discussions of the PV algorithm refer to the modified PV algorithm.

Let J be an n-dimensional real cube with midpoint m = m(J) and side length w = w(J). The explicit formula for the C 0 test appears in the proof of Corollary 14. Since the C 1 test is based on the function g whose domain is 2n-dimensional and the square J × J has midpoint (m, m), side length w, the C 1 test simplifies, in terms of multi-index notation, to

|α|+|β|≥1 n i=1 ∂ α+e i f (m)∂ β+e i f (m) ∇ f (m) 2 α!β! w 2 |α|+|β| < 1, (C 1 )
where e i is the i-th standard basis vector. We additionally note that since f and J are real, ∇ f (m) 2 = g(m, m). Additionally, for future reference, we collect and adapt the statement of Corollary 16 to the case of g on the region J × J in the following corollary:

Corollary 19. Let f ∈ R[x 1 , . . . , x n ] and define g ∈ [x 1 , . . . , x n , y 1 , . . . , y n ] as g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) . Let J ⊆ R n , and suppose that there is a point (a, b) ∈ J × J such that w ≤ 2 ln 1 + 2 2-4n dist C ((a, b), g) 2 2n+1 (deg( f ) -1) + √ 2n ln 1 + 2 2-4n
.

Then, C 1 (J) is true.

We end this section by collecting a corollary of Corollary 19 which resembles a diameterdistance test and will be used in the next section:

Corollary 20. Let f ∈ R[x 1 , . . . ,
x n ] and define g ∈ [x 1 , . . . , x n , y 1 , . . . , y n ] as g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) . Let J ⊆ R n , and suppose that there is a point x ∈ J such that

Diam(J) ≤ 2 √ n ln 1 + 2 2-4n dist C ((x, x), g) 2 2n+1 (deg( f ) -1) + √ 2n ln 1 + 2 2-4n
.

Then, C 1 (J) is true.

Worst-Case Bounds

In this section, we provide worst-case complexity bounds for the modified PV algorithm. We bound both the number of regions produced by the subdivision as well as the overall bitcomplexity of the algorithm. In the next section, we give examples which show that these bounds are tight in the worst case.

Non-adaptive Bounds

In this section, we use Proposition 3 to bound the number of regions produced by the PV algorithm. We assume9 that f ∈ Z[x 1 , . . . , x n ] and the fixed initial input region I has corners in Z n . Suppose that we can find a δ so that

0 < δ ≤ min x∈I max {dist C (x, f ), dist C ((x, x), g)} ,
and define

K = min          2 √ n ln 1 + 2 2-2n 2 n deg( f ) + √ n ln 1 + 2 2-2n , 2 √ n ln 1 + 2 2-4n 2 2n+1 (deg( f ) -1) + √ 2n ln 1 + 2 2-4n         
.

We observe that the terms in the definition of K are the coefficients in Corollaries 16 and 20. With a slight modification to Proposition 3, we can substitute K and δ from above into Equation (1) to get the following corollary:

Corollary 21. Let f ∈ R[x 1 , . . . , x n ].
Then, the PV algorithm performs at most max 1, 2 Diam(I) Kδ n subdivisions.

We spend the remainder of this section computing a lower bound for δ. We begin by observing that f is a polynomial in n variables and g is a polynomial in 2n variables. In other words, the varieties V C ( f ) and V C (g) are embedded in two different spaces. It becomes easier to study and compare the varieties if they are subsets of the same space; therefore, we consider the image of V C ( f ) in the diagonal of a 2n-dimensional space. In particular, let the variables of C 2n be {x 1 , . . . , x n , y 1 , . . . , y n }. The diagonal ∆ consists of all points of the form x i = y i ; then, ∆ is ndimensional and we can identify C n with ∆. In our case, we write f ∆ for the polynomial system f (x 1 , . . . , x n ) and x iy i for all i. We note that for

x ∈ C n , dist C ((x, x), f ∆ ) = √ 2 dist C (x, f
). Therefore, we are interested in computing a lower bound for min

x∈I max 1 √ 2 dist C ((x, x), f ∆ ), dist C ((x, x), g) ≥ 1 √ 2 min x∈I max dist C ((x, x), f ∆ ), dist C ((x, x), g) . ( 11
)
We now focus on computing a lower bound on the RHS of Inequality ( 11). First we introduce some notation. Let I ∆ be the image of I in ∆, i.e.,

I ∆ = {(x, x) ∈ ∆ : x ∈ I}. Moreover, let C ε = ([-ε, ε] × [-iε, iε]
) 2n be the cube of side length 2ε centered at the origin in C 2n . Then, we write I ∆ ε = I ∆ ⊕ C ε , where ⊕ denotes the Minkowski sum. We observe that for all (x, y) ∈ I ∆ ε , the distance from (x, y) to I ∆ is at most 2 √ nε since that is the largest distance from a point in C ε to the origin. Similarly, if (x, y) ∈ C 2n is not in I ∆ ε , then the distance from (x, y) to I ∆ is more than ε since C ε contains the closed ball of radius ε centered at the origin.

Suppose that we can find a positive integer k so that for any (x, x) ∈ V C ( f ∆ , g), the distance between (x, x) and I ∆ is at least

√ n 2 k-1
. Then, we may use a bound of [START_REF] Jeronimo | On the minimum of a polynomial function on a basic closed semialgebraic set and applications[END_REF] to find a lower bound for the RHS of Inequality (11) as follows:

Proposition 22. Let f ∈ Z[x 1 , . . . , x n ] be of degree d, and define g ∈ R[x 1 , . . . , x n , y 1 , . . . , y n ] as g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) . Suppose that I ⊆ R n is an axisaligned n-dimensional cube whose corners have integral coordinates. Let H be the maximum absolute value of the coefficients of f and coordinates of the corners of I. Suppose that f ∆ = { f, x iy i } is the polynomial system corresponding to the image of V C ( f ) in the diagonal of C 2n and I ∆ = {(x, x) : x ∈ I} is the image of I in the diagonal of C 2n . Let k be a positive integer so that for any (x, x) ∈ V C ( f ∆ , g), the distance between (x, x) and I ∆ is more than

√ n 2 k-1 . Then, min x∈I max dist C ((x, x), f ∆ ), dist C ((x, x), g) ≥ 1 2 k+1 2 4-4n max 2 (2d-2)(k+1) nd 2 H 2 , 60n + 8 (2d -2) 8n -4n2 8n (2d-2) 8n . (12)
Proof. If d = 1, then g is a nonzero constant, so the bound holds trivially. Therefore, we assume that d ≥ 2. We observe that g has degree 2d -2, and, since the maximum absolute value of the coefficients of ∇ f is dH, the maximum absolute value of the coefficients of g is nd 2 H 2 . Let ε = 1 2 k ; by the assumption on k, it follows that V C ( f ∆ , g) ∩ I ∆ ε is empty. We proceed by applying a homothety centered at the origin by a factor of 2 k in C 2n . Therefore, if we let I ∆ be the image of I ∆ after applying the homothety, then applying the homothety to I ∆ ε results in I ∆ 1 . Let f and g be the images of f and g under the homothety and a suitable scaling to restore integer coefficients. Then, the maximum absolute value of the coefficients of f is 2 dk H and the maximum absolute value of the corners of I ∆ is 2 k H. Additionally, the linear terms of the form x iy i are unchanged and the maximum absolute value of the coefficients of g is 2 (2d-2)k nd 2 H 2 .

Next, we identify C 2n with R 4n by decomposing each complex variable into two real variables. This doubles the number of polynomials and scales the maximum absolute value of the coefficients by binomial coefficients, which can be trivially bounded by 2 d for the polynomials coming from f and 2 2d-2 for the polynomials coming from g. Hence, the maximum absolute value of the coefficients coming from f is at most 2 d(k+1) H, and the maximum absolute value of the coefficients coming from g is at most 2 (2d-2)(k+1) nd 2 H 2 .

We observe that if I = [a i , b i ], then I ∆ 1 can be defined by the inequalities:

2 k a i -1 ≤ (x i ), (y i ) ≤ 2 k b i + 1 -1 ≤ (x i ), (y i ) ≤ 1 -2 ≤ (x i ) -(y i ) ≤ 2.
This system accounts for 10n inequalities with largest absolute value of the coefficients at most 2 k H + 1. Moreover, f ∆ corresponds to 2n + 2 equalities while g corresponds to 2 equalities. By applying Jeronimo et al. (2013, Theorem 1.2), we get that the distance between V C ( f ∆ ) and

V C ( g) within I ∆ 1 is at least 2 4-4n max 2 (2d-2)(k+1) nd 2 H 2 , 60n + 8 (2d -2) 8n -4n2 8n (2d-2) 8n .
By scaling this by 1 2 k to remove the homothety and appealing to the triangle inequality, we get the desired result.

In the remainder of this section, we find an upper bound for k. We find this bound by computing a separation bound between I ∆ and V C ( f ∆ , g).

Adaptive Bounds

In this section, we use continuous amortization to adaptively compute the number of boxes created by the PV algorithm. We follow the formulation of continuous amortization in Theorem 5. While Corollary 17 shows that the C 0 test can be substituted directly into the integral of Proposition 7, the C 1 test is slightly more challenging to use, even with Corollary 20 in hand, since it involves both n-dimensional and 2n-dimensional spaces. We, therefore, return to the original formulation of continuous amortization in Theorem 5. We observe that Corollary 16 can be reformulated into a local size bound since the volume of an n-dimensional cube is the width of the cube to the n th power, namely,

G 0 (x) =         2 ln 1 + 2 2-2n dist C (x, f ) 2 n deg( f ) + √ n ln 1 + 2 2-2n        
n is a local size bound for the C 0 test. Similarly, Corollary 19 can be reformulated into a local size bound as follows:

G 1 (x) =         2 ln 1 + 2 2-4n dist C ((x, x), g) 2 2n+1 (deg( f ) -1) + √ 2n ln 1 + 2 2-4n         n .
In this case, even though the test in Corollary 19 uses points (a, b) ∈ J × J, since the statement is existential, the upper bound only gets smaller when restricted to the points in J × J on the diagonal ∆. Applying these local size bounds to Theorem 5 gives the following result:

Theorem 25. Let f ∈ R[x 1 , . . . ,
x n ], and define g ∈ R[x 1 , . . . , x n , y 1 , . . . , y n ] as g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) . Suppose that I ⊆ R n is an axis-aligned n-dimensional cube. The number of (terminal) regions after the subdivision performed by the PV algorithm (before balancing) is bounded above by the maximum of 1 and

2 n I min                  2 n deg( f ) + √ n ln 1 + 2 2-2n 2 ln 1 + 2 2-2n dist C (x, f )         n ,         2 2n+1 (deg( f ) -1) + √ 2n ln 1 + 2 2-4n 2 ln 1 + 2 2-4n dist C ((x, x), g)         n          dV n
where dV n is the n-dimensional volume form. Moreover, the algorithm does not terminate if and only if the integral diverges.

Proof. This is a straight-forward application of continuous amortization from Theorem 5 with ε 1 = 2 -n . The only statement left to prove is that if the integral diverges, then the algorithm does not terminate. The integral diverges if and only if there exists a point x ∈ I so that dist C (x, f ) = 0 and dist C ((x, x), g) = 0. This, however, only happens when f has a real singularity, and regions containing real singularities never pass either of the C 0 or C 1 tests.

This integral provides a more adaptive and accurate estimate on the complexity than the worst-case a priori bounds based on the size of the input because it does not assume that the worst case occurs at every point (or even at any point). Moreover, this integral can be evaluated even when the input polynomial has complex (but not real) singularities. Additionally, this integral applies even when f does not have integral coefficients.

Overall Bit-complexity Bound

In this section, we extend the results of Theorems 24 and 25 to bound the bit-complexity of the PV algorithm using both adaptive and non-adaptive approaches. We begin by bounding the cost for evaluating each of the tests C 0 and C 1 on an arbitrary n-dimensional cube. In this section, we use O(•) and O B (•) to denote the arithmetic complexity and bit-complexity, respectively. The soft-O notation, O(•) and O B (•), mean that we are ignoring logarithmic factors of the dominant term.

A closer look at the predicates C 0 and C 1 and the centered form (see Section 3) reveals that each step of the PV algorithm consists of a multivariate Taylor shift. In particular, given a polynomial F ∈ Z[x 1 , . . . , x n ] and dyadic rational numbers a 1 , . . . , a n , we recursively compute the coefficients of F(x 1 + a 1 , . . . , x n + a n ), cf [START_REF] Mantzaflaris | On continued fraction expansion of real roots of polynomial systems, complexity and condition numbers[END_REF].

Lemma 26. Consider a polynomial F ∈ Z[x 1 , • • • , x n ]
of total degree d and whose coefficients have maximum bit-size τ, and integers a 1 , . . . , a n of bit-size at most . The Taylor shift F(x 1 + a 1 , . . . ,

x n + a n ) costs O B (d n+1 + d n τ).
Proof. We begin the proof with two observations: The maximum degree of any polynomial appearing in this proof is d and the logarithm of the bit-size of the coefficients is O(d + τ), see, e.g., von zur Gathen and Gerhard (1997, Lemma 2.1). We prove this lemma by induction; when n = 1, this is a univariate Taylor shift, whose complexity is O B (d 2 + dτ) by von zur Gathen and Gerhard (1997, Theorem 2.4).

For the inductive step, we assume that d+1 is a power of 2. We begin by calculating (x n +a n ) 2 i for i = 0, . . . , lg d. Since each of these polynomials has coefficients of maximum bit-size O B (d ), and these can be computed through successive squaring, the total cost is O B (d 2 ). We now write

F(x 1 + a 1 , . . . , x n + a n ) = F 0 (x 1 + a 1 , • • • , x n + a n ) + (x n + a n ) d/2 F 1 (x 1 + a 1 , • • • , x n + a n )
where in each F i , the degree in x n is at most d/2. The cost to compute the product (

x n + a n ) d/2 F 1 (x 1 + a 1 , • • • , x n + a n ) is O B (d n (d + τ))
. By continuing this computation recursively, we see that the number of polynomials doubles each time and the maximum degree of x n halves each time, so the total cost of multiplication remains O B (d n (d + τ)) at every step. The recursion has depth lg(d + 1), and the final step of the recursion requires (d + 1) Taylor shifts on (n -1) variables. The result then follows from the inductive hypothesis.

Using Theorem 24 and Lemma 26, we can calculate the overall bit-complexity of the PV algorithm.

Theorem 27. Let f ∈ Z[x 1 , . . . , x n ] be smooth and of degree d and I ⊆ R n be an axis-aligned n-dimensional cube whose corners have integral coordinates. Let τ = lg H be the maximum bitsize of the coefficients of f and the corners of I. The overall bit-complexity of the PV algorithm is

2 O(n 3 2 24n d 12n+1 (d+τ+n lg d)) O B (2 26n d 14n+2 (d + τ)).
Proof. We observe that, after each subdivision in the PV algorithm, the bit-size of the center of the Taylor shift increases by at most 1. To simplify the calculation, we charge each n-dimensional cube in the final partition for all intermediate n-dimensional cubes that contain it, proportionally to their relative areas. Following the approach of Burr (2016, Section 7.1), it follows that the total complexity cost of the PV algorithm is at most twice the cost incurred by the terminal regions themselves.

We observe that the maximum bit-size of a Taylor shift is O(-lg δ) from Theorem 24, so we replace in Lemma 26 by the bound from this theorem. We also recall, from Proposition 22, that g is a polynomial of degree 2d -2 in 2n variables whose coefficients have maximum bit-size O(τ + d + lg n). By substituting these values into Lemma 26 and multiplying by the maximum number of regions, we arrive at the overall bit-complexity of

2 O(n 3 2 24n d 12n+1 (d+τ+n lg d)) O B ((2d -2) 2n+1 (n 2 2 24n d 12n+1 (d + τ + n lg d)) + (2d -2) 2n (τ + d + lg n)),
which simplifies to the desired expression.

We observe that in the 2-dimensional case that frequently occurs in applications, the overall bit-complexity of the PV algorithm is as follows:

Corollary 28. The bit-complexity of the PV algorithm for curves is

2 O(d 25 (d+τ+lg d)) O B (d 30 (d + τ)).
We may also use Theorem 5 to find an adaptive bound for the bit complexity. To be able to use this Theorem, we need to define the appropriate functions h 0 and h 1 that compute the charges to the terminal regions depending on the C 0 and C 1 tests. The main complexity costs in the C 0 and C 1 tests are the costs for the Taylor shifts. Therefore, we use Lemma 26 to derive appropriate cost functions. We observe that for an n-dimensional cube J, the bit-size of the appropriate Taylor shift is at most (lg w(I)lg w(J)). By the discussion above, since the complexity cost of the PV algorithm is at most twice the complexity cost of the terminal regions, we may focus on terminal regions.

If J passes C 0 , since f is a degree d polynomial in n variables whose coefficients have maximum bit-size τ, it follows that the charge associated to J is O B (d n+1 lg w(I)d n+1 lg w(J) + d n τ). Since the functions in Theorem 5 are based on the measure of J and not its width, we define the function

h 0 (y) = d n+1 lg w(I) - d n+1 n lg y + d n τ k 0 (d, τ, n)
where k 0 (d, τ, n) is the maximum value over I of the suppressed terms in the O B . We observe that h 0 (µ(J)) is an upper bound on the bit-cost to compute the Taylor shift for the C 0 test for J.

On the other hand, if J passes C 1 , since g is a degree 2d -2 polynomial in n variables whose coefficients have maximum bit-size O(τ + d + lg n), it follows that the charge associated to J for the C 1 test is O B (2 2n d 2n+1 lg w(I) -2 n d 2n+1 lg w(J) + 2 2n d 2n (τ + d + lg n)), which simplifies to O B (2 2n d 2n+1 lg w(I) -2 2n d 2n+1 lg w(J) + 2 2n d 2n τ). As above, we define the function

h 1 (y) = 2 2n d 2n+1 lg w(I) - 2 2n d 2n+1 n lg y + 2 2n d 2n τ k 1 (d, τ, n)
where k 1 (d, τ, n) is the maximum value over I of the suppressed terms in the O B . We observe that h 1 (µ(J)) is an upper bound on the bit-cost to compute the Taylor shift for the C 1 test for J.

We use these two functions along with G 0 and G 1 as defined in Section 5.2 to develop adaptive bounds on the bit-complexity of the PV algorithm as follows:

Theorem 29. Let f ∈ Z[x 1 , . . . , x n ] be smooth and of degree d and I ⊆ R n be an axis-aligned n-dimensional cube whose corners have integral coordinates. Let τ = lg H be the maximum bitsize of the coefficients of f and the corners of I. The overall bit-complexity of the PV algorithm is the maximum of h 0 (w(I) n ), h 1 (w(I) n ), and

2 n I min h 0 (2 -n G 0 (x)) G 0 (x) , h 1 (2 -n G 1 (x)) G 1 (x)
dV n .

Examples

The bounds in Theorems 24 and 25 are both exponential with respect to the degree of the polynomial and the number of variables. They remain exponential even if we assume that the number of variables is constant. In [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], the authors show that for several examples the computation time is efficient in practice. The following examples illustrate that:

• The exponential behavior is optimal, up to constants in the exponents and • In particular cases, the complexity is provably better than the worst-case. ) where w = w(I) is the width of he initial region. (b) The output of the PV algorithm on 1000x 4 y 4 -1. We observe that the near-singularity at (0, 0) does not cause exponentially many subdivisions. Instead, the pair of curves with the same asymptote contribute to this behavior since the width of boxes along the horizontal (vertical) axis must be less than the vertical (horizontal) distance between the two branches.

Lemma 30. The bound of Theorem 24 is asymptotically tight.

Proof. Following the construction in [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF], consider the Mignotte polynomial P(x) = x d -2(ax -1) 2 and the related polynomial P 2 (x) = x d -(ax -1) 2 where a is a sufficiently large positive integer. The product P(x)P 2 (x) is of degree 2d and the largest coefficient is of size 2a 4 . In [START_REF] Eigenwillig | Almost tight recursion tree bounds for the descartes method[END_REF], it is shown that the product P(x)P 2 (x) has (at least) three roots in the interval (a -1h, a -1 + h) where h = a -d/2-1 , see Figure 3(a). Treating P(x)P 2 (x) as a polynomial in n variables, we see that the PV algorithm to approximate the variety in an n-dimensional cube I of side length w = w(I) requires subdividing until the regions have side length at most 2h to separate the three vertical hyperplanes in the interval (a -1h, a -1 + h). Since this occurs along an entire hyperplane of the input region, the number of small boxes is, at least,

w n-1 2h = 1 2 w(I) n-1 a d/2+1
, which is exponential in both the size of the input region and the size of the coefficients of the polynomial.

The previous example, while illustrating that the bounds are tight, raises the question of whether exponential behavior is due to the fact that the example is a one-dimensional problem lifted to higher dimensions. We now show provide an example that shows that this exponential behavior can be observed for a curve involving both x and y in two dimensions. In particular, in Lemma 30, the exponential behavior in two-dimensions was caused by two curves which were close together, but had a curve of critical points between them. We can mimic that behavior for a curve in two-dimensions by considering a situation where two local components of the curve share an asymptote.

Example 31. Fix ε > 0 and consider f (x 1 , x 2 ) of one of the following forms:

• f (x 1 , x 2 ) = x a 1 1 x a 2 2 -ε a 1 +a 2 where a 1 and a 2 are both positive integers and at least one is even or

• f (x 1 , x 2 ) = x a 1
1 x a 2 2 + ε a 1 +a 2 where a 1 and a 2 are both positive integers and exactly one is even.

In either of these cases, the PV algorithm produces exponentially many regions in the size of the input box and the size of the coefficients of the polynomial, see Figure 3(b).

Since all of the cases are similar, we focus on the case where f (x 1 , x 2 ) = x a 1 1 x a 2 2 -ε a 1 +a 2 and a 2 is even. In this case, we show that the number of regions which intersect the positive x-axis is exponential in the size of the input. Since ∇ f is zero on the positive x-axis, any box which is terminal and intersects the positive x-axis must satisfy Condition C 0 . For any positive x,

       x, ± ε a 1 +a 2 x a 1 1/a 2       
are points on the variety V R ( f ). Therefore, any region which is terminal and contains (x, 0) must have width at most 2 ε a 1 +a 2 x a 1 1/a 2 (13) since, otherwise, the region would contain a point of V R ( f ) and could not satisfy Condition C 0 . Let J be a terminal region which intersects the positive x-axis and let [s 1 , s 2 ] be the intersection of J with the positive x-axis. Then, consider the integral

1 2 s 2 s 1 x a 1 ε a 1 +a 2 1/a 2 dx ≤ w(J) 2 s a 1 2 ε a 1 +a 2 1/a 2 , ( 14
)
where the inequality follows since the integrand is increasing. Since (s 2 , 0) ∈ J, by the bound in Expression (13), it follows that Expression ( 14) is at most 1.

Suppose that the intersection of the initial region I with the positive x-axis is [r 1 , r 2 ]. Then, by the bound on Integral ( 14) from above, it follows that 1 2

r 2 r 1 x a 1 ε a 1 +a 2 1/a 2 dx = a 2 2(a 1 + a 2 )       r 2 ε a 1 +a 2 a 2 - r 1 ε a 1 +a 2 a 2      
is a lower bound on the number of regions formed by the PV algorithm along the positive x-axis. This region count is exponential in both the size of the input region and the size of the coefficients of the polynomial.

We remark that the example above is intrinsically hard for the algorithm and it can be adapted to higher dimensions and applies even under a change of coordinates. We also note that the exponential behavior does not come from the near singularity at (0, 0), but from the curves sharing asymptotes. For the centered form, see Section 3, the complex portions of the curve also affect subdivisions, so, when using the centered form for the tests C 0 and C 1 , the exponential behavior from the analysis above can be extended for all positive integers a 1 and a 2 such that a 1 + a 2 > 2. Even though our bounds are optimal, in practice, these are often quite pessimistic, as the actual separation bounds do not follow the worst case behavior. We illustrate this better behavior in the following two examples:

Example 32. Fix ε > 0 and consider f (x 1 , x 2 ) = x 2 1 + x 2 2 + ε 2 . Then,

dist C ((x 1 , x 2 ), f ) = x 2 1 + x 2 2 2 + ε 2
and dist C ((x 1 , x 2 , x 1 , x 2 ), g) = x 2 1 + x 2 2 .

Let I be the initial input square where w = w(I) is the width of I. By substituting these bounds into Theorem 25, we find that the number of regions constructed by the PV algorithm is O(lg(w)lg(ε)), see Figure 4(a).

Example 33. Fix ε > 0 and consider f (x 1 , x 2 ) = x 2 1 + x 2 2 -ε 2 . Then,

dist C ((x 1 , x 2 ), f ) =            x 2 1 + x 2 2 -ε x 2 1 + x 2 2 ≤ 4ε 2 x 2 1 +x 2 2 2 -ε 2 x 2 1 + x 2 2 > 4ε 2
and dist C ((x 1 , x 2 , x 1 , x 2 ), g) = x 2 1 + x 2 2 . Let I be the initial input square where w = w(I) is the width of I. By substituting these bounds into Theorem 25, we find that the number of regions constructed by the PV algorithm is O(lg(w)lg(ε)), see Figure 4(b).

Moreover, we observe that for each of these examples, the minimum distance between V C ( f ∆ ) and V C (g) is at most ε. Therefore, a bound coming from Theorem 24 would be much larger than the bound continuous amortization provides.

It remains an open question to deduce adaptive complexity bounds for the PV algorithms from Theorem 25 based on geometric and a priori parameters. We observe that since the complexity of the algorithm can be exponential in the inputs, the integral must be described in terms of more parameters than the degree of f and the size of the coefficients of f .

Figure 3 :

 3 Figure3: (a) The output of the PV algorithm for f (x, y) = (x d -2(ax -1) 2 )(x d -(ax -1) 2 ). The solutions to f (x, y) = 0 are close vertical lines (the illustrated case is when d = 3 and a = 3). The width of boxes between vertical lines is at most2a -d2 -1 and they extend the entire length of the initial region. The number of regions is bounded from above by Ω(wad 2 +1) where w = w(I) is the width of he initial region. (b) The output of the PV algorithm on 1000x 4 y 4 -1. We observe that the near-singularity at (0, 0) does not cause exponentially many subdivisions. Instead, the pair of curves with the same asymptote contribute to this behavior since the width of boxes along the horizontal (vertical) axis must be less than the vertical (horizontal) distance between the two branches.

Figure 4 :

 4 Figure 4: (a) The output of the PV algorithm for f (x, y) = x 2 + y 2 -ε 2 . The number of regions is bounded by O(lg(w)lg(ε)). (b) The output of the PV algorithm for f (x, y) = x 2 +y 2 +ε 2 . The number of regions is bounded by O(lg(w)-lg(ε)).

Our approach applies to similar subdivision-based methods for approximating curves such as[START_REF] Yap | Adaptive isotopic approximation of nonsingular curves: the parameterizability and nonlocal isotopy approach[END_REF]. The final complexity results are similar, and we leave the details to the interested reader.

The correctness depends on the curve being bounded, but the termination of the algorithm depends only on the smoothness. SeeBurr et al. (2012) for an extension of this algorithm which includes a correctness statements for unbounded curves.

The theory continues to apply even with weaker conditions, but there are a few technicalities that arise. For most applications, this assumption is not an additional constraint. For weaker conditions, we leave the details to the interested reader.

We note that the expression for cos(nx) is a multiple of of the Chebyshev polynomials of the first kind.

This analysis can be extended to the case of a region which is not an n-dimensional cube, by considering the smallest n-dimensional cube containing the n-dimensional box. We leave the details to the interested reader.

w. By substituting this into Expression (10), the desired result follows.By combining Corollary 14 with Lemma 15, we explicitly show that the predicate C is a distance-diameter test:

The argument in this section can be directly generalized for f ∈ Q[x 1 , . . . , x n ] and I whose corners are in Q n . We leave the details to the interested reader.

Acknowledgments.

The authors are grateful to Martin Sombra for his comments on the previous versions of this paper. The authors are also grateful to Hoon Hong for his observations on how to generalize the initial results. The authors thank Bernard Mourrain for his discussions about the difficult cases for the algorithm.

1 Partially supported by a grant from the Simons Foundation (#282399 to Michael Burr) and National Science Foundation Grant CCF-1527193.

2 Partially supported by the National Science Foundation under Grants CCF-1407623, DMS-1403062 and DMS-1547399.

3 Partially supported by the ANR JCJC "GALOP" (ANR-17-CE40-0009) and the PGMO grant ALMA.

Proposition 23. Let f ∈ Z[x 1 , . . . , x n ] be smooth and of degree d, and define g ∈ R[x 1 , . . . , x n , y 1 , . . . , y n ] as g(x 1 , . . . , x n , y 1 , . . . , y n ) = ∇ f (x 1 , . . . , x n ), ∇ f (y 1 , . . . , y n ) . Suppose that I ⊆ R n is an axisaligned n-dimensional cube whose corners have integral coordinates. Let H be the maximum absolute value of the coefficients of f and coordinates of the corners of I. Suppose that f ∆ = { f, x iy i } is the polynomial system corresponding to the image of V C ( f ) in the diagonal of C 2n and I ∆ = {(x, x) : x ∈ I} is the image of I in the diagonal of C 2n . Let (x, x) ∈ V C ( f ∆ , g), then the distance between (x, x) and I ∆ is at least

.

Proof. If d = 1, then g is a nonzero constant, so the bound holds vacuously. Therefore, we assume that d ≥ 2. Throughout this proof, we restrict our attention to I ∆ 1 and we observe that if

1 is bounded away from the real points in the diagonal, i.e., R 2n ∩ ∆. Moreover, since V C ( f ∆ , g) contains no real points, it follows that 2 (x i ) 2 is bounded away from zero for all (x, x)

1 and the sum is a lower bound on the square of the distance to I 1 ∆ . We now proceed to find a lower bound on this sum.

As in Proposition 22, we identify C 2n with R 4n . Since no homotheties are required, I ∆ 1 corresponds to 10n inequalities with maximum coefficient size H + 1, f ∆ corresponds to 2n + 2 equalities with coefficient size at most 2 d H, and g corresponds to 2 equalities with coefficient size at most 2 (2d-2) nd 2 H 2 . Finally, the sum of interest is of degree 2 with maximum coefficient size of 2. By applying Jeronimo et al. (2013, Theorem 1.1), we get that on

Since the sum 2 (x i ) 2 is the square of the distance from (x, x) to R 2n , which contains I ∆ 1 , so, by taking the square root, the result follows.

By combining Corollary 21 with Propositions 22 and 23, we obtain an explicit bound for the number of terminal regions produced by the PV algorithm.

Theorem 24. Let f ∈ Z[x 1 , . . . , x n ] be smooth and of degree d and I ⊆ R n be an axis-aligned n-dimensional cube whose corners have integral coordinates. Let H be the maximum absolute value of the coefficients of f and coordinates of the corners of I. The number of (terminal) regions produced by the PV algorithm is 2 O(n 3 2 24n d 12n+1 (d+lg H+n lg d)) .

Proof. We observe that for any positive constant a, lg(ln

Next, since k is an integer and an the exponent of 2, k can be chosen to be within 1 of the base 2 logarithm of the bound in Proposition 23. Therefore, k = O(n2 8n d 4n (d + lg H + n lg d)). Substituting this into the bound for δ, we find thatlg δ = O(n 2 2 24n d 12n+1 (d + lg H + n lg d)). Substituting the bounds into the expression in Corollary 21 results in the stated complexity.