Thibault Tricard

Semyon Efremov

Cédric Zanni

Fabrice Neyret

Sylvain Lefebvre

Procedural Phasor Noise

Keywords: CCS Concepts:, Computing methodologies → Texturing procedural, textures, noise, Gabor, pattern, texture synthesis

oscillating field with prescribed main frequencies and preserved contrast oscillations. In addition, the profile of each oscillation is directly controllable (e.g. sine wave, sawtooth, rectangular or any 1D profile). Our technique builds upon a reformulation of Gabor noise in terms of a phasor field that affords for a clear separation between local intensity and phase. Applications range from texturing to modeling surface displacements, as well as multi-material microstructures in the context of additive manufacturing.

Fig. 1. High-contrast patterns produced by our approach. Note how the profile of the oscillations smoothly transition from a rectangular wave (20% black), to a square wave, to a triangular profile and finally a sine wave. At the same time, the orientation of the waves changes from left to right. The field visualized here is purely procedural. It is obtained by feeding our phasor noise into periodic profile functions (shown in blue), that are interpolated from left to right.

Procedural pattern synthesis is a fundamental tool of Computer Graphics, ubiquitous in games and special effects. By calling a single procedure in every pixel -or voxel -large quantities of details are generated at low cost, enhancing textures, producing complex structures within and along surfaces. Such procedures are typically implemented as pixel shaders. We propose a novel procedural pattern synthesis technique that exhibits desirable properties for modeling highly contrasted patterns, that are especially well suited to produce surface and microstructure details. In particular, our synthesizer affords for a precise control over the profile, orientation and distribution of the produced stochastic patterns, while allowing to grade all these parameters spatially. Our technique defines a stochastic smooth phase field -a phasor noisethat is then fed into a periodic function (e.g. a sine wave), producing an

INTRODUCTION

The Computer Graphics community is in a perpetual quest to extend the range of patterns that can be generated at low memory and computational costs, from a single procedure. Most techniques rely on procedural noises [START_REF] Ares Lagae | A Survey of Procedural Noise Functions[END_REF], which generate scalar fields with prescribed frequency content. These base noises are then combined through various functions to produce interesting patterns [START_REF] David S Ebert | Texturing & modeling: a procedural approach[END_REF].

A key limitation of standard procedural noises is the lack of direct control over the local characteristics of the patterns, such as contrast, 57:2 • T. Tricard and S. Efremov, C. Zanni, F. Neyret, J. Martínez, S. Lefebvre. sharpness, as well as the range and histogram of values within a neighborhood. A degree of control is obtained though thresholding a base noise, or through lookup tables (color maps) reshaping the global histogram. Unfortunately, noise patterns generally exhibit variations such as local loss of contrast and undesired variations in scales and orientations. These variations make it very difficult to obtain an overall consistent appearance through lookup tables, and it is rarely possible to find parameters producing a consistent, stable aspect, such as the waves shown in Figure 1. Lookup tables also cannot distinguish between the 'sides' of an oscillation, preventing the synthesis of patterns with asymmetric aspects, such as the sawtooth wave in Figure 1, middle.

In this paper we propose a procedural synthesizer that has the ability to generate patterns exhibiting precisely controllable oscillations, in terms of their profile shape, minimum/maximum values, orientation and frequencies (scale). The synthesizer can freely grade the pattern parameters, within planes (2D) and volumes (3D). Our noise also works along surfaces -see surface noise in [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]] -which affords for the synthesis of patterns along objects, changing their appearance and surface roughness profiles. We explore the various possibilities in Section 4. Thanks to its procedural nature, our approach is implemented as an efficient OpenGL shader.

Contrary to most existing techniques, our procedure does not directly generate a scalar field of values. Instead it generates a phase field designed to be subsequently modulated through a periodic function. This phase field -our phasor noise -is such that, once seen through the periodic function, it produces a pattern with the specified frequency and orientation. In addition, the shape of the oscillation profile remains that of the periodic function without any undesirable fluctuations. This affords for a precise control over the shape of the produced stochastic structures, illustrated in Figure 1.

We achieve this by reformulating Gabor noise [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]] as the explicit product between an intensity field, and a modulated sine wave. The phasor noise is obtained as the instantenous phase of the Gabor oscillations. Its name originates from our reformulation, which exposes a sum of phasors 1 within Gabor noise.

Feeding a phasor noise into a sine wave produces a signal that ressembles Gabor noise, has a similar spectrum but oscillates with a corrected, uniform local contrast (see Section 3). This can serve as a fundamental harmonic to reproduce any 1D periodic profile by an approach akin to inverse Fourier transform (i.e. additive synthesis). Profiles are synthesized as weighted sums of integer harmonics of the base pattern, see Figure 13.

Being able to orient synthesized patterns while controlling their profile and distribution is especially well suited for the definition of patterns such as hatches, stripes, cracks, ridges, scales and ripples. Our technique also finds applications in modeling for additive manufacturing, for sculpting grip patterns or for precise control of proportions in multi-material mixtures. Contributions. In summary, our contributions are:

• The definition of phasor noises, enabling pattern synthesizers with precisely controlled oscillation's profiles. • An analytical separation of Gabor noise into an intensity and a modulated sine wave component.

1 See e.g. https://en.wikipedia.org/wiki/Phasor.

• A numerical analysis of the resulting signals, in terms of their spectrum and spectrum of variance. • A complete implementation with available source code.

PREVIOUS WORK

Procedural noises [START_REF] Ares Lagae | A Survey of Procedural Noise Functions[END_REF]] have a long history in Computer Graphics and are widely used in applications ranging from games to movie special effects. They allow for arbitrary amounts of content to be synthesized at low memory and computational costs. Procedural noise methods focus on producing seemingly random patterns in the spatial domain, while providing an accurate control over their frequency content. A variety of 2D patterns -texturesare obtained by combining layers of noises having different spectral contents [START_REF] David S Ebert | Texturing & modeling: a procedural approach[END_REF].

One of the most recent development in this area are techniques based on the so-called Gabor noise [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]. A Gabor noise is an instance of a sparse convolution noise [START_REF] Lewis | Texture Synthesis for Digital Painting[END_REF][START_REF] Lewis | Methods for Stochastic Spectral Synthesis[END_REF], where a kernel is convolved in the spatial domain with a random impulse noise. An important property is that the spectrum of the obtained noise is that of the kernel. Thus, precisely controlling the kernel spectral content translates to a similar control over the noise. The Gabor kernel offers such a precise, elementary control.

Several improvements and variants of Gabor noise have been proposed, for isotropic noises [START_REF] Lagae | Filtering Solid Gabor Noise[END_REF], filtering of solid (3D) noises [START_REF] Lagae | Filtering Solid Gabor Noise[END_REF], fitting Gabor noises to examples [START_REF] David S Ebert | Texturing & modeling: a procedural approach[END_REF][START_REF] Gilet | Procedural Descriptions of Anisotropic Noisy Textures by Example[END_REF], as well as modifying the underlying formulation for faster synthesis [START_REF] Tavernier | Making Gabor Noise Fast and Normalized[END_REF] and explicit manipulation of phase [START_REF] Gilet | Local Random-phase Noise for Procedural Texturing[END_REF].

The original Gabor noise [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]] supports two important features in our context. First, the ability to spatially grade the noise content, freely changing frequencies and orientations. Second, the basic formulation can be applied in 2D and 3D, but also along any surface equipped with a tangent field. This is achieved without the need for a global parameterization (texture coordinates). Preserving these properties while manipulating phase is especially difficult. For instance, the formulation by [START_REF] Gilet | Local Random-phase Noise for Procedural Texturing[END_REF] localizes kernels at the nodes of a grid to synchronize their phases. While this provides excellent results in the context of 2D texture synthesis, the phase synchronization cannot be preserved through an orientable field or along a surface.

loss of contrast interferences

When it comes to synthesizing highly contrasted patterns, Gabor noise suffers from two interrelated defects: local loss of contrast and interferences, as illustrated in the Figure inset. To our knowledge, these have been identified and studied for the first time by [START_REF] Neyret | Understanding and controlling contrast oscillations in stochastic texture algorithms using Spectrum of Variance[END_REF], who provide deep insights regarding these defects and focus on analyzing and filtering them. In particular, they show that the spectrum of variance of the noise -the Fourier transform of the squared signal -characterizes the local loss of contrast. We follow their insight and rely on the spectrum of variance to analyze the properties of our noise. We provide additional discussion later in the paper.

Procedural textures in additive manufacturing. The procedural texture synthesizer we propose is general purpose, but one application we aim at is additive manufacturing, and more specifically for the synthesis of multi-material patterns.

There are several major challenges when considering pattern synthesis for AM. First, the patterns should afford for a wide range of effects that the user should be able to freely choose from. In addition, it is expected that the patterns can be progressively spatially graded, triggering different effects in different regions of space.

Second, the patterns should be very efficient to compute. In most cases, thousands of high-resolution slices have to be generated for a single part and sent to the printer. It is also important for the user to be able to preview the result beforehand.

To address these challenges, one line of research advocates for the use of procedural approaches, inspired from procedural texturing in Computer Graphics [START_REF] David S Ebert | Texturing & modeling: a procedural approach[END_REF][START_REF] Livesu | From 3D Models to 3D Prints: An Overview of the Processing Pipeline[END_REF]. This is especially interesting to grade materials having anisotropies: the absence of underlying grid affords for complete freedom in orientation [START_REF] Martínez | Orthotropic K-nearest Foams for Additive Manufacturing[END_REF]]. Several frameworks have shown how to exploit procedural synthesis to achieve efficient process planning in AM [START_REF] Vidimce | Foundry: Hierarchical material design for multi-material fabrication[END_REF][START_REF] Vidimče | OpenFab: a programmable pipeline for multi-material fabrication[END_REF].

Our work fits well within this trend. However, rather than proposing a specific type of structure, we seek to develop a novel type of procedural texture generator well suited for pattern synthesis in modeling for additive manufacturing. In particular, we seek to design multi-material patterns akin to laminates, with precise control over the material mixtures. We also consider surface displacements, with the ability to produce different roughness profiles.

PHASOR NOISE

We now describe our approach in details, and analyze its properties. Applications and results are described in Section 4.

The crux of our technique is to rewrite a Gabor noise as a single sine wave, with a clear separation of intensity and modulation:

Gabor Noise (x) = I (x) sin (ϕ(x)) (1)
where x is the point of evaluation. For conciseness we denote Gabor noise as G(x) in the remainder of this paper. The phasor noise is the instantaneous phase ϕ of the rewritten Gabor noise:

Phasor Noise (x) = ϕ(x) (2)
As we will see, the phasor noise is the argument (angle) of a complex valued version of Gabor noise, which is obtained by summing phasors 2 representing the sine waves of the Gabor kernels. This change of formulation is defined in the most general case of Gabor noise, where kernels of different frequencies, orientation, bandwidth and amplitudes may be summed and where all these parameters can be spatially graded. Thus, phasor noise works in 2D, 3D and for surface noise [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]. It is a pure procedural function that can be implemented in a GLSL pixel shader.

Having clearly separated the instantaneous phase from the noise intensity has two direct benefits. First, it affords for the definition of a noise that perfectly oscillates without local loss of constrast, akin to a form of normalized Gabor noise:

PhasorSinewave(x) = sin(Phasor Noise(x)) (3)
We call the result a phasor sine wave and denote it by P(x) for conciseness. We further analyze its properties in the following Sections.

The second benefit is to afford for controlling the shape (profile) of the noise oscillations themselves. A direct way to achieve this is to replace the sine wave by another periodic function:

Pro f iledNoise(x) = f (Phasor Noise(x)) (4)
where f is a function of period 2π . Alternatively, phasor noise can be used to generate harmonics of a base phasor sine wave, allowing to produce different profiles through an inverse Fourier transform (additive synthesis):

Noise (x) = K k =1 w k • sin k • Phasor Noise(x) + ϕ k
This is a weighted sum of K > 0 integer harmonics with synchronized phases. The weights w k and phase shifts ϕ k control the produced profile, allowing to synthesize any noise profile, e.g. rectangular, sawtooth, or any other user controlled profile. While similar to specifying a profile function, this offers novel possibilities such as optimizing the weights -profile synthesis -or dynamically adjusting them for filtering purposes (e.g. anti-aliasing).

To understand our approach and its properties we will study noises of increasing complexities, starting by the simplest phasor noise using a single main frequency and orientation, to richer noises.

For convenience we summarize the main notations introduced throughout the paper in Table 1.

The bi-lobe case

The simplest form of phasor noise is obtained from a Gabor noise using a single kernel. As a reminder, a Gabor noise is obtained by convolving a Gabor kernel with a random impulse process. The spectrum of the noise is then the spectrum of the kernel [START_REF] Ares Lagae | Procedural Noise using Sparse Gabor Convolution[END_REF]].

In the frequency domain, a Gabor kernel is two symmetric Gaussians, produced by the convolution of two delta functions and a Gaussian. This is illustrated in Figure 2. The delta functions come from the sine wave in the spatial domain, and are precisely positioned in the spectrum by changing the frequency and orientation of the wave. The Gaussian defines the bandwidth of the noise, e.g. its frequency spread around the perfect sine wave. As in [START_REF] Neyret | Understanding and controlling contrast oscillations in stochastic texture algorithms using Spectrum of Variance[END_REF] we call this simplest form of noise the bi-lobe case.

The main observation making our work possible is that, in the bi-lobe case, the noise is obtained by summing in the spatial domain a large number of sine waves of same frequency weighted by Gaussians. This can be rewritten as a single weighted sine wave, with an explicit separation of intensity and phase. This operation is akin to phasor addition in signal processing.

We start from a bi-lobe Gabor noise, defined in 2D (the same derivations apply in 3D). It is parameterized by a main frequency F 2π , along direction u. In every position x, using the simplified expression of Gabor noise with sine waves and no random weights [START_REF] Tavernier | Making Gabor Noise Fast and Normalized[END_REF] 3 . It evaluates as:

G (x) = n j=0 a x -x j sin F • x -x j • u (5)
where a is a centered Gaussian of bandwidth b:

a (x) = e -π b 2 | |x| | 2 (6)
Through phasor addition, we rewrite the noise as:

G (x) = I (x) sin (F • x • u + φ (x)) ( 7
)
where I is the intensity and φ the phase. This rewritten Gabor noise is strictly equivalent to the initial one. In this case the phasor noise is the instantaneous phase field ϕ

(x) = F • x • u + φ (x)
, and we denote the bi-lobe phasor sine wave as P (x) = sin(ϕ (x)).

Our claim is that, in the bi-lobe case, removing the intensity Ikeeping only the sine wave -creates a noise oscillating with perfect local contrast, i.e. there are no low frequencies in the spectrum of variance. Yet, the local shape of the sine wave is preserved, and the overall frequency remains F . This means that harmonics of the phasor sine wave are preserved as well, and thus the characteristics of the profile function -that can be seen as a Fourier combination of harmonics of the base phasor sine wave -are preserved. While this might seem obvious when looking at Equation 7, the difficulty comes from the fact that the new phase φ (x) does depend on x. By varying it could reshape the sine wave arbitrarily, significantly changing the frequency content of the result.

Before entering the details, let us give the full expression of I and φ obtained from the sum of sine waves. This is a direct application of the phasor addition formulas 4 . The intensity is obtained as:

I (x) = j a j (x) sin (φ j ) 2 + j a j (x) cos (φ j ) 2 (8)
and the phase is obtained as:

φ (x) = atan2 j a j (x) sin (φ j ), j a j (x) cos (φ j ) (9)
where a j (x) = a xx j and φ j = -x j • u • F , and where atan2(y, x) returns the polar angle of the vector (x, y).

It is worth noting that atan2 is only defined if the vector norm is non zero (I (x) 0). We later discuss singularities occurring when I vanishes (Section 3.1.4).

In the next sections we will verify that by omitting the intensity we eliminate the local loss of contrast in the bi-lobe case of the noise. At the same time, the overall frequency content is preserved: the sine wave still oscillates at the same main frequency.

3.1.1 Phasor field. The rewrite of Gabor noise in terms of a single sine wave exposes an underlying complex valued field. Indeed, I (x) and φ (x) -Equations 8 and 9 -can be directly interpreted as the modulus and argument of a sum of complex numbers:

n j=1 a j (x) cos (φ j ) + i sin (φ j ) = n j=1 a j (x) e iφ j
where e iφ j are random complex numbers of unit modulus. This is the phasor field, illustrated in Figure 3.

To understand its effect, let us consider the complex version of Gabor noise underlying the phasor addition:

G (x) = e i F •x•u n j=1 a j (x) e iφ j (10)
The real valued Gabor noise is obtained as

G (x) = Im(G (x)).
G is the product of two terms: a perfectly oscillating complex wave (e i F •x•u ) and the phasor field. Thus, in a bi-lobe Gabor noise, the phasor field is the source for both the local loss of contrastvia its varying modulus I (x) -and the phase shifts φ (x) -via its argument. A key question is whether the phase shifts φ (x) could be strong enough to introduce frequencies lower than F .

Interestingly, the phasor field is a sparse convolution obtained by summing spatially located Gaussians (a j ) multiplied by random, unit complex numbers. Therefore, the produced random field has the bandwidth b of the Gaussians, and the speed at which φ changes is driven by b. When b is much smaller that F , the variations in phase due to φ are much slower than that of the base frequency F . This implies that φ has little impact on the local profile of the sine function, while its impact on the overall frequency is limited by b.

This hints at a first result: the spectrum of the phasor sine wave is similar to that of the Gabor noise -since φ has little impact on F -while its spectrum of variance exhibit no low frequency (no local loss of contrast in the signal) -since φ cannot reshape the sine wave profile significantly. Note that this is only true for b smaller than F .

The phasor field also reveals that the phase will exhibit singularities. Indeed, in some singular locations the modulus of the phasor field vanishes: this is where the sum of complex numbers becomes null. Around these singularities φ rotates abruptly. We discuss this effect in Section 3.1.4.

3.1.2 Spectrum of the phasor sine wave. In this section and the following we only outline the main results of our analysis of phasor noise. Please refer to the supplemental material for full details and derivations. We perform the analysis on periodic instances of the phasor sinewave. This is achieved by sampling kernels in a periodic domain. In this case, the Fourier transform of P is :

F x [P (x)] (ω) = ω 0 +∞ k =-∞ F x [B (x) P (x)](kω 0 ) δ kω 0 (ω) (11)
with B a windowing function verifying the partition of unity property and where ω 0 = 1 T 0 with T 0 the analysis period. The amplitude of the deltas δ in F x [P (x)] is driven by the spectrum of F x [B (x) P (x)] which is:

F x [B (x) P (x)](ω) = 1 2 F x [B (x) (sin (φ (x)) -i cos (φ (x)))] ω - F 2π u + 1 2 F x [B (x) (sin (φ (x)) + i cos (φ (x)))] ω + F 2π u
The two terms are functions convoluted with two symmetric deltas located at the main frequency of the noise (ω ± F 2π ). These two functions are smooth and their Fourier transforms are rapidly decaying (holds for one dimension, verified experimentally in 2D and 3D due to singularities -please refer to the supplemental material for details). This explains the presence of two lobes in the spectrum.

3.1.3 Spectrum of variance the phasor sine wave. We use the same approach to derive the spectrum of variance, noting that if P is periodic, then P 2 is also periodic. The spectrum of variance is the spectrum of the squared signal, F x [B (x) P 2 (x)]. In this case, the amplitude of the deltas in the Fourier transform are obtained as (coloring the main terms to help the description):

F x [B (x) P 2 (x)] (ω) = 1 2 δ 0 (ω) + 1 4 F x [B (x) (1 -2 cos 2 (φ (x)))] ω - 2F 2π u - i 2 F x [B (x) sin (φ (x)) cos (φ (x))] ω - 2F 2π u + 1 4 F x [B (x) (1 -2 cos 2 (φ (x)))] ω + 2F 2π u + i 2 F x [B (x) sin (φ (x)) cos (φ (x))] ω + 2F 2π u
(12) where the terms in orange and green are again two functions that are convolved with symmetric deltas located at twice the main frequency of the noise (± 2F 2π u). These terms correspond to the lobes shown in Figure 4. The central dirac δ 0 corresponds to the average value of the variance, which is non zero. The most important observation is the absence of any central lobe, which indicates that the phasor sine wave suffers no local loss of contrast.

The non-centered lobes could however contribute low frequencies in the spectrum of variance, depending on the ratio between F and b (the bandwidth of the Gaussians a j ). In this case removing the intensity terms from the Gabor noise does not remove all low frequencies from the spectrum of variance. Going back to the interpretation of the bi-lobe noise in terms of underlying random phasor field (Section 3.1), this is an instance where the phasor field varies more rapidly than the base sine wave frequency, thereby changing its profile. Such a case is illustrated in Figure 5.

3.1.4 When the phasor modulus goes to zero. The intensity I (x) of the Gabor noise goes to zero in a few singular locations. Note that we do not suffer from numerical issues since we define the phasor sine wave P by omitting to multiply by I . The only potential source of trouble (atan2) is well identified. In the implementation null phasors never truly occur thanks to the small rounding errors.

In such locations, we can observe singularities in the phase φ: abrupt rotations around a singular point or abrupt variations across ridges (see Figure 6). By going to zero I hides these occurrences in the original Gabor noise. When removing I we bring these singularities into the signal. This creates sharp points and edges, which correspond to locations in the Gabor noise where the sine wave oscillations are in phase opposition.

There are two main types of issues (see Figure 6). First, those which occur at a point where one sine wave band disappears. Second, those which occur along abrupt variations in the field, where the wave fronts on both sides oppose. These correspond to curved ridges that are connecting two point singularities.

The singularities only add high frequencies in the direction orthogonal to the waves. As they are relatively few and localized, their influence on the overall signal is limited. The number of singularities is impacted by the choice of Gabor bandwidth b -the smaller the bandwidth the closer to an ideal sine wave.

It is worth noting that [START_REF] Neyret | Understanding and controlling contrast oscillations in stochastic texture algorithms using Spectrum of Variance[END_REF] proposed an iterated scheme to filter away these issues -paving the way for potential solutions. However, the several iterations of direct and inverse Fourier transforms cannot be done in a procedural manner. Also, such a filtering approach does not give access to an explicit separation between a base sine wave and its argument, precluding the ability to control the oscillation's profiles.

While we currently do not address these singularities, this is a clear venue of future work. We believe the phasor field interpretation to open interesting possibilities by directly exposing the source of these defects.

General case

We now present our approach to produce general phasor noise patterns with richer frequency content.

As with Gabor noise, kernels of different orientations and frequencies are combined together. We then extract the instantaneous phase modulation ϕ(x) and evaluate a single sine wave at the end. Thus, in the most general case we define the noise as:

G (x) = I (x) sin (ϕ (x)) (13)
It contains Gabor kernels at arbitrary frequencies F j and direction u j (unit vectors), that is:

G(x) = I (x) sin (ϕ(x)) = n j=1 a j (x) sin (ϕ j (x)) (14)
where

a j (x) = a ||x -x j || and ϕ j (x) = F j • (x -x j ) • u j (15)
and the phasors are a j (x) • cos (ϕ j (x)) , sin (ϕ j (x)) . Through phasor addition, we obtain I and ϕ as: Equivalently, I and ϕ are the modulus and argument of a complex version of the noise G, defined similarly as in Section 3.1.1. Note that Equation 13 remains strictly equivalent to the original Gabor noise, but rewritten as a single sine wave. In the following we will study the properties of the phasor sine wave P = sin(ϕ(x)) obtained by omitting I .

I (x) =
To understand the behavior in the general case, we start by discussing the case of two bi-lobes in Section 3.2.1. We then proceed to the general case combining kernels of different frequencies and orientations in Section 3.2.2. It can be understood as a generalization of the two bi-lobe case.

Two bi-lobes.

Let us consider a Gabor noise combining two bi-lobe noises of same frequency F but with different orientations u k and u m . Applying the phasor addition we obtain:

G (x) = G u k (x) + G u m (x) = I u k (x) sin (ϕ u k (x)) + I u m (x) sin (ϕ u m (x)) = I (x) sin(ϕ (x))
with:

I (x) = (I u k cos ϕ u k + I u m cos ϕ u m ) 2 + (I u k sin ϕ u k + I u m sin ϕ u m ) 2 ϕ (x) = atan2( I u k sin ϕ u k + I u m sin ϕ u m , I u k cos ϕ u k + I u m cos ϕ u m )
This is the sum of two phasors, one for each initial Gabor noise G u k and G u m . These phasors vary in space and have different behaviors along different directions.

In the following we discuss the spectrum of variance and spectrum of the two bi-lobe phasor sine wave. Detailed derivations, obtained from the complex version of the noises, are provided in the accompanying supplemental material.

Please note that the angle θ between the two bi-lobes, that is between directions u k and u m , plays a role in the analysis.

Spectrum of variance: two bi-lobe Gabor noises. In total, the spectrum of variance for Gabor noise in the two bi-lobes case exhibits Gaussians in nine locations, shown in Figure 8, left. Five of these Gaussians correspond to the squared bi-lobe noises: centered Gaussians and Gaussians at ±2u k and ±2u m coming from respectively G 2

u k and G 2 u m (recall that the spectrum of variance is the spectrum of the squared noise). The four other Gaussians emerge from interactions between the two bi-lobes. They appear when squaring the Gabor noise:

G 2 (x) = G 2 u k (x) + G 2 u m (x) + 2G u k G u m .
The last product term leads to the remaining four Gaussians located at ±(u k -u m ) and ±(u k + u m ): the product in the spatial domain is a convolution in the spectral domain, resulting in the diagonal lobes.

Spectrum of variance: two bi-lobe phasor sine waves. We now consider the spectrum of variance of the phasor sine wave resulting from omitting I , see Figure 7 middle. We observe the absence of centered Gaussians -implying there is no loss of local contrast -and the presence of lobes at ±2u k and ±2u m . The Gaussians at ±(u k -u m ) have disappeared ; however, those at ±(u k +u m ) remain. This is only partially satisfying as the Gaussians at ±(u k + u m ) may fold back onto low frequencies depending on the angle θ between u k , u m . It is possible to obtain the symmetric effect, ensuring that Gaussians at ±(u k + u m ) are removed in place of those that ±(u ku m ), see Figure 7, right. This is achieved by inverting the direction of one of the base noises. This may be used to select the approach depending on the angle θ between u k , u m , avoiding low frequencies to be introduced by the diagonal bi-lobes (e.g. in Figure 7 selecting the middle case instead of the right one where the lobes get closer to the center). The switch should occur when u k • u m changes sign, that is when the angle θ between both directions exceeds π 2 . As a side note, we can also choose to cancel only the centered Gaussians in the Gabor spectrum of variance. To achieve this, we keep the noise intensity I but divide it by I 2 u k + I 2 u m . The full noise expression becomes

I (x) I 2 u k (x)+I 2 um (x)
sin ϕ (x). An example is shown in Figure 8, right. The rationale is that the central lobe originates from that of both individual Gabor noises. Hence, by rescaling the noise with

I (x) I 2 u k (x)+I 2 um (x)
we cancel the center lobes but preserve the others. However, this renormalization is different from our initial intent with phasor noise: it only applies to a sine wave profile.

Phase singularities. As with the bi-lobe case, phase singularities appear where intensity vanishes. However, another type of phase defect appears, which depends on the angle θ between the two bi-lobe directions.

To understand this effect let us consider -without loss of generality -a single direction. For simplicity we select u k . Along this direction the phasors take different frequencies, respectively F and F • u m • u k (this stems from the dot product with u j s in Equation 15). Now, depending on the angle θ between u m and u k , the two phasors will interfere with one another in different manners. The different cases we are about to discuss are illustrated in Figure 9.

If u m = u k the pattern is exactly similar to the bi-lobe case. For small angles θ , we obtain a smooth pattern with gentle orientation changes. For larger angles however, interferences start to introduce sharp turns in the phase. This is most visible when the angle θ between the two bi-lobes is π 2 . Visually this results in an abrupt interruption of the waves ('chopping'), shown in Figure 9, rightmost. High frequencies appear in the spectrum reflecting the perceivable edges. This defect is easily explained by considering the phasor diagram of two sine waves in the π 2 case, as detailed in Figure 10. [START_REF] Lagae | Filtering Solid Gabor Noise[END_REF]6,32) with orientations regularly distributed over π . Bottom row: Corresponding phasor sine wave. Both the spectrum and spectrum of variance (outlined in blue) are shown. The achieved power spectrum are similar. The phasor noise has some additional energy in higher frequencies due to singularities. The phasor sine wave eliminates a large part of the low frequency disc in the spectrum of variance. The parts that remain are due to the twosided normalization observed with the two bi-lobes case, here generalized for large number of bi-lobes. We discuss this further in Section 3.2.2.

Whether the visual artifact is problematic depends on the target application. However, these defects appear progressively. This hints at a possibility to reduce their impact by locally avoiding mixing patterns with large differences, which we discuss in Section 3.2.2.

Mixed kernels.

The most general case involves summing kernels of randomly sampled orientations and/or frequencies [Lagae Figure 11 illustrates the patterns obtained when mixing an increasing number of bi-lobes. Each additional bi-lobe introduces new pairwise interactions, which are as many additional Gaussians in the spectrum of variance (the number of Gaussians doubles for each new bi-lobe).

We observe a generalization of the behavior of two bi-lobes: A large portion of the low frequency disk in the spectrum of variance of the phasor sine wave is removed, albeit in a two-sided manner. The remaining low frequencies are due to pairs of bi-lobes separated by an angle θ above π 2 . However, if the angle between all bi-lobe pairs remains below π 2 , all low frequencies are suppressed. We exploit this to formulate an isotropic phasor noise that strongly reduces the mixing of incompatible bi-lobes.

Locally coherent isotropic noise. In the standard approach, the isotropic noise is produced by randomly sampling orientations in every kernel, i.e. u k = random_unit_vector (x k ). To better control which kernels mix together in the noise we modify the sampling of directions. Instead of using random independent sampling, we rely on a Gaussian sparse convolution noise with controlled bandwidth:

u j = n j=1 д x -x j r j (ω f ) || n j=1 д x -x j r j (ω f )|| where д(x) = e -π •b 2 f • | |x| | 2
and r j (ω f ) are random unit vectors sampled in the range of angles [-ω f , ω f ]. The impulses x j are independent from those of the phasor noise. By changing b f we control how locally coherent the direction field is. By increasing ω f we smoothly go from the bi-lobe case to arcs, and ultimately to a full ring. To obtain an isotropic noise we set ω f = 2.4 and b f = b 2 with b the bandwidth of the phasor noise. This does not remove all areas of abrupt changes, but strongly reduces the spurious interferences. As shown in Figure 12 the ring in the spectrum of variance is now close to perfect. Of course, while the spectral content is now satisfactory, this significantly changes the aspect of the field. Nevertheless, this provides an additional control to trade-off between the more standard isotropic aspect and an aspect which, while being locally similar to the bi-lobe case, is globally isotropic.

Finally, let us note that this incurs a performance penalty, as the direction field has to be sampled at every kernel.

Oscillation's profiles

Phasor noise provides an instantaneous phase modulation that affords for the synthesis of noises with arbitrarily shaped oscillation profiles.

This can be achieved in two ways. The most direct is to feed the phasor noise into a periodic function having a 2π period. This generates a field where the main oscillations follow that of a phasor sine wave, but where the oscillations are shaped by the new periodic function. In addition, the periodic function itself may vary spatially, offering a unique control over the patterns, as shown in Figure 1.

The second way is to perform weighted sums of integer harmonics of a base phasor sine wave, as illustrated in Figure 13. This is a form of inverse Fourier transform, using the phasor sine wave as the fundamental harmonic. The linear summation has the advantage of being analytic and offering a direct control over introduced frequencies. It also paves the way for optimizing the profiles or dynamically adjusting them during rendering (e.g. anti-aliasing).

Figure 23 shows various profile shapes and how they modify the spectrum and spectrum of variance. We further explore profile reshaping in Section 4.

RESULTS

Our technique is especially well suited to produce highly contrasted patterns and strongly oriented, laminate-like stochastic structures. Direct control over the profile enables further spatial gradations.

We start by showing the high degree of control made possible by phasor noise in Section 4.1. We continue with applications in modeling surface displacements and textures in Section 4.2. We finally illustrate potential applications for synthesizing multi-material laminates, in the context of additive manufacturing, in Section 4.3.

Throughout the section, please keep in mind that the shown effects are all purely procedural [START_REF] Ares Lagae | A Survey of Procedural Noise Functions[END_REF], and computed in a few milliseconds from a pixel shader.

Controls

Figure 15 illustrates the variety of controls enabled by phasor noise. In this particular case the profile being used is a rectangular signal which 'up' time in each period can be modified -visually modifying the white/black ratio and producing thicker or thinner lines. We refer to this profile as a Pulse Width Modulation (PWM) profile.

As can be seen in Figure 15 all parameters can be freely spatially graded, including the profile parameters themselves. Figure 14 shows a similar control using a grayscale image as a source. Reshaping the profiles is a powerful tool, especially when targeting a precise histogram and homogeneous appearance. For instance, Figure 16 compares the synthesis of a noise with a fixed, chosen width (25% high). With phasor noise we immediately get the corresponding profiled pattern, which exhibits the chosen average histogram. On the contrary, thresholding a Gabor noise involves manually tuning the threshold, and the achieved result is not spatially homogeneous.

Figure 1 illustrates how different oscillation profiles may be blended together while controlling other parameters. The profile functions are linearly interpolated based on the abscissa: the two profiles enclosing the lookup point are accessed and the result is linearly interpolated.

Surface displacements and texturing

Another promising application of phasor noise is its use for modeling surface displacements. The profile, frequency and orientation controls allow to produce different displacements along the parts, changing significantly the appearance of surfaces, and possibly their friction and tactile properties. Figure 17 shows a lava sphere which surface is modeled by an isotropic phasor noise with sawtooth profile, as well as fingerprints generated by a phasor sine wave from a user specified direction field. Changing the random seed generates new sets of unique fingerprints. Figure 18 shows a complex 3D bracelet which outer crust is modeled with a phasor noise displacement.

Phasor noise also provides novel possibilities in terms of procedural texturing. For instance, Figure 19 shows cracked tree bark textures that answer different requirements along cone cross-sections: the bark stripes are controlled independently in their width and number. Such precise controls would not be possible through simple lookup tables applied to existing noise patterns. These effects being produced from a pixel shader they could also easily evolve with a dynamic geometry. This provides a fully procedural alternative to optimization-based parameterization methods for generating stripe patterns [START_REF] Knöppel | Stripe Patterns on Surfaces[END_REF].

Figure 21 illustrates how the oscillation's profile can directly mimic some naturally occurring phenomena such as the asymmetric sand ripples along desert dunes.

Finally, Figure 20 shows a variety of textures such as cracked paints and tree barks applied to tori. These texture combine several layers of phasor noises using rectangular and sawtooth profiles.

Multi-material laminates

Multi-material laminates are typically made of two sandwiched materials, one more rigid than the other. Such laminates are of high interest for additive manufacturing, as they produce orthotropic elastic behaviors [START_REF] Martínez | Orthotropic K-nearest Foams for Additive Manufacturing[END_REF]]: the part is more flexible orthogonally to the laminate than in the other two directions. Ideally, for maximal design flexibility, the laminates should be freely orientable in space, and the proportions between both materials should be precisely controlled.

Phasor noise provides such flexibility through PWM profiles. Materials A and B are directly mapped to the up/down parts of the profile. This affords for a direct control of the orientation and material proportions. Figure 22 illustrates a cross shaped plate where we can observe that the laminate pattern triggers different bending behaviors in the extremities. The Figure also shows a miniature 3D shoe insole with a spatially varying 3D laminate pattern. Both models are driven by control fields allowing to paint orientation, mixture ratios and degree of isotropy.

Note that this could easily accommodate for more than two materials by creating a profile with multiple steps.

LIMITATIONS

In the bi-lobe case the noise is overall well behaved. Point singularities are acceptable, however we would like to reduce the impact of ridges between them. These produce waves in opposition and 'chopping' of the sine waves. Now that we have exposed the phasor field underlying the noise, it would be interesting to consider how to filter singularities. However, doing this in a procedural manner seems difficult. Some undesirable interactions exist when mixing kernels using very different directions, with sine waves being abruptly cut off. Controlling the underlying direction field to obtain a locally coherent noise helps greatly and in practice, most interesting profiles have themselves high frequencies and mask these issues. However, ultimately we would like to reduce or remove these defects.

CONCLUDING REMARKS

There are links between phasor noise and the work of [START_REF] Neyret | Understanding and controlling contrast oscillations in stochastic texture algorithms using Spectrum of Variance[END_REF]. In particular, they introduce the notion of a normalized Gabor noise, which is filtered to suffer no local loss of contrast. This actually defines a full family of noises, as the normalization embeds a notion of window size (the smallest window resulting in a zero-thresholded noise -a noise normalized at all scales and hence with a square profile). Our phasor noise is normalized but has no notion of normalization scale. The normalization is an intrinsic property of its definition. Having a window size to chose from offers additional controls to the user, which is often desirable for texturing applications. One benefit of our approach, however, is that no a-posteriori filtering is required, while the filtering approaches in [START_REF] Neyret | Understanding and controlling contrast oscillations in stochastic texture algorithms using Spectrum of Variance[END_REF] involve either a local averaging through super-sampling, or an iterative Fast Fourier Transform algorithm. Finally, the ability to reshape the noise oscillation's profiles is quite unique to our method.

In a sense, the phasor sine wave is a generalized sine wave, which scale and orientations can be freely controlled in space. This leads to the possibility of performing additive synthesis to generate various profiles. We believe this is a novel tool in procedural synthesis that will lead to further developments. In particular, procedural noise is typically to add fine scale unstructured details [START_REF] Guingo | Bi-Layer textures: a Model for Synthesis and Deformation of Composite Textures[END_REF]. Phasor noise paves the way to the ability to synthesize increasingly structured patterns from purely procedural synthesizers.

A remarkable property of phasor noise is that it is purely analytical. Even in the case of using a PWM there is no thresholding: the result is obtained as a weighted sum of harmonics. This opens interesting venues of future work, such as optimizing the pattern profile directly under specific objectives (elasticity, appearance).

ACKNOWLEDGMENTS

This work was supported partly by the french PIA project « Lorraine Université d'Excellence » (ANR-15-IDEX-04-LUE), by the ANR MuFFin (ANR-17-CE10-0002), and the project IMPRIMA(ANR-18-CE46-0004). We thank Samuel Hornus for proofreading.

  Fig. 2. A bi-lobe Gabor noise and its Fourier transform. The orientation is θ , the frequency of the sine wave F and the bandwith b. They are all clearly visible in the spectrum. Note that the image of the Fourier transform has been cropped for readability.

Fig. 3 .

 3 Fig.3. From left to right: a phasor sine wave, its underlying phasor field (directions as hue, modulus as luminance) and associated phase field φ.

Fig. 4 .

 4 Fig. 4. Noise (left) and spectrum of variance (right). First row: Gabor noise. Second row: Phasor sine wave. The colors are the same as the ones used for equation 12. Note the absence of low frequency -no center lobe -in the spectrum of variance of the phasor sine wave.

Fig. 5 .

 5 Fig. 5. Noise and spectrum of variance. (Left: Gabor noise, Right: phasor sine wave).In this case the ratio between the bandwith of the gaussian and the frequency is large enough for φ to significantly reshape the main sine wave.

Fig. 6 .

 6 Fig. 6. Left: Intensity field. Middle: Color coded phase field (angled mapped on hue). Right: Phasor sine wave. Two types of artifacts appear in the phase field. The first are point singularities (in orange) around which the phase rotates rapidly. It corresponds to the disappearance of a wave front. When two such singularities are close to one another in the direction orthogonal to the noise (in cyan), a ridge appears in the phase with opposing wave fronts across.

  x) sin (ϕ j (x)), n j=1 a j (x) cos (ϕ j (x))

Fig. 7 .

 7 Fig. 7. Spectrum of variance of two bilobes. Left: Gabor noise exhibiting local loss of contrast, resulting in a central lobe in the spectrum of variance. Middle: Phasor sine wave with the two bi-lobes separated by π 4 . The local loss of contrast is completely eliminated: no central lobe. Note the twosided spectrum of variance. Right: Same phasor sine wave using the dual formulation to eliminate the other lobes instead.

Fig. 8 .

 8 Fig. 8. Two bi-lobes noise. Left: Gabor noise and its spectrum of variance.The nine Gaussian locations are clearly visible; in yellow the center Gaussians, in blue and green the bi-lobes Gaussians at ±2u k and ±2u m and outlined by orange and purple arrows the interaction Gaussians at (respectively) ±(u k + u m ) and ±(u ku m ). Right: Normalized noise where the center lobe is removed. This is not a phasor sine wave as an intensity term remains, please refer to the text.

Fig. 9 .

 9 Fig. 9. Phasor sine wave in the two bi-lobes case, with separation angles θ of respectively (left to right) π 10 , π 6 , π 4 , π 2 . The insets at the bottom right show the power spectrum of the patterns. As the angle θ increases, high frequencies due to interferences between the two bi-lobes appear.

Fig. 11 .

 11 Fig. 11. Top row: Gabor noises mixing an increasing number of bi-lobes (3,[START_REF] Lagae | Filtering Solid Gabor Noise[END_REF] 6, 32) with orientations regularly distributed over π . Bottom row: Corresponding phasor sine wave. Both the spectrum and spectrum of variance (outlined in blue) are shown. The achieved power spectrum are similar. The phasor noise has some additional energy in higher frequencies due to singularities. The phasor sine wave eliminates a large part of the low frequency disc in the spectrum of variance. The parts that remain are due to the twosided normalization observed with the two bi-lobes case, here generalized for large number of bi-lobes. We discuss this further in Section 3.2.2.

Fig. 12 .

 12 Fig. 12. Left: Phasor sine wave with randomly sampled orientations in [-π 2 , π 2 ]. The spectrum of variance has a two-sided aspect. Middle: Phasor noise with locally coherent random direction field. The spectrum of variance (bottom right) is now close to an ideal ring, with little residual energy in the low frequencies. The top right image shows the direction field with angle mapped on hue.

Fig. 13 .

 13 Fig. 13. Additive synthesis of a sawtooth profile. Top left: Base phasor sine wave. Top row: integer harmonics (without their weight for clarity). Middle row : Progressive summation of the weighted harmonics. Bottom: Obtained pattern with sawtooth profile ; a few more harmonics would be required to properly converge.

Fig. 14 .

 14 Fig. 14. Pattern produced from phasor noise, controlling a PWM profile width from the image gray level and the orientations from the gradients.

  57:10 • T. Tricard and S. Efremov, C. Zanni, F. Neyret, J. Martínez, S. Lefebvre.

Fig. 15 .

 15 Fig. 15. Columns, from left to right: Control fields, frequency control, PWM 'width' control, orientation control, all together.

Fig. 16 .

 16 Fig.16. Left: phasor noise through a PWM profile having 25% white (global average verified at 25.007 %). Right: Gabor noise thresholded to obtain an average of 25% (final at 0.25902%) the threshold was manually found as it depends on the noise parameters.

Fig. 17 .

 17 Fig. 17. Left: A lava sphere which surface is synthesized by an isotropic phasor noise with sawtooth profile. Right: Fingerprints modeled with a phasor sine wave following a user specified direction field.

Fig. 18 .

 18 Fig.18. Bracelet modeled with a phasor noise displacement. By varying the profile according to the local section of the torus, a tearing effect is achieved.

Fig. 19 .

 19 Fig. 19. Cones covered with a cracked tree bark pattern, using a PWM profile. Each cone achieves a different effect through different controls. Left: constant width and constant frequency, resulting in less stripes at the top. Middle: varying width and constant frequency, preserves the amount of white along cross-sections. Right: constant width and varying frequency.

Fig. 20 .

 20 Fig. 20. Top: Left torus is covered by ridges produced with a locally coherent isoptropic phasor noise and a triangular profile; right torus is covered by a grip pattern, obtained as the product of two correlated (same kernels distribution), orthogonal phasor sine waves. Middle: Cracked patterns. The left texture is obtained by crossing two phasor noises using PWM patterns of varying width. The right texture is obtained from a single triangular profile, clamped at a varying height. Bottom: Bark patterns. These barks combine several phasor noises with triangular and sawtooth profiles to achieve the final appearance, with varying degrees of isotropy.

  57:12 • T. Tricard and S. Efremov, C. Zanni, F. Neyret, J. Martínez, S. Lefebvre.

Fig. 21 .

 21 Fig.21. A desert scene. The sand surface pattern is a phasor noise fed into a sawtooth profile. A second higher frequency phasor noise adds details to the surface. The color is also modulated by the profile.

Fig. 22 .

 22 Fig. 22. Left: Cross shaped plate patterned with a phasor noise through a PWM profile. Printed using two materials (one flexible, one rigid). The user controlled fields are shown to the side: a ratios of material, b orientation and c isotropy. The patterns vary across the plate. Right: Shoe insole using a multimaterial pattern procedural synthesized with phasor noise. Both objects are printed on an Ultimaker 3.

Fig. 23 .

 23 Fig. 23. Profile synthesis using phasor noise, and effect on the spectrum and spectrum of variance. Top row: Noise patterns. Middle row: Noise spectrum. Bottom row: Noise spectrum of variance. From Left to Rigth: Gabor noise followed by phasor noises with sine, square, triangular and sawtooth profiles.

Table 1 .

 1 Summary of the main notations.

	2 See e.g. https://en.wikipedia.org/wiki/Phasor.

ACM Trans. Graph., Vol. 38,[START_REF] Lagae | Filtering Solid Gabor Noise[END_REF], Article 57. Publication date: July 2019.

This is not mandatory for our reformulation to be well defined, we adopt this point of view for the sake of simplicity.

see e.g. https://ccrma.stanford.edu/~jos/filters/Proof_Using_Trigonometry.html ACM Trans. Graph., Vol. 38, No. 4, Article

Publication date: July 2019.


Cyann Donnot

Antoine Genitrini
email: antoine.genitrini@lip6.fr

Yassine Herida

Unranking Combinations Lexicographically: an efficient new strategy compared with others

Keywords: Unranking Algorithm, Combination, Lexicographic Order, Complexity Analysis

We propose a comparative study of four algorithms dedicated to the lexicographic unranking of combinations. The three first ones are algorithms from the literature. We give some key ideas how to analyze them in average and describe their strengths and weaknesses. We then introduce a new algorithm based on a new strategy inside the classical factorial numeral system (or factoradics). We conclude the paper with an experimental study, in particular when n = 10000; k = 5000, our algorithm is 36 (resp. 53) times as fast as our C++ implementation of Matlab's (resp. Sagemath's) algorithm.

One of the most fundamental combinatorial object is called combination. It consists of a selection of items from a collection. In many enumerating problems it appears either as the main combinatorial structure, or as a core fundamental block because of its simplicity and counting characteristics.

In the 60s while resolving some optimization problem about scheduling Lehmer rediscover an important property linking natural numbers and a mixed radius numeral system based on combinations. This relationship gives him the possibility to exhibit some greedy approach for a ranking algorithm that transforms (bijectively) a combination into an integer. This numeral system is now commonly called combinatorial numbers system or combinadics. It is often used for the reverse of Lehmer's problem: generating the u-th combination (for a given order about the combinations). This approach is substituted to the exhaustive generation once the latter is not possible anymore due to the combinatorial explosion of the number of objects when their sizes increase. In the context of combination, the explosion appears really quickly. This generation strategy of a single element is classically called an unranking method. It is today often used as a basic brick in scheduling problems [START_REF] Tamada | Parallel algorithm for learning optimal bayesian network structure[END_REF] but also e.g. in software testing [START_REF] Myers | k-out-of-n:g system reliability with imperfect fault coverage[END_REF].

In order to unrank elements one must first define an order over the elements. The one that is usually used is the lexicographic one that can been seen as a generalization of the alphabetical order. The lexicographic order is humanly easy to handle with, and thus has been extensively studied. But, as Ruskey [17, p. 59] mentions, the lexicographic generation is usually not the most efficient, thus a particular care must be taken while unranking in this order.

The classical approach for the construction of combinatorial structures presenting a recursive decomposition schema consists in taking advantage of this decomposition in order to build a bigger object from a smaller one. The method has been extensively detailed in the famous book of Nijenhuis and Wilf [START_REF] Nijenhuis | Combinatorial algorithms[END_REF]. They did not directly study the example of combinations, but the case of integer compositions is proposed and is not so far from our problem. The method has been then applied generically on decomposable objects in the sense of analytic combinatorics, first in the context of recursive generation [START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF], and then in the context of unranking approaches [START_REF] Martínez | A generic approach for the unranking of labeled combinatorial classes[END_REF]. Fig. 1: Time (in ms) for unranking a combination, with n = 10000 and k = 0..n Aside such generic approaches there are several ad hoc algorithms but to the best of our knowledge no practical experiments does exist to determine which one should be used. We thus will start with te recall of three classical algorithms, supported by a revisited average complexity analysis. Then we will describe our new strategy and finally we will compare them in some experiments. In Fig. 1 we represent the time efficiency of some improved versions of 3 of these algorithms, our being depicted in green. The figure will be described later. But as a foretaste, when n = 10000; k = n/2, our algorithm is experimentally 36 (resp. 53) times as fast as our C++ implementation of Matlab's (resp. Sagemath's) actual algorithm. Along the paper, we represent combinations as follows.

Definition 1. Let n and k be two integers with 0 ≤ k ≤ n. We represent a combination of k elements among n denoted by {0, 1, . . . , n-1} as a tuple containing k distinct elements increasingly sorted from left to right.

For example, let n and k be respectively 5 and 3. The tuples (0, 1, 2) and (0, 2, 4) are combinations of k among n, but (0, 2, 1) and (0, 1, 2, 3) are not. There is another representation by using a 2-letters alphabet, but we do not deal with it through this paper. However we could rewrite the paper in this context also.

There are several orders for comparing combinations. In the following we restrict our attention to orders comparing combinations of the same length, i.e. the same number of elements.

Definition 2. Let A = (a 0 , a 1 , . . . , a k-1 ) and B = (b 0 , b 1 , . . . , b k-1 ) be two dis- tinct combinations.
-In the lexicographic order, we say that A is smaller than B if and only if both combination have the same prefix (eventually empty) such that (a 0 , . . . , a p-1 ) = (b 0 , . . . , b p-1 ) and a p < b p . -In the co-lexicographic order, we say that A is smaller than B if and only if the tuple (a k-1 , . . . , a 0 ) is smaller than (b k-1 , . . . , b 0 ) for the lexicographic order. -An order being given such that A is smaller than B, then, for the reverse order, B is smaller that A.

Definition 3. Let n, k define combinations and A be a combination. For a given order, the rank u of A belongs to {0, 1, . . . n k -1} and is such that A is the u-th smallest element.

With these definitions, we can enter the core of the paper organized as follows. Section 1 describes a recursive algorithms solving the lexicographic combination unranking problem. Despite its simplicity due to the approach, the fact it is recursive is problematic to handle very big combinations. Then Section 2 describes two classical algorithms based on combinadics. It seems that it is the first time both are compared and the reason why one is better is explained. In Section 3 we enter the context of factoradics, and then describe our new unranking algorithm in this numeral system. Finally in Section 4 we compare all algorithms experimentally and give some possibility for really improve the time complexity of all approaches.

Unranking through recursion

We deal with a combinatorial structure: combinations, that is very well understood from a combinatorial sense. Thus while trying to develop an unranking algorithm the first idea consists in developing a recursive algorithm based on the classical recursive generation presented in [START_REF] Nijenhuis | Combinatorial algorithms[END_REF]. Let us directly introduce such an algorithm. L := Rec_Generation(n, k, u) 3:

Algorithm 1 Recursive Unranking

return (n -1 -L[k -1 -i] for i from 0 to k -1) binomial(n, k) computes the value of n k ; append(A, a): appends element a in A; For i ≥ 0, L[i] is the element of index i in L.
In line 3 and 5, the outputs are tuples built by comprehension. The simple one in line 5 is (0, 1, . . . , k -1).

1: function Rec_Generation(n, k, u) 2:

if k = 0 then 3:

return () 4:

if n = k then 5:

return (i for i from 0 to k -1)

6: b := binomial(n -1, k -1) 7: if u < b then 8: R := Rec_Generation (n -1, k -1, u) 9:
append(R, n -1) 10:

return R 11:

else 12: return Rec_Generation (n -1, k, u -b)
Proposition 1. The function Rec_Generation(n, k, •) computes the combinations for n, k in the reverse co-lexicographic order.

Corollary 1. The function Unranking_Rec(n, k, •) computes the combinations for n, k in the lexicographic order.

The proposition is proved by induction and the corollary is a direct observation given in [11, p. 47].

In order to analyze the algorithm in details, we are interested in the average number of calls to the function binomial, when u describes the whole range of integers from 0 to n k -1. Ruskey justifies such a measure by supposing the table of all binomial coefficients precomputed, thus each call is equivalent. Later, in Section 4 we will discuss this measure.

Let us introduce the sequence u n,k computing the cumulative number of calls for the whole range for u.

Lemma 1. Let u n,k be the cumulative number of calls to binomial while unranking all possible u from 0 to n k -1. The sequence satisfies: u n,0 = 0 and u n,n+i = 0 for all n and i ≥ 0 and otherwise

u n,k = n k + u n-1,k-1 + u n-1,k .
Proof. The algorithm relies on the recursion underlying the binomial coefficients, with the additive cost 1 for each call to the function Rec_Generation(n, k, •) (in line 6). Furthermore, by calling the function once with the parameters n, k and u and a second time with n, k and v (such that u = v) then the following recursive calls will also be with distinct triplets. There is thus no multi-counting.

Theorem 1. Let U (z, y) be the ordinary generating function associated to This sequence is a shifted version of the sequence stored under the reference OEIS A0597973 . We thus can complete the properties in OEIS using our results.

(u n,k ), such that U (z, y) = n≥0 n k=0 u n,k y k z n . Then U (z, y) = 1 1 -z -zy 1 1 -z -zy - 1 1 -z - zy 1 -zy , thus, u n,k = n k k n + 1 k + 1 - 1 n -k + 1 .
Due to the values of the extreme cases when k = 0 and k = n and the symmetry in the recurrence we obviously obtain the fact that u n,k = u n,n-k , reflecting the symmetry of the binomial coefficients.

u n,k n k = n→∞ k=αn n + 2 - 1 α(1 -α) + O 1 n .
This average value is interesting to be noted, because it indicates that obtaining an extreme case (for the recurrence) during the recursive calls is obtained in average after (n + O(1)) recursive calls.

In the latter strategy, the recursive approach is a drawback for some programming languages that do not handle recursion efficiently (due to the depth of the stack). In fact, today the computer are able to handle combinations for very big values of n and k thus the recursive approach should be avoided. Naturally other strategies have been suggested in the literature.

Unranking through combinadics

In 1887, E. Pascal [START_REF] Pascal | Sopra una formula numerica[END_REF] and later D. H. Lehmer (detailed in the book [1, p. 27]) presented an interesting way to decompose a natural number, in what we call today a mixed radix numeral system, in their case it is the combinatorial number system, or combinadics. The decomposition relies on combinations.

Fact 2 Let n ≥ k be two positive numbers. For all integers u, with 0 ≤ u < n k , there exist a unique sequence

0 ≤ c 1 < c 2 < • • • < c k < n such that 4 u = c 1 1 + c 2 2 + • • • + c k-1 k -1 + c k k .
The writing (c 1 , . . . , c k ) is called the combinatoric of u.

For example when n = 5 and k = 3, the number 8 is represented as

1 1 + 3 2 + 4 3 , thus the combinadic of 8 is (1, 3, 4).
In 2004, using this representation, McCaffrey exhibited in the MSDN article [START_REF] Mccaffrey | Generating the mth Lexicographical Element of a Mathematical Combination[END_REF], an algorithm to build the u-th element (in lexicographic order) of the combinations of k elements among n. But in fact, this algorithm was already published in [11, p. 47] and can also be seen as an extension of the work of Lehmer. This algorithm is interesting in the sense it corresponds to the implementation used in the mathematics software Sagemath [19] 5 . In Fig. 3, on the left handside we present the algorithm used in Sagemath and on the right handside we exhibit the combinadics of some ranks of combinations of 2 elements among 6. More precisely we have chosen the reverse of the ranks to be able to compare them with their combinadics and the corresponding combination (due to line 3 in the algorithm). The following algorithm is also close to Er's algorithm [START_REF] Er | Lexicographic ordering, ranking and unranking of combinations[END_REF] whose representation for combinations is distinct but the computations are analogous; furthermore in his paper, Theorem 2 corresponds exactly to the combinadic decomposition.

The function Unranking_Combinadic(n, k, •) computes the combinations for n, k in the lexicographic order but the core of the algorithm is reverse colexicographic (in fact, the tuple L represents the reverse of the combinadic for efficiency reasons). The correction of the algorithm is presented in [START_REF] Kreher | Combinatorial Algorithms: generation, enumeration, and search[END_REF]. 4 We extend the definition of binomial coefficients with r-1 r = 0. 5 Unranking algorithm from Sagemath is stored in the Sofware Heritage Archive swh:1:cnt:c60366bc03936eede6509b23307321faf1035e23;lines=473-537

Algorithm 2 Unranking a combination

1: function Unranking_Combinadic(n, k, u) 2: L := () 3: u := binomial(n, k) -1 -u 4: v := n 5: for i from 0 to k -1 do 6: v := v -1 7: b := binomial(v, k -length(L)) 8: while u < b do 9: v := v -1 10: b := binomial(v, k -length(L)) 11: u := u -b 12:
append(L, v) 13:

return (n -1 -for in L)

We take n = 6 and k = 2.

rev. rank combinadic combination 14 (4, 5) (0, 1) 13

(3, 5) (0, 2) 12

(2, 5) (0, 3) 11

(1, 5) (0, 4) 10 (0, 5) (0, 5)

9 (3, 4) (1, 2) . . . . . . . . . 4 (1, 3) (2, 4) 3 (0, 3) (2, 5) 2 (1, 2) (3, 4) 1 (0, 2) (3, 5) 0 (0, 1) (4, 5)
Fig. 3: Algorithm in Sagemath; Examples of combinadics and their combination A first step in the complexity analysis of this algorithm consists in calculating the number of calls to the function binomial in Algorithm 2. The values n and k being given, the worst cases are obtained when v is as small as possible at then end of the loop, thus for all u whose combinadic satisfy c 1 = 0. To obtained a detailed analysis, we are also interested in the average number of calls to binomial, when u describes the whole range from 0 to n k -1. To reach this goal, let us introduce the function u n,k6 computing the cumulative number of calls for the whole range for u.

Lemma 2. Let u n,k be the cumulative numbers of calls to binomial while unranking all possible u from 0 to n k -1. The sequence satisfies: u n,0 = 1 and u n,n+i = 0 for all n and i > 0 and otherwise The single difference with Lemma 1 relies on the extreme cases. The key-ideas for the proof are given in Appendix A.2.

u n,k = n k + u n-1,k-1 + u n-1,k .
The leftmost part in Fig. 4 corresponds to the sequence given in Lemma 2. Our sequence is a shifted version of the sequence OEIS A264751. Both combinatorial objects can be put in bijection, and thus some conjectures stated there, are solved in the following.

Note, in this case the cumulative numbers are not symmetrical u n,k = u n,n-k . In fact the computation of the combinadics is not symmetrical. Theorem 3. Let U (z, y) be the generating function associated to (u n,k ). Then

U (z, y) = 1 1 -z -zy 1 1 -z -zy - z 1 -z ; u n,k = n k n + 1 - n -k k + 1 .
Corollary 3. The average number of calls to binomial in Algorithm 2 for n being large and k being of the form αn for α ∈]0, 1[ is

u n,k n k = n→∞ k=αn n + 2 - 1 α + O 1 n .
Algorithm 3 Unranking a combination faster

1: function Unranking_Combinadic2(n, k, u) 2: L := (0) 3:
r := 0 4:

for i from 0 to k -2 do 5:

while u >= r do 6:

L[i] := L[i] + 1 7: b := binomial(n -L[i], k -i -1) 8: r := r + b 9: append(L, L[i]) 10: L[i] := L[i] -1 11: r := r -b 12: append(L, L[i] + u + 1 -r) 13:
return L

In the literature there is another algorithm based on the combinadics. But there, in the computation of the combinadic for a given rank, the coefficients are computed from the smallest one, c 1 , to the second largest one, c k-1 , and finally the value for c k is directly deduced with no need of further trials. In this algorithm, the variable L contains the combinadic (not its reverse). First the computations that are needed are based on some smaller values for the binomial coefficients: in fact, the trials start with small parameters for the binomials and further they are increasing until we reach the solution, thus the biggest (non necessary) trials are saved. And second the last coefficient is directly deduced. This algorithm should be faster in practice. Such a modification in the order of computations can be put in parallel in the context of recursive generation or unranking, with the Boustrophedon traversal presented in [START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF] whose improvement relies also in the order of the computations.

Unranking_Combinadic2(n, k, •) is a lexicographic unranking for combinations. This algorithm is presented by Buckles and Lybanon [START_REF] Buckles | Algorithm 515: Generation of a Vector from the Lexicographical Index [G6][END_REF]. The correction of the algorithm is presented in [START_REF] Kreher | Combinatorial Algorithms: generation, enumeration, and search[END_REF]. Just as a comment, for the calculation of the next value L[i + 1], we start with the preceding value L[i] + 1 (cf. line 9). Finally note it is approximately the implementation in Matlab [START_REF] Butler | Function kSubsetLexUnrank, MATLAB central file exchange[END_REF] (the latter is a little bit less efficient: line 12 is omitted and some more trials are done in the loop). Lemma 3. Let u n,k be the cumulative numbers of calls to binomial while unranking all possible u from 0 to n k -1. The sequence satisfies: u n,k = 0 for all n, k = 1, 2 or k > n and otherwise

u n,k = n k + u n-1,k-1 + u n-1,k .
The result is proved in an analogous way as the one in Lemma 2, summing over c k-1 instead of c 1 . In Fig. 4 (right) we computed the first values of (u n,k ). We note the first values are smaller than the previous ones, but what about their asymptotic behavior?

Theorem 4. Let U (z, y) be the generating functions associated to (u n,k ). Then

U (z, y) = z 2 y 2 (1 -z) 2 (1 -z -zy) 2 ; u n,0 = 0; ∀k > 0, u n,k = n k k -1 k + 1 (n + 1).
Corollary 4. The average number of calls to binomial in Algorithm 3 for n being large and k being of the form αn for α ∈]0, 1[ is

u n,k n k = n→∞ k=αn n + 1 - 2 α + O 1 n .
The improvement for the efficiency of the algorithm cannot be deduced directly from the latter result in comparison to the one of Corollary 3.

In [17, p. 65] Ruskey presents an alternative implementation based on combinadics, in co-lexicographical order. But the first coefficients are computed starting by the largest ones, like in Algorithm 2.

Unranking through factoradics: a new strategy

Aside the combinatorial number system, there is another classical mixed radix numeral system: the factorial number system or factoradics. This decomposition relies on factorial numbers.

Fact 5 For all positive integers u, we define n such that (n -1)! < u ≤ n!, there exists a unique sequence (f ) ∈{0,...,n-1} , with 0 ≤ f ≤ for all such that

u = f 0 • 0! + f 1 • 1! + • • • + f n-2 • (n -2)! + f n-1 • (n -1)!.
The writing (f 0 , f 1 , . . . , f n-1 ) is called the factoradic decomposition of u.

First remark that f 0 = 0 for all u. Take the number u = 2020 as an example, we obtain its decomposition:

0 • 0! + 0 • 1! + 2 • 2! + 0 • 3! + 4 • 4! + 4 • 5! + 2 • 6!, thus its factoradic is (0, 0, 2, 0, 4, 4, 2).
Definition 4. Let n be a positive integer. A n-permutation is an ordering of the elements from {0, 1, . . . , n -1}.

We represent a n-permutation as a tuple of n components. For example the tuple (2, 4, 0, 3, 1) is a 5-permutation.

The factorial number system is used to define a one-to-one correspondence between integers and permutations: it is thus an unranking method for permutation. The algorithm implemented in the function Unranking_Permutation is an obvious adaptation of Fisher-Yates random sampler for permutations [START_REF] Fisher | Statistical tables for biological, agricultural and medical research[END_REF].

Algorithm 4 Unranking a permutation

1: function Unranking_Permutation(u, n) 2:

F := factoradic(u) 3:

while length(F ) < n do 4:

append(F, 0) 5:

return Extract(F, n, n)

1: function Extract(F, n, k) 2: P := (0, 1, . . . , n -1) 3:

L := () 4:

for i from 0 to k -1 do 5:

append(L, P [F [n -1 -i]]) 6: remove(P, F [n -1 -i]) 7:
return L Since the factoradic (with 8 components) of 2020 is (0, 0, 2, 0, 4, 4, 2, 0), the 8-permutation of rank 2020 is (0, 3, 6, 7, 1, 5, 2, 4).

Note, in the function Extract, the way the data structure P is handled in memory is very important for efficiency reasons. The best way is probably to build a dynamic balanced tree as presented in [START_REF] Bonet | Efficient algorithms to rank and unrank permutations in lexicographic order[END_REF], or a multiset (whose elements have weight 1 or 0) as presented in the appendix of [START_REF] Bodini | A Quantitative Study of Pure Parallel Processes[END_REF]. It seems there does not exist any algorithm based on some swap operation giving an in-place shuffle to unrank permutation in the lexicographic order: in fact Durstenfeld's algorithm [START_REF] Durstenfeld | Algorithm 235: Random permutation[END_REF] cannot be easily adapted for the lexicographic order.

Let us now turn to the unranking of a combination through factoradics. The basic ideas driving our algorithm are the following:

1. we define a bijection among the combinations of k elements among n and a subset of the permutations of n elements; 2. we transform the combination rank u into the rank u of the appropriate permutation; 3. we build the k-prefix of the permutation of rank u by using the Algorithm 4. Definition 5. Let n and k be two integers with 0 ≤ k ≤ n. We define the application P such that to the combination ( 0 , 1 , . . . , k-1 ) it associates the n-permutation

( 0 , 1 , . . . , k-1 , d k , . . . d n-1 ) such that d k < d k+1 < • • • < d n-1 .
Thus, by definition, for n = 5 and k = 3, the permutations associated to the combinations (0, 1, 2) and (0, 2, 4) are respectively (0, 1, 2, 3, 4) and (0, 2, 4, 1, 3). Proposition 2. The integers n and k being given, the application P transforming a combination into a permutation is one-to-one.

Remark that the 5-permutation (0, 1, 2, 3, 4) is the permutation associated to both combinations (0, 1) and (0, 1, 2) but in these examples the values of k are distinct. In fact, there are exactly 6 combinations associated to the latter permutation: P depends in n and k, it is a bijection when n and k are given. The proof is obvious.

The algorithm implemented in the function Rank_Conversion (see Algorithm 5) transforms the rank of a combination for n, k into the rank of its associated permutation through the computation of its factoradic. The main idea of the rank conversion is similar to the one for determining the combinadic through trials. It consists in building the (reverse of the) factoradic by detailing each component by subtracting the permutation ranks that must be omitted, when the component is still too small. Proposition 3. The function Unranking_Factoradic(n, k, •) computes the combinations for n, k in the lexicographic order. During the conversion from the rank of the combination to the one of the associated permutation, the coefficients are obtained via trials for f n-1 to f n-k , remarking that through our bijection P the latter sequence is weakly increasing. Thus the worst cases are obtained when the value f n-k is as large as possible, that is n -k. Thus for such elements, the number of calls to binomial is n. For the average number of calls to the function binomial, unranking all combinations u when it describes the whole range from 0 to n k -1 introduces again the following cumulative sequence. Lemma 4. Let u n,k be the cumulative numbers of calls to binomial while unranking all possible u from 0 to n k -1. The sequence satisfies: u n,k = 0 and u n,n+i = 0 for all n and i > 0 and otherwise

u n,k = n k + u n-1,k-1 + u n-1,k .
Aside from the extreme cases, the recursive formula is identical to the previous ones. Here again the proof is obtained through the same kind of proof than the ones of Section 2. See the Appendix A.3 for key ideas for the proof.

In Fig. 6 we exhibit the first values for u n,k . Here again we obtain a sequence stored under the reference OEIS A127717. The bijection between both structures is direct, and thus we have new information about this sequence in the following.

Theorem 7. Let U (z, y) be the generating functions associated to (u n,k ). Then

U (z, y) = 1 1 -z -zy 1 1 -z -zy - 1 1 -z ; u n,k = n k n - n -k k + 1 .
The proof is similar to the one of Theorem 1.

Corollary 5. The average number of calls to binomial in Algorithm 5 for n being large and k being of the form αn for α ∈]0, 1[ is

u n,k n k = n→∞ k=αn n + 1 - 1 α + O 1 n .
Comparing the latter result with the asymptotic behavior for the second algorithm based on combinadics, even if the dominant parts are equal, our approach seems not as efficient. But again, our complexity measure assumes the binomial coefficients being pre-computed.

Experimental comparison and algorithm improvements

For all algorithms, n and k being given, we proved that the average number of calls to the function binomial is of the same order, equivalent to n. In Fig. 7 we have represented the following statistics 7 . We take n = 10000 and in the abscissa-axis we let the value k ranges from 25 to 9975 with an iteration step of 50. In the ordinate-axis we have represented the average time in ms to unrank one combination (for each step we sampled, uniformly at random, 100 distinct combinations and we computed the average time for one). We have implemented all algorithms in C++ using the classical GMP library for big integers. The time needed by the algorithms presented in the previous sections are depicted in dashed lines: the red color is for the recursive Algorithm 1, blue corresponds to Algorithm 2, brown corresponds to its improvement, Algorithm 3, and finally green is for our factoradic Algorithm 5. We first remark that Algorithm 2 is worse than the 3 others whose computation times are very close. Further, asymptotically, the average number of calls to binomial is linear in n, whatever the value of k, this characteristics is not reflected in the experiments.

Usually the way for evaluating n k is the following left to right computation: n

•(n-1)/2•(n-2)/3 • • • (n-k +1)/k. Thus it needs Θ(k) arithmetic operations.
It is done this way in GMP. There is a classical improvement used for iterated binomial evaluations. We compute the first coefficient we need, and then we correct this value to obtain the second one that is then modified to obtain the third one and so on. By studying all the previous algorithms we remark that the successive binomial coefficients are such that their parameters are only shifted by ±1. Thus obtaining the next binomial coefficient is directly obtained with a multiplication and a division.

b • (n -m + 1 -k + i)/(n -m) if n = m else 1.
Doing these improvements for all algorithms, we finally obtained the solid curves in Fig. 7. Obviously all algorithms are faster than their first version, but Algorithm 2 is the less improved (in particular when k > n/2). The order of the computations is really penalizing in this approach. By focusing only on the new versions of Algorithms 1, 3 and 5, we obtain Fig. 1. We thus exhibit that our Algorithm 5, in green, is the most efficient while k ranges the interval. where L(n) is a function in log 2 n depending on the algorithmic used for the integer multiplication.

The theorem is a by product of Lemma 5 saying that the average number of calls to binomial is equivalent to n that is amended to obtain the average number of multiplications or divisions being equivalent to Θ(n) through the previous improvements and the result about the bit complexity for the computation of n! (cf. [4, p. 270]). In Fig. 8 we merge the time complexity of our algorithm (in green) with a normalized curve (in black) from Theorem 8. and f n-k . In this case (in front of Lemma 2) we need weak compositions since the sequence f n-k , f n-k+1 , . . . , f n-1 is weakly decreasing. Furthermore for each weak composition determining a combination, we need (f n-k + 1 + k) calls to binomial to determine the combination. We thus obtain:

u n,k = n-k f n-k =0 (f n-k + k) • f n-k + k -1 k -1 .
An induction completes the proof.

A.4 Optimized factoradics algorithm

In this following new version of Algorithm 5 we present 4 improvements:

-The algorithm is not symmetrical anymore. In fact, we observe that when k > 2 then the previous version is faster than when k < 2. Thus when unranking in the range 0..n/2, with some suitable substitutions we perform the unranking in the other range. As a consequence we get a small stall in this curve around k = n/2 in Fig. 1.

Let C be a combination for k < n/2, then its complement C is somehow computed, but we store C and not C (lines 5-9 and 23-26 and 39-44). -Only the first binomial coefficient is computed. Then it is adjusted to obtain the right value (lines 20 and 36). -The last component is computed without any trial, just as a consequence of the other components (line 49). -We avoid the call to Extract. In fact, when unranking the first part of a combination, the sequence of values we extract is an increasing sequence, we thus can completely avoid the expensive call to Extract just be saving how many values have already been extracted (lines 28 and 46 when k > n/2 and the more complicated cases lines 23-25 and lines 39-44).

1 :

 1 function Unranking_Rec(n, k, u) 2:

Fig. 2 :

 2 Fig. 2: First values of u n,k for n = 1..6 and k = 0..n The proof of Theorem 1 is given in Appendix A.1This sequence is a shifted version of the sequence stored under the reference OEIS A0597973 . We thus can complete the properties in OEIS using our results.Due to the values of the extreme cases when k = 0 and k = n and the symmetry in the recurrence we obviously obtain the fact that u n,k = u n,n-k , reflecting the symmetry of the binomial coefficients.

Corollary 2 .

 2 The function Unranking_Rec(n, k, •) needs in average u n,k / n k calls to the function binomial. For n being large and k being of the form α • n for α ∈]0, 1[, we get

Fig. 4 :

 4 Fig. 4: First values of u n,k for n = 1..6 and k = 0..n (left: Algorithm 2; right: Algorithm 3)

Algorithm 5 a combination 1 : 1 :

 511 length(F ): computes the number of components in F ; remove(F, i): removes from F the element at index i; composition(F ): computes the integer whose factoradic is the reverse of F . Unranking function Unranking_Factoradic(n, k, u) 2: u := Rank_Conversion(n, k, u) 3: F := Unranking_Permutation(u , n) 4: return Extract(F, n, k) function Rank_Conversion(n, k, u) 2: D := (0, . . . , 0) n components in D 3: m := 0; i := 0 4: while i < k do 5: b := binomial(n -1 -m, k -1 -i) 6: if b > u then 7: D[i] = m 8: i := i + 1 9: n := n -1 10: else 11: u := u -b 12: m := m + 1 13: return composition(D) D contains the reverse factoradic tuple

Fig. 5 :Fact 6

 56 Fig. 5: Algorithms for the unranking of combination through factoradics

Fig. 6 :

 6 Fig. 6: First values of u n,k for n = 1..6 and k = 0..n

Fig. 7 :

 7 Fig. 7: Time (in ms) for unranking a combination, with n = 10000 and k = 0..n

  For example, in Algorithm 5, we insert between line 3 and 4 the instruction b := binomial(n, k), then between line 9 and 10 we put b := b • (k -i)/(n -m) if n = m else 1, and finally, between line 12 and 13 we put b :=

Fig. 8 :Theorem 8 .

 88 Fig. 8: Merge of Algorithm 5 and theoretical complexity

Throughout this paper, a reference OEIS A• • • points to Sloane's Online Encyclopedia of Integer Sequences www.oeis.org.

We use several times the same writings u n,k and U (z, y) for distinct algorithms, but they rely on the same definition.

The experiments have been driven through a standard laptop PC with an I7-8665U CPU, 32Gb RAM running Ubuntu Linux.

A Appendix

A.1 Appendix dedicated to Section 1 Proof (of Theorem1). The first step of the proof consists in exhibiting the ordinary generating function associated to U (z, y). In order to obtain the equation for U , we start from the result presented in Lemma 1. The extreme cases areu n,0 = 0 and u n,n+i = 0 for all n and i ≥ 0. And the recursive equation is

We remark the constant n k in the equation. We thus need the bivariate generating function for binomial coefficient. Le us denote it by B(z, y); it satisfies:

In order to follow the extreme cases, we must remove the first column k = 0 and the diagonal k = n:

By summing the recursive equation by taking care of the extreme cases we have:

We thus deduce

The second step in the proof consists in extracting the coefficient u n,k .

By extraction the coefficient in front of z n :

The latter result correspond to the distribution when k ranges from 0 to n. But we can further extract the coefficient of z n y k :

Using Hockey-Stick identity we obtain

Thus we conclude

The result is proved.

A.2 Appendix dedicated to Section 2

Proof (ideas of Lemma 2). Once c 1 is given, all possible compositions of (n -c 1 ) in k parts are possible for the values for c 2 , . . . , c k , and for each composition, we need (n -c 1 ) calls to binomial. Furthermore, for each rank we computes its reverse, with one call to binomial. Thus

Using the latter equation, the recursive equation is directly proved by induction.

A.3 Appendix dedicated to Section 3

To prove the correction of our Algorithm 5 and also its average complexity, we rely on the following lemma.

Lemma 5. The factoradics of the combinations k elements among n are all the tuples (0, . . . , 0,

From the later lemma we deduce the proof of Proposition 3 and Lemma 4. In fact, the rank conversion consists in enumerating the factoradics corresponding to combinations, thus avoiding permutations that do not represent combinations.

In fact we can directly omit them by using their numbers given by binomial coefficients.

Proof (ideas for Lemma 4). Here the last coefficient to be determined is f n-k , it is the greatest one. Once it is given, all other values f n-k+1 , . . . , f n-1 are just determined through a weak composition of the available integers between 0 if B > r then 16:

if i < k -2 then 17:

if n -1 -m = 0 then 18:

B := 1 19: else 20:

if inverse then 23:

for e from m2 to m + i do 24:

D return D