Luc Forget

Gauthier Harnisch

Ronan Keryell

Amd Adaptive

Florent De Dinechin

A single-source C++20 HLS flow for function evaluation on FPGA and beyond

This paper presents a framework to reuse the intelligence of RTL generators in a single-source HLS setting. This framework is illustrated by a C++ fixed-point library to generate mathematical function evaluator. A compiler flow from C++20 to Vivado IPs has been developed to make the library usable with Vitis HLS. This flow is demonstrated on two applications: an adder for the logarithmic number system, and additive sound synthesis. These experiments show that the approach allows to easily tune the precision of the types used in the application. They also demonstrate the ability to generate arbitrary function evaluator at the required precision.

INTRODUCTION

High Level Synthesis (HLS) has been able to improve FPGA programming productivity compared to traditional register transfer level (RTL) flows. This higher productivity is obtained by describing hardware algorithm with software programming languages that offers high level abstractions. Compiling as-is a program targeting CPU with an HLS compiler can however lead to suboptimal results. Indeed, FPGAs expose much more parallelism and flexibility than CPUs, and a CPU program must often be adapted to exploit this potential. In order to use the parallelism, the algorithm implementation can be pipelined, loops can be unrolled, and buffer elements can be streamed in a dataflow fashion.

Most of the recent HLS tools offer language extensions to this purpose.

A second way to exploit the flexibility of FPGAs is to tailor the numeric formats used in the computation to fit the application accuracy requirement. Not every application has its optimal numeric formats covered by the few standard formats supported by usual software programming languages. Recent work on machine learning or signal processing applications have for instance highlighted the interest for reduced-precision floating point formats [START_REF] Mueller | Dlfloat: A 16-b floating point format designed for deep learning training and inference[END_REF][START_REF] Burgess | Bfloat16 processing for Neural Networks[END_REF][START_REF] Bita | Pushing the Limits of Narrow Precision Inferencing at Cloud Scale with Microsoft Floating Point[END_REF][START_REF] Henry | Leveraging the Bfloat16 Artificial Intelligence Datatype For Higher-Precision Computations[END_REF][START_REF]NVIDIA A100 Tensor Core GPU Architecture[END_REF], or even radically different numerical representations [START_REF] Syed Asad Alam | Low-precision Logarithmic Number Systems: Beyond Base-2[END_REF][START_REF] Johnson | Rethinking floating point for deep learning[END_REF][START_REF] Lu | Evaluations on Deep Neural Networks Training Using Posit Number System[END_REF]. To address this, most vendor HLS tools provides support for custom-sized fixed point formats, and external HLS libraries offer basic support for custom-sized floating points or even more exotic formats [START_REF] Thomas | Templatised Soft Floating-Point for High-Level Synthesis[END_REF][START_REF] Uguen | Evaluating the Hardware Cost of the Posit Number System[END_REF]. Such HLS libraries can benefit from all the optimizations performed by the HLS tools [START_REF] Thomas | Templatised Soft Floating-Point for High-Level Synthesis[END_REF][START_REF] Uguen | Application-Specific Arithmetic in High-Level Synthesis Tools[END_REF]. However, these libraries only provide support for basic arithmetic operations. They usually miss, for instance, the elementary functions (exponential, logarithm, trigonometric) available only for floating-point formats.

Indeed, a third way to exploit the FPGA flexibility is to build ad-hoc architectures for numerical functions appearing in application code, at the precision required by the application. These functions may be elementary functions, but also composite ones such as log 2 (1 + 2 ) or sin( + ) .

Publication rights licensed to ACM. ACM acknowledges that this contribution was authored or co-authored by an employee, contractor or affiliate of a national government. As such, the Government retains a nonexclusive, royalty-free right to publish or reproduce this article, or to allow others to do so, for Government purposes only.

Tools exist [START_REF] Florent De Dinechin | Automatic generation of polynomial-based hardware architectures for function evaluation[END_REF][START_REF] Hsiao | Table Size Reduction Methods for Faithfully Rounded Lookup-Table-Based Multiplierless Function Evaluation[END_REF][START_REF] Lee | Hierarchical Segmentation Schemes for Function Evaluation[END_REF][START_REF] Thomas | A general-purpose method for faithfully rounded floating-point function approximation in FPGAs[END_REF] that can generate an evaluator for an arbitrary function, usually provided as an arbitrary mathematical expression. These tools are complex program that analyze the function in order to use various approximation techniques before producing RTL code. This generated RTL IP may be used in an HLS tool, but it comes with some drawbacks:

• There is a loss in productivity and portability of the solution, since inclusion of RTL IPs in an HLS design is a tool-specific process and requires extra plumbing; • HLS designs are pipelined at the HLS stage, which manipulates high level representations. HLS tools are not able to blend the RTL IPs in this pipeline process. Due to this, the RTL IPs behaves like optimization barriers for the HLS stage. Instead of optimizing the whole component, HLS tools will only be able to optimize separately the sub-components above and below the black box, which may lead to suboptimal synthesized designs.

Filling this gap requires to fully integrate in the HLS flow the complex process of generating an optimized architecture from an arbitrary function specification. This paper presents an HLS C++ library that is a first step towards this goal. It provides a single-source way to generate an evaluator for a function given as a mathematical expression, on custom-sized fixed-point inputs. This library is described in section 2.

Section 3 presents the compiler flow, which involves two compilation stages, and its interaction with external tools that perform the mathematical approximation of the function. In its current state, the library offers plain tabulation and bipartite function approximation [START_REF] Sarma | Faithful bipartite ROM reciprocal tables[END_REF][START_REF] Sunderland | CMOS/SOS frequency synthesizer LSI circuit for spread spectrum communications[END_REF] introduced in more details in Section 5.

In order to use modern C++ constructs on AMD FPGAs, an experimental compiler targeting Vitis IP has been developed from clang++. This is necessary because AMD Vitis toolchain only supports C++14. This compiler is described in section 4.

A few applications have been developed in order to demonstrate the library capabilities. These applications are described in section 6.

Finally, improvement possibilities are discussed in section 7.

LIBRARY PRESENTATION

The library1 allows the user to define custom fixed point formats and computations on value of these formats, as well as conversions between these formats and standard C++ types.

In all this paper, a fixed point format is defined by the weights msb_weight and lsb_weight associated respectively to its most-significant and least-significant bits, and its signedness signedness. Signed values are represented in 2's complement. The width of a format (in bits) is width = msb_weight -lsb_weight + 1. The representable values in such a format are scaled integers of the form × 2 lsb_weight , where is an integer such that

∈ 0, 2 width -1 if the format is unsigned, or ∈ -2 width-1 , 2 width-1 -1 if the format is signed.

C++ types for fixed-point number

In the library, a format is represented by the template FixedFormat. For instance, the type FixedFormat<31, 0, signed> represents a fixed-point format which is equivalent to a 32 bit signed integer.

A FixedFormat is used to parameterize the template class FixedNumber which represents a value of a given format.

using format = FixedFormat<4, -3, unsigned>; FixedNumber<format> X{0b00000010};

Listing 1. Construction of a FixedNumber from its representation.

A FixedNumber can be constructed from the representation of its value. For instance, in Listing 1, the variable X is initialized with the value 2 -2 because only its second bit (of weight -2) is set.

A FixedNumber may also be constructed from the value of a standard arithmetic type, thanks to the static method get_from_value. This conversion returns the value of the target format closest to the source value. In case of a tie, the greatest of the two possible values is returned.

In addition, a user-defined literal (UDL) allows to create a FixedNumber from a value without having to specify the format. The narrower format that can represent exactly the expression will be automatically deduced from the operator.

As decimal floating-point literals may have non-finite representation (for instance 0.1 10 = 0.000110011001100... 2 ), this UDL only accepts hexadecimal floating-point literals.

Finally, conversion from FixedNumber to standard arithmetic types are expressed with the get_as template method from FixedNumber. This method takes a target arithmetic type as template parameter and returns the value of this type which is the closest to the represented value (with the same tie-breaking rule as the other conversions). Listing 2 presents an example of computation performed with the library demonstrating many of its capabilities.

Classical arithmetic computations

The proposed library defines standard arithmetic operations between FixedNumber of different FixedFormat, in a way very similar to mainstream HLS fixed point libraries such as ap_fixed or ac_fixed.

Addition, subtraction and multiplication are computed exactly, the FixedFormat of the result being deduced by the library to be wide enough to hold the result. Rounding, overflow or saturation must be managed explicitely:

FixedNumber provides additional methods to extract a part of it, round to a specific precision, or expand the number to a wider format.

using radius_fmt = FixedFormat<3, -4, unsigned>; auto circumference(FixedNumber<radius_format> radius) { // For simplicity, FixedNumber<m, l, s> replaces FixedNumber<FixedFormat<m,l,s>> // Generate a constant of value two constexpr auto two = 2._fixed; // two is the same as FixedNumber<1, 1, unsigned>{1}; // Generate an approximation of pi rounded to 2^(-4) constexpr auto pi = FixedNumber<1, -4, unsigned>::get_from_value(PI); // Diameter is a FixedNumber<4, -3, unsigned>; auto diameter = two * radius; // circumference is a FixedNumber<6, -7, unsigned> auto circumference = diameter * pi; // Eventually we get rid of low bits // Return type is FixedNumber<6, -4, unsigned> return circumference.round_to<-4>(); } The proper way to handle division and remainder is still under study.

The main novelty introduced by the proposed library is the presence of an evaluator generator for more sophisticated mathematical expression. The mechanism behind this is introduced in the next section.

Defining and evaluating a mathematical function of a FixedNumber argument

The library provides a standard set of operators and elementary functions. As this set is limited, applications often require the evaluation of mathematical functions which are themselves defined as compositions of operations and functions from this standard library, for instance ( ) = log(1 + ). Such a function can be evaluated at run-time by a composition of library components (here an exponential, an addition and a logarithm). It would however often be more efficient to directly compile this whole expression into a single component, as illustrated by Figure 1. For this purpose, the proposed library also provides a mechanism to describe a function of one variable, then generate an accuracy-constrained fixed-point evaluator for this function.

In the proposed approach, such a composite function may be defined as a C++ expression with a specific property: it has a free variable defined by the FreeVariable constructor. More precisely, this expresssion must be a tree of mathematical operations or elementary functions, and the leaves are either free variables or constants. A valid expression, for the purpose of building an evaluator, should have at most one free variable. Then the expression tree defines a function, and the FixedFormat associated to the free variable defines the domain of this function.

Listing 3 shows how a simple expression tree can be constructed. At line 3 of the listing, a FreeVariable node is created from a fixed-point variable. The FreeVariable gets associated to the format of the FixedNumber, so the input domain of the resulting function consists of the set of values representable on FixedFormat<5, -12, unsigned>, or [0, 2 5 -2 -12 ]. At line 4, a Constant leaf is created using a user-defined literal 2 . Similarly to what happen with the _fixed literal, the FixedFormat to associate to the constant will be deduced to fit exactly the constant storage requirement, but it is less relevant here as the constant will usually not directly appear in the final architecture. The library also has special constants to represent numbers like , which are kept symbolic (infinitely accurate) in the expression tree defining a function. The f variable at the end of the listing represents the function : → sin( + 1.5) Once the function is properly defined, it is possible to generate an accuracy constrained evaluator for it. For a given accuracy constraint > 0, an accuracy-constrained evaluator ˜ for a given function over a domain is a function such that ∀ ∈ , ( ) -˜ ( ) < This is done using the evaluate library function. This function is templated by an integer lsb_out. The result to this call is a fixed-point number having lsb_outfor lsb. The distance between the real value of the function and the evaluate function result is strictly lower than 2 lsb_out . This accuracy constraint is often called "faithful rounding" in the literature.

over the input domain [0, 2 5 -2 -12 ].
For instance, in listing 3, the variable functionOutput will contain an approximation of sin(1.5 + functionInput) faithful to 2 -17 . The output domain of the evaluator is FixedFormat<1, -17, signed>. This format is deduced by the library from the input domain and expression tree.

This evaluator can resort to plain tabulation for small precisions, as described in the next section, or to more sophisticated approximation methods including arbitrarily complex optimization techniques, an example of which is described in Section 5.

INTERACTION WITH EXTERNAL TOOLS IN A SINGLE-SOURCE SETTING

While building the evaluators in a completely metaprogrammed way would have been ideal to reach the goal of the library, it has two important drawbacks. First, the evaluation of C++ functions at compilation-time is way slower than running their compiled equivalent, which can be problematic as the exploration space for evaluator architectures is important. Secondly, depending on the evaluation technique used, it requires either to be able to compute accurately the function on all the inputs or to do some analysis of the expression to evaluate. This requires specific tools, that to the authors knowledge are not provided by any existing C++ header-only library. The development of such a library would be a huge task, therefore the presented library opted for a hybrid approach. The complete flow for using the library is described in figure 2. The source file, containing all the expression tree of the function and the associated calls to evaluate, has first to be compiled for CPU.

app
The produced CPU program links at runtime to an architecture generation library, which is responsible to determine an efficient hardware implementation for the problem that cannot be solved by metaprogrammatically. In the case of evaluator generation, it first builds a runtime representation of each expression tree for which an evaluation was required. The resulting runtime expression description is used to build a description understandable by the Sollya library [START_REF] Chevillard | Sollya: An Environment for the Development of Numerical Codes[END_REF]. Sollya is used to build the highly accurate reference table for the function. This table is then processed to find a good bi-partite decomposition. If no good decomposition is found, then the tool falls back to a plain tabulation.

The chosen algorithm and associated tables are then saved in the specialization header. With this specialization header, it is now possible to compile the component with an HLS compiler. However the library types and the specialization mechanism make use of multiple C++20 constructs. To our knowledge no vendor compiler supports C++20, therefore we had to design a custom compiler. This custom compiler interact with the Vitis toolchain to build Vitis IP. It is described in the next section.

COMPILER FOR RECENT C++ STANDARD HLS ON AMD FPGAS

This section describes the custom C++ compiler that has been developed for enabling C++20 constructs on HLS targeting AMD FPGAs. The Vitis frontend, which has recently been open sourced, consists in a derivative from Clang/LLVM [START_REF] Llvm | The LLVM Compiler Infrastructure[END_REF],

with added support for custom pragmas and some specific optimization pass. The LLVM version at which the tool is frozen is LLVM 7. Moving the Vitis frontend back from the past to align it with the current LLVM 15 seems to represent an enormous amount of work. Instead of doing so, the strategy has been to start from a modern Clang++ 15 compiler, and use its output to feed the Vitis frontend with LLVM IR. This modern clang++ is able to support C++ 20 and some of C23 constructs that the generator needs, such as concepts. However, the LLVM IR version produced by the recent Clang/LLVM is not directly compatible with Vitis for two reasons:

• first, Vitis understands only LLVM 7 IR, and LLVM 15 IR has some constructs invalid in LLVM 7 IR;

• secondly, Vitis HLS expects some IR patterns that are no longer emitted by clang++. Without them, the quality of result can degrade abruptly. For instance, the canonical way to move objects in recent LLVM is via a memcpy instructions, while Vitis works better with loads and stores.

So a few LLVM passes that downgrade modern LLVM IR to Vitis understandable IR has been developed. This compiler also offers a few custom C++ decorators that replace Vitis HLS pragmas. Most of this work is shared with the experimental porting of SYCL3 on AMD FPGAs [START_REF] Amd | SYCL for Vitis: Experimental fusion of triSYCL with Intel SYCL oneAPI DPC++ up-streaming effort into Clang/LLVM[END_REF]. To allow the usage of the tool outside the SYCL environment, new Vitis-IP target has been developed and added too this compiler. With this target, the downgrading passes are run and vitis_hls is invoked on the result, to produce an IP package that can be imported in a Vivado block design.

A command to launch an IP synthesis is then similar to a clang++ invocation such as clang++ --target=vitis_ip-xilinx vitis_ip.cpp -std=c++20 -o ip.zip --vitis-ip-part=xc7vx330t-ffg1157-1 ↩→

A NON-TRIVAL FUNCTION APPROXIMATION METHOD : THE BIPARTITE APPROXIMATION

The flow described in figure 2 introduces the opportunity to run arbitrarily complex programs to optimize the implementation of the function. To illustrate this opportunity, we implemented a classical but non-trivial function approximation technique: the bipartite table method [START_REF] Sarma | Faithful bipartite ROM reciprocal tables[END_REF][START_REF] Sunderland | CMOS/SOS frequency synthesizer LSI circuit for spread spectrum communications[END_REF]. As illustrated by Figure 3, this method is based on a piecewise linear approximation. The input domain is partitionned in 2 sub-intervals by splitting the input as = + 2 -, where (the interval index) consists of the leading bits of and (the index within one interval) consists of the least significant bits. The linear approximation on each subdomain indexed by is ( ) ≈ ( ) + 2 -× ( ) where ( ) is the slope of the approximation segment. The idea of the bipartite method is to tabulate the multiplications × ( ).

To reduce the cost of this tabulation, the slopes are shared between 2 consecutive intervals, so instead of × ( ) it is possible to tabulate × ( ) where consists of theleading bits of . Many evolutions to this method have been developed (see [START_REF] Hsiao | Table Size Reduction Methods for Faithfully Rounded Lookup-Table-Based Multiplierless Function Evaluation[END_REF] and references therein). Due to lack of time, the current tool flow only provides a basic bipartite approximation: its purpose is to demonstrate that the approach is not restricted to simple tabulation. The current implementation computes the parameters and , and fills the corresponding tables for a faithful approximator. Plain tabulation is used when this bipartite approximation fails (for instance when the function does not offer sufficient linearity).

// lns_t is a specialization of FixedNumber lns_t LNSAdder(lns_t op1, lns_t op2) { auto min_exp = min(op1, op2); auto max_exp = max(op1, op2); min_exp -= max_exp; auto diff_fv = FreeVariable{min_exp}; auto f = log2(0x1.p0_cst + pow(0x2.p0_cst, diff_fv)); auto rounded_f = f + lns_t::rounding_constant; auto result_exp_diff = evaluate<lns_t::add_eval_precision>(expr_rounded); max_exp += result_exp_diff; return max_exp; } Listing 6. A simplified LNS adder

APPLICATION STUDY

This section presents applications that illustrates the usage of the library. The first application, a LNS adder, consists mostly in one function evaluation call. The next one, an additive synthesizer, illustrates the usage of the the evaluator integrated in a more complex computation. It also demonstrate the compromises on QOR that are made possible by the combination of variable precision fixed-point and HLS. Finally, the LNS adder component is used to build a SYCL kernel, to illustrate the advantages of a single-source system.

Narrow LNS Adder

The Logarithm Number System (or LNS) is an alternative to floating-point. Instead of representing a number by an exponent and a fraction, only the exponent is kept, but it is kept with a fraction part. In other words, a number is represented by a sign bit and a fixed-point 2's complement base logarithm , coded in some FixedFormat < , , signed> format with < 0. In other words the real value represented by the and fields is = (-1) × .

(

The main advantage of LNS is the simplicity of multiplication, division and square root, implemented respectively by addition, subtraction and right shift of the logarithms of the operands. Besides, multiplication and division are exact (they incur no rounding error).

The main drawback of LNS is that addition and subtraction are much more complicated than in floating-point. In the case of the sum, if and are two positive numbers such that > , we have

+ = log ( + ) = log (1 + -) = + ⊕ ( -), with ⊕ ( ) = log (1 + ),
The function ⊕ ( ) is a good example of composite functions that can be implemented using the techniques described in the present article, for a simple implementation of LNS. With this approach the constant parameter can be changed easily [START_REF] Syed Asad Alam | Low-precision Logarithmic Number Systems: Beyond Base-2[END_REF].

In order to demonstrate the capabilities of the library, a simple LNS adder has been developed. The simplified adder is presented in listing 6.

The LNS Adder has been synthesized for an lns_t having an exponent of FixedFormat<5,-6, unsigned>. This requires to evaluate ⊕ with an accuracy constraint of 2 -8 . The report after synthesis and place-and-route shows that the design uses two BRAM, that matches what was expected.

The next section describes a more complex application that was developed after having checked that this simple experiment gives a design with sound resource consumption.

Additive sound synthesis

Additive synthesis is a sound synthesis technique typically used in electronic music to create timbres by adding sine waves together [START_REF] Lockhart | The electrical production of music[END_REF] as in the following computation :

( ) = =1 ( ) sin(2 + ) .
It allows very rich sounds when a lot of sines are used, but then it is very compute-intensive. When used in an interactive configuration with hardware in the loop to simulate real instruments interacting in real-time with the real world physics, it requires very low latency.

In order to check the impact of the library, an additive synthesis kernel was developed both with FixedNumber and float computation. Both version have an underlying collection of 256 oscillators, having a frequency = * 1000 256 and phases = 0. The amplitudes is an input of the kernel. In both versions, the oscillators are instantiated using template unrolling, ensuring that they are inlined.

In the FixedNumber implementation, the are fixed-point number of format FixedFormat<-1, -8, unsigned>.

The sine is computed from a wave table of 2 12 entries generated by the library. The function thus evaluated is ( ) = (2 ) for in [0; 1). Two experiments, FixedNumber 8 and FixedNumber 16 have been performed with an evaluator accuracy constraint of 2 -8 and 2 -16 respectively.

In the float implementation, the std::sin function from the mathematical library is used. All the computations are done in float.

Function value over one period for all the formats and difference with float for FixedNumber 8 and FixedNumber 16 is plotted on Figure 4. The graphs for all the format coincide perfectly, and as expected, the absolute difference to float is smaller for FixedNumber 16 than for FixedNumber 8 .

The evaluation is performed using our compiler and targeting a Virtex ultrascale plus FPGA with part number xcvu13p-fhga2104-3-e. Both pipelined and unpipelined solution have been generated for the three formats. Results are given after place and route. The design metrics are are reported in table 1.

These results shows that fixed-point is definitely a good option for this application. Three factors can explain this better quality of results. The first one is a question of bitwidth of the argument that have to be moved in the design.

Indeed, the 256 8 bits amplitude coefficients input of the FixedNumber case can be implemented as a very wide input register while it is not possible for the 256 32bits coefficients of the input. Secondly, the intrinsic arithmetic of floating point is more complex than fixed-point on the operation involved in this application. Lastly, the fact that floating-point arithmetic is based on external IP integrated by the design can hinder some optimization opportunities.

Single-source end-to-end LNS with SYCL

As our compiler also offers an experimental support for compiling SYCL programs targeting AMD FPGAs, we tried to use the library in a full end-to-end single-source application. In the SYCL compilation flow, the computation kernels are outlined by the compiler, and compiled separately for the hardware targets. Instead of producing a Vitis IP, it generates a x86 fat binary that embeds the hardware kernel binary, and contains all the glue to invoke the kernel execution on the FPGA. Using the previously described LNS adder and the expression evaluator generator framework, we built a small SYCL application containing a addition between two lns_t having an exponent of format FixedFormat<7,-8> running on FPGA. While being able to run it and launch it on an alveo u200 board, the process was far from straightforward.

Due to the shell of the board that is used and the M_AXI controller synthesized to move inputs and outputs between the host and the kernel, it is difficult to comment on the QOR of the design. It however validates that it is possible to integrate our expression evaluator generator framework into a different compilation flow.

FUTURE WORK

A limitation of the current evaluator architecture is its restriction on function of only one free variable. Accuracy constrained efficient evaluation of function of many variables is in the general case is an open problem. However, for some specific functions like composition of sums and products, it is possible to build evaluators that have better performance than the execution of all the individual operations. Specialized efficient evaluators include for instance truncated tiled multipliers or bit heaps [START_REF] Brunie | Arithmetic core generation using bit heaps[END_REF]. An other issue is the lack of mechanism to represent multiple-output functions. While these functions can be computed with one sub-function for each input, this can shadow some sharing optimization. An example of such sharing opportunities can be seen with multiple constant multipliers [START_REF] Gustafsson | A Difference Based Adder Graph Heuristic for Multiple Constant Multiplication Problems[END_REF][START_REF] Voronenko | Multiplierless Multiple Constant Multiplication[END_REF].

Representing expressions as trees can also hinder optimization. While this is not an issue for the currently supported class of functions, an expression would be more adequately represented as a directed acyclic graph. Indeed, the same variable can be reused twice. For instance, the expression x + x can be optimized as 2 * x if we are able to understand that the two x variables represent the same value. This problem is more fundamental with the current architecture. As the expression representation is ultimately depending on the type of the nodes, being able to identify the reuse of a variable would require each variable to have a unique type. There is no mechanism that allows to do this automatically within the current C++ standard, so the user would need to add manually a unique ID in the declaration of each variable.

This way of doing would be far from convenient and error-prone.

An alternate approach that can lift this limitation would be to move inside the compiler all the work that is currently done when running the specialization header generation executable. The compiler can identify the shared leaves from the invocation context.

If going in this direction, it can also be beneficial to move toward using the MLIR framework [START_REF] Lattner | MLIR: Scaling Compiler Infrastructure for Domain Specific Computation[END_REF] to describe and optimize expression DAG instead of relying on a custom expression tree runtime class.

Finally, some improvement could be done on FixedNumber. One of them would be to parameterize the rounding and overflow behavior of the class, to give more freedom of usage. The second one would be to associate bounds with the types. This would allow better automatic format deduction for operand result. As of now for instance, adding together 3 FixedNumber<0,0, ...> will produce an output of the deduced type FixedNumber<2, 0, ...>, while FixedNumber<1, 0, ...> would be sufficient.

CONCLUSION

This paper presents a framework to perform high-level optimization in a single source manner at HLS stage. This framework usage is illustrated by a generator of accuracy-constrained evaluator for mathematical functions. In order to evaluate this generator, a fixed-point library and a custom C++ compiler have been developed. This compiler allows to run Vitis HLS on C++20 programs. The system is shown to be able to build efficient (while simple) evaluators for functions of one free variable. Limitations of the current infrastructure have been identified as well as directions to lift them.

Listing 2 .FixedFormat<5Fig. 1 .

 21 Fig. 1. Example of function evaluator (here correspondig to listing 3).

FixedNumber<FixedFormat<5, - 12 ,

 12 unsigned>> functionInput; functionInput = (...) // HLS code computing functionInput FreeVariable x{functionInput}; // defining a free variable auto c = 1.5_cst; auto f = sin(x + c); auto functionOutput = f.evaluate<-17>(); Listing 3. Construction of the expression tree corresponding to ( ) = sin( + 1.5) (lines 3-5), and construction of a fixed-point evaluator for this function (line 6).

Fig. 2 .

 2 Fig. 2. The library usage flow. File extensions describe the type of each state outcome.

Fig. 3 .

 3 Fig.3. Example bipartite approximation architecture, here replacing a table of 26 entries with two tables of 2 4 entries.

  This header contains template specialization for the required Evaluator, which is the object produced under the hood by the calls to evaluate. Listing 4 shows the outline of the generic Evaluator class, while listing 5 illustrates what a specialization looks like. This listing is a simplification, and of course there is no need to use the evaluator mechanism for a simple addition.

	template<ExpressionType ET, prec_t OutputPrec>
	struct Evaluator {
	Evaluator(){ /* Find evaluation plan and dump it to specialization header */ }
	auto evaluate(ET expr) {/* Return marker type indicating unevaluated context */}
	};

Listing 4. Outline of the generic Evaluator class. template<> struct Evaluator<SumExpr<FreeVariable<...>, Constant<1, ...>, /*OutputPrec=*/-17> {, .. auto evaluate(SumExpr<FreeVariable<...>, Constant<1, ...>> & expr) { return expr.variable + 1; } }; Listing 5. One Evaluator specialization for + 1 at with a target accuracy of 2 -17 .

Table 1 .

 1 Fig.4. Output of an additive synthesis with two frequencies of equal amplitudes and the computational error (as absolute difference on a logarithmic scale). Area and timing metrics comparison between float and various FixedNumber for an additive synthesizer of 256 oscillators. II stands for Initiation Interval, the number of clock cycles that pass before the pipeline is ready to be fed a new input.

								2 -8
	sin(4πt) sin(2πt) +	2	-0.75 -0.50 -0.25 0.00 0.25 0.50 0.75	float FixedNumber 16 FixedNumber 8		2 -14 2 -11 2 -17 2 -20 2 -23		k = 16 k = 8	10 -5 10 -4 10 -3 10 -6 10 -7	k | FixedNumber |float -
			0.00	0.25	0.50	0.75	1.00		0.00	0.25	0.50	0.75	1.00
					t						t
			Experiment			Area				Timing
			Pipeline	Format		LUT	Flip-flop DSP BRAM Latency Critical path (ns)	II
				float 12389	15222	175	0	653	2.787	-
				FixedNumber 6	7932	6985	256	129	257	2.907	-
		Unpipelined	FixedNumber 8	8882	7284	256	257	297	2.901	-
				FixedNumber 16	5777	6917	256	513	284	2.995	-
				FixedNumber 22	6551	7304	256	769	270	2.945	-
				float 141491 166723 1304	0	226	2.995	128
				FixedNumber 6	7442	4977	256	128	5	2.803	1
			Pipelined	FixedNumber 8	7568	6500	256	256	5	2.701	1
				FixedNumber 16	4513	8806	256	512	8	2.806	1
				FixedNumber 22	5627	9385	256	768	8	3.074	1

The library is available on github at this address https://github.com/lforg37/marto, as subpart of the Marto library.

It is also possible to create it by specifying the exact type.

SYCL is a programming standard from the Khronos Group[START_REF] Khronos® | SYCL™ 2020 Specification[END_REF] to program applications using heterogeneous accelerators with a modern single-source C++ framework

Acknowledgements

This work was partly funded by the Imprenum project of Agence Nationale de la Recherche.


J-P Bongiraud

G Cauffet
email: cauffet@leg.ensieg.inpg.fr

C Jeandey

Ph Le Thiec

An Optimization Method for Magnetic Field Generator

Keywords: Magnetic Environment Simulator, Air coils Optimization Noncircular coils field Expansion, Tchebychev polynomials

We propose an optimization method to design an elongated three-axes magnetic field generator with given criteria specified over a large volume. The approach is based on the field expansion in Spherical Harmonics and Tchebychev polynomials, for noncircular symmetrical coils arrangement.

We developed a specific tool, to get a "flat" or a given "equal-ripple" solution, over a chosen length. The parameters to be defined are: dimensions, coil positions and Amp-turns, associated with the different axes. Once these parameters have been computed, a program predicts the field for the whole structure.

A major interest of using such a method lies on the fact that, once the true optimal solution is found, any deviation from theoretical results (for instance building inaccuracy) can be compensated by adjustment of any other design parameters (i.e. current), restoring the initial homogeneity.

This method has been successfully applied to the simulator of the Magnetic Metrology Laboratory for Low Field (Laboratoire de Métrologie Magnétique en Champ Faible in French). The experimental results correspond to the theoretical computation.

Introduction

The Magnetic Metrology Laboratory for Low Field -MMLLF -( LMMCF in French) is an experimental facility conducting research and measurements in the area of very low magnetic fields (typically under 1nT of noise). For this purpose, we need a magnetic environment simulator able to compensate local earth field and create any field between ± 50 000 nT. The homogeneity and accuracy must be better than 10 -3 on the largest usable volume within the dimensions of the building (27m*9m*9m). To reach this goal, we decided to design a tri-axial set of coils, respectively:

• longitudinal simulator or L coils(Y axis, N-S) • vertical simulator or V coils (Z axis)

• transversal simulator or T coils (X axis, W-E).

The theoretical homogeneity must be close to 10 -4 over the largest volume.

To define and optimize the field uniformity of such set of coils, two basic approaches can be considered:

• The first one uses interactive computer programs to improve an initial coil configuration by trials and minimization of errors. Here, the objective function must take in account the desired homogeneity on the volume and geometrical constraints due to the building. For example it can be the sum of the squares of the field deviations on particular points selected by the designer [START_REF] Lugansky | On optimal synthesis of magnetic fields[END_REF]. This method's success bases strongly on the experience and the intuition of the designer, and the choice of the starting point ("the seed") is very important. A bad choice can drive the optimization function toward a local minimum which presents no interest. • The second approach uses analytical methods. One usual way consists in a field expansion of spherical harmonics. This design method developed in the 1950's was widely used in NMR experiments and became standard in the shimming of MRI scanner [START_REF] Garret | Axially symmetric systems for generating an measuring magnetic fields Part I[END_REF]; it gives a "flat" response by cancellation of successive derivatives at the center of symmetry [3][4]. However, for a prolate volume where a small "ripple" is allowed, this method is less efficient and Tchebychev polynomials expansion is preferentially used because they approximate a function over the greatest length with a minimum pk to pk error or "equal-ripple". This technique was first introduced by CARTER [START_REF] Carter | Coil System Design for production of uniform magnetic fields[END_REF] for the design of coils and recently shown again by M. LEIFER [START_REF] Leifer | R F solenoid with extended equiripple field profile[END_REF]. CARTER mentioned also that spherical harmonics solution could be obtained from Tchebychev expansion, by letting the specified range converge toward zero.

The longitudinal field simulator

When starting the project, we decided to use hexagonal coils for the L simulator. This geometry is closed to the circular one that gives the lower ripple for equivalent surface and space between coils [START_REF] Bongiraud J-P | The Magnetic Metrology Laboratory for Low Field in Grenoble[END_REF]. Besides it is also well integrated within the roof and the bottom of the building especially designed. In order to keep clear the median vertical plane, interesting for sensors location, we also chose an even number of coils, arranged symmetrically with respect to the origin. Thus, we have a set of equal coaxial coils to provide a field with a given uniformity on the maximum length. The parameters of this system are adjusted according to the method.

The number of coils is determined with respect to the minimal specified pk to pk error. If the expected homogeneity is not reached after optimization, coils number is increased.

2-1 Tchebychev polynomials expansion

The field generated on its axis by a regular polygonal coil of n sides, inscribed in a circle of radius a, at a distance d and supplied by a current I, is: The addition of several coils gives rise to a ripple on the axis. The expression (1) may be expanded in series of equal-ripple functions as Tchebychev polynomials. The position and current may be adjusted to cancel the successive coefficients. So, the axial field will be homogeneous within the high-order terms, which are negligible, provided they decrease in magnitude. Tchebychev polynomials are given by: T n (x)=cos(n*arc cos(x)) for -1 ≤ x ≤ 1 [START_REF] Garret | Axially symmetric systems for generating an measuring magnetic fields Part I[END_REF] They are orthogonal when integrated over [-1;1] with a weight function

( ) ( ) 2
(1-x 2 ) -1/2 : () ( )      = = ≠ = ≠ = - ∫ - 0 n m 0 n m n m 2 / 0 dx x 1 ) x ( T * x T 1 1 2 m n π π (3)
As a consequence of this property, it is possible to expand a function f(x) within the range [-1;1] as a serie of Tchebychev polynomials:

f(x)= t 0 +t 1 *T 1 (x)+t 2 *T 2 (x)+…. ( 4
)
Where:

∫ - - 1 1 2 0 dx ) x 1 ( ) x ( f 1 π = t (5) dx ) x 1 ( ) x ( f * ) x ( Tn 2 t 1 1 2 n ∫ - - = π for n≠0 (6)
For an hexagonal coil pair at y = ± d and for a= 1, the field is: [START_REF] Bongiraud J-P | The Magnetic Metrology Laboratory for Low Field in Grenoble[END_REF] If L is the half-length where the ripple has to be minimized, we make the substitution y = L*x and cosθ = y/L to expand Hy(y,d) over [-L;L]; expressions (5) and (6) become:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )         + + + + + - + - + = 2 2
() θ θ π π d ) d , cos L ( H 1 d t 0 y 0 ∫ = (8) () θ θ θ π π d ) n cos( * ) d , cos L ( H 2 d t 0 y n ∫ = (9)
Let's take a double pair of hexagonal coils as an example (figure 2). By symmetry of the system with respect to the origin, the odd coefficients are cancelled. By choosing a unity current (I 1= 1) in the middle pair of coils, only three parameters must be adjusted: I 2 , d 1 and d 2 . It means that we will be able to minimize three even coefficients of the expansion (terms of 2 nd , 4 th and 6 th degree).

B 2 =1*t 2 (d 1 )+I 2 *t 2 (d 2 ) B 4 =1*t 4 (d 1 )+I 2 *t 4 (d 2 ) B 6 =1*t 6 (d 1 )+I 2 *t 6 (d 2 )
B 2 , B 4 and B 6 represent the residual error on the field. We chose to minimize the function f, sum of the quadratic deviations associated with the Tchebychev polynomials of order 2,4, and 6.

f(I 2 ,d 1 ,d 2 )=B 2 2 +B 4 2 +B 6 2 (10)

2-2 Application to the L coils

To obtain the equal-ripple allowed on the length expected, the above mentioned method leads to install 7 pairs of coils (see further down). It means that 13 parameters have to be adjusted: 6 currents (I 2 ,…I 7 ) and 7 distances (d 1 ,d 2 ,…d 7 ). We must obtain the simultaneous cancellation of the 13 coefficients, B 2 to B 26 : To find the optimum, we have to minimize the "target" function f: A direct use of traditional optimization algorithms for a system of nonlinear equations with 13 unknown is not obvious and it takes too much computation time. Moreover the convergence zone is very restricted for such a set of parameters.

B 2 =1*t 2 (d 1 )+I 2 *
f(I 2 ,
To improve the solution, we will proceed in two steps. In a first step, starting from an initial set of parameters (d 1 ,……d 7 ) we can notice that half of the coefficients can be determined by solving in a linear way the system: (13) We obtain the optimal values for (d 1 ,...d 7 ) by minimizing expression (13) with available standard algorithm (optimization toolbox MatLab). Furthermore, the t 2n (d) coefficients are more efficiently calculated by using a Gauss-Legendre quadrature approximation instead of integrating expression (9) [START_REF] Press | Numerical Recipies[END_REF].

B 2 =B 4 =B 6 =B 8 =B 10 =B 12 =0 (12) i. e.            - = + + + + + - = + + + + + - = + + + + + - = + + + + + - = + + + + + - = + + + + + ) d ( t )
There is only one optimal solution; we call "the canonical solution". This solution is the best one because any deviation from the theoretical results does not change in a significant way the field homogeneity.

As mentioned above, in the first step, we have to start with an initial set for (d 1 …d n ). With a low number of coils, the choice for this initial set is not critical and we have a large length range available that gives us the canonical solution. When the coil number is increasing, some local minimum due to a bad set could give unstable solutions and it is not realistic to directly optimize a system with 7 pairs of coils. We first apply this method to a double pair of coils whose the results are used as "initial set" for 3 pairs of coils, and so on, until the imposed ripple is reached.

The progression toward the result is summarized in table 1. With 7 pair of coils, we got the expected uniformity (<2 10 -4 ) on the maximum length (L=2.65). The distances d 1 to d 7 are in normalized units, compared to the radius of the circle inscribed within the hexagon (a = 1). Fig. 3 represents these results and shows how it is possible, step by step to found the canonical solutions when increasing the number of coils. The parameters (7 distances and 6 currents) of the LMMCF simulator corresponding to the optimal (canonical) solution are given in table 2 (normalized units). L is the half-optimization length from the origin. Vertical axis represents the relative homogeneity variation from the origin.

Number of pair of coils

N.B For a set of n coils, the canonical solution leads to n field oscillations, regularly decreasing over the optimization length. This characteristic belongs only to the canonical solution making it a good way to discriminate from sub-optimal solutions.

2-3 Final design

For simplicity and stability, all the coils are connected in series and supplied by a bipolar generator such as only one current has to be controlled. The six inside pairs of coils carry the same number of Amp-turns ; the two end coils are identical but with a different number of Amp-turns. Considering this configuration, only 9 parameters must be adjusted -2 currents and 7 positions-. Table 3 gives the new parameters of the optimal solution for the final realization: These results are slightly different from canonical solution and the field homogeneity is compared in Fig. 5. We can notice that the number of oscillations is only six instead of seven, which confirm that this optimal solution is not the canonical solution (previous N.B). However, as we started from the canonical set for (d 1 ,…d 7 ) defined on table 2, the field homogeneity is slightly decreased but the ripple specification is still respected over the required length.

Finally, the longitudinal simulator consists of 14 coils connected in series with the positions and Amp-turns ratio defined in table 3 (standardized units).

Vertical and transversal simulator

Once the main longitudinal simulator has been determined, we must define the coils set for the two perpendicular directions (vertical and transversal). The basic shape will be the same and consists of several rectangular loops whose length and width are close to those of the building. The structure for V or T simulator is the same, except it's rotated by 90°, so this enables us to study only one set.

While for the L coil set, the ripple is only optimized along the Y axis, (the offaxis homogeneity then depends of the outer circle radius), for V or T we have to take into account two directions for each one:

along the axis normal to the loop (Z for V and X for T) which define an homogeneity in the vertical median plan parallel to the plane of the loop, on the Y axis.

3-1 Homogeneity in the vertical plane

In a first step, let us consider a system with infinite wires. We have to determine the number of infinite wire pairs, giving us the homogeneity required over the maximum length along the axis in a vertical plane perpendicular to these wires.

Using spherical harmonics [START_REF] Hoult | Magnetic field profiling, analyzing and correcting coil design[END_REF], we express the field of an infinite wire as a development of the nth derivative in Taylor series about the origin. This development allows optimization of positions and currents for infinite wires in a cross section, by cancellation of even successive derivatives.

Results are given in Fig. 6 and table 3. T 2 , T 4 or T 6 type means that derivatives are cancelled up to 2 nd , 4 th or 6 th order.

Table 4: Optimal conductor location and Amp-turns vs number of infinite wire pairs

The T 4 design has been chosen because it's the only one witch allows to use the 45° support for both, vertical and transversal simulator. It needs a ratio 2 between the mid coil current, and the upper and lower ones.

As mentioned earlier, these results come from the optimization using infinite straight conductors. They concern only the homogeneity along the normal axis of the loops and lead to the simplest structure with 3 rectangular coils, as shown in fig. 7.

This elementary layout gives the expected results along the Z axis, however on longitudinal direction Y, the terminal segment effects limit the homogeneous length and they must be taken in account. One way to reduce these effects is to split the end segments of median coil into two parts and push back those of side coils. The modified main structure of V simulator with horse saddle shape for each coil, is represented in Fig. 8 B' C'

D' A' D" A" C" B" D A B C Me d i a n V c o i l Lower coi l Upper coi l D' A' D" A" B' C' C" B"
x y z Fig. 8: Modified V simulator structure with horse saddle.

Choosing a horse saddle shape for the coils, brings up several advantages: access in the simulator is much easier than with a straight connection in the horizontal plane, especially for median coil. field homogeneity on median longitudinal axis is also increased. moreover, with hexagonal ends, the same frame as the longitudinal simulator is usable and both coil sets (V ant T) will be fitted with L coils.

3-2 Longitudinal Y axis optimization

Although horse saddle shape with hexagonal ends provides a better homogeneity than rectangular loops, this partial improvement is not sufficient to meet the specifications along the Y axis.

To increase significantly this homogeneity length, we install additional coils similar to the main one, as shown in Fig. 9. These additional coils produce a ripple along the Y axis. By adjusting the positions and the Amp-turns of each coil with the Tchebychev polynomial expansion method, we minimize this ripple.

The system, with 5 unknown parameters -2 currents and 3 positions-is solved in the same way as for the longitudinal simulator. Only expression [START_REF] Bongiraud J-P | The Magnetic Metrology Laboratory for Low Field in Grenoble[END_REF] for the field is modified to the right one corresponding to a multi segments coil [START_REF] Durand | Magnétostatique[END_REF].

The results of this optimization are given in table 5: simulator (L=2). The geometry of the upper and lower ends is not the same as the median one, so we obtain two sets of positions and currents. There is only a slight difference between the two sets. In order to simplify the construction, we took the average of the 2 values, after verifying it had no effect on the overall homogeneity.

Theoretical and experimental results.

The components of the local earth field are: H L :20700nT; H V : 40400nT and H T <100nT

For hexagonal coils with 4m sides (a = 4) and to create a field of ± 50 000 nT, the final parameters for the field simulator are summarized on the tables 6 and 7, respectively for L and V coils. For the V simulator, the upper and lower coils (or side coils for T) are located at 45°, on the side of the hexagon, so the distance to be considered is 3,69m (a' = 3,69)

4-1 Field verification

Once the geometry and the Amp-turns of the simulator have been determined, an another program using the Biot and Savart law, computes the field created by each axis (or by the complete three-axes simulator).

As an example, Fig. 10 and 11 give the results for the field along Y axis, for L and V simulators. The residual ripple meets our requirements, i. e. lower than 2.10 -4 .

The designed homogeneity along all the axes of this three-axial simulator (9 is summarized in table 8. Table 8: Field Homogeneity on the 3 axis.

4-2 Measurements

The magnetic simulator has been built in accordance with the geometry defined thanks to this optimization and is now in use.

The first experimental measurements gave excellent results, very close to the predicted one (ripple <2 10 -4 on a 15m length )

As an example; Fig. 12 shows the measured field created by the horizontal simulator(L).along the Y axis This measurement has been done with a fluxgate sensor; fixed on a cart moving along Y axis. We present only a length 5.5 m but the homogeneity is the same over 20 m.

The field created by L coils was close to 20700 nT, such as to compensate horizontal component of local earth field. The measured homogeneity reach 2.3 10 -4 (4 nT/20700 nT).

Conclusion

An analytical approach to optimize the design of a three-axes magnetic field simulator has been described. This method, using Spherical harmonics and Tchebychev polynomial expansion gives us a flat or an equal-ripple solution we called "canonical solution".

When increasing the coils number, the program developed is really efficient to minimize this ripple, in order to satisfy given specifications.

Moreover, the complementary program predicting the air coils field in all the space is very useful to check the theoretical results given by the optimization routine.

By applying these tools to our simulator, we obtained the required homogeneity, i.e. lower than 2.10 -4 , in a large volume (1,6m*1,6m*15m).

The experimental measurements have validated the method and the theoretical results.

We have now a flexible tool, enable us to design any air coils structure, in order to satisfy given criteria.

As a final remark, Tchebychev polynomials method could be taken into account more than it is today for magnetic system design.

6 Fig. 1 :

 61 Fig. 1: Hexagonal Coil.

2 :

 2 System with two pairs of coils.

CoilFig. 3 :

 3 Fig. 3: Coils position (normalized unit).

Current 1 .

 1 000 0.981 0.959 0.942 0.950 1.098 2.378

Fig. 4 Fig. 4 :

 44 Fig.4 gives the relative homogeneity of the field along the Y axis.-4∆H/H ( * 10 -4 )

Fig. 5 :

 5 Fig.5: Field homogeneity for the two solutions.

Fig. 6 :

 6 Fig. 6: Disposition of infinite wires in a cross section for vertical field. Homogeneity Number of planes

Fig. 7 :

 7 Fig. 7: Elementary structure for the vertical or transversal (rotated by 90°) simulator.One way to reduce these effects is to split the end segments of median coil into two parts and push back those of side coils. The modified main structure of V simulator with horse saddle shape for each coil, is represented in Fig.8

Fig. 9 :

 9 Fig. 9: 1/4 of the complete V simulator.

Fig. 11 :

 11 Fig. 10: Field along the Y axis for the L simulator.B(%)

Fig. 12 :

 12 Fig. 12: Axial homogeneity for L Coils.

  t 2 (d 2 )+I 3 *t 2 (d 3 )+I 4 *t 2 (d 4 )+I 5 *t 2 (d 5 ) +I 6 *t 2 (d 6 )+I 7 *t 2 (d 7 ) M B 26 =1*t 26 (d 1 )+I 2 *t 26 (d 2 )+I 3 *t 26 (d 3 )+I 4 *t 26 (d 4 ) +I 5 *t 26 (d 5 )+I 6 *t 26 (d 6 )+I 7 *t 26 (d 7 )

Number of pair of coils -L=2.65

	2	3	4	5	6	7
	D 5			3.0075 2.3793 2.0495
	D 6				3.0763 2.511
	D 7					3.176

D 1 0.8678 0.58 0.4352 0.3482 0.2782 0.2371 D 2 2.4133 1.665 1.2731 1.0277 0.8273 0.7061 D 3 2.6965 2.0408 1.6643 1.3572 1.1637 D 4 2.8754 2.2667 1.8657 1.61

Table 1 :

 1 Results from 2 to 7 pair of coils.

Table 2 :

 2 Canonical Solution for LMMCF

(7 currents and 7 distances, L=2.65).

Table 3 :

 3

Optimal Solution for final realization (2 currents and 7 positions) L=2.65

Table 5 :

 5 Optimal solutions for vertical

Table 6 :

 6 L simulator (real units with a = 4).

	Pair N°	1	2	3	4	5	6	7
	Distance (in meter)	0.977 2.930 4.878 6.818 8.740 10.632 12.977
	Amp Turn 14	14	14	14	14	14	30

Table 7 :

 7 V simulator(real units for a'= 3,69)

Acknowledgment

This works and this project has been supported by the DGA/DCE/GESMA. B. Barreyres was very helpful in translating the original H.P. Basic software into an up to date Matlab version.