# Variable Conversion and Degree Reduction The Amplify SDK allows you to create models that include real and integer variables and polynomials of any degree. On the other hand, combinatorial optimization solvers are generally limited in the types of variables and degrees of polynomials they can handle. This section describes the variable conversion and degree reduction methods that the Amplify SDK implements. These methods match the model with the variable types and the polynomial degree that the solver can handle. ## Retrieving an intermediate model Using the {py:meth}`~amplify.Model.to_intermediate_model` method on a model returns an intermediate model transformed to match the specified variable type and polynomial degree. When it builds the intermediate model, the {py:func}`~amplify.solve` function uses this method. So it is helpful to see how the Amplify SDK performs variable conversion and degree reduction. ```{attention} Usually, there is no need to call the {py:meth}`~amplify.Model.to_intermediate_model` method. Use it to see how the Amplify SDK performs variable transformations and order reduction within the {py:func}`~amplify.solve` function. ``` You can specify the intermediate model's variable type and polynomial degree in the first argument of the {py:meth}`~amplify.Model.to_intermediate_model` method. If you want to use the degree that the solver can handle, pass the {py:attr}`~amplify.AmplifyAEClient.acceptable_degrees` property as follows. ```{testcode} from amplify import VariableGenerator, Model, AmplifyAEClient, AcceptableDegrees gen = VariableGenerator() q = gen.array("Binary", 3) model = Model(q[0] * q[1] * q[2]) client = AmplifyAEClient() client.parameters.time_limit_ms = 1000 client.solver = "Qubo" im, mapping = model.to_intermediate_model(client.acceptable_degrees) ``` The {py:meth}`~amplify.Model.to_intermediate_model` method returns an intermediate model and a map of variable conversions from the input model to the intermediate model. For example, QUBO mode of {py:class}`~amplify.AmplifyAEClient` can handle a quadratic polynomial of binary variables as an objective function. See "[](#create-intermediate-model)." In the above example, the input model is a cubic polynomial of binary variables. However, the {py:meth}`~amplify.Model.to_intermediate_model` method converts the intermediate model's objective function to a quadratic polynomial of binary variables. The method imposes a [{py:class}`~amplify.QuadratizationMethod.Substitute` constraint](#substitution-multiplier) on the intermediate model. ```{doctest} >>> print(im) minimize: q_2 q'_0 subject to: q_0 q_1 - q'_0 == 0 (weight: 1) ``` ```{note} The variables in the intermediate model use a different {py:class}`~amplify.VariableGenerator` than the input variables. In other words, it is impossible to operate between variables in the input model and the intermediate model. ``` Pass an instance of the {py:class}`~amplify.AcceptableDegrees` class to the {py:meth}`~amplify.Model.to_intermediate_model` method. You can then freely specify the type and degree of the variable to be converted. ```{testcode} bq = AcceptableDegrees(objective={"Binary": "Quadratic"}) im, mapping = model.to_intermediate_model(bq) ``` In addition to specifying the variable type and degree, the {py:meth}`~amplify.Model.to_intermediate_model` method can specify variable conversion and degree reduction algorithms. The following explains the variable conversion and order reduction methods implemented in the Amplify SDK by looking at the {py:meth}`~amplify.Model.to_intermediate_model` method results. ## Variable conversion The Amplify SDK converts a variable in an input model to another variable type. This conversion transforms the model into an intermediate model with variable types that the solver can handle. Currently, the following variable conversions are implemented. * Conversion from integer to binary variables * Conversion from real to binary variables * Conversion from binary to Ising variables and vice versa (encode-integer)= ### Conversion from integer to binary variables For an integer variable $n$ with a lower limit of $l$ and an upper limit of $u$, $k$ binary variables $q_0, q_1, ..., q_{k-1}$ are newly issued. Now, $n$ is converted to a polynomial: $$ a_0 q_0 + a_1 q_1 + \cdots + a_{k-1} q_{k-1} + l. $$ Here, the integer sequence $a_0, a_1, \ldots, a_{k-1}$ is determined so that the possible range of this polynomial matches the set of all integers between $l$ and $u$. Four different algorithms are implemented to determine the number of variables $k$ and the sequence of coefficients $a_0, a_1, \ldots, a_{k-1}$ of the polynomial. As an example, consider a model whose objective function is the integer variable `n`{l=python}, which takes integers between -10 and 10. Suppose the intermediate model is composed of binary variables as follows. Let us see how the integer variable `n`{l=python} is transformed into a polynomial of binary variables in each algorithm. ```{testcode} gen = VariableGenerator() n = gen.scalar("Integer", bounds=(-10, 10)) # Issue an integer variable model = Model(n) bq = AcceptableDegrees(objective={"Binary": "Quadratic"}) ``` The conversion algorithm from integer variables to binary variables is specified in the `integer_encoding_method` keyword argument of the {py:meth}`~amplify.Model.to_intermediate_model` method and the {py:func}`~amplify.solve` function. The default is {py:class}`~amplify.IntegerEncodingMethod.Default`. {py:class}`~amplify.IntegerEncodingMethod.Unary` : Sets the number of variables $k = u - l$, and $a_0 = a_1 = \cdots = a_{k-1} = 1$. ```{testcode} :hide: gen = VariableGenerator(); n = gen.scalar("Integer", bounds=(-10, 10)); model = Model(n) ``` ```{testcode} im, mapping = model.to_intermediate_model(bq, integer_encoding_method="Unary") ``` ```{doctest} >>> print(mapping[n]) # doctest: +NORMALIZE_WHITESPACE q_0 + q_1 + q_2 + q_3 + q_4 + q_5 + q_6 + q_7 + q_8 + q_9 + q_{10} + q_{11} + q_{12} + q_{13} + q_{14} + q_{15} + q_{16} + q_{17} + q_{18} + q_{19} - 10 ``` {py:class}`~amplify.IntegerEncodingMethod.Linear` : Assigns $(a_0, a_1, a_2, \ldots, a_{k-2}, a_{k-1}, a_k) = (1, 2, 3, \ldots, k-1, k, m)$. $k$ is the maximum integer satisfying $1 + 2 + 3 + \cdots + k \leq u - l$, and $m$ are the fraction. ```{testcode} :hide: gen = VariableGenerator(); n = gen.scalar("Integer", bounds=(-10, 10)); model = Model(n) ``` ```{testcode} im, mapping = model.to_intermediate_model(bq, integer_encoding_method="Linear") ``` ```{doctest} >>> print(mapping[n]) q_0 + 2 q_1 + 3 q_2 + 4 q_3 + 5 q_4 + 5 q_5 - 10 ``` {py:class}`~amplify.IntegerEncodingMethod.Binary` : Assigns $(a_0, a_1, a_2, a_3, \ldots, a_{k-1}, a_k) = (1, 2, 4, 8, \ldots, 2^{k-1}, m)$. $k$ is the maximum integer satisfying $1 + 2 + 4 + \cdots + 2^{k-1} \leq u - l$, and $m$ is the fraction. ```{testcode} :hide: gen = VariableGenerator(); n = gen.scalar("Integer", bounds=(-10, 10)); model = Model(n) ``` ```{testcode} im, mapping = model.to_intermediate_model(bq, integer_encoding_method="Binary") ``` ```{doctest} >>> print(mapping[n]) q_0 + 2 q_1 + 4 q_2 + 8 q_3 + 5 q_4 - 10 ``` {py:class}`~amplify.IntegerEncodingMethod.Default` (Default) : Of the above three algorithms, the one that uses the least number of variables $k$ is used. ```{note} In {py:class}`~amplify.IntegerEncodingMethod.Linear` and {py:class}`~amplify.IntegerEncodingMethod.Binary`, if the fractions $m$, $m_1$, $m_2$ become $0$, the corresponding auxiliary variable is omitted from being issued. ``` (encode-real)= ### Conversion from real to binary variables For a real variable $x$ with a lower limit of $l$ and an upper limit of $u$, newly create $k$ binary variables $q_0, q_1, ..., q_{k-1}$. Then, $x$ is converted to a polynomial: $$ a_0 q_0 + a_1 q_1 + \cdots + a_{k-1} q_{k-1} + l, $$ where the coefficients $a_0, a_1, \ldots, a_{k-1}$ are positive real numbers satisfying $a_0 + a_1 + \cdots a_{k-1} = u - l$. A Randomized algorithm [^5] is implemented to determine the coefficients of the polynomial $a_0, a_1, \ldots, a_{k-1}$. For example, consider the model with the real variable `x`{l=python} as the objective function, taking reals between -10 and 10. It is converted into a polynomial of binary variables as follows. Let us see how the real variable `x`{l=python} is transformed into a polynomial of binary variables in each algorithm. ```{testcode} gen = VariableGenerator() x = gen.scalar("Real", bounds=(-10, 10)) model = Model(x) bq = AcceptableDegrees(objective={"Binary": "Quadratic"}) ``` The conversion algorithm from real variables to binary variables is specified in the `real_encoding_method` keyword argument of the {py:meth}`~amplify.Model.to_intermediate_model` method and the {py:func}`~amplify.solve` function. The default is {py:class}`~amplify.RealEncodingMethod.Random16`. {py:class}`~amplify.RealEncodingMethod.Random` : Set $(a_0, a_1, a_2, a_3, \ldots, a_{k-1})$ by uniform random numbers. The number of binary variables $k$ is one of (4, 8, 16, 32), which is extracted from the last digits of the conversion algorithm name. The following is an example for $k = 16$. ```{testcode} :hide: import amplify gen = VariableGenerator(); x = gen.scalar("Real", bounds=(-10, 10)); model = Model(x) ``` ```{testcode} amplify.set_seed(0) im, mapping = model.to_intermediate_model(bq, real_encoding_method="Random16") ``` ```{doctest} >>> print(mapping[x]) 1.99333941505556 q_0 + 1.05184485013684 q_1 + 1.18505654902564 q_2 + 1.49764241671521 q_3 + 1.76252706549575 q_4 + 0.841695701394766 q_5 + 0.21960146400976 q_6 + 0.635503335788004 q_7 + 1.80253453941945 q_8 + 0.173438994097798 q_9 + 2.04763527252344 q_{10} + 0.611762636976138 q_{11} + 1.95487069012893 q_{12} + 1.62228108547751 q_{13} + 1.9691590736731 q_{14} + 0.631106910082103 q_{15} - 10 ``` ```{note} You can set the seed of random numbers used for variable conversion with the {py:func}`~amplify.set_seed` function. ``` ```{note} In {py:class}`~amplify.Poly`, terms with coefficients less than `1e-10`{l=python} are ignored to deal with numerical errors in symbolic processing. Therefore, when the difference between the upper and lower limits $u - l$ is very small, some terms of the transformed binary variable may be lost. ``` ### Conversion from Ising to binary variables For an Ising variable $s$, issue a new binary variable $q$ to convert $s$ to $2 q -1$. ```{testcode} gen = VariableGenerator() s = gen.scalar("Ising") model = Model(s) bq = AcceptableDegrees(objective={"Binary": "Quadratic"}) im, mapping = model.to_intermediate_model(bq) ``` ```{doctest} >>> print(mapping[s]) 2 q_0 - 1 ``` ### Conversion from binary to Ising variables For a binary variable $q$, issue a new Ising variable $s$ to convert $q$ to $(s + 1) / 2$. ```{testcode} gen = VariableGenerator() q = gen.scalar("Binary") model = Model(q) iq = AcceptableDegrees(objective={"Ising": "Quadratic"}) im, mapping = model.to_intermediate_model(iq) ``` ```{doctest} >>> print(mapping[q]) 0.5 s_0 + 0.5 ``` ## Degree reduction of binary variable polynomials The Amplify SDK provides functionality to reduce the degrees of polynomials of binary or Ising variables to meet solver requirements. By applying variable conversion and degree reduction, the Amplify SDK can transform polynomials involving integer variables into lower-order polynomials with binary variables. (quadratization-method)= ### Degree reduction algorithm Two algorithms are implemented: {py:class}`~amplify.QuadratizationMethod.Substitute` [^1] and {py:class}`~amplify.QuadratizationMethod.IshikawaKZFD` [^2][^3][^4]. The degree reduction algorithm can be specified with the `quadratization_method` keyword argument of the {py:func}`~amplify.solve` function and the {py:meth}`~amplify.Model.to_intermediate_model` method. The default is {py:class}`~amplify.QuadratizationMethod.Substitute`. To see how the intermediate model is constructed as the degree is reduced in each algorithm, we create a model. Its objective function is the following third-order binary variable polynomial. The intermediate model is also specified so that the objective function is binary quadratic and unconstrained. ```{testcode} gen = VariableGenerator() q = gen.array("Binary", 4) model = Model(q[0] * q[2] * q[3] - q[1] * q[2] * q[3]) bq = AcceptableDegrees(objective={"Binary": "Quadratic"}) ``` (substitution-multiplier)= {py:class}`~amplify.QuadratizationMethod.Substitute` (Default) : For a product $q_a q_b$ of two binary variables that appears in a product of third order or higher, issue the auxiliary variable $x$. Replace all $q_a q_b$ in each term of the polynomial expression with $x$. The following constraints are then added to the intermediate model. $$ \begin{aligned} \text{subject to:} \quad & q_a q_b = x \\ \text{penalty function:} \quad & q_a q_b - 2 q_a x - 2 q_b x + 3 x \end{aligned} $$ The above operations are repeated until the polynomial becomes the desired order. ```{testcode} :hide: gen = VariableGenerator(); q = gen.array("Binary", 4); model = Model(q[0] * q[2] * q[3] - q[1] * q[2] * q[3]) ``` ```{testcode} im, mapping = model.to_intermediate_model( bq, quadratization_method="Substitute" ) ``` ```{doctest} >>> print(im) minimize: q_0 q'_0 - q_1 q'_0 subject to: q_2 q_3 - q'_0 == 0 (weight: 2) ``` In the above, the auxiliary variable `q'_0` replaces `q_0 q_1`. The order reduction adds constraints. The weights of their penalty function are initially set to the sum of the absolute values of the coefficients of the substituted terms. This is because the weights must be large enough compared to the objective function. The initial values of the weights are sometimes too large. Then the `substitution_multiplier` keyword argument can specify a coefficient for these initial values. ```{testcode} :hide: gen = VariableGenerator(); q = gen.array("Binary", 4); model = Model(q[0] * q[2] * q[3] - q[1] * q[2] * q[3]) ``` ```{testcode} im, mapping = model.to_intermediate_model(bq, quadratization_method="Substitute", substitution_multiplier=0.5 ) ``` ```{doctest} >>> print(im) minimize: q_0 q'_0 - q_1 q'_0 subject to: q_2 q_3 - q'_0 == 0 (weight: 1) ``` ```{note} See "[](#penalty-weight)" for information on adjusting constraint weights. ``` {py:class}`~amplify.QuadratizationMethod.IshikawaKZFD` : For each higher-order term of the polynomial, an auxiliary variable is issued. The degree of a higher-order term is reduced to a quadratic polynomial. The value does not change when the auxiliary variable takes the optimal value. The following conversions are performed for third-order or higher terms with binary variables $k q_1 q_2 \cdots q_n$. * If the coefficient $k$ of the term is negative, an auxiliary variable $x$ is issued and the term $k q_1 q_2 \cdots q_n$ is converted to: $$ kx (q_1 + q_2 + \cdots + q_n + 1 - n). $$ * If the coefficient $k$ of the term is positive and the degree $n$ is even, issue $(n-2)/2$ auxiliary variables $x_1, x_2, \ldots, x_{\left(n-2 \right)/2}$. These convert the term $k q_1 q_2 \cdots q_n$ to: $$ \begin{aligned} & k \left(\frac{1}{2} S\left(S-1\right) - \sum_{i = 1}^{\frac{n-2}{2}} x_i \left(2 \left(S - 2i\right) + 1\right) \right) \\ & \text{where} \quad S = q_1 + q_2 + \cdots + q_n. \end{aligned} $$ * If the coefficient of the term $k$ is positive and the degree $n$ is odd, issue $(n-1)/2$ auxiliary variables $x_1, x_2, \ldots, x_{\left(n-1\right)/2}$. These convert the term $k q_1 q_2 \cdots q_n$ to: $$ \begin{aligned} & k k \left(\frac{1}{2} S\left(S-1\right) - \sum_{i = 1}^{\frac{n-1}{2} - 1} x_i \left(2 \left(S - 2i\right) + 1\right) - x_{\frac{n-1}{2}} \left(S - n + 2\right) \right) \\ & \text{where} \quad S = q_1 + q_2 + \cdots + q_n. \end{aligned} $$ ```{testcode} :hide: gen = VariableGenerator(); q = gen.array("Binary", 4); model = Model(q[0] * q[2] * q[3] - q[1] * q[2] * q[3]) ``` ```{testcode} im, mapping = model.to_intermediate_model( bq, quadratization_method="IshikawaKZFD" ) ``` ```{doctest} >>> print(im) # doctest: +NORMALIZE_WHITESPACE minimize: q_0 q_2 + q_0 q_3 - q_0 q'_0 - q_1 q'_1 + q_2 q_3 - q_2 q'_0 - q_2 q'_1 - q_3 q'_0 - q_3 q'_1 + q'_0 + 2 q'_1 ``` Here, `q'_0` and `q'_1` are the auxiliary variables issued for each term. ```{attention} In the {py:class}`~amplify.QuadratizationMethod.IshikawaKZFD` algorithm, the degree of the constraints cannot be reduced. If you want to pass a constraint with a third-order constraint expression to a solver that accepts a second-order constraint expression, use the {py:class}`~amplify.QuadratizationMethod.Substitute` algorithm. ``` ```{dropdown} For Ising variables The following conversions are performed for third-order or higher terms with Ising variables $k s_1 s_2 \cdots s_n$. * If the coefficient of the term $k$ is positive and the degree $n$ is odd, or coefficient $k$ is negative and the degree $n$ is even, then $\mathrm{floor} \left( n / 2 \right)$ number auxiliary variables $x_1, x_2, \ldots, x_{\mathrm{floor} \left( n / 2 \right)}$ are issued to convert the term $k s_1 s_2 \cdots s_n$ to: $$ \begin{aligned} & \left| k \right| \left( 2 S^2 - 8 \sum_{i = 1}^{\mathrm{floor} \left( n / 2 \right)}{ \left( \frac{x_i + 1}{2} \right) \left( S - 2 i + 1 \right) } - 1 \right) \\ & \text{where} \quad S = \sum_{i =1}^{n}{\frac{s_i + 1}{2}}. \end{aligned} $$ * If the coefficient of the term $k$ is positive and the degree $n$ is even, or coefficient $k$ is negative and the degree $n$ is odd, then $\mathrm{floor} \left( \left(n-1\right) / 2 \right)$ number auxiliary variables $x_1, x_2, \ldots, x_{\mathrm{floor} \left( \left(n-1\right) / 2 \right)}$ are issued to convert the term $k s_1 s_2 \cdots s_n$ to: $$ \begin{aligned} & \left| k \right| \left( 2 \left(S-1\right)^2 - 8 \sum_{i = 1}^{\mathrm{floor} \left( \left(n-1\right) / 2 \right)}{ \left( \frac{x_i + 1}{2} \right) \left( S - 2 i \right) } - 1 \right) \\ & \text{where} \quad S = \sum_{i =1}^{n}{\frac{s_i + 1}{2}}. \end{aligned} $$ ``` ### The characteristics of each algorithm Each degree reduction algorithm has the following characteristics. Depending on which algorithm is used, the number of variables in the intermediate model can vary greatly, which can affect the solver's feasibility and results. ```{list-table} :width: 100% :header-rows: 1 :stub-columns: 1 :widths: 1 2 2 * - Algorithm - {py:class}`~amplify.QuadratizationMethod.Substitute` - {py:class}`~amplify.QuadratizationMethod.IshikawaKZFD` * - Variables of the polynomial reducible to lower degree - Binary variables only - Binary or Ising variables * - Degree of Polynomial After Reduction - Arbitrary degree, although quadratic equality constraints or penalties may be added - Always quadratic or lower * - Reduction of Constraint Degree - Possible - Not possible; penalties must be reduced instead * - Number of Auxiliary Variables - The number of auxiliary variables can be smaller compared to Ishikawa-KZFD when cross-term substitution is effective - Roughly proportional to half the degree multiplied by the number of terms (proportional to the number of terms if polynomial coefficients are negative) ``` (intermediate-result)= ## Solver execution result and intermediate model The result of the {py:func}`~amplify.solve` function execution can also provide information about the intermediate model and the result of solving for the intermediate model. You can obtain this information from the {py:attr}`~amplify.Result.intermediate` attribute of the {py:class}`~amplify.Result` class (the {py:class}`~amplify.Result.ModelConversion` class). The {py:class}`~amplify.Result.ModelConversion` class has the following attributes. ```{list-table} :width: 100% :header-rows: 1 * - Attribute - Data type - Details * - {py:attr}`~amplify.Result.ModelConversion.model` - {py:class}`~amplify.Model` - Intermetiate model * - {py:attr}`~amplify.Result.ModelConversion.mapping` - {py:class}`~amplify.Result.ModelConversion.IntermediateMapping` - Variable transformation map from input model variables to intermediate model variables * - {py:attr}`~amplify.Result.ModelConversion.num_variables` - {py:class}`int` - Number of variables used in the intermediate model's objective function and constraints (including penalty functions, if generated) * - {py:attr}`~amplify.Result.ModelConversion.values_list` - {py:class}`~amplify.Result.ValuesList` - Result of the intermediate model ``` As an example, let's execute {py:func}`~amplify.solve` with {py:class}`~amplify.AmplifyAEClient` for the following model with integer variables in the objective function. ```python from amplify import VariableGenerator, Model, AmplifyAEClient, solve gen = VariableGenerator() n = gen.scalar("Integer", bounds=(1, 3)) # Issue an integer variable model = Model(n) client = AmplifyAEClient() # client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" client.parameters.time_limit_ms = 1000 result = solve(model, client) ``` ```{testcode} :hide: import os from amplify import VariableGenerator, Model, AmplifyAEClient, solve gen = VariableGenerator() n = gen.scalar("Integer", bounds=(1, 3)) # Issue an integer variable model = Model(n) client = AmplifyAEClient() client.url = os.getenv("TEST_AE_URL", "") client.token = os.getenv("TEST_AE_TOKEN", "") client.parameters.time_limit_ms = 1000 result = solve(model, client) ``` You can obtain information about the intermediate model from the {py:attr}`~amplify.Result.intermediate` attribute of the result. ```{doctest} >>> print(result.intermediate.model) minimize: q_0 + q_1 + 1 ``` You can see how the integer variable `n`{l=python} is transformed in {py:attr}`~amplify.Result.ModelConversion.mapping` as follows. ```{doctest} >>> print(result.intermediate.mapping[n]) q_0 + q_1 + 1 ``` You can check the solution of the intermediate model with {py:attr}`~amplify.Result.ModelConversion.values_list`. {py:attr}`~amplify.Result.ModelConversion.values_list` can contain multiple solutions, so we obtain the first solution here. ```{doctest} >>> result.intermediate.values_list[0] Values({Poly(q_0): 0, Poly(q_1): 0}) ``` The solution for the input model is obtained from the conversion of the solution of the intermediate model based on {py:attr}`~amplify.Result.ModelConversion.mapping`. Although a bit complicated, all intermediate model solutions can be transformed into input model solutions as follows. ```{doctest} >>> im = result.intermediate >>> { ... var: im_mapped.evaluate(im_values) ... for im_values in im.values_list ... for var, im_mapped in im.mapping.items() ... } {Poly(n_0): 1.0} ``` We can see that this is consistent with the {py:attr}`~amplify.Result.solutions` that are the output of {py:func}`~amplify.solve`. [^1]: I. G. Rosenberg, “Reduction of Bivalent Maximization to the Quadratic Case,” Cahiers du Centre d’Etudes de Recherche Operationnelle 17, 71–74 (1975). [^2]: [V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?” IEEE Transactions on Pattern Analysis and Machine Intelligence 26, 147–159 (2004).](https://doi.org/10.1109/TPAMI.2004.1262177) [^3]: [D. Freedman and P. Drineas, “Energy Minimization via Graph Cuts: Settling What is Possible,” in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), Vol. 2 (IEEE, 2005) pp. 939–946.](https://doi.org/10.1109/CVPR.2005.143) [^4]: [H Ishikawa, “Transformation of General Binary MRF Minimization to the First-Order Case,” IEEE Transactions on Pattern Analysis and Machine Intelligence 33, 1234–1249 (2011).](https://doi.org/10.1109/TPAMI.2010.91) [^5]: [K. Endo, Y. Matsuda, S. Tanaka, and M. Muramatsu, “Novel real number representations in Ising machines and performance evaluation: Combinatorial random number sum and constant division,” PLOS ONE 19, e0304594 (2024).](https://doi.org/10.1371/journal.pone.0304594)