Gilles Barthe
email: gilles.barthe@imdea.org

César Kunz
email: cesar.kunz@inria.fr

David Pichardie
email: david.pichardie@irisa.fr

Julián Samborski-Forlese
email: julian.sf@imdea.org

Preservation of proof obligations for hybrid verification methods

Program verification environments increasingly rely on hybrid methods that combine static analyses and verification condition generation. While such verification environments operate on source programs, it is often preferable to achieve guarantees about executable code. We show that, for a hybrid verification method based on numerical static analysis and verification condition generation, compilation preserves proof obligations and therefore it is possible to transfer evidence from source to compiled programs. Our result relies on the preservation of the solutions of analysis by compilation; this is achieved by relying on a bytecode analysis that performs symbolic execution of stack expressions in order to overcome the loss of precision incurred by performing static analyses on compiled (rather than source) code. Finally, we show that hybrid verification methods are sound by proving that every program provable by hybrid methods is also provable (at a higher cost) by standard methods.

Introduction

Program verification techniques are widely used to reason about the correctness of applications and play an important role in fields such as mobile code and embedded systems where strong guarantees are required. However, program verification, and in particular deductive program verification, is traditionally applied to source code, whereas guarantees are often required about executable code. This discrepancy is particularly acute in the context of mobile code, where code consumers do not trust code producers and may not have access to the source code. Therefore, it is of interest to develop methods to transfer evidence from This work is partially supported by the EU project MOBIUS. source code verification to the code consumers. In earlier work [START_REF] Barthe | Preservation of proof obligations for Java[END_REF], we focussed on transferring proofs from source to compiled programs in the context of verification methods based on verification condition generators, that are commonly used in Proof Carrying Code [START_REF] Necula | Proof-carrying code[END_REF] and program verification environments. We have shown that non-optimizing compilation preserves proof obligations, and therefore that the certificates of source programs can be reused directly to validate compiled programs. However, state-of-the-art verification tools do not use plain verification condition generation; instead, they rely on hybrid methods, that combine static analyses and verification condition generation.

The objective of the paper is to extend preservation of proof obligations to hybrid verification methods. For concreteness, we consider a small imperative language with arrays, and we focus on a hybrid method based on a generic numerical analysis [START_REF] Miné | Weakly Relational Numerical Abstract Domains[END_REF][START_REF] Besson | Result certification for relational program analysis[END_REF] and that can be instantiated to several numeric domains, including polyhedra.

We first define a hybrid verification method in which programs are subjected to static analysis, and then to verification condition generation. The VCgen exploits the information of the analysis in two useful ways: on the one hand, verification conditions that originate from spurious edges in the control-flow graph are discarded: more precisely, the VCgen ignores the case of out-of-bound accesses whenever the analysis ensures that accesses are within bounds. This leads to fewer, smaller verification conditions. Furthermore, the VCgen adds the results of the analysis as additional assumptions to help the user prove the verification conditions. This is particularly useful for the relational analyses considered as they can provide part of the invariants required to prove programs correct.

Then, we prove preservation of proof obligations using the techniques of Barthe et al. [START_REF] Barthe | Preservation of proof obligations for Java[END_REF]. The proof relies on knowing that the solutions of the analysis are preserved by compilation. Although analyzing compiled programs is known to be less precise than analyzing source programs, as stated by Logozzo and Fähndrich [START_REF] Logozzo | On the relative completeness of bytecode analysis versus source code analysis[END_REF], we achieve preservation of solutions by defining at bytecode level an analysis that performs a symbolic execution of stacks [START_REF] Xi | Towards array bound check elimination in java tm virtual machine language[END_REF][START_REF] Wildmoser | Prototyping proof carrying code[END_REF][START_REF] Besson | Result certification for relational program analysis[END_REF].

Finally, we relate hybrid verification to standard verification. We show that programs that are provably correct using our hybrid method, remain provably correct using standard verification condition generation; to this end, we define a compiler that transforms a hybrid specification (combing logical assertions and analysis results) into a logical one by giving a logical interpretation of the analysis results.

Proof Carrying Code and Hybrid Methods

The general goal of our work [START_REF] Barthe | MO-BIUS: Mobility, ubiquity, security. objectives and progress report[END_REF] is to develop efficient, scalable and trustworthy security architectures for mobile code. We adopt the view of Proof Carrying Code (PCC) [START_REF] Necula | Proof-carrying code[END_REF], which emphasizes security through verifiable evidence, and require that mobile code components come equipped with a certificate that can be checked efficiently and independently by the code consumer to ensure that the downloaded components issued by the producer respect its policy.

There are two main flavours of PCC, that respectively rely on type systems and program verification. We briefly outline the issues of certificate generation and certificate checking for both of them.

In the type-based approach, certificates are generated automatically by an analyzer and packaged with the code; on the consumer side, the certificate is verified efficiently by a checker that is tighly coupled with the analyzer (in the spirit of lightweight bytecode verification). In the logic-based approach, certificates may be generated automatically by certifying compilers, provided the policies are sufficiently simple. However, certificates must be constructed interactively once policies become more complex; preservation of proof obligations allows the construction of the certificate to arise at source level. The certificate (that may include some program annotations) is packaged with the program. On the consumer side, the protocol proceeds in three steps: first, annotations are checked against the policy; second, verification condition generation is used to generate proof obligations; third, certificate checking is performed using a proof checker, that verifies whether the certificates prove the proof obligations.

Hybrid methods aim to provide a tight integration of type systems and logical methods. There are two approaches to hybrid program verification. In the explicit approach, the user provides safety annotations that are used by the verification condition generator, and checked by a annotation checker. In contrast, the implicit approach advocates that the annotations are inferred by a static analyzer, and then used by verification condition generation. Both approaches are used (sometimes in conjunction) in deductive program verification, as well as in type-based analyses.

For the hybrid method developed in this paper, program verification proceeds as follows: first, the annotated program is subjected to the static analysis, and a solution is computed (note the analysis does not take advantage of the annotations). Then, the verification condition generator uses the solution to compute a smaller set of proof obligations, which must then be discharged interactively. On the consumer side, the program arrives packaged with the solution of the analysis (or as in lightweight bytecode verification a partial solution that contains sufficient information to recompute the solution efficiently), the program annotations, and the certificate, and the checking proceeds in four steps; first, one checks the correctness of the analysis. The remaining three steps are as in logic-based PCC.

Setting

This section introduces the source language (an imperative language with arrays of integers), the target language (a stack-based language with jumps), and the compiler.

We assume given two disjoint sets V s of scalar variables and V a of array variables, and let V denote V s + V a . Each variable in V a has an associated size. Furthermore, we assume given two sets V old s and V old a in 1-1 correspondence with V s and V a , which are used to store initial values. We also consider a special variable res, which is used to represent the value of the program result. Finally, we assume given a set Lab ⊂ N of labels.

Source Language

Programs are defined as commands, and are decorated with labels in order to express analysis results:

e ::= e op e | n | x | a[e] c ::= Skip | [x:=e] k | [a[e]:=e] k | c; c | [return e] k | if [e ⋊ ⋉ e] k then c else c | while [e ⋊ ⋉ e] k do c
where x, a, n and k respectively range over V s , V a , Z and Lab, op ranges over (binary) arithmetic operations, and ⋊ ⋉ over arithmetic comparisons. We assume that labels occur at most once in commands. The semantics of source programs is formalized by a small-step transition relation between states. States may be intermediate, in which case they consist of a statement and of a memory, or final, in which case they consist of a memory, and possibly a tag to denote abnormal termination. Memories are modeled as pairs of mappings respectively from variables to values and from arrays to indices to values. We assume that each array a comes equipped with its size |a| and define the semantic domains of the source language as follows:

V M em = V s → Z AM em = Πa ∈ V a . {i | 1 ≤ i ≤ |a|} → Z M em = V M em × AM em State s I = Stmt × V M em × AM em State s F = V M em × AM em × (Z + {AOB}) State s = State s I + State s F
The operational semantics of programs is standard and, thus, omitted. (See the next subsection for the semantics of instructions that manipulate arrays.)

Bytecode Language

A bytecode program is defined as a list of instructions. Instructions either manipulate the memory that stores the values of variables and the contents of arrays, or manipulate the operand stack, or perform a conditional or unconditional jump. The set of instructions is defined by the grammar

ins ::= prim op | push v | load x | store x | return | aload a | astore a | cjmp ⋊ ⋉ l | jmp l | nop
We denote by ṗ[l] the instruction at position l of a bytecode program ṗ. The semantics of bytecode programs is formalized using a transition relation between states. States may either be intermediate or final; intermediate states consist of a program counter, an operand stack, that stores the results of intermediate computations, and a memory. The semantic domains of the bytecode language are defined as follows, where we implicitly assume that the program counter is within the bounds of programs:

Stack = Z ⋆ State b I = N × V M em × AM em × Stack State b F = V M em × AM em × (Z + {AOB}) State b = State b I + State b F
The operational semantics of programs is standard. We only provide the operational semantics of the instructions aload and astore; these instructions may cause abrupt termination if array accesses are out-of-bound. The rules are given in Figure 2, where we use the notation [f | s → r] to refer the function that is identical to f everywhere except in r that returns s, for any sets R and S and any function f : R → S.

Compiler

The compiler is standard, and defined in Figure 1; we use the function init : Stm → Lab to associate to each statement its initial label. We assume label compatibility, i.e. that the label of a source statement is the same as the label of the program point for its compilation. Throughout the rest of the paper we let P be a source program, and the bytecode program ṗ the result of the compilation of program P .

; s2 = s 1 ; s 2 [return e] k = k : e e; return [Skip] k = k : nop if [e 1 ⋊ ⋉ e 2 ] k then s 1 else s 2 =
k : e 2 e; e 1 e; cjmp ⋊ ⋉ k 1 ; k 2 : s 2 ; jmp l; k 1 : s 1

where

k 1 = init(s 1 ) = k 2 + | s 2 | + 1 k 2 = init(s 2 ) = k + | e 2 e; e 1 e| + 1 l = k 1 + | s 1 | while [e 1 ⋊ ⋉ e 2 ] k do s = k : e 2 e; e 1 e; cjmp ⋊ ⋉ k 1 ; jmp l; k 1 : s ; jmp k where k 1 = k + | e 2 e| + | e 1 e| + 2 l = k 1 + | s | + 1 Figure 1. Compiler P [i] = aload a 0 ≤ n < |a| i, ρv, ρa, n :: s i + 1, ρv, ρa, ρa a n :: s P [i] = aload a ¬ 0 ≤ n < |a| i, ρv, ρa, n :: s EX ρv, ρa, AOB P [i] = astore a 0 ≤ n < |a| i, ρv, ρa, n :: v :: s i+1, ρv, [ρa | a → [ρa a | n → v]], s P [i] = astore a ¬ 0 ≤ n < |a| i,

Preservation of solutions

It is folklore that compilation potentially yields a loss of precision for relational analyses. The purpose of this section is to show that solutions of abstract interpretations are preserved by compilation, provided one uses symbolic expressions, as done in [START_REF] Xi | Towards array bound check elimination in java tm virtual machine language[END_REF][START_REF] Wildmoser | Prototyping proof carrying code[END_REF][START_REF] Besson | Result certification for relational program analysis[END_REF], to mitigate the presence of the operand stack and to recover the loss of precision incurred by compilation.

Symbolic Expressions

Expressions and guards serve as the interface with the numerical relational domain in the analysis for bytecode. Be-

γ(d 1 ⊓ d 2 ) ⊇ γ(d 1 ) ∩ γ(d 2 ) γ(d 1 ⊔ d 2 ) ⊇ γ(d 1 ) ∪ γ(d 2 ) γ( x:=e ♯ (d)) ⊇ { ρv[x → v], ρa | ρ ∈ γ(d) ∧ v ∈ e ρ} γ( x[e 1 ]:=e 2 ♯ (d)) ⊇ { ρv, [ρa | x → [ρa a | v 1 → v 2 ]] | ρ ∈ γ(d) ∧ v 1 ∈ e 1 ρ ∧ v 2 ∈ e 2 ρ} {ρ | t ρ} ⊆ γ(assume ♯ (t))
where ρ = ρv, ρa . Expr ∋ e ::

= n | x | x[e] | ? | ?[e] | e op e x ∈ V Guard ∋ t ::= e ⋊ ⋉ e
The expression ? represents an unknown value; therefore, expressions are interpreted as sets of possible values. Formally, the semantics e ρ and t ρ of expressions with respect to an environment ρ = ρ v , ρ a are defined by the clauses:

n ρ = {n} x ρ = ρ v x ? ρ = Z ?[e] ρ = Z x[e] ρ = {ρ a x v | v ∈ e ρ } e 1 op e 2 ρ = {n 1 op n 2 | n1 ∈ e 1 ρ , n 2 ∈ e 2 ρ } e 1 ⋊ ⋉ e 2 ρ ⇐⇒ ∃ n 1 ∈ e 1 ρ , n 2 ∈ e 2 ρ • n 1 ⋊ ⋉ n 2
Note that the expression ? is not required for analyzing bytecode programs that are achieved by compilation of the source program, since the stack is empty after storing a value in an array. However, it provides more precision when dealing with programs that are not obtained by compilation.

Abstract domain

Following Miné [START_REF] Miné | Weakly Relational Numerical Abstract Domains[END_REF], we assume given an abstract numerical domain interface, which can be instantiated with standard relational abstract domains. The interface consists of a domain D equipped with a partial order ⊑ ⊆ D × D, meet and join operators ⊓, ⊔ : D × D → D, a least element ⊥ and a greater element ⊤. We also assume given abstract assignment functions x:=e ♯ , x[e 1 ]:=e 2 ♯ : D → D, and a function assume ♯ that maps guards to abstract elements.

Finally, we assume given a monotone concretization function γ : D → P(V M em × AM em) mapping abstract elements to sets of environments in V M em × AM em, and satisfying the properties in Figure 3.

We define the abstract test

t ♯ : D → D of a guard t ∈ Guard by t ♯ (l ♯ ) = assume ♯ (t) ⊓ l ♯ .

Source Code Analysis

The source code analysis is specified by abstract transfer functions that map elements of the abstract domain into elements of the abstract domain.

Definition 1 (Abstract Domain for High-Level) A result of the analysis for the source program P is described by a mapping Loc in the lattice State ♯ = Lab → D.

Definition 2 (Solution) A mapping Loc for the source program P is a solution of the analysis if it verifies the constraint system defined in Figure 4, i.e. Loc ⊢ P holds.

Byte Code Analysis

As for the source code analysis, the bytecode analysis is defined by abstract transfer functions that map abstract states into abstract states. In this case, the abstract states are pairs of the form (s ♯ , l ♯ ) where l ♯ is an element of the abstract domain, and the list of symbolic expressions s ♯ abstracts the operand stack. The symbolic abstract domain for stacks is Expr ⋆ , where for any set A, A ⋆ denotes the domain of lists with elements in A. The set of variables considered by the bytecode analysis is the same as in the source code analysis. Definition 3 (Bytecode Abstract Domain) A result of the analysis for ṗ is described by a mapping l ȯc in the lattice

ṡtate ♯ = Lab → (Expr ⋆ L × D)
An analysis result is a solution of the analysis if it satisfies the constraint system associated to each program. The constraint system is defined in Figure 5. For instructions other than branching or return instructions, the constraint is defined by partial transfer functions in Expr ⋆ × D ⇀ (Expr ⋆ × D), most of them defined as a symbolic execution affecting the abstract representation of the operand stack.

Definition 4 (Solution) A mapping l ȯc for the bytecode program ṗ is a solution of the analysis if it satisfies the constraint system of Figure 5, i.e. if l ȯc ⊢ ṗ holds.

Preservation of Solutions

We define first the compilation of a source code analysis solution and then show that it is a solution for the byte code analysis. For notational convenience, we denote by

ḟs1;••• ;sn (s ♯ , l ♯ ) the composition ḟsn (• • • ( ḟs1 (s ♯ , l ♯ )) • • • ),
where s 1 ; • • • ; s n is a sequence of byte code instructions.

Let succ(l) denotes the set of successors of a label l, e.g. succ(l) = ∅ and succ(l

) = {l + 1, l ′ } respectively for ṗ[l] = return and ṗ[l] = cjmp ⋊ ⋉ l ′ . The set pred(l) is defined as {l ′ | l ∈ succ(l ′ )}. stm ∈ {if t then s 1 else s 2 , while t do c, s 1 ; s 2 , return e} Loc ⊢ {Loc(i)} [stm] i {Fstm(Loc(i))} Loc ⊢ {Loc(i)} [return e] i {⊥} Loc ⊢ { t ♯ (Loc(i))} s 1 {l ♯ 1 } Loc ⊢ { ¬t ♯ (Loc(i))} s 2 {l ♯ 2 } Loc ⊢ {Loc(i)} if [t] i then s 1 else s 2 {l ♯ 1 ⊔ l ♯ 2 } Loc ⊢ { t ♯ (Loc(i))} s {l ♯ } l ♯ ⊑ Loc(i) Loc ⊢ {Loc(i)} while [t] i do s { ¬t ♯ (Loc(i))} Loc ⊢ {l ♯ } s 1 {l ♯ 1 } Loc ⊢ {l ♯ 1 } s 2 {l ♯ 2 } Loc ⊢ {l ♯ } s 1 ; s 2 {l ♯ 2 } Loc ⊢ {⊤} P {l ♯ } Loc ⊢ P where Fstm(l ♯ ) = 8 < : l ♯ if stm = Skip x:=e ♯ (l ♯ ) if stm = x:=e a[e 1 ]:=e 2 ♯ (l ♯ ) if stm = a[e 1 ]:=e 2
s ♯ , l ♯ ) push n (s ♯ , l ♯ ) → (n :: s ♯ , l ♯ ) load r (s ♯ , l ♯ ) → ( r :: s ♯ , l ♯ ) store r (e :: s ♯ , l ♯ ) → (s ♯ [?/r], r:=e ♯ (l ♯ )) aload a (e :: s ♯ , l ♯ ) → ( a[e] :: s ♯ , l ♯ )
astore a (e 1 :: e 2 ::

s ♯ , l ♯ ) → (s ♯ [?/a], a[e 1 ]:=e 2 ♯ (l ♯ )) nop (s ♯ , l ♯ ) → (s ♯ , l ♯ ) Instr ∈ { jmp i ′ , cjmp ⋊ ⋉ i ′ , return} ḟinstr (l ȯc(i)) ⊑ l ȯc(i + 1)
l ȯc ⊢ i : Instr l ȯc ⊢ i : return l ȯc(i) = (e 1 :: e 2 ::

s ♯ , l ♯ ) (s ♯ , ¬(e 1 ⋊ ⋉ e 2 ) ♯ (l ♯ )) ⊑ l ȯc(i + 1) (s ♯ , e 1 ⋊ ⋉ e 2 ♯ (l ♯ )) ⊑ l ȯc(j) l ȯc ⊢ i : cjmp ⋊ ⋉ j l ȯc(i) ⊑ l ȯc(j) l ȯc ⊢ i : jmp j ⊤ ⊑ l ȯc(0) ∀i ∈ dom( ṗ) • l ȯc ⊢ i : ṗ[i] l ȯc ⊢ ṗ Figure 5
. Definition of the constraint system for the byte code analysis.

Remark 5 For each byte code program ṗ, we can extract from the previous constraint system a set of transfer functions ( ġi,j ) (i,j)∈Lab 2 such that l ȯc ⊢ ṗ if and only if

k ′ ∈pred(k) ġk ′ ,k (l ȯc(k ′ )) ⊑ l ȯc(k) for all k ∈ dom( ṗ).
We can extend a partial function loc partial ∈ ṡtate ♯ to a total function l ȯc on dom( ṗ) if we set

l ȯc(k) = loc partial (k) if k ∈ dom( loc partial ) k ′ ∈pred(k) ġk ′ ,k (l ȯc(k ′ )) otherwise
This definition only makes sense if, by considering the control flow graph of ṗ whose edges are {(i, j) | i ∈ dom( ṗ) ∧ j ∈ succ(i)}, every loop contain a label in dom( loc partial ). We refer to the function l ȯc as the completion of the partial function loc partial . This definition can be shown to be well defined from the facts that Loc annotates every loop in P and each loop in the control flow graph of ṗ contains a label of a loop in P .

Lemma 7 Let ṗ1 , ṗ2 and e s.t. ṗ = ṗ1 :: l : e e :: l ′ : ṗ2 . Then, l ȯc(l ′ ) = f i1;...;i k (s ♯ , l ♯ ) = (e :: s ♯ , l ♯ ) where (s ♯ , l ♯ ) = l ȯc(l) and [i 1 ; . . . ; i k ] = e e .

PROOF. We prove the lemma by structural induction over expression e.

The following lemma states the main result of this section: compilation preserves analysis solutions.

Lemma 8 If Loc is s.t. Loc ⊢ P , then the analysis result l ȯc compiled from Loc is s.t. l ȯc ⊢ ṗ, i.e. it is a solution of the bytecode analysis. PROOF. By structural induction over s we can prove that if ṗ = ṗ1 :: k 1 : s :: k 2 : ṗ2 and exists l ♯ such that

Loc ⊢ {Loc(k 1 )} s {l ♯ } and l ♯ ⊑ l ȯc(k 2 ) then we have: ∀k ∈ [k 1 , k 2 ), l ȯc ⊢ k : ṗ[k].

Preservation of proof obligations

In this section we define two verification frameworks, respectively for source programs and for unstructured bytecode of previous sections. As a specification language we consider first order formulae, namely the domain of assertions A. The validity of an assertions in a particular execution state η ∈ State s is standard. In particular, an assertion that contains the expression a[e] is invalid in those execution states in which e is out of the bounds of the array a.

We consider as a program specification a tuple (pre, annot, post, χ), where the assertion pre is a precondition, post and post are respectively normal and abnormal postconditions, and the partial function annot : Lab ⇀ A maps program labels to internal points specifications. The special variable res may only occur in post, and pre only refers to variables from V . When specifying a bytecode program, assertions may refer to the special variable s representing the operand stack.

We say that a program satisfies the specification (pre, annot, post, χ), if every execution starting in a state that satisfies pre only reaches normal final states satisfying post or abnormal states satisfying χ, and only reaches intermediate l-labeled points satisfying annot(l). Given a program specification (pre, annot, post, χ), a verification condition generator (VCgen) framework provides a set of sufficient proof obligations that ensures that the program satisfies the specification.

The VCgens defined in this section are hybrid in the sense that they take advantage of a previously computed analysis to reduce the size of proof obligations. We assume that the result of a relational analysis (Loc and l ȯc respectively for source and bytecode programs) is given as input to the VCgen. A further improvement over standard VCgens consists of reusing the result of the analysis to strengthen loop invariants. This technique helps reducing the size of annotations and the burden of interactive specification. To that end, we assume a concretization function γ a : D → A to interpret abstract elements d ∈ D as assertions.

VCgen for Source Programs

Consider a specification (pre, annot, post, χ) for the source program P . In this section, we assume that annot sufficiently annotates the program P , that is, for every subprogram while [t] l do c of P , we have that l ∈ dom(annot).

A VCgen for source programs is defined by the set of proof obligations PO defined as {pre ⇒ φ[ V / V old ]} ∪ θ, where φ, θ = WP(P, post), φ[ V / V old ] represents the result of substituting in φ any array or scalar variable x old in V old s + V old a by x, and the function WP is defined in Figure 6. In the figure, the assertion inB(e) stands for the condition that must satisfy an execution state to ensure that every array access in e is within bounds. For instance, if e does not contain array expressions inB(e) is defined as true and inB(a[e]) as 0 ≤ e < |a|. We follow the simplifying assumption that expressions contain no more than one array access. For any array variable a and expressions e 1 and e 2 , upd(a, e 1 , e 2 ) is interpreted as the array a ′ such that a ′ [e] is evaluated to e 2 if e 1 = e and to a[e] otherwise. To simplify the presentation of examples, proof obligations for while statements are split into two assertions corresponding to the true and false branches.

The function WP considers the result of the analysis Loc to reduce the size of proof obligations. That is, if the abstract value Loc(l) associated to the program point under consideration indicates that any array access in the statement is within bounds, the returned predicate is simplified by omitting the exceptional postcondition. Consider the program of Figure 7. If the analysis is able to compute at label k 1 an abstract value d such that d |= 0 ≤ i < |A|, the WP function will return the assertion

upd(A, i, A[0])[i + 1 -1] = A[0]
, which together with the loop invariant at label k yields the proof obligation:

A[i -1] = 0 ∧ 0 ≤ i ≤ |A| ⇒ i < |A| ⇒ upd(A, i, A[0])[i + 1 -1] = A[0]
where the boxed assertion 0 ≤ i ≤ |A| represents the result of the analysis at the loop entry point.

In contrast, if we do not take advantage of the result of the analysis we must prove the bigger formula:

A[i -1] = 0 ∧ 0 ≤ i ≤ |A| ⇒ i < |A| ⇒ (0 ≤ i < |A| ⇒ upd(A, i, A[0])[i + 1 -1] = A[0] ∧¬(0 ≤ i < |A|) ⇒ false)
As can be seen from the definition of WP, proof obligations are of the form φ 1 ∧ γ a (d) ⇒ φ 2 , whereas a standard VCgen outputs the stronger proof obligation φ 1 ⇒ φ 2 . In consequence, one can provide the code with a weaker invariant φ 1 as long as the analyzer is able to eventually infer the missing information γ a (d). For instance, for the simple program of Figure 7, a standard VCgen will return the invalid proof obligation

A[i -1] = A[0] ⇒ ¬(i < |A|) ⇒ A[|A| -1] = A[0]
for the path that does not enter the loop. It is sufficient to provide a stronger invariant, i.e. in conjunction with the condition i ≤ |A|, to prove the program correct. However, as an alternative to increasing the size of the program annotations, assuming the condition i ≤ |A| is inferred by the analysis, the hybrid VCgen generates the weaker (and valid) proof obligation:

A[i -1] = A[0] ∧ 0 ≤ i ≤ |A| ⇒ ¬(i < |A|) ⇒ A[|A| -1] = A[0]
VCgen for Bytecode Programs Let (pre, annot, post, χ) be a specification for the bytecode program ṗ. As with the VCgen for source programs defined above, the precondition pre and the internal annotations annot(l) are strengthened with the result of the analysis. To that end, we interpret the result of the analysis with the aid of the concretization functions γ a : D → A and γa : (Expr ⋆ × D) → A. A VCgen for bytecode is defined by extracting the set of proof obligations:

po = {pre ⇒ wpi(0)[ V / V old ]}∪ {annot(l) ∧ γa (l ȯc(l)) ⇒ wpi(l) | l ∈ dom(annot)}
where the predicate transformer wp is shown in Figure 8.

If the program point is annotated, the function wp returns annot(l). Otherwise it applies the weakest precondition transformer wpi, defined in terms of the instruction at program point l, taking as parameters the annotations computed for the successor program points. The definition of wp and wpi is done by induction along the control flow paths of the program. A program ṗ is sufficiently annotated if the control flow graph of the program ṗ does not contain unannotated loops. The induction principle following from the definition of sufficiently annotated programs is sufficient to ensure that wp and wpi are well defined For a list s, s[0] and s [1] represent the first and second element of s, and ↑ s denotes the result of removing the first element from s.

Preservation of Proof Obligations

Consider the specification (pre, annot, post, χ) for source program P , and assume that annot is a sufficient annotation for P , i.e. every loop is annotated. Let (pre, annot, post, χ) define as well the specification for the bytecode program ṗ.

From previous results, we know that if annot is a sufficient annotation for P then it is also a sufficient annotation for the result of the compilation ṗ. Let Loc be a solution of the analysis for the source program P , and l ȯc a solution of the analysis for the bytecode program ṗ, compiled from Loc as described in Section 4.

We assume that the concretization functions satisfy the property γa ([], d) = γ a (d), so that the interpretation of abstract analysis results in the source and bytecode sides coincides (recall that by definition l ȯc(l) = ([ ], Loc(l)) for every l in dom(Loc).) In addition, for any expression e and any d ∈ D, if e does not contain array expressions, i.e. inB(e) = true, then d |= inB(e).

The following result about the compilation of expressions is helpful to prove preservation of proof obligations: Lemma 9 Let ṗ be equal to ṗ1 :: l 1 : e e :: l 2 : ṗ2 . Then wpi(l 1 ) is equal to wpi(l 2 )[ e::s / s ] if l ȯc(l 1 ) |= inB(e) and equal to inB(e) ⇒ wpi(l 2 )[ e::s / s ] ∧ ¬inB(e) ⇒ χ otherwise.

The coincidence of the sets of proof obligations PO and po is stated in the following lemma, provided the bytecode program ṗ is the result of compiling the source program P .

Proposition 10 For every subprogram c of P , proof obligations corresponding to c are equal to the proof obligations in ṗ that correspond to the subsequence c .

Consider, the bytecode program of Figure 9 compiled from the example in Figure 7. One can see that the proof obligation at label k is

A[i -1] = A[0] ∧ 0 ≤ i ≤ |A| ⇒ (i < |A| ⇒ (A[i -1] = A[0])[ upd(A,i,A[0]),i+1 / A,i ]) ∧ (¬(i < |A|) ⇒ A[|A| -1] = A[0])
which is equal to the proof obligation at label k for the source program of Figure 7.

From hybrid VCgen to VCgen

In this section we show a correspondence between the hybrid VCgen for bytecode of previous sections with a standard VCgen that does not take advantage of the result of the analysis. More precisely, interpreting the abstract result as logical formulae, we show an equivalence between the proof obligations of both VCgen's. Assuming that the relation |= satisfies a correctness condition, soundness of Given a specification (pre, ânnot, post, χ) for the bytecode program ṗ, a non-hybrid VCgen extracts the set of proof obligations po ∪ {pre ⇒ ŵpi(0)[ V / V old ]}, where ŵpi and po are defined in Figure 10. To avoid ambiguity, in the sequel we make explicit some parameters needed in the definition of wpi, wp, ŵpi and ŵp. We write for instance ŵpi(l, ânnot, post, χ) instead of simply ŵpi(l).

Let l ȯc be a result of the analysis for the bytecode program ṗ. Consider the specifications (pre, annot, post, χ) and (pre, ânnot, post, χ) for program ṗ, such that for all l in dom(annot), ânnot(l) is defined as annot(l) ∧ γa (l ȯc(l)).

(l) = wp(l + 1)[ s[0] op s[1]::↑ 2 s /s] ṗ[l] = prim op wpi(l) = wp(l + 1)[ v::s /s] ṗ[l] = push v wpi(l) = wp(l + 1)[ s[0],↑s /x,s] ṗ[l] = store x wpi(l) = wp(l + 1)[ x::s /s] ṗ[l] = load x wpi(l) = ckB(wp(l + 1)[ upd(a,s[0],s[1]),↑ 2 s /a,s]) ṗ[l] = astore a wpi(l) = ckB(wp(l + 1)[ a[s[0]]::↑s /s]) ṗ[l] = aload a wpi(l) = s[0] ⋊ ⋉ s[1] ⇒ wp(l ′ )[ ↑ 2 s /s] ṗ[l] = cjmp ⋊ ⋉ l ′ ∧¬(s[0] ⋊ ⋉ s[1]) ⇒ wp(l + 1)[ ↑ 2 s /s] wpi(l) = wp(l ′ ) ṗ[l] = jmp l ′ wpi(l) = wp(l + 1) ṗ[l] = nop wpi(l) = post[ s[0] /res] ṗ[l] = return where ckB(ψ) stands for ψ if l ȯc(l) |= inB(x[s[0]]) and inB(x[s[0]]) ⇒ ψ ∧ ¬inB(x[s[0]]) ⇒ χ otherwise.
[i -1] = A[0] and A[|A| -1] = A[0]
in k and k ′′ respectively. Let ânnot be defined by strengthening annot with the result of the analysis, i.e. ânnot(k) = annot(k) ∧ 0 ≤ i ≤ |A| (we can let ânnot(k ′′ ) = annot(k ′′ )). Let Ψ be the weakest precondition computed by the non hybrid VCgen at label k 1 :

0 ≤ i < |A| ⇒ (upd(A, i, A[0])[i + 1 -1] = A[0] ∧ 0 ≤ i + 1 ≤ |A|) ∧¬(0 ≤ i < |A|) ⇒ false
which, from Lemma 11 is implied by the hybrid wp and the

k 0 : push 1 store i k : jmp k ′ k 1 : push 0 aload A load i astore A k 2 : push 1 load i prim + store i k ′ : push |A| load i cjmp < k 1 k ′′ : . . . Figure 9. Program example ŵpi(l) = ckB( ŵp(l + 1)[ upd(x,s[0],s[1]),↑ 2 s /x,s]) ṗ[l] = astore x ŵpi(l) = ckB( ŵp(l + 1)[ x[s[0]]::↑s /s]) ṗ[l] = aload x ŵpi(l) = wpi(l) otherwise where ckB(ψ) stands for inB(x[s[0]]) ⇒ ψ ∧ ¬inB(x[s[0]]) ⇒ χ re- gardless of whether l ȯc(l) |= inB(x[s[0]]) is satisfied. ŵp(l) =  annot(l) if l ∈ dom(annot) ŵpi(l) otherwise po = { ânnot(l) ⇒ ŵpi(l) | l ∈ dom( ânnot)}
[i + 1 -1] = A[0] ∧ 0 ≤ i < |A|
As stated in Proposition 12, if the proof obligations returned by the hybrid VCgen are valid, and assuming the analysis is verifiable, we have that

A[i -1] = A[0] ∧ 0 ≤ i ≤ |A| ⇒ i < |A| ⇒ upd(A, i, A[0])[i + 1 -1] = A[0] and 0 ≤ i ≤ |A| ⇒ i < |A| ⇒ 0 ≤ i < |A| are provable.
Then, it follows that the verification condition

A[i -1] = A[0] ∧ 0 ≤ i ≤ |A| ⇒ i < |A| ⇒ Ψ
returned by the standard VCgen is provable. The above results establish that hybrid verification methods can be mapped to standard verification methods. In the context of PCC, one would like to establish the stronger result that hybrid certificates can be compiled into standard certificates. It is in fact possible to prove such a result, using the framework of [START_REF] Barthe | Certificate translation in abstract interpretation[END_REF]. However, the compilation of hybrid certificates into standard certificates requires using a certifying analyzer, that generates automatically logical proofs of correctness of the results of the analysis. While it is possible to avoid hybrid methods altogether, e.g. to rely on standard PCC architectures, hybrid methods are beneficial both for the code producer because they significantly reduce the number of proof obligations required to certify code, and for the code consumer, because they yield certificates that are more compact and more efficient to check. Cetificate Translation from a hybrid to a standard VCgen is interesting for PCC scenarios in which original proof obligations are generated by a hybrid VCgen, but in which the targeted Trusted Computed Base has no support for hybrid certificates.

Related work

Some authors have formalized and proved the soundness of hybrid verification methods. For example, Wildmoser, Chaieb and Nipkow [START_REF] Wildmoser | Bytecode analysis for proof carrying code[END_REF] have used Isabelle/HOL to prove the soundness of hybrid methods for Java bytecode; they rely on interval analyses to detect arrays out of bounds, and implement a proof-producing version of the analysis that generates proofs that the results of the analysis is correct. More recently, Grégoire and Sacchini [START_REF] Grégoire | Combining a verification condition generator for a bytecode language with static analyses[END_REF] have formalized in Coq a hybrid verification method for JVM programs. They focus on a null-pointer analysis; although the analyzer is not formalized in Coq, Hubert and co-workers [START_REF] Hubert | Semantic foundations and inference of non-null annotations[END_REF] provides a good starting point for carrying such an implementation. The method of Grégoire and Sacchini [START_REF] Grégoire | Combining a verification condition generator for a bytecode language with static analyses[END_REF] supports bidirectional interaction between the analysis and verification condition generation, as the static analysis can extract useful information from the program annotations, at the same time as the results of the analysis are exploited by the verification condition generator to reduce the number of proof obligations (although we have not done so, it is relatively easy to integrate such bidirectional interaction in our work).

It is folklore that deductive verification methods can be viewed as abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF][START_REF] Cousot | Semantic foundations of program analysis[END_REF]. Logical abstract interpretation [START_REF] Gulwani | Combining abstract interpreters[END_REF] explores the interaction between analysis and verification from the perspective of using theorem proving to improve the precision of abstract interpretations, and combinations of them.

Finally, the paper is closely related to previous works on proof-transforming compilation, and proof-producing program analyses. Saabas and Uustalu [START_REF] Saabas | A compositional natural semantics and Hoare logic for low-level languages[END_REF] provide an algorithm to transform proofs in Hoare logic into proofs in compositional proof systems for assembly programs. Moving to more realistic languages, Müller et al. define prooftransforming compilation for Java and Eiffel [1,[START_REF] Müller | Proof-transforming compilation of programs with abrupt termination[END_REF][START_REF] Nordio | Proof-transforming compilation of eiffel programs[END_REF],

The aforementioned works, as the current paper, focus on non-optimizing compilation. Compiling proofs along program optimizations require using proof-producing analyses, that produce formal proofs of their correctness. Such analyses are also required to extend the results of Section 3 to certificates. Proof-producing analyses have been studied by several authors, including Wildmoser, Chaieb and Nip-kow [START_REF] Wildmoser | Bytecode analysis for proof carrying code[END_REF] in the context of verification condition generation for a bytecode language and Seo, Yang and Yi [START_REF] Seo | Automatic Construction of Hoare Proofs from Abstract Interpretation Results[END_REF] in the context of a Hoare logic for a simple imperative language.

Conclusion

Program verification environments increasingly rely on hybrid methods to prove software correctness. In this paper, we have shown that hybrid verification methods at source and bytecode levels are tightly related, and that hybrid verification methods can be "compiled" into standard verification methods, which ensure that hybrid methods are sound.

Our next goal is to extend our results to more realistic languages and analyses. Modern languages include features, e.g. exceptions, that can potentially yield very large control flow graphs, making hybrid certificates particularly necessary; we expect that our results on preservation of proof obligations for Java [START_REF] Barthe | Preservation of proof obligations for Java[END_REF] will scale to hybrid methods.

Besides, it would be beneficial to allow hybrid methods to rely on more advanced static analyses that provide valuable information for proving properties of programs. We intend to focus on the recent analysis of [10], and to develop a hybrid verification method based on this analysis.

  n e = push n x e = load x x[e] e = e e; aload x e 1 op e 2 e = e 2 e; e 1 e; prim op [x:=e] k = k : e e; store x [a[e 1 ]:=e 2 ] k = k : e 2 e; e 1 e; astore a s1

Figure 2 .

 2 Figure 2. Semantics of bytecode (excerpts)

Figure 3 .

 3 Figure 3. Requirements over function γ

Figure 4 .

 4 Figure 4. Definition of the constraint system for the source code analysis.

Definition 6 (

 6 Compiled analysis results) Given an analysis result Loc for the program P , an analysis result com-piled from Loc is the completion of the function loc partial defined on each k ∈ dom(Loc) by loc partial (k) = ([ ], Loc(k)).

  For the abstract domain D, we consider a relation |= ⊆ D × A such that for any guard b and any d ∈ D, d |= b indicates that the interpretation of the abstract element d ensures the validity of the condition b. For example, when accessing an array in the expression a[x] we shall check that the value of the variable x is within the bounds of the array a. If we instantiate D with the domain of convex polyhedra, each element d ∈ D represents a set of linear constraints from which we can discover whether the condition 0 ≤ x < |a| is satisfied.

Figure 7 .

 7 Figure 7. Program example

  We say that the relation |= ⊆ D × Guard is valid if for every abstract element d ∈ D and b ∈ Guard we have that d |= b implies the universal validity of γ a (d) ⇒ b. The result of the analysis l ȯc is said verifiable if the set of proof obligations po(true, γa • l ȯc, true, true) are provable. Lemma 11 For every label l in the program ṗ: wpi(l, annot, post, χ)∧ γa (l ȯc(l)) ⇒ ŵpi(l, ânnot, post, χ) provided the relation |= ⊆ D × Guard is valid, and the analysis l ȯc is verifiable.

  wpi

Figure 8 .

 8 Figure 8. VCgen for bytecode programs

Figure 10 .

 10 Figure 10. Non-hybrid bytecode VCgen


Alberto Luis

Naomi Escudero

G Feldheim

Jos É Luis Ünther Koliander

Romero

EFFICIENT COMPUTATION OF THE ZEROS OF THE BARGMANN TRANSFORM UNDER ADDITIVE WHITE NOISE

Keywords: 2020 Mathematics Subject Classification. 65R10, 62M30, 60G70, 60G55, 60G15, 30H20 Bargmann transform, random analytic function, short-time Fourier transform, zero set, computation, Wasserstein metric. L

We study the computation of the zero set of the Bargmann transform of a signal contaminated with complex white noise, or, equivalently, the computation of the zeros of its short-time Fourier transform with Gaussian window. We introduce the adaptive minimal grid neighbors algorithm (AMN), a variant of a method that has recently appeared in the signal processing literature, and prove that with high probability it computes the desired zero set. More precisely, given samples of the Bargmann transform of a signal on a finite grid with spacing δ, AMN is shown to compute the desired zero set up to a factor of δ in the Wasserstein error metric, with failure probability O(δ 4 log 2 (1/δ)). We also provide numerical tests and comparison with other algorithms.

1. Introduction 1.1. The Bargmann transform and its zeros. The Bargmann transform of a real variable function f ∈ L 2 (R) is the entire function

F (z) = 2 π e -z 2 /2 R f (t)e -t 2 +2tz dt, z ∈ C. (1.1)
Originally introduced as a link between configuration and phase space in quantum mechanics [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform[END_REF], the Bargmann transform was later recognized as a powerful tool in signal analysis [START_REF] Daubechies | Frames in the Bargmann space of entire functions[END_REF] because it encodes the correlations between the signal f and the time-frequency shifts of the Gaussian function g(t) = 2 In the jargon of time-frequency analysis, the right-hand side of (1.2) is known as the short-time Fourier transform of f with Gaussian window, and measures the contribution to f (t) of the frequency y near t = x.

In practice, the values of the short-time Fourier transform of a signal f are only available on (a finite subset of) a grid {(δk, δj) : k, j ∈ Z}, δ > 0. (1.3) The goal of Gabor analysis is to extract useful information about f from such limited measurements. Equivalently, by (1.2), the task is to capture the analytic function F given a limited number of its samples on a grid. This second point of view led to the most conclusive results in Gabor theory, such as the complete description of all grids (1.3) for which the Gabor transform fully retains the original analog signal f [START_REF] Daubechies | Frames in the Bargmann space of entire functions[END_REF][START_REF] Lyubarskiȋ | Frames in the Bargmann space of entire functions[END_REF][START_REF] Seip | Density theorems for sampling and interpolation in the Bargmann-Fock space[END_REF][START_REF] Seip | Density theorems for sampling and interpolation in the Bargmann-Fock space[END_REF].

While Gabor signal analysis has traditionally focused on large values of the short-time Fourier transform (1.2), recent work has brought to the foreground the rich information stored in its zeros, especially when the signal is contaminated with noise. Heuristically, the zeros of the Bargmann transform of noise exhibit a rather rigid random pattern with predictable statistics, from which the presence of a deterministic signal can be recognized as a salient local perturbation [START_REF] Gardner | Sparse time-frequency representations[END_REF][START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF][START_REF] Flandrin | The sound of silence: Recovering signals from time-frequency zeros[END_REF]. Remarkably, the Bargmann transform of white noise has been identified as a certain Gaussian analytic random function [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF][START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF], and consequently the wellresearched statistics of their zero sets [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF][START_REF] Nazarov | What is. . .a Gaussian entire function?[END_REF] can be leveraged in practice [START_REF] Flandrin | Explorations in time-frequency analysis[END_REF]Chapters 13 and 15], [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF]. The particular structure observed in the zeros of the Bargmann transform under even a moderate amount of white noise has also been invoked as explanation for the sparsity resulting from certain non-linear procedures to sharpen spectrograms [START_REF] Gardner | Sparse time-frequency representations[END_REF], as the zeros of the Bargmann transform are repellers of the reassignment vector field [START_REF] Flandrin | Explorations in time-frequency analysis[END_REF]Chapter 12]. The practical exploitation of such insights requires an effective computational link between finitely given data on the one hand and zeros of Bargmann transforms of an analog signal on the other. 1.2. Computation of zero sets. Suppose that the values of the Bargmann transform F of a signal f are given on a grid Λ = {δk + iδj : k, j ∈ Z}, δ > 0, (1.4) and we wish to compute an approximation of {F = 0}, the zero set of F , within the square Ω L = {x + iy : |x|, |y| ≤ L}. (1.5) More realistically, we only have access to samples of F on those grid points near the computation domain, e.g., on Λ L = {δk + iδj : k, j ∈ Z, |δk|, |δj| ≤ L}. (1.6) The inverse of the spacing of the grid, 1/δ, will be called the resolution of the data. Thresholding. The most naive approach to compute {F = 0} is thresholding: one selects all grid points λ such that |F (λ)| is below a certain threshold ε > 0:

e -1 2 |λ| 2 |F (λ)| < ε. (1.7)
The normalizing weight e -1 2 |λ| 2 is motivated by (1.2), as the short-time Fourier transform of a typical signal can be expected to be bounded. One disadvantage of this approach is that it requires an educated choice for the threshold ε. Moreover, computations with various reasonable choices of thresholds, such as quantiles of e -1 2 |λ| 2 |F (λ)| calculated over all grid points λ, either fail to compute many of the zeros or capture too many points (see Figure 1).

Extrapolation. One may consider using the samples of F on the finite grid (1.6) to reconstruct the signal f , resample F at arbitrarily high density, and thus calculate more easily the zero set {F = 0}. However, computation of zeros through extrapolation may be inaccurate: while the samples of F on the infinite grid (1.4) determine F as soon as √ π • δ < 1 [START_REF] Daubechies | Frames in the Bargmann space of entire functions[END_REF][START_REF] Lyubarskiȋ | Frames in the Bargmann space of entire functions[END_REF][START_REF] Seip | Density theorems for sampling and interpolation in the Bargmann-Fock space[END_REF], the truncation errors involved in the approximation of F near Ω L from finite data (1.6) can only be neglected at very high resolution 1/δ. Even if the values of F are successfully extrapolated to a higher resolution grid, the remaining computation is still not trivial, as, for example, simple thresholding may fail even at high resolution (see Figure 4 and Section 5).

Minimal Grid Neighbors. A greatly effective numerical recipe for the computation of zeros of the Bargmann transform can be found in the code accompanying [START_REF] Flandrin | Time-frequency filtering based on spectrogram zeros[END_REF] -although not explicitly described in the text. A grid point λ is selected as a numerical approximation of a zero if e -

|F (λ)| ≤ e -1 2 |µ| 2 |F (µ)|, |λ -µ| ∞ = δ, (1.8)
where |z| ∞ = max{|x|, |y|}. The subset of points that pass the test furnish the computation of {F = 0}. This method, which we call minimal grid neighbors (MGN), performs impressively as long as the grid resolution is moderately high. Indeed, we understand that the method is behind the simulations in [START_REF] Flandrin | Explorations in time-frequency analysis[END_REF]Chapter 15] which quite faithfully reproduce the statistics of the zeros of the Bargmann transform of complex white noise (that are known analytically [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF][START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]). The MGN algorithm was also used to produce the plots in [START_REF] Bardenet | On the zeros of the spectrogram of white noise[END_REF], as pointed out in [5, Section 5.1. 1.]; see also [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF]Section 5], [20, Section IV], and [START_REF] Abreu | Filtering with wavelet zeros and gaussian analytic functions[END_REF]. Heuristically, the test (1.8) succeeds in identifying zeros due to the analyticity of F , which implies that |F (z)|e -1 2 |z| 2 does not have non-zero local minima. Remarkably, (1.8) is also effective even if the comparison involves only neighboring grid points.

The MGN algorithm performs equally well when calculating the zeros of the Bargmann transform of a signal

f = f 1 + σ • N (1.9)
composed of a deterministic real-variable function f 1 plus complex white noise with variance σ 2 > 0. The presence of a certain amount of randomness must be behind the success of the algorithm, as, for σ = 0, the method cannot be expected to succeed. Indeed, one can engineer a deterministic signal f where the detection of zeros fails, as the value of its Bargmann transform F can be freely prescribed on any given finite subset of the computation domain [START_REF] Seip | Density theorems for sampling and interpolation in the Bargmann-Fock space[END_REF]. We are unaware of performance guarantees for MGN. 1.3. Contribution. In this article we introduce a variant of MGN, called Adaptive Minimal Grid Neighbors (AMN). The algorithm is based on a comparison similar to (1.8) but incorporates an adaptive decision margin, that depends on the particular realization of F . While AMN has the same mild computational complexity and similar practical effectiveness as MGN, we are able to estimate the accuracy and confidence of the computation with AMN in terms of the grid resolution. In this way, we show that AMN is probably approximately correct for the signal model (1.9), in the sense that it computes the zero set with high probability up to the resolution of the data.

On the one hand, we present what to the best of our knowledge are the first formal guarantees for the approximate computation of zero sets of analytic functions from grid values. In fact, besides its main purpose of computation with specific data, the AMN algorithm offers a computationally attractive and provably correct method to simulate zero sets of the Gaussian entire function (2.10), by running the procedure with simulated inputs. On the other hand, our analysis is a first step towards understanding the performance of MGN.

Main Result

2.1. The adaptive minimal neighbors algorithm. We now introduce a new algorithm to compute zero sets of Bargmann transforms. Suppose again that samples of a function F : C → C are given on those points of the grid (1.4) that are near the computation domain (1.5), say, on

Λ L+2δ = {δk + iδj : k, j ∈ Z, |δk|, |δj| ≤ L + 2δ}. (2.1)
For each grid point λ strictly inside the computation domain, λ ∈ Λ L = Λ∩Ω L , we use the neighboring sample at λ + δ to compute the following comparison margin:

η λ = e -1 2 |λ| 2 max |F (λ)| , 3 4 e -δ λF (λ + δ) -F (λ) . (2.2)
The margin η λ therefore depends on the particular realization of F . To motivate the definition, note that, when λ = 0, the maximum in (2.2) is taken over |F (0)| and the absolute value of an incremental approximation of F ′ (0). The comparison margin thus incorporates the size and oscillation of F at z = 0. In general, η λ has a similar interpretation with respect to the covariant derivative

∂ * F (z) = z F (z) -∂F (z), (2.3)
and, indeed, η λ is defined so that

η λ ≈ e -1 2 |λ| 2 max |F (λ)|, 3 4 ∂ * F (λ) δ . (2.4)
The differential operator ∂ * F plays a distinguished role in the analysis of vanishing orders of Bargmann transforms [START_REF] Brekke | Density theorems for sampling and interpolation in the Bargmann-Fock space[END_REF][START_REF] Escudero | Multiple sampling and interpolation in weighted Fock spaces of entire functions[END_REF] because it commutes with the translational symmetries of the space that they generate (the Bargmann-Fock shifts defined in Section 3.2).

The first step of the algorithm selects all grid points λ ∈ Λ L that pass the following comparison test:

e -1 2 |µ| 2 |F (µ)| ≥ e -1 2 |λ| 2 |F (λ)| + η λ , whenever |λ -µ| ∞ = 2δ, µ ∈ Λ. (2.5)
In contrast to (1.8), the comparison in (2.5) does not involve the immediate grid neighbors of λ but rather those points lying on the square centered at λ with half-side-length 2δ; see Figure 2. (In particular, the test only involves grid points µ ∈ Λ L+2δ .) Intuitively, the larger distance between λ and µ permits neglecting the error in the differential approximation (2.4). Selection step: For each grid point λ ∈ Λ L = Λ ∩ Ω L inside the target domain Ω L , we define the following comparison margin:

η λ = e -1 2 |λ| 2 max |F (λ)| , 3 4 e -δ λF (λ + δ) -F (λ) . (2.6)
The grid point λ is then selected if the following test is satisfied:

e -1 2 |µ| 2 |F (µ)| ≥ e -1 2 |λ| 2 |F (λ)| + η λ , whenever |λ -µ| ∞ = 2δ, µ ∈ Λ. (2.7)
(Note that the test (2.7) only involves grid points µ ∈ Λ L+2δ .) Let Z 1 be the set of all selected grid points.

Sieving step: Use an off-the-shelf clustering algorithm to select a subset Z ⊂ Z 1 that is 5δ separated:

inf |λ -µ| ∞ : λ, µ ∈ Z, λ ̸ = µ ≥ 5δ (2.8)
and maximal with respect to this property, i.e., no proper superset satisfies (2.8).

(One concrete implementation of the sieving step is described in Section 5.2.1.)

Output: The set Z.

Performance guarantees for AMN.

To study the performance of the AMN algorithm we introduce the following input model, which, as we will explain, corresponds to the Bargmann transform of an arbitrary signal contaminated with complex white noise of an arbitrary intensity.

Input model. We consider a random entire function on the complex plane

F = F 1 + σ • F 0 , (2.9)
where F 1 : C → C is a deterministic entire function, F 0 is a (zero mean) Gaussian analytic function with correlation kernel:

E F 0 (z) • F 0 (w) = e z w,
z, w ∈ C, (2.10) and σ > 0 is the noise level. We assume that the deterministic function F 1 satisfies the quadratic exponential growth estimate

|F 1 (z)| ≤ A • e 1 2 |z| 2 , z ∈ C, (2.11)
for some constant A ≥ 0.

As for F 0 , the assumption means that for each z 1 , . . . , z n ∈ C, (F 0 (z 1 ), . . . , F 0 (z n )) is a normally distributed (circularly symmetric) complex random vector, with mean zero and covariance matrix e z k z ℓ k,ℓ . Alternatively, F 0 can be described as

F 0 (z) = n≥0 ξn √ n! z n , (2.12)
where (ξ n ) n≥0 are independent standard complex random variables [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF].

Discussion of the model. The random function F 0 is the Bargmann transform of standard complex white noise N . As each realization of complex white noise is a tempered distribution, the computation of its Bargmann transform (1.1) is also to be understood in the sense of distributions, as in [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform. Part II. A family of related function spaces. Application to distribution theory[END_REF]. (See [START_REF] Bardenet | Time-frequency transforms of white noises and Gaussian analytic functions[END_REF] and [START_REF] Haimi | Zeros of Gaussian Weyl-Heisenberg functions and hyperuniformity of charge[END_REF]Section 5.1] for a detailed discussion on this, and alternative approaches.)

We can similarly interpret F 1 as the Bargmann transform of a distribution f 1 on the real line [START_REF] Bargmann | On a Hilbert space of analytic functions and an associated integral transform. Part II. A family of related function spaces. Application to distribution theory[END_REF]. The assumption (2.11) means precisely that f 1 belongs to the modulation space M ∞ (R) consisting of distributions with Bargmann transform bounded with respect to the standard Gaussian weight -or, equivalently, with bounded short-time Fourier transform [START_REF] Bényi | Modulation Spaces: With Applications to Pseudodifferential Operators and Nonlinear Schrödinger Equations[END_REF]. The modulation space M ∞ (R) includes all square-integrable functions f 1 ∈ L 2 (R) and also many of the standard distributions used in signal processing.

In summary, the input model (2.9) corresponds exactly to the Bargmann transform of a random signal

f = f 1 + σ • N , (2.13)
where f 1 ∈ M ∞ (R) and σ • N is complex white noise with standard deviation σ.

Performance analysis. We now present the following performance guarantees, pertaining to the computation domain (1.5) and the acquisition grid (2.1). To avoid immaterial technicalities, we assume that the corners of the computation domain lie on the acquisition grid.

Theorem 2.1. Fix a domain width L ≥ 1, a noise level σ > 0, and a grid spacing δ > 0 such that L/δ ∈ N. Let a realization of a random function F as in (2.9) with (2.10) and (2.11) be observed on Λ L+2δ , and let Z be the output of the AMN algorithm.

There exists an absolute constant C such that, with probability at least

1 -C L 2 exp A 2 8σ 2 max 1, log 2 (1/δ) δ 4 , (2.14)
there is an injective map Φ : {F = 0} ∩ Ω L → Z with the following properties:

• (Each zero is mapped into a near-by numerical zero)

|Φ(ζ) -ζ| ∞ ≤ 2δ, ζ ∈ {F = 0} ∩ Ω L . (2.15)
• (Each numerical zero that is away from the boundary arises in this form)

For each λ ∈ Z ∩ Ω L-2δ there exists ζ ∈ {F = 0} ∩ Ω L such that λ = Φ(ζ).
A proof of Theorem 2.1 is presented in Section 4. We remark some aspects of the result.

• The AMN algorithm does not require knowledge of the noise level σ and is homogeneous in the sense that F and cF , with c ∈ C \ {0} produce the same output.

• Within the estimated success probability, the computation is accurate up to a factor of the grid spacing.

• The analysis concerns an arbitrary deterministic signal impacted by noise and is uniform over the class (2.11). As usual in such smoothed analysis, the success probability grows as the signal to noise ratio A/σ decreases, because randomness helps preclude the very untypical features that could cause the algorithm to fail [START_REF] Spielman | Smoothed analysis: an attempt to explain the behavior of algorithms in practice[END_REF].

• Up to a small boundary effect, the guarantees in Theorem 2.1 comprise an estimate on the Wasserstein distance between the atomic measures supported on {F = 0} ∩ Ω L and on the computed set Z. More precisely, for a tolerance θ > 0 let us define the boundary-corrected Wasserstein pseudo-distance between two sets U, V ⊆ C as

W L,θ (U, V ) = inf Φ max z∈U |Φ(z) -z| ∞ ,
where the infimum is taken over all injective maps Φ :

U → V such that V ∩ Ω L-θ ⊆ Φ(U )
. (The definition is not symmetric on U and V , but this is not important for our purpose.) Then Theorem 2.1 reads

P W L,2δ ({F = 0} ∩ Ω L , Z) > 2δ ≤ CL 2 exp A 2 8σ 2 max 1, log 2 (1/δ) δ 4 .
• The presented analysis concerns a signal contaminated with complex white noise. This is a mathematical simplification; we believe that with more technical arguments a similar result can be derived for real white noise. The case of colored noise seems more challenging and will be the object of future work.

Numerical experiments.

In Section 5, we report on numerical experiments that compare the AMN and MGN algorithms. We also include a modified version of thresholding (ST), that uses a thereshold proportional to the grid spacing AMN and MGN ST and incorporates a sieving step as in AMN (while standard thresholding without sieving performs extremely poorly, as seen in Figure 1).

(a) f 1 (t) = exp(-t 2 ), A = 1 (b) f 1 (t) = t exp(-t 2 ), A = 100
The performance of AMN, MGN, and ST is first tested indirectly, by using these algorithms to simulate the zero sets of the random functions in the input model (2.9). We then compare theoretically derived statistics of the zeros of (2.9) to empirical statistics obtained from the output of AMN and MGN under various simulated realizations of (2.9).

Second, we perform a consistency experiment that aims at estimating the probability of computing a low-distortion parametrization of the zero set of F , as in Theorem 2.1. Specifically, we simulate a realization of the random input F sampled at high-resolution and use the output of AMN or MGN as a proxy for the ground truth {F = 0}. We then test the extent to which this set is captured by the output of AMN, MGN, or ST from lower resolution subsets of the same simulated data.

The performance of AMN and MGN is almost identical, although the minimal resolution at which MGN starts to perform well is slightly lower than that for AMN. (This is to be expected, as the constants 2 and 5 used in (2.7) and (2.8) are not adequate for low resolutions.) Both AMN and MGN significantly outperform ST. See also Figure 4 for an illustration.

The favorable performance of AMN is interesting also when the input is just noise, as it gives a fast and provably accurate method to simulate the zeros of the Gaussian entire function (2.10). (The simulations that we present in Section 5 use certain heuristic shortcuts to accelerate the simulation of the input (2.10) -see Section 5.1; although we formally do not analyze these, they are implicitly validated, as the simulated point process reproduces the expected theoretical statistics.)

All numerical experiments can be reproduced with openly-accessible software and our code is available at https://github.com/laescudero/discretezeros. Our implementation of the Bargmann transform uses [START_REF] Bailey | The fractional Fourier transform and applications[END_REF].

2.4. Organization. Section 3 introduces the notation and basic technical tools about analytic functions, Bargmann-Fock shifts, and their applications to random functions and their zeros. Theorem 2.1 is proved in Section 4, while numerical experiments are presented in detail in Section 5.

Preliminaries

3.1. Notation. For a complex number z = x + iy, we use the notation |z| ∞ = max{|x|, |y|}, while |z| denotes the usual absolute value. The zero set of F is denoted by {F = 0}. The differential of the (Lebesgue) area measure on the plane will be denoted for short dm, while the measure of a set E is |E|. With a slight abuse of notation, we also denote the cardinality of a finite set Z by |Z|. Squares on the complex plane are denoted by Q r (z) = {w ∈ C : |z -w| ∞ ≤ r}. For two non-negative functions f, g, we write f ≲ g if there exists an absolute constant C such that f (x) ≤ Cg(x), for all x. We write f ≍ g if f ≲ g and g ≲ f .

A Gaussian entire function (see [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]Ch. 2] and [START_REF] Nazarov | What is. . .a Gaussian entire function?[END_REF]) is a random function F : C → C that is almost surely entire, and such that for every z 1 , . . . , z n ∈ C, F (z 1 ), . . . , F (z n ) is a circularly symmetric complex normal vector. We will be only concerned with the random function F given in (2.9). We also use the notation (1.4), (1.5), (1.6), possibly for distinct values of L.

3.2.

Bargmann-Fock shifts and stationarity of amplitudes. The analysis of the AMN algorithm is more transparent when formulated in terms of the Bargmann-Fock shifts. For a function F : C → C we let

F w (z) = e -1 2 |w| 2 -zw F (z + w). (3.1)
The amplitude of an entire function F is defined as the weighted magnitude

G(z) = e -1 2 |z| 2 |F (z)|,
and satisfies

G w (z) := G(z + w) = e -1 2 |z| 2 |F w (z)|. (3.2)
The comparison margin of the AMN algorithm (2.6) can be expressed in terms of Bargmann-Fock shifts as

η λ = max |F λ (0)| , 3 4 F λ (δ) -F λ (0) , (3.3)
and leads to the approximation (2.4) because

|F ′ λ (0)| = e -1 2 |λ| 2 ∂ * F (λ) .
(Here and throughout we write F ′ λ (0) for [F λ ] ′ (0).) Similarly, in terms of amplitudes, the test (2.7) reads

G(µ) = e -1 2 |µ-λ| 2 |F λ (µ -λ)| ≥ |F λ (0)| + η λ = G(λ) + η λ , µ ∈ Λ, |µ -λ| ∞ = 2δ. (3.4)
With respect to the input model (2.9) we note that, if F 0 is the (zero mean) Gaussian entire function with correlation kernel (2.10), then the Bargmann-Fock shifts F 0 → F 0 w preserve the stochastics of F 0 , as they leave its covariance kernel invariant. As a consequence, for any w ∈ C, F 0 w (0), F 0 w ′ (0) are independent standard complex normal random variables (with zero mean and variance 1).

∆ log H(z) = ∆ |z| 2 2 -log |F (z)| = 2.
Hence, the maximum principle for subharmonic functions together with (3.6) implies that log H(z) and therefore |F (z)|e -1 2 |z| 2 is constant on D. For z ∈ D, we compute

0 = ∂ z |F (z)| 2 e -|z| 2 = ∂ z F (z) e -|z| 2 F (z) = F ′ (z) -zF (z) F (z) e -|z| 2 .
As F is non-vanishing on D, it follows that F ′ -zF = 0 on D, and therefore Lemma 3.2. Let F be as in (2.9). Then there exists an absolute constant C > 0 such that for all L ≥ 1 and t ≥ A,

0 = ∂ z [F ′ -zF ] = -F, on D.
P sup w∈Ω L ,|z|≤10 |z| -2 F w (z) -F w (0) + F ′ w (0)z > t ≤ CL 2 e -(t-A) 2 /(8σ 2 ) , P sup w∈Ω L ,|z|≤10 |z| -2 F w (z)e -1 2 |z| 2 -F w (0) + F ′ w (0)z > t ≤ CL 2 e -(t-A) 2 /(8σ 2 ) .
Proof. We consider the Taylor expansion of F :

F (z) = F (0) + F ′ (0)z + E 2 (z)z 2 ,
where

sup |z|≤10 |E 2 (z)| ≲ |ζ|≤12 |F (ζ)| dm(ζ) ≍ |ζ|≤12 |G(ζ)| dm(ζ).
We also note that for |z| ≤ 10,

|F (z) -F (z)e -1 2 |z| 2 | = |F (z)| 1 -e -1 2 |z| 2 ≲ |z| 2 |ζ|≤12 |F (ζ)| dm(ζ) ≍ |z| 2 |ζ|≤12 |G(ζ)| dm(ζ).
Hence, for |z| ≤ 10,

F (z) -F (0) + F ′ (0)z ≲ |z| 2 |ζ|≤12 |G(ζ)| dm(ζ), F (z)e -1 2 |z| 2 -F (0) + F ′ (0)z ≲ |z| 2 |ζ|≤12 |G(ζ)| dm(ζ).
We apply the previous bounds to F w , note that, by (3.2), G w (z) = G(z + w), and obtain that for |z| ≤ 10,

A w (z) := |z| -2 F w (z) -F w (0) + F ′ w (0)z ≲ |ζ-w|≤12 |G(ζ)| dm(ζ), B w (z) := |z| -2 F w (z)e -1 2 |z| 2 -F w (0) + F ′ w (0)z ≲ |ζ-w|≤12 |G(ζ)| dm(ζ). Hence sup |z|≤10,w∈Ω L A w (z) + B w (z) ≲ sup w∈Ω L |ζ-w|≤12 |G(ζ)| dm(ζ) ≲ sup |ζ|≤L+12 |G(ζ)|.
Let G 0 and G 1 be the weighted magnitudes corresponding to F 0 and F 1 , respectively. Then by (2.11),

|G(ζ)| ≤ σ • G 0 (ζ) + |G 1 (ζ)| ≤ A + σ • |G 0 (ζ)|, ζ ∈ C.
Hence,

P sup |ζ|≤L |G(ζ)| > t ≤ P sup |ζ|≤L G 0 (ζ) > t -A σ .
To conclude, we claim that the following excursion bound holds:

P sup |ζ|≤L |G 0 (ζ)| > t ≤ CL 2 e -t 2 /8 , t ≥ 0,
where C > 0 is an absolute constant. For L ≤ 1/4 this follows for example from [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]Lemma 2.4.4]. In general, we cover the domain with ≲ L 2 squares of the form w + [-1/4, 1/4] 2 , apply the previously mentioned bound to G 0 w (z) = G 0 (z + w), and use a union bound. This completes the proof. □ 3.5. Almost multiple zeros. It is easy to see that, almost surely, the random function (2.9) has no multiple zeros. In the analysis of the AMN algorithm, we will also need to control the occurrence of zeros that are multiple up to a certain numerical precision, in the sense that F and its derivative are simultaneously small. The following lemma is a first step in that direction, as it controls the probability of finding a grid point that is an almost multiple zero.

Lemma 3.3. Let F be as in (2.9) and α, β > 0. Then the probability that for some grid point λ ∈ Λ L the following occurs:

|F λ (0)| ≤ α, and |F ′ λ (0)| ≤ β (3.7) is at most CL 2 α 2 β 2 δ -2 σ -4
, where C is an absolute constant.

Proof. For each grid point λ ∈ Λ L , F λ (0) and F ′ λ (0) are independent complex normal variables with possibly non-zero means µ 1 , µ 2 and variance σ 2 . Therefore,

P |F λ (0)| ≤ α = 1 πσ 2 |ζ|≤α e -1 σ 2 |ζ-µ 1 | 2 dm(ζ), (3.8) P |F ′ λ (0)| ≤ β = 1 πσ 2 |ζ|≤β e -1 σ 2 |ζ-µ 2 | 2
dm(ζ). (3.9) By Anderson's lemma [START_REF] Anderson | The integral of a symmetric unimodal function over a symmetric convex set and some probability inequalities[END_REF], the right-hand sides of (3.8) and (3.9) are maximal when µ 1 = 0 and µ 2 = 0, respectively. Direct computation in those cases yields

P (|F λ (0)| ≤ α) ≲ α 2 σ -2 and P (|F ′ λ (0)| ≤ β) ≲ β 2 σ -2
. By independence, the probability of (3.7) is ≲ α 2 β 2 σ -4 . On the other hand, there are ≲ L 2 δ -2 grid points under consideration, so the conclusion follows from the union bound. □ 3.6. First intensity of zeros. The following proposition is not used in the proof of Theorem 2.1, but rather as a benchmark in the numerical experiments (Section 5). Proposition 3.4. Let F be as in (2.9). Then for every Borel set

B ⊂ C, E[|{z ∈ B : F (z) = 0}|] = B ρ 1 (ζ)dm(ζ)
where

(3.10) ρ 1 (ζ) = 1 π e -1 σ 2 |F 1 (ζ)| 2 e -|ζ| 2 1 + e -|ζ| 2 σ 2 ∂ ζ F 1 (ζ) -ζF 1 (ζ) 2 .
Proof. The set of zeros and thus ρ 1 (ζ) does not change if we scale F by a fixed constant. Hence, by considering the function 1 σ F in place of F , we can assume that σ = 1. The expected number of points {z ∈ B : F (z) = 0} of a Gaussian random field F is given by Kac-Rice's formula:

(3.11) E[|{z ∈ B : F (z) = 0}|] = B E |det DF (ζ)| F (ζ) = 0 p F (ζ) (0) dm(ζ),
see, e.g., [3, Th. 6.2]. Since F is analytic, the determinant in question can easily be seen to simplify to |det

DF (ζ)| = |∂ ζ F (ζ)| 2 . The joint vector (F (z), ∂ z F (z)) has mean (F 1 (z), ∂ z F 1 (z)) and covariance (3.12) Cov[(F (z), ∂ z F (z))] = e |z| 2 ze |z| 2 ze |z| 2 (1 + |z| 2 )e |z| 2 .
Following a Gaussian regression approach, see, e.g., [3, Prop.

E |det DF (ζ)| F (ζ) = 0 = e |ζ| 2 + ∂ ζ F 1 (ζ) -ζF 1 (ζ) 2 .
Inserting into (3.11) while noting that p

F (ζ) (0) = 1 π e -|ζ| 2 e -|F 1 (ζ)| 2 e -|ζ| 2
yields (3.10). □

Proof of Theorem 2.1

We present the proof of Theorem 2.1 in several steps. Assume all the hypothesis, and denote by Z 1 the set produced by the AMN algorithm after the selection step. Recall that L ≥ 1. By choosing a sufficiently large constant in (2.14), we can assume that δ ≤ 1/5; otherwise, the success probability would be trivial. For the same reason, we can assume that

δ 4 exp A 2 8σ 2 ≤ 1. (4.1)
4.1. Excluding bad events. We let γ = 8σ and wish to apply Lemma 3.2 with t = γ log(1/δ). By (4.1),

δ 4 exp A 2 16σ 2 ≤ δ 4 exp A 2 8σ 2 ≤ 1.
Hence t ≥ A and we can apply Lemma 3.2 to conclude that

F w (z) -F w (0) + F ′ w (0)z ≤ γ • log(1/δ) • |z| 2 ≤ 2γ • log(1/δ) • |z| 2 ∞ , (4.2)
F w (z)e -1 2 |z| 2 -F w (0) + F ′ w (0)z ≤ γ • log(1/δ) • |z| 2 ≤ 2γ • log(1/δ) • |z| 2 ∞ , (4.3)
for all w ∈ Ω L+1 and |z| ≤ 10, except for an event of probability at most

CL 2 exp -(t -A) 2 /(8σ 2 ) , where C is an absolute constant. Since (t -A) 2 ≥ t 2
2 -A 2 , we further have

CL 2 exp - (t -A) 2 8σ 2 ≤ CL 2 exp A 2 8σ 2 δ γ 2 16σ 2 = CL 2 exp A 2 8σ 2 δ 4 .
Second, we select a large absolute constant κ > 1 to be specified later, and use Lemma 3.3 with

α = κγ • log(1/δ) • δ 2 , β = 2κγ • log(1/δ) • δ, to conclude that, for each grid point λ ∈ Λ L+2δ , either |F λ (0)| > α, or |F ′ λ (0)| > β, ( possibly
both), (4.4) except for an event with probability at most ≲ L 2 log 2 (1/δ)δ 4 .

Overall we have excluded events with total probability

≲ L 2 exp A 2 8σ 2 log 2 (1/δ)δ 4 .
In what follows, we show that under the complementary events the conclusions of Theorem 2.1 hold.

4.2.

The true zeros are adequately separated. We claim that, by taking κ sufficiently large, the set {F = 0} ∩ Ω L+2δ satisfies:

inf |ζ -ζ ′ | ∞ : ζ, ζ ′ ∈ {F = 0} ∩ Ω L+2δ , ζ ̸ = ζ ′ > 7δ. (4.5) Suppose that ζ, ζ ′ ∈ {F = 0} ∩ Ω L+2δ are such that 0 < |ζ -ζ ′ | ∞ ≤ 7δ. Since L/δ ∈ N,
F ζ (0) = 0 while |ζ -ζ ′ | ≤ √ 2|ζ -ζ ′ | ∞ ≤ 7 √ 2δ ≤ 10 to obtain: F ′ ζ (0) |ζ ′ -ζ| ≤ γ • log(1/δ) • |ζ ′ -ζ| 2 . Since ζ ̸ = ζ ′ , we conclude: F ′ ζ (0) ≤ γ • log(1/δ) • |ζ ′ -ζ|. (4.6)
Second, we similarly apply (4.3) with w = ζ and z = λ -ζ, to obtain

F ζ (λ -ζ) • e -1 2 |λ-ζ| 2 -F ′ ζ (0) • (λ -ζ) ≤ γ • log(1/δ) • |λ -ζ| 2 .
Combining the last equation with (4.6) yields

F ζ (λ -ζ) • e -1 2 |λ-ζ| 2 ≤ γ • log(1/δ) • |λ -ζ| 2 + |ζ ′ -ζ| • |λ -ζ| . (4.7)
Third, we apply (4.2) with w = λ and z = ζ -λ to obtain

|F λ (0) + F ′ λ (0) • (ζ -λ)| ≤ γ • log(1/δ) • |ζ -λ| 2 . (4.8) Note that |F λ (0)| = |F ζ (λ -ζ)| • e -1 2 |λ-ζ| 2 .
Hence, combining (4.7) and (4.8) we obtain:

|F λ (0)| ≤ γ • log(1/δ) • |λ -ζ| 2 + |ζ ′ -ζ| • |λ -ζ| , |F ′ λ (0)| • |λ -ζ| ≤ γ • log(1/δ) • 2|λ -ζ| 2 + |ζ ′ -ζ| • |λ -ζ| . Since 0 < |λ -ζ| ≤ δ and |ζ -ζ ′ | ∞ ≤ 7δ, we conclude that |F λ (0)| ≤ γ • log(1/δ) • |λ -ζ| 2 + |ζ ′ -ζ| • |λ -ζ| ≤ 11 • γ • δ 2 • log(1/δ), |F ′ λ (0)| ≤ γ • log(1/δ) • (2|λ -ζ| + |ζ ′ -ζ|) ≤ 12 • γ • δ • log(1/δ
). Assuming as we may that κ > 11, this contradicts (4.4). Thus, (4.5) must indeed hold.

4.3.

Linearization holds with estimated slopes. For each λ ∈ Λ L , we use the notation

τ λ = F λ (δ) -F λ (0) δ ,
and observe that, by (4.2),

τ λ -F ′ λ (0) ≤ γ • log(1/δ) • δ. (4.9)
Combining this with (4.3), we conclude that for |z| ∞ ≤ 2δ and λ ∈ Λ L ,

F λ (z)e -1 2 |z| 2 -F λ (0) + τ λ z ≤ F λ (z)e -1 2 |z| 2 -F λ (0) + F ′ λ (0)z + |z| F ′ λ (0) -τ λ ≤ 2γ • log(1/δ) • |z| 2 ∞ + |z| • γ • log(1/δ) • δ (4.10) ≤ (8 + 2 √ 2) • γ • δ 2 • log(1/δ). ≤ 11 • γ • δ 2 • log(1/δ). (4.11)
4.4. After the selection step, each true zero is close to a computed zero. We show that

({F = 0} ∩ Ω L ) ⊆ Z 1 + Q δ/2 (0). (4.12) Let ζ ∈ Ω L be a zero of F . Since L/δ ∈ N, we can find λ ∈ Λ L such that |ζ -λ| ∞ ≤ δ/2. We show that λ ∈ Z 1 .
Let us first prove that

|τ λ | ≥ κγ log(1/δ) • δ. (4.13)
Suppose to the contrary that |τ λ | < κγ log(1/δ) • δ. We will show that this contradicts (4.4). Assuming as we may that κ ≥ 1, by (4.9),

|F ′ λ (0)| ≤ κγ + γ log(1/δ) • δ ≤ β, while, by (4.2), |F λ (0)| = |F λ (0) -F λ (ζ -λ)| ≤ F λ (ζ -λ) -F λ (0) + F ′ λ (0)(ζ -λ) + F ′ λ (0)(ζ -λ) ≤ γ log(1/δ)|ζ -λ| 2 + |F ′ λ (0)||ζ -λ| ≤ γ log(1/δ) √ 2δ 2 2 + |F ′ λ (0)| √ 2δ 2 ≤ γ 2 log(1/δ)δ 2 + √ 2 2 (κγ + γ) log(1/δ)δ 2 = 1 + √ 2 2 + √ 2 2 κ γ log(1/δ)δ 2 ≤ α, provided κ ≥ 1+ √ 2 2- √ 2 .
This indeed contradicts (4.4). We conclude that (4.13) holds.

Second, we show that λ ∈ Z 1 by showing that the test (2.7) is satisfied. By (4.10), and since |ζ -λ| ∞ ≤ δ/2, we have

|F λ (0)| ≤ F λ (ζ -λ)e -1 2 |ζ-λ| 2 -(F λ (0) + τ λ (ζ -λ)) + |τ λ ||ζ -λ| ≤ 2γ log(1/δ)δ 2 + √ 2 2 |τ λ |δ. (4.14) Choosing κ ≥ 8 3-2 √
2 , (4.14) and (4.13) further imply

|F λ (0)| ≤ 3 -2 √ 2 4 + √ 2 2 |τ λ | δ = 3 4 |τ λ | δ. (4.15)
Hence,

η λ = 3 4 |τ λ |δ. (4.16)
Let µ ∈ Λ be an arbitrary lattice point with |µ -λ| ∞ = 2δ. By (4.11),

|F λ (µ -λ)|e -1 2 |µ-λ| 2 = |F λ (µ -λ)e -1 2 |µ-λ| 2 -F λ (ζ -λ)e -1 2 |ζ-λ| 2 | = F λ (µ -λ)e -1 2 |µ-λ| 2 -(F λ (0) + τ λ (µ -λ)) -F λ (ζ -λ)e -1 2 |ζ-λ| 2 + (F λ (0) + τ λ (ζ -λ)) + τ λ (µ -ζ) ≥ |τ λ ||µ -ζ| -(11 + 2)γ log(1/δ)δ 2 ≥ |τ λ ||µ -ζ| ∞ -13γ log(1/δ)δ 2 ≥ |τ λ |(|µ -λ| ∞ -|ζ -λ| ∞ ) -13γ log(1/δ)δ 2 ≥ 3 2 |τ λ |δ -13γ log(1/δ)δ 2 .
Together with (4.14), this implies

|F λ (µ -λ)|e -1 2 |µ-λ| 2 ≥ |F λ (0)| + 3 - √ 2 2 |τ λ |δ -15γ log(1/δ)δ 2 .
Finally, we use (4.13) to analyze the obtained comparison margin against (4.16):

3 - √ 2 2 |τ λ |δ -15γ log(1/δ)δ 2 ≥ 3 4 |τ λ |δ + 3 -2 √ 2 4 κγ log(1/δ)δ 2 -15γ log(1/δ)δ 2 = 3 4 |τ λ |δ + 3 -2 √ 2 4 κ -15 γ log(1/δ)δ 2 ≥ η λ ,
where we fixed the value of κ so that 3-2

√ 2 4
κ -15 ≥ 0. Therefore, the point λ passes the selection test (2.7) (as formulated in (3.4)), i.e., λ ∈ Z 1 , as claimed. 4.5. After the selection step, each computed zero is close to a true zero. We show that

Z 1 ⊆ {F = 0} + Q 2δ (0). (4.17)
Let λ ∈ Z 1 be a computed zero, and let us find a zero z of F with |λ -z| ∞ ≤ 2δ. In terms of the Fock shift F λ , the success of the test (2.7) reads,

|F λ (µ)|e -1 2 |µ| 2 ≥ |F λ (0)| + η λ , for all µ ∈ Λ such that |µ| ∞ = 2δ; (4.18)
see (3.4). For an arbitrary z ∈ C with |z| ∞ = 2δ, we can find a lattice point µ ∈ Λ with |µ| ∞ = 2δ such that |z -µ| = |z -µ| ∞ ≤ δ/2. Hence, by (4.11),

|F λ (z)e -1 2 |z| 2 -F λ (µ)e -1 2 |µ| 2 | ≤ F λ (z)e -1 2 |z| 2 -F λ (0) + τ λ z + F λ (µ)e -1 2 |µ| 2 -F λ (0) + τ λ µ + |τ λ | |z -µ| ≤ 1 2 |τ λ |δ + 22γ • log(1/δ)δ 2 .
By (4.4) and (4.9), either

|τ λ |δ ≥ κγ log(1/δ) • δ 2 or |F λ (0)| ≥ κγ log(1/δ) • δ 2 ≥ |τ λ |δ. Choosing κ ≥ 88 ensures in the first case that |F λ (z)e -1 2 |z| 2 -F λ (µ)e -1 2 |µ| 2 | ≤ 3 4 |τ λ |δ ≤ η λ
and in the second case

|F λ (z)e -1 2 |z| 2 -F λ (µ)e -1 2 |µ| 2 | ≤ 3 4 |F λ (0)| ≤ η λ .
Combining this with (4.18), we conclude that By (4.12), there exists µ ∈ Z 1 such that |ζ -µ| ∞ ≤ δ/2. By (4.20), Z ⊊ Z ∪ {µ}. We claim that Z ∪ {µ} is 5δ-separated. For this, it suffices to check that

|F λ (z)|e -1 2 |z| 2 ≥ |F λ (0)|, for all z ∈ C such that |z| ∞ = 2δ. (4.19) By Lemma 3.1, there exists w λ ∈ C with |w λ | ≤ 2δ such that F λ (w λ ) = 0. This means that z λ := w λ + λ is a zero of F that satisfies |z λ -λ| ∞ ≤ 2δ, as desired.
|µ -λ| ∞ > 4δ, λ ∈ Z.
If λ ∈ Z, by (4.17), there exist

ζ ′ ∈ {F = 0} such that |ζ ′ -λ| ∞ ≤ 2δ. If ζ ′ = ζ, then |ζ -λ| ∞ ≤ 2δ, contradicting (4.20). Thus ζ ̸ = ζ ′ , while, ζ ′ ∈ Z+Q 2δ ⊂ Ω L+2δ .
Hence, we use (4.5) to conclude that

|µ -λ| ∞ ≥ |ζ -ζ ′ | ∞ -|µ -ζ| ∞ -|λ -ζ ′ | ∞ ≥ 7δ -δ/2 -2δ > 4δ.
Thus, the set is 5δ-separated:

inf |λ -λ ′ | ∞ : λ, λ ′ ∈ Z ∪ {µ}, λ ̸ = λ ′ ≥ 5δ,
contradicting the maximality of Z. It follows that a point λ ∈ Z such that |ζ -λ| ∞ ≤ 2δ must exist. We choose any such point, and define Φ(ζ) = λ.

4.7.

Verification of the properties of Φ. By construction, the map Φ satisfies (2.15). We now show the remaining properties. To show that Φ is injective, assume that Φ(ζ) = Φ(ζ ′ ). Then, by (2.15),

|ζ -ζ ′ | ∞ ≤ |Φ(ζ) -ζ| ∞ + |Φ(ζ ′ ) -ζ ′ | ∞ ≤ 4δ.
Hence, by (4.5), we must have

ζ = ζ ′ .
Finally, assume that λ ∈ Z ∩ Ω L-2δ and use (4.17) to select a zero ζ ∈ {F = 0} such that |ζ -λ| ∞ ≤ 2δ. Then ζ ∈ Ω L , and, by (2.15),

|Φ(ζ) -λ| ∞ ≤ |Φ(ζ) -ζ| ∞ + |ζ -λ| ∞ ≤ 4δ. (4.21)
As λ, Φ(ζ) ∈ Z and Z is 5δ-separated (see (2.8)) we conclude that λ = Φ(ζ), as claimed.

This concludes the proof of Theorem 2.1. □

Numerical Experiments

In this section we perform a series of tests of the AMN algorithm and compare its performance with MGN and thresholding supplemented with a sieving step (ST). 5.1. Simulation. We first discuss how to simulate samples from the input model (2.9). To make simulations tractable, we introduce a fast method to draw samples of the Gaussian entire function F 0 given by (2.10) on the finite grid (1.6). The method is based on the relation between the Bargmann transform and the shorttime Fourier transform (1.2) and amounts to discretizing the underlying signal f . We fix L > 0, T > 0, and δ > 0. For convenience, we further let σ = 1 and assume that T δ -1 is an integer. Recall that we also assumed that Lδ -1 is an integer.

To model a discretization of N , we take i.i.d. noise samples in the interval [-T -L, T +L] ⊆ R spaced by a distance δ. More specifically, we consider a random vector w = (w -(T +L)δ -1 , . . . , w (T +L)δ -1 ), where the elements w s ∼ N C (0, δ π 2 ) are independent, i.e., E[w s w s ] = δ π 2 , and E[w s w s ′ ] = 0. Here, w s can be interpreted as an integration of N over the interval [δs, δ(s + 1)].

Let φ = g| [-T,T ] be the restriction of g(t) = ( 2π )

1 2 e -t 2 to the compact support [-T, T ] and define

H(k + ij) := T δ -1 +k s=-T δ -1 +k w s + δf 1 (δs) φ(δ(s -k))e -2isjδ 2 , (5.1)
for k, j ∈ {-Lδ -1 , . . . , Lδ -1 }. The mean of H is given by

E[ H(k + ij)] = δ T δ -1 +k s=-T δ -1 +k f 1 (δs)φ (δ(s -k))e -2isjδ 2 , k, j ∈ {-Lδ -1 , . . . , Lδ -1 },
and approximates the integral

∞ -∞ f 1 (t)g(t -
x)e -2iyt dt = e -ixy e -1 2 (x 2 +y 2 ) F 1 (z), with x = δk and y = δj. Furthermore, the covariance of H is

Cov H(k + ij), H(k ′ + ij ′ ) = δ π 2 T δ -1 +k ′ s=-T δ -1 +k ′ φ (δ(s -k ′ )) φ (δ(s -k))e -2is(j-j ′ )δ 2 .
For small δ and sufficiently large T , this is an approximation of the integral

π 2 ∞ -∞
g t -δk ′ g t -δk e -2it(j-j ′ )δ dt = e -u 2 +v 2 +x 2 +y 2 2 e i(uv-xy) e (x-iy)(u+iv) ,

(5.2) with x = δk, y = δj, u = δk ′ , and v = δj ′ . Therefore, if we take T large enough so that we can ignore the numerical error introduced by the truncation of the normalized Gaussian window g, we obtain in (5.1) a random Gaussian vector whose covariance structure approximates (2.10) on a grid Λ L , provided that δ is small.

By conjugating z and multiplying by the deterministic factor e -ixy , we thus obtain an approximate sampling of (2.9) with weight e -1 2 |z| 2 :

(5.

3) e -1 2 |z| 2 F (z) ≈ e -ixy H(z)

for z = δk + iδj. We carry out all computations with the weighted function (5.3), as the unweighted version can lead to floating point arithmetic problems. Note that, for a grid point λ, the comparison margin (2.6) can be expressed in terms of

e -1 2 | • | 2 F ( • ) as η λ = max e -1 2 |λ| 2 |F (λ)| , 3 4 e 1 2 δ(2i Im(λ)+δ) e -1 2 |λ+δ| 2 F (λ + δ) -e -1 2 |λ| 2 F (λ) .
5.2. Specifications for the experiments.

5.2.1.

Implementation of the sieving step in AMN. In order to fully specify the AMN algorithm we need to fix an implementation of the sieving step, which provides a subset Z ⊆ Z 1 satisfying (2.8), and such that no proper superset Z 1 ⊇ Z ⊋ Z satisfies (2.8). We choose an implementation that uses knowledge of the input F to decide which points are to be discarded. We assume that Z 1 is non-empty, otherwise Z is trivial.

Algorithm S1: Obtain a maximal subset that is 5δ separated.

Input: Values of a function F on a grid Λ L . A discrete non-empty set Z 1 ⊆ Λ L .

Step 1: Copy the set Z 1 to Z aux 1 .

Step 2: Consider the (pre)ordered set (Z aux 1 , ≼), where

λ ≼ µ ⇐⇒ e -1 2 |λ| 2 |F (λ)| ≤ e -1 2 |µ| 2 |F (µ)| , λ, µ ∈ Λ L . (5.4)
Step 3: Choose a minimal point λ ∈ (Z aux 1 , ≼).

Step 4: Add λ to Z.

Step 5: Remove all µ ∈ Z aux 1 such that (5.5) 0 ≤ |λ -µ| ∞ ≤ 4δ.

Step 6: If the set Z aux 1 is not empty, repeat Steps 3-6. If the set Z aux 1 is empty, then the algorithm ends.

Output: The set Z.

The resulting set Z ⊆ Z 1 always satisfies (2.8). Moreover, any superset Z ⊋ Z included in Z 1 must contain some of the discarded points µ ∈ Z 1 , which by construction satisfy (5.5) for some λ ∈ Z, and therefore Z is not 5δ separated. Thus, Z is indeed maximal with respect to (2.8).

The choice of λ ∈ Z aux 1 in Step 3 of S1 is not essential. Our particular choice is motivated by finding the zeros of F ; however, we did not observe any significant performance difference when using other algorithms than S1 as the sieving step of AMN.

Specification of the compared algorithms.

Given the values of a function F : C → C on the grid Λ L , we consider the following three algorithms to compute an approximation of {F = 0} ∩ Λ L-1 .

• AMN: the AMN algorithm run with domain length L -1 and with sieving step S1 implemented as described in Section 5.2.1,

• MGN: outputs the set of all grid points λ ∈ Λ L-1 such that

(5.6) e -1 2 |λ| 2 |F (λ)| ≤ e -1 2 |µ| 2 |F (µ)|, |λ -µ| ∞ = δ.
• ST: outputs the set of grid points λ ∈ Λ L-1 obtained as the result of applying the sieving algorithm S1 to

λ ∈ Λ L-1 : e -1 2 |λ| 2 |F (λ)| ≤ 2δ .
Note that each of the algorithms relies only on the samples of F on Λ L+2δ-1 . The use of a common input grid Λ L simplifies the notation when considering various grid spacing parameters δ.

5.2.3.

Varying the grid resolution. In the numerical experiments, we start with a small minimal spacing value δ = δ Hi , that provides a high resolution approximation in (5.1), and simulate F as in Section 5.1. We then incrementally double δ to produce coarser grid resolutions and subsample F accordingly. More precisely, each element of the grid Λ L can be written as

λ k,l = (-L + kδ) + i(-L + lδ), 0 ≤ k ≤ M , 0 ≤ l ≤ N, (5.7)
for adequate M , N > 0. If F is given on Λ L , we subsample it by setting S(F )(λ k,l ) := F λ 2k+i2l , (5.8) for values (k, l) such that the indices 2k + i2l are valid. 5.3. Faithfulness of simulation of zero sets. As a first test, we simulate random inputs from the model (2.9), as specified in Section 5.1, apply the abovedescribed three different algorithms, and test whether this process faithfully simulates the zero sets of the random function (2.9). To this end, we estimate first or second order statistics on the computed zero sets by averaging over several realizations of (2.9), and compare them to the corresponding expected values concerning the zero sets of (2.9).

5.3.1.

No deterministic signal. We first consider the case F 1 ≡ 0 and σ = 1 in (2.9). Let F δ Hi 1 , . . . , F δ Hi R be R independent realizations of samples of (2.9) on a grid Λ L with resolution δ = δ Hi , simulated as in Section 5.1. These are then subsampled with (5.8) yielding F δ k r = S (k) (F δ Hi r ) and used as input for AMN, MGN, and ST, as specified in Section 5.2.2. The corresponding output sets are denoted Z δ r where we omit the dependence on the method to simplify the notation. These sets should approximately correspond to {F r = 0} ∩ Λ L-1 , for R independent realizations of (2.9). We now put that statement to test.

[|{F = 0} ∩ Θ|] = Θ 1 π dm(ζ) = |Θ| π , (5.9)
see, e.g., [START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF]Section 2.4]. We define the following empirical estimator for the first intensity ρ 1 = 1/π:

ρ(Θ, r, δ) = | Z δ r ∩ Θ| |Θ| . (5.10)
If the computed set Z δ r is replaced by {F = 0} in (5.10), the estimator would be unbiased. The mean of the estimation error ρ(Θ, r, δ) -1/π thus measures the quality of the algorithm used to compute Z δ r , as it should be close to zero when the algorithm is faithful. In Table 1, we present the empirical means and the empirical standard deviations of the estimation error over R = 1000 independent realizations F δ r for L = 7, Θ = Ω L-1 , T = 6, and various grid sizes δ. To derive a benchmark for the empirical standard deviation of ρ(Θ, r, δ) -1/π, we express the variance of |{F = 0} ∩ Θ|/|Θ| in terms of the second intensity function Table 2. Functions f 1 and their Bargmann transforms F 1 = B(f ).

f 1 F 1 f 1 (t) = e -t 2 F 1 (ζ) = 1 f 1 (t) = e -t 2 2t F 1 (ζ) = ζ
A formula for ρ 2 (ζ, ζ ′ ) is provided in [START_REF] Hannay | Chaotic analytic zero points: exact statistics for those of a random spin state[END_REF] and numerical integration over Θ = Ω L-1 results in Var[|{F = 0} ∩ Θ|/|Θ|] ≈ 0.01165. We see in Table 1 that the methods AMN and MGN almost perfectly match the expected mean and standard deviation while ST does not.

5.3.2. Deterministic signal plus noise. We now consider the input model (2.9) with F 1 ̸ = 0 and σ = 1. We choose F 1 from Table 2 and rescale it so that A = sup ζ∈C e - We only test first order statistics of the computed zero sets. The benchmark is provided by Proposition 3.4: the expected number of zeros of F in Θ is (5.12)

E[|{F = 0} ∩ Θ|] = Θ ρ 1 (ζ) dm(ζ),
where ρ 1 is given by (3.10) (with σ = 1). For each of the tested algorithms, we define an estimator for the error resulting from replacing {F = 0} in (5.12) by the computed set Z δ r (for 1 ≤ r ≤ R):

(5.13)

β(Θ, r, δ) = | Z δ r ∩ Θ| -Θ ρ 1 (ζ) dm(ζ) |Θ| .
As before, we simulate R = 100 realizations of F = F 0 + F 1 on a grid with a certain spacing δ. The empirical average of β(Θ, r, δ) over all realizations is denoted β R (Θ, δ). As ρ 1 is not constant when F 1 ̸ = 0, this time we calculate β R (Θ, δ) on Θ = Ω L 1 for several values of L 1 .

The results for δ = 2 -9 are depicted in Figure 5. We see that the performance of AMN and MGN is indistinguishable, while ST may perform poorly even at such high resolution. Lower grid resolutions yield similar results. 5.4. Failure probabilities and consistency as resolution decreases. Having tested the statistical properties of the computed zero sets under the input model (2.9) we now look into the accuracy of the computation for an individual realization F . We aim to test the existence of a map as in Theorem 2.1, that assigns true zeros to computed ones with small distortion and almost bijectively. As a proxy for the (unavailable) ground truth {F = 0} we will use the output of AMN from data at very high resolution (computations with MGN yield indistinguishable results). We thus conduct a consistency experiment, where the zero set of the same realization of F is computed from samples on grids of different resolution, and the existence of a map as in Theorem 2.1 between both outputs is put to test.

Suppose that samples of a function F are simulated on a high-resolution grid Λ L with spacing δ = δ Hi and restricted to the low-resolution grid Λ L with spacing δ = δ Lo by subsampling. We compute Z δ Hi ⊆ Ω L-1 from the high-resolution data using AMN, and Z δ Lo ⊆ Ω L-1 from the low-resolution data, using one of the algorithms described in Section 5.2.2.

Second we construct a set U ⊆ Z δ Hi and a map ϕ : U → Z δ Lo r with the following greedy procedure:

π 1 2

 1 e -t 2 :e -ixy e -1 2 (x 2 +y 2 ) F (x -iy) = R f (t)g(t -x)e 2ity dt, x, y ∈ R. (1.2)

Figure 1 .

 1 Figure 1. Calculation of zero sets by thresholding: Values below the same threshold (marked by circles) fail to detect some zeros and at the same time cannot clearly separate other zeros.

  (a) One of the immediate neighbors is used to calculate a comparison margin. (b) The weighted values of F are compared against points in the larger box.

Figure 2 .

 2 Figure 2. The selection step of AMN.

  (a) A point passes the selection test. (b) The same zero causes a second detection.

Figure 3 .

 3 Figure 3. Sliding the test box.

Figure 4 .

 4 Figure 4. A realization of e -1 2 |z| 2 |(F 0 (z) + F 1 (z))| with F 1 the Bargmann transform of f 1 . The deterministic functions are scaled to obtain the prescribed A. Zeros computed with AMN, MGN, and ST are calculated from grid samples with δ = 2 -9 . Zeros from AMN and MGN coincide (circle), while ST (cross) fails either by detecting false zeros (left) or by not capturing all of them (right).

3. 3 . 1 2 |z| 2 |F

 312 Minimum principle for weighted amplitudes. The following weighted version of the minimum principle is at the core of the success of MGN and AMN. Lemma 3.1. Let F : C → C be entire, r > 0, and assume that |F (0)| ≤ |F (z)|e -1 2 |z| 2 , for all z ∈ C such that |z| ∞ = r. (3.5) Then there exists z ∈ C with |z| ∞ ≤ r such that F (z) = 0. Proof. Let D := {z ∈ C : |z| ∞ < r} and suppose that F does not vanish on D. Then the function H(z) = e (z)| is well defined on D and satisfies log H(0) ≥ log H(z), z ∈ ∂D. (3.6) By the analyticity of F , ∆ log |F | = 0 and thus

  1.2], the conditional expectation of |∂ ζ F (ζ)| 2 given F (ζ) = 0 is the same as the expectation of |W | 2where W is a shifted circularly symmetric complex Gaussian random variable with mean ∂ ζ F 1 (ζ) -ζF 1 (ζ) and variance e |ζ| 2 . Thus,(3.13)

  we can select a lattice point λ ∈ Λ L+2δ such that 0 < |λ -ζ| ≤ δ. We now use repeatedly (4.2) and (4.3). First, we use (4.2) with w = ζ and z = ζ ′ -ζ, and note that F ζ (ζ ′ -ζ) = 0 and

4. 6 .

 6 Definition of the map Φ. We now look into the sieving step of the AMN algorithm and analyze the final output set Z. Given ζ ∈ {F = 0}∩Ω L we claim that there exists λ ∈ Z such that |ζ -λ| ∞ ≤ 2δ. Suppose to the contrary that |ζ -λ| ∞ > 2δ, λ ∈ Z. (4.20)

2 ( 5

 25 ρ 2 (ζ, ζ ′ ) of {F = 0} as follows:E |{F = 0} ∩ Θ| -|Θ| π .11) = E |{F = 0} ∩ Θ| • (|{F = 0} ∩ Θ| -1) -|Θ| 2 π 2 + |Θ| π = Θ Θ ρ 2 (ζ, ζ ′ ) dm(ζ) dm(ζ ′ ) -|Θ| 2 π 2 + |Θ| π .

1 2 |ζ| 2

 2 |F 1 (ζ)| holds for the signal intensities A = 1 and 100.

1 , 2 - 9 ) f 1 1 , 2 - 9 ) f 1 Figure 5 .

 129112915 Figure 5. Empirical mean of β(Θ, r, δ) for different choices of f 1 and A, increasing domain Θ = Ω L 1 for L 1 < L, and the three methods. Note the different scale in the bottom right plot illustrating a systematic error in the ST method.

  This contradiction shows that F must vanish on D. □3.4. Linearization.In what follows, we derive basic facts about the input model (2.9), and always assume that (2.10) and (2.11) hold.

The following is a strengthened version of

[START_REF] Hough | Zeros of Gaussian analytic functions and determinantal point processes[END_REF] Lemma 2.4.4]

.

Table 1 .

 1 Empirical means ± standard deviations of the estimation errors ρ(Θ, r, δ) -1/π for Θ = Ω L-1 , L = 7, and 1000 independent realizations. Benchmark values for a faithful computation are 0 for the mean and 0.01165 for the standard deviation.

	δ	AMN	MGN	ST
	2 -4 -0.00120 ± 0.01171 -0.00048 ± 0.01150 +0.01868 ± 0.02858
	2 -5 -0.00062 ± 0.01164 -0.00057 ± 0.01162 +0.02189 ± 0.04047
	2 -6 -0.00065 ± 0.01156 -0.00064 ± 0.01155 +0.02280 ± 0.05391
	2 -7 -0.00068 ± 0.01153 -0.00068 ± 0.01153 +0.02354 ± 0.06774
	2 -8 -0.00062 ± 0.01155 -0.00062 ± 0.01155 +0.02424 ± 0.07429
	2 -9 -0.00067 ± 0.01158 -0.00067 ± 0.01158 +0.02390 ± 0.07237
	The expected number of zeros of the random function F on a Borel set Θ ⊆ C
	is
		E

and J. L. R. gratefully acknowledge support from the Austrian Science Fund (FWF): Y 1199 and P 29462. N. F. gratefully acknowledges support from the Israel Science Foundation (ISF) grant no. 1327/19.

Construction of U and ϕ

Input: Two subsets of Ω L-1 : Z δ Hi and Z δ Lo .

Step 1: Choose a total order on Z δ Hi . Let U and U ′ be empty sets. If Z δ Hi is empty, output U = ∅ and ϕ = ∅. Otherwise proceed to Step 2.

Step 2: Let λ be the first element of Z δ Hi \ (U ∪ U ′ ).

Step 3:

If Φ(λ) is non-empty, add λ to U , and choose ϕ(λ) ∈ Φ(λ) such that

Step 4:

Output: The set U and the map ϕ.

The resulting function ϕ is injective and satisfies

We say that the computation of Z δ Lo r was certified to be accurate if

∩ Ω (L-1)-2δ Lo ⊆ ϕ(U ). (5.14) In this case, the map ϕ satisfies properties analogous to the ones in Theorem 2.1. Conceivably, other such maps may exist even if the one constructed in the greedy fashion fails to satisfy (5.14). We define the following computation certificate:

The experiment to estimate failure probabilities as a function of the grid resolutions is fully specified as follows. We consider the input model (2.9) with σ = 1. We choose F 1 from Table 2 and rescale it so that A = sup ζ∈C e -1 2 |ζ| 2 |F 1 (ζ)| holds for the signal intensities A = 1 and 100. We fix L > 0 and δ Hi > 0 and let F δ Hi 1 , . . . , F δ Hi R be R independent realizations of samples of (2.9) on a grid Λ L with resolution δ = δ Hi , simulated as in Section 5.1. These are then subsampled j times with (5.8) yielding F δ k r = S (k) (F δ Hi r ), 1 ≤ k ≤ j. We use AMN with input F δ Hi r to obtain a set Z δ Hi r . Further, for each 1 ≤ k ≤ j, we use each of the algorithms M = AMN, MGN, or ST with input F δ k r to obtain sets Z δ k r,M . Finally, we compute all the certificates M( Z δ Hi r , Z δ k r,M ) and average them over all realizations to obtain the following estimated upper bound for the failure probability of the method M with grid spacing δ = δ k :

(5.15)

We present in Table 3 values obtained for p(δ k , M ) for a resolution starting as high as δ Hi = 2 -9 , with a truncation of the window g at T = 6, in the target domain Ω L-1 for L = 7, and R = 1000 realizations of a zero-mean F . We also present the results for F 1 as in Table 2, rescaled to achieve a signal intensity A = 1 or A = 100. We see that both AMN and MGN deliver very low failure probabilities (with MGN slightly outperforming AMN at lower resolutions). In contrast, ST delivers large failure probabilities even at high resolution.

Table 3. Estimation of the failure probability p(δ k , M ) in the sense of Theorem 2.1, in the domain Ω L-1 with parameters δ Hi = 2 -9 , T = 6, and L = 7. Averages are computed over R = 1000 and R = 100 realizations for the pure noise and signal f 1 plus noise cases, respectively.

2 -4 0.082 0.001 0.665 0.07 0.00 0.67 0.13 0.00 0.87 0.07 0.00 0.64 0.18 0.00 1.00 2 -5 0.007 0.000 0.536 0.00 0.00 0.50 0.00 0.00 0.75 0.00 0.00 0.52 0.01 0.00 1.00 2 -6 0.001 0.000 0.419 0.00 0.00 0.41 0.00 0.00 0.70 0.00 0.00 0.42 0.00 0.00 1.00 2 -7 0.000 0.000 0.389 0.00 0.00 0.32 0.00 0.00 0.72 0.00 0.00 0.34 0.00 0.00 1.00 2 -8 0.000 0.000 0.369 0.00 0.00 0.29 0.00 0.00 0.65 0.00 0.00 0.33 0.00 0.00 1.00 2 -9 0.000 0.000 0.359 0.00 0.00 0.31 0.00 0.00 0.71 0.00 0.00 0.32 0.00 0.00 1.00