EVALUATING QUANTUM COMPUTING CIRCUITS IN VIEW OF THE RESOURCE COSTS OF A QUANTUM ALGORITHM

Abstract

Methods for evaluating quantum computing circuits in view of the resource costs of a quantum algorithm are described. A processor-implemented method for performing an evaluation of a polynomial corresponding to an input is provided. The method includes determining a polynomial interpolation for a set of sub-intervals corresponding to the input. The method further includes constructing a quantum circuit for performing, in parallel, polynomial evaluation corresponding to each of the set of sub-intervals.

Description

BACKGROUND

Quantum computing may solve certain problems much faster than classical devices. Examples of such problems include Shor's algorithm, unstructured search problems, and simulation of quantum mechanical systems. The crossover point at which the choice of quantum algorithm versus classical algorithm changes is based on the implementation details of the functions being implemented by the quantum circuits corresponding to the quantum algorithm.

SUMMARY

In one example, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include determining a polynomial interpolation for a set of sub-intervals corresponding to the input. The method may further include constructing a quantum circuit for performing, in parallel, polynomial evaluation corresponding to each of the set of sub-intervals.

In another example, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include determining a polynomial interpolation of an initial degree d for an initial set of sub-intervals corresponding to an input. The method may further include determining whether a required target precision is achievable based on the initial degree d and the set of sub-intervals. The method may further include if the target precision is achievable, then constructing a quantum circuit for performing, in parallel, polynomial evaluation corresponding to each of the set of sub-intervals. The method may further include if the target precision is not achievable, then either incrementing the initial degree d or subdividing the initial set of subintervals into additional intervals until the target precision is achievable and then constructing the quantum circuit.

In yet another example, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include for each interval of a set of sub-intervals corresponding to the input, determining whether to perform the evaluation of the polynomial in a first mode or a second mode. The method may further include, in the first mode, performing the evaluation of the polynomial by: determining a polynomial interpolation of at most degree d for a first subset of the set of sub-intervals corresponding to the input, and constructing a quantum circuit for performing polynomial evaluation corresponding to each of the first subset of the set of sub-intervals. The method may further include, in the second mode, performing the evaluation of the polynomial using Newton's method for a second subset of the set of sub-intervals corresponding to the input, wherein the second subset is different from the first subset.

This Summary is provided to introduce a selection of concepts in a simplified form that are further described below in the Detailed Description. This Summary is not intended to identify key features or essential features of the claimed subject matter, nor is it intended to be used to limit the scope of the claimed subject matter. Furthermore, the claimed subject matter is not limited to implementations that solve any or all disadvantages noted in any part of this disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is illustrated by way of example and is not limited by the accompanying figures, in which like references indicate similar elements. Elements in the figures are illustrated for simplicity and clarity and have not necessarily been drawn to scale.

FIG. 1 shows a parallel polynomial evaluation circuit in accordance with one example;

FIG. 2 shows a circuit for computing the label register l in accordance with one example;

FIG. 3 shows a flow chart of a method for performing an evaluation of a polynomial corresponding to an input in accordance with one example;

FIG. 4 shows a flow chart of another method for performing an evaluation of a polynomial corresponding to an input in accordance with one example;

FIG. 5 shows a flow chart of another method for performing an evaluation of a polynomial corresponding to an input in accordance with one example;

FIG. 6 shows a flow chart of another method for performing an evaluation of a polynomial corresponding to an input in accordance with one example;

FIG. 7 shows an example system environment for implementing the methods described in the present disclosure;

FIG. 8 shows absolute errors of an inverse square root function before tuning the constants in accordance with one example;

FIG. 9 shows absolute errors of an inverse square root function after tuning the constants in accordance with one example;

FIG. 10 shows absolute errors of a reversible implementation of the arcsine function in accordance with one example;

FIG. 11 shows a circuit for the n-th Newton iteration of computing of the inverse square root of a variable in accordance with one example; and

FIG. 12 shows absolute errors of a reversible implementation of the square root function for a certain number of Newton iterations in accordance with one example.

DETAILED DESCRIPTION

Examples described in this disclosure relate to evaluating quantum computing circuits in view of the resource costs of a quantum algorithm. Certain examples relate to the problem of implementing numerical operations on a quantum computer. Such operations arise, for instance, in the implementation of classical functions on a superposition of inputs and are omnipresent in the quantum algorithm literature. One example case is the computation of amplitudes from a given subroutine specification. In this case, a need for computing arcsine functions or arcsine functions composed with a square root function arises. In order to facilitate implementation of numerical operations, their testing, and estimation of resource requirements, a quantum software module is implemented. This module allows one to automatically generate quantum circuits for evaluating piecewise smooth functions of inputs that are given in the computational basis. The present disclosure provides details of this module and provides resource estimates for some often-encountered classical functions such as Gaussians, tanh (hyperbolic tangent), and sine/cosine. In addition, the present disclosure describes implementations of inverse square root and arcsine, which may also be used in quantum chemistry and machine learning applications. The provided cost estimates might help to identify and address bottlenecks in future quantum applications. One technique relates to triaging inputs into subintervals, which are evaluated using a circuit for parallel polynomial evaluation.

Certain examples address the implementation of numerical functions such as 1/x, 1/sqrt(x), arc sin(x), sin(x), tanh(x), and in general the implementation of smooth functions via a piece-wise polynomial approximation, followed by functions that are used in specific applications. Some examples combine the evaluation scheme for smooth functions with the inverse square root, which is for instance, applicable to the case of implementation of arc sin(x). Having reversible implementations of these functions available enables more detailed cost analyses of various quantum algorithms such as the HHL algorithm for matrix inversion. Other use cases arise in Quantum Metropolis sampling, Gibbs state preparation and in the widely applicable framework of Quantum Rejection Sampling to transform one or more samples of a given quantum state into a quantum state with potentially different amplitudes, while maintaining relative phases.

Further applications of numerical functions can be anticipated in the quantum machine learning, e.g., threshold functions that need to be evaluated on a superposition of values, e.g., tanh(x) or sigmoid functions. In addition, 1/√{square root over (x)} can be used, e.g., for (re-)normalization of intermediate results. In quantum algorithms for chemistry further examples for numerical functions arise, e.g., for on-the-fly computation of the one- and two-body integrals. There, 1/√{square root over (x)} as well as the evaluation of smooth functions such as Gaussians is needed, and on-the-fly computation of finite element matrix elements often involves the evaluation of smooth functions such as sin(x) and cos(x).

A basic scheme to evaluate a single polynomial on a quantum computer in the computational basis is the classical Horner scheme, which evaluates:

$P\left(x\right)=\sum _{i=0}^d{a}_i{x}^i$

by iteratively performing a multiplication by x, followed by an addition of a_i for i ∈ {d, d−1, . . . ,0}. This may include performing the following operations:

a _d x+a _d−1 a _d x ² +a _d−1 x+a _d−2 . . . a _d x ^d + . . . a _o.

A reversible implementation of this scheme may simply store all intermediate results. At iteration i, the last iterate y_i−1 is multiplied by x into a new register y_i, followed by an addition by the constant a_i. Due to the linear dependence of successive iterates, a pebbling strategy can be employed in order to optimize the space/time trade-offs according to some chosen metric.

Oftentimes, the degree d of the minimax approximation over a domain Ω must be chosen to be very high in order to achieve a certain L_∞(Ω)-error. In such cases, it makes sense to partition Ω, i.e., find Ω_i such that

$\begin{array}{cc}\Omega =\underset{i=0}{\bigcup ^M}{\Omega }_i,& {\Omega }_i\bigcap {\Omega }_j=\varnothing \forall i\ne j,\end{array}$

and to then perform a case distinction for each input, evaluating a different polynomial for x ∈ Ω_i than for y ∈ Ω_j if i≠j. A straight-forward generalization of this approach to the realm of quantum computing would loop over all subdomains Ω_i and, conditioned on a case-distinction or label register l, evaluate the corresponding polynomial. Thus, the cost of this inefficient approach grows linearly with the number of subdomains.

In order to improve upon this approach, one can parallelize the polynomial evaluation if the degree d is constant over the entire domain Ω. FIG. 1 shows a parallel polynomial evaluation circuit 100 in accordance with one example. NEXT_a 112 changes the register to hold the next set of coefficients (in superposition) Σ_l|l|a_{l, i−1}Σ_l|l|a_l,i. MUL 110 and ADD 114 perform a multiplication and an addition, respectively. The small triangle indicates the output of the ADD 114 and MUL 110 gates. In this example, merely adding the label register l mentioned above and performing:

|_l,i−₁x|0|l|y_l,i−₁x|a_l,i|l  (1)

|y_l,i−₁x+a_l,i|a_l,i|l  (2)

|y_l,i|0|l  (3)

enables the evaluation of multiple polynomials in parallel. The impact on the circuit size is minor. This is because, in this example, the depth of the circuit remains unaltered, since the initialization step (1) can be performed while multiplying the previous iterate y_i−1 by x; see FIG. 1. An illustration of the circuit for computing the label register l can be found in FIG. 2. The LABEL gate 210 (additional details in 220) initializes the label register |l which consists of ┌log₂(M)┐ qubits, to indicate the subdomain Ω_l to which x belongs. P_i (e.g., P₁ 230, P₂ 232, and P_M 234) computes the predicate indicating whether x ∈ Ω_i into the ancilla qubit. Conditioned on this result, the label is then initialized to the value chosen to represent the i-th interval. The purpose of the SET blocks 222, 224, 226 is to create in the label register a particular bit-pattern that corresponds to the specific label that identifies the subdomain Ω_l to which x belongs. To this end, the SET blocks 222, 224, 226 consist of bit-flip operations that flip the bits of the bit-pattern corresponding to the subdomain Ω_l when starting from the all-zero-bit pattern. A slight drawback of this parallel evaluation is that it requires one extra ancilla register for the last iteration, since the in-place addition circuit can no longer be used. The small overhead of using many intervals results in good approximations already for low-degree polynomials (and thus using few qubit registers).

Using reversible pebble games, it is possible to trade the number of registers needed to store the iterates with the depth of the resulting circuit. The parameters are: the number n of bits per register, the total number m of these n-qubit registers, the number r of Horner iterations, and the depth d of the resulting circuit. The trade-space we consider involves m, r, and d. In particular, the optimal circuit depth for a fixed number m of registers and a fixed number r of iterations is determined. In one example, dynamic programming is used to construct the optimal strategies as the dependency graph is just a line which is due to the sequential nature of Horner's method. The optimal number of pebbling steps as a function of m and r can be found in Table 1 below.

TABLE 1
m\r	1	2	3	4	5	6	7	8	16	32	64
1	1	∞	∞	∞	∞	∞	∞	∞	∞	∞	∞
2	1	3	∞	∞	∞	∞	∞	∞	∞	∞	∞
3	1	3	5	9	∞	∞	∞	∞	∞	∞	∞
4	1	3	5	7	11	15	19	25	∞	∞	∞
5	1	3	5	7	9	13	17	21	71	∞	∞
6	1	3	5	7	9	11	15	19	51	193	∞
7	1	3	5	7	9	11	13	17	49	145	531
8	1	3	5	7	9	11	13	15	47	117	369

In order to enable automatic compilation of a classical oracle used in a quantum algorithm, certain examples use the Remez algorithm as a preferred embodiment to obtain a piece-wise polynomial approximation, which can be implemented using a quantum circuit for parallel polynomial evaluation. This allows one to triage inputs into different sub-intervals and to apply the parallel polynomial evaluation for different polynomials, depending on the particular subinterval the input x falls into. This method works for superposition of input x also.

FIG. 3 shows a flow chart for implementing an evaluation of a polynomial on a quantum computer with a given set of intervals and a given degree. Step 310 may include receiving an input, including a function f, an interval subdivision of input, and a maximal degree d. A system (e.g., quantum computing system 710 of FIG. 7) may include access to at least one classical algorithm for polynomial interpolation and quantum circuit libraries for integer arithmetic.

Step 320 may include determining the best polynomial of degree at most d on each interval using an algorithm for polynomial evaluation. In one example, the Remez algorithm may be used to determine a piecewise polynomial approximation. In order to enable automatic compilation of a classical oracle used in a quantum algorithm, the Remez algorithm can be used in a subroutine to determine a piece-wise polynomial approximation, which can then be implemented using the circuit described later.

In particular, the oracle may be implemented with a given precision, accuracy, and number of available quantum registers (or, equivalently, the polynomial degree d) over a user-specified interval Ω=[a, a+L). The example algorithm proceeds as follows: in a first step, run the Remez algorithm which, given a function ƒ(x) over a domain Ω ⊂ and a polynomial degree d, finds the polynomial P(x) which approximates ƒ(x) with minimal L_∞(Ω)-error, and check whether the achieved error is low enough. If it is too large, reduce the size of the domain Ω₁:=[a, a+L/2) and check again.

Repeating this procedure and carrying out binary search on the right interval border will eventually lead to the first subdomain Ω₁=[a, b₁) with the corresponding degree d polynomial achieving the desired accuracy on Ω₁. Next, the next subdomain Ω₂=[b₁, b₂) is determined using the same procedure. This is iterated until b_i≥b, meaning that all required subdomains and their corresponding polynomials have been determined and ƒ(x) can be implemented using a parallel polynomial evaluation circuit. This algorithm was implemented and then run for various functions, target accuracies, and polynomial degrees in order to determine approximate resource estimates for these parameters.

In one example, the Remez algorithm may be implemented using the circuit in FIG. 1. As explained earlier, the circuit in FIG. 1 may perform parallel evaluation of a polynomial. The MUL box 110 and the ADD box 114 may perform a multiplication and an addition, respectively. The small black triangles in FIG. 1 indicate the output of the ADD and MUL gates, respectively. Referring still to FIG. 1, the NEXT_a block 112 changes the register to hold the next set of coefficients (in superposition). FIG. 2 may be used to label the coefficients to indicate the subdomains in which an input x belongs. In this example, the subdomain may be a superset of the intervals.

Table 2 below provides example code corresponding to a function to implement the Remez's algorithm to find a polynomial for a given interval. The code also includes a function to subdivide an interval into subintervals.

TABLE 2
import math
import numpy as np
from matplotlib import pyplot as plt
# returns n chebyshev nodes on the interval (a,b)
def getChebyshevNodes(n, a, b):
nodes = [.5*(a+b) + .5*(b−a)*math.cos((2*k+1)/(2.*n)*math.pi) for k in
range (n)]
return nodes
# returns the error on given nodes of a polynomial with coefficients
polycoeff
# approximating the function with function values exactvals (on these
nodes).
def getErrors(exactvals, polycoeff, nodes):
ys = np.polyval(polycoeff, nodes)
for i in range(len(ys)):
ys[i] = abs(ys[i]−exactvals[i])
return ys
# returns the coefficients of a polynomial of degree d approximating the
function fun over the interval (a,b)
# actually, cn2 should be chosen adaptively but for most functions this
works just fine.
def RunRemez(fun, a, b, d = 5, odd_or even=False, even = False, tol =
1.e−13):
finished = False
cn = getChebyshevNodes(d+2, a, b)
cn2 = getChebyshevNodes(50*d, a, b)
it = 0
while not finished and len(cn)==d+2 and it < 50:
it += 1
b = np.array([fun(c) for c in cn])
A = np.matrix(np.zeros([d+2,d+2]))
for i in range(d+2):
x = 1.
if odd_or_even and not even:
x *= cn[i]
for j in range(d+2):
A[i, j] = x
x *= cn[i]
if odd_or_even:
x *= cn[i]
A[i, −1] = (−1)**(i+1)
res = np.linalg.solve(A, b)
revlist = list(reversed(res[0:−1]))
sccoeff = [ ]
for c in revlist:
sccoeff.append(c)
if odd_or_even:
sccoeff.append(0)
if even:
sccoeff = sccoeff[0:−1]
errs = getErrors([fun(c) for c in cn2], sccoeff, cn2)
maximum_indices = [ ]
if errs[0]>errs[1]:
maximum_indices.append(0)
for i in range(1,len(errs)−1):
if errs[i]>errs[i−1] and errs[i]>errs[i+1]:
maximum_indices.append(i)
if errs[−1]>errs[−2]:
maximum_indices.append(−1)
finished = True
for idx in maximum_indices[1:]:
if abs(errs[idx]−errs[maximum_indices[0]]) > tol:
finished = False
cn = [cn2[i] for i in maximum_indices]
return (max(abs(errs)),sccoeff)
# Partitions the interval [a,b] and approximates the function f
# using piecewise minimax polynomial approximation of degree d and
max. error eps
# The function roundIntervalBorder can modify the right border of the
interval
# e.g., round it to a specific #bits
def RunPartitioning(f, a, b, d, eps, roundIntervalBorder, Odd, Even):
A = [roundIntervalBorder(a)]
coeffs = [ ]
numIntervals = 0
while A[−1] < roundIntervalBorder(b)−1.e−14:
a = A[numIntervals]
deltaL = (roundIntervalBorder(b) − a)*.5
L = roundIntervalBorder(b) − a
while deltaL > 1.e−12:
if roundIntervalBorder(a + L) − a < 1.e−14:
raise RuntimeError(“roundIntervalBorder rounded interval
length to zero! Please decrease epsilon or increase the resolution of the
interval borders.”)
err, coeff = RunRemez(f, a, roundIntervalBorder(a + L), d, Odd
or Even, Even, 1.e−14)
if err > eps:
L −= deltaL
elif a+L < roundIntervalBorder(b):
L += deltaL
else:
L = roundIntervalBorder(b)−a
deltaL *= 0.5
coeffs.append(coeff)
A.append(roundIntervalBorder(a + L))
numIntervals += 1
return A[0:numIntervals+1], coeffs
# Function to plot the resulting partitioning
def PlotApproximation(f, A, coeffs):
numSteps = 100
nodes = [ ]
approx_y = [ ]
for i in range(len(A)−1):
curvalues = [A[i](A[i+1]−A[i])*k/numSteps for k in
range(numSteps)]
approx_y.extend(np.polyval(coeffs[i], curvalues))
nodes.extend(curvalues)
ys = [f(x) for x in nodes] # exact values
yerrs = [abs(ys[i]−approx_y[i]) for i in range(len(ys))]
plt.plot(nodes, yerrs)
plt.yscale(‘log’)
plt.show( )
if_name_ == “_main_”:
#################################
# Let's test the implementation:#
#################################
def f(x):
return math.exp(−x)
err = 1.e−7
Odd = False
Even = False
a = 0.0
b = −math.log(err) #math.pi/2
A, coeffs = RunPartitioning(f, a, b, 5, err, lambda x: x, Odd, Even)
print(“Segments: { }”.format(len(A)))
for i in range(len(A)−1):
print(“{ },{ }):\t{ }”.format(A[i], A[i+1], coeffs[i]))
PlotApproximation(f, A, coeffs)

Step 330 may include determining the overall precision that is achievable for any input that falls into any interval.

Step 340 may include constructing a quantum circuit for parallel polynomial evaluation of polynomials on each subinterval. After all of the subdomains and the corresponding polynomials have been determined, ƒ(x) may be implemented using the parallel polynomial evaluation circuit shown in FIG. 2.

Step 350 may include generating an output, including a quantum circuit for function evaluation of f.

FIG. 4 shows a flow chart for implementing an evaluation of a polynomial on a quantum computer with a target precision epsilon and an interval subdivision of the input. Thus, unlike the method described with respect to FIG. 3, in this method the maximal degree of precision is not specified a priori. Thus, as part of this example method, the degree of the polynomial is incremented by 1 each time until the target precision is achieved. Step 410 may include receiving an input, including a function f, a target precision epsilon, and an interval subdivision of input. A system (e.g., quantum computing system 710 of FIG. 7) may include access to at least one classical algorithm for polynomial interpolation and quantum circuit libraries for integer arithmetic. Step 420 may include starting with a degree d=0 of the polynomial.

Next, step 430 may include determining the best polynomial interpolation of degree at most d on each interval using an algorithm for polynomial evaluation. In one example, the Remez algorithm may be used to determine a piecewise polynomial. In one example, the Remez algorithm may be implemented using the circuit in FIG. 1.

Next, step 440 may include determining whether a target precision is achievable. If not, then the degree of the polynomial may be incremented by 1, as shown in step 450. If yes, then step 460 may include constructing a quantum circuit for polynomial evaluation of each polynomial on each subinterval. After all of the subdomains and the corresponding polynomials have been determined, ƒ(x) may be implemented using the parallel polynomial evaluation circuit shown in FIG. 1.

Next, step 470 may include generating an output, including a quantum circuit for evaluation of the function f.

FIG. 5 shows a flow chart for implementing an evaluation of a polynomial on a quantum computer with a target precision epsilon and a maximal degree d. Thus, unlike the method described with respect to FIG. 1, in this method the sub intervals are not known. Thus, as part of this example method, the intervals are subdivided until the target precision is achieved. Step 510 may include receiving an input, including a function f, a target precision epsilon, and an interval subdivision of input. A system (e.g., e.g., quantum computing system 710 of FIG. 7) may include access to at least one classical algorithm for polynomial interpolation and quantum circuit libraries for integer arithmetic. Step 520 may include starting with a degree d=0 of the polynomial.

Next, step 530 may include determining the best polynomial interpolation of degree at most d on each interval using an algorithm for polynomial evaluation. In one example, the Remez algorithm may be used to determine a piecewise polynomial approximation. In one example, the Remez algorithm may be implemented using the circuit in FIG. 1. The circuit in FIG. 1 may perform parallel evaluation of a polynomial. As described earlier, the MUL box 110 and the ADD box 114 may perform a multiplication and an addition, respectively. The small black triangles in FIG. 1 indicate the output of the ADD and MUL gates, respectively. Referring still to FIG.1, the NEXT_a block 112 changes the register to hold the next set of coefficients (in superposition). As explained earlier, the circuit for labeling shown in FIG. 2 may be used to label the coefficients to indicate the subdomains in which an input x belongs. In this example, the subdomain may be a superset of the intervals.

Next, step 540 may include determining whether a target precision is achievable. If not, then the intervals may be further subdivided, as shown in step 550. Table 3 below provides an example of code for interval partitioning functions, such as sin(x), arc sin(x), tanh(x), exp(x) etc.

TABLE 3
import math
import numpy as np
from matplotlib import pyplot as plt
from IntervalPartitioner import *
def run_gaussian(eps, deg):
def f(x):
return math.exp(−x*x)
Odd = False
Even = True
a = 0
b = 5
A, coeffs = run_fun(f, a, b, deg, eps, Odd, Even)
return len(A)
def run_tanh(eps, deg):
def f(x):
return math.tanh(x)
Odd = True
Even = False
a = 0
if eps > 3.e−6:
b = 7
elif eps > 3.e−8:
b = 9
else:
b = 10
A, coeffs = run_fun(f, a, b, deg, eps, Odd, Even)
return len(A)
def run_sin(eps, deg):
def f(x):
return math.sin(x)
Odd = True
Even = False
a = 0
b = math.pi/2.
A, coeffs = run_fun(f, a, b, deg, eps, Odd, Even)
return len(A)
def run_fun(f, a, b, deg, eps, Odd, Even):
A, coeffs = RunPartitioning(f, a, b, deg, eps, lambda x: x, Odd, Even)
#PlotApproximation(f, A, coeffs)
return A, coeffs
def run_arcsin(eps, deg):
def f(x):
return math.asin(x)
#eps = 1.e−7
a = 0
b = 0.5
#deg = 4 # polynomial degree to use
p = math.log(b,2) # point position
A, coeffs = run_fun(f, a, b, deg, eps, True, False)
return len(A)
def run_exp(eps, deg):
def f(x):
return math.exp(−x)
#eps = 1.e−7
a = 0
b = −math.log(eps)
#deg = 4 # polynomial degree to use
p = math.log(b,2) # point position
A, coeffs = run_fun(f, a, b, deg, eps, False, False)
return len(A)
PlotApproximation(f, A, coeffs)
print(“Can do it using { }intervals:”.format(len(A)))
print(A)
print(“\nCoefficients:”)
print(coeffs)
num = len(A) # number of intervals
n = int(math.ceil(math.log(1./eps, 2.)) + p + 8)
qubits = int((deg+1)*n + math.ceil(math.log(num,2.)))
d = deg+1
toffolis = int(math.ceil(1.5*n**2*d + 3*n*p*d + 3.5*n*d − 3*p**2*d
+3*p*d − d + 2*num*d*(4*math.ceil(math.log(num,2))−8) +4*num*n))
print(“\n\n***Estimates:***\nToffolis: { }\nQubits:
{ }”.format(toffolis,qubits))
if_name_ == “_main_”:
#run_exp( )
#exit( )
def cols(num_at):
return “ &”*(6-num_at) + “\\\\”
epsilons = [10.**(−5), 10.**(−7), 10.**(−9)]
tex_epsilons = [‘$10{circumflex over ( )}{−5}$’, ‘$10{circumflex over ( )}{−7}$’, ‘$10{circumflex over ( )}{−9}$’]
degrees = [3, 4, 5, 6]
print(“Func\tEpsilon\tdeg\tnum_intervals”)
functions = [run_tanh, run_gaussian, run_sin, run_exp, run_arcsin]
functionname = [‘$\\tanh(x)$’, ‘$\\exp(−x{circumflex over ( )}2)$’, ‘$\\sin(x)$’, ‘$\\exp(−
x)$’,‘$\\arcsin(x)$’]
pointpositions = [8, 7, 3, 0, 0]
for i in range(len(functions)):
print(“{ }{ }”.format(functionname[i], cols(1)))
fun = functions[i]
for e in range(len(epsilons)):
eps = epsilons[e]
if i <= 2:
p = pointpositions[i]
elif i == 3:
p = math.log(−math.log(eps), 2)
elif i == 4:
assert functionname[i] == ‘$\\arcsin(x)$’
p = 1
print(“& {}{}”.format(tex_epsilons[e], cols(2)))
for deg in degrees:
num = fun(eps, deg)
n = int(math.ceil(math.log(1./eps, 2.)) + p + 8)
qubits = int((deg+1)*n + math.ceil(math.log(num,2.)))
d = deg+1
toffolis = int(math.ceil(1.5*n**2*d + 3*n*p*d + 3.5*n*d −
3*p**2*d + 3*p*d − d + 2*num*d*(4*math.ceil(math.log(num,2))−8)
+4*num*n))
print(“& & { } & { } & { } & { }\\\\”.format(deg, num, qubits, toffolis))

If yes, then step 560 may include constructing a quantum circuit for polynomial evaluation of each polynomial on each subinterval. After all of the subdomains and the corresponding polynomials have been determined, ƒ(x) may be implemented using the parallel polynomial evaluation circuit shown in FIG. 2.

Next, step 570 may include generating an output, including a quantum circuit for evaluation of the function f.

While FIGS. 3-5 describe methods for polynomial interpolation of certain quantum algorithms or functions, there are other functions that may require implementation using a Newton iteration. As an example, functions such as inverse square root may be better evaluated using a method like the Newton-Raphson method. In addition, for certain functions the evaluation may be better performed using the polynomial interpolation described with respect to FIG. 3 for certain intervals; but for other intervals the evaluation may be better performed using the Newton-Raphson method.

FIG. 6 shows a flow chart or implementing evaluation of a polynomial on a quantum computer based on given intervals, using polynomial interpolation or Newton's method depending upon the interval. In step 610, an input, including a function f, target precision epsilon, and intervals subdividing the input may be received. A classical algorithm for polynomial interpolation as well as a quantum circuit library to implement basic arithmetic operations such as ADD and MUL can be provided in step 615.

In step 620, for each interval, a system (e.g., quantum computing system 710 of FIG. 7) may determine whether to interpolate or apply Newton's method.

In step 630, the system may evaluate whether with respect to the next interval, interpolation or Newton's method should be used. In the case of interpolation, as part of step 640, a quantum circuit for the polynomial evaluation may be determined. Any of the techniques described earlier may be used as part of this step.

In the case of application of the Newton-Raphson method, in step 650, the number of rounds of Newton's method may be determined. Next, in step 660, an initial constant for the Newton's method may be determined. Next, in step 670, a space-optimized circuit for the Newton's method may be constructed. This step may include access to an application of a pebbling strategy to reversibly compute programs with linear dependency graphs 665. In one example, as part of the Newton-Raphson method, the pebbling technique may be used to conserve space while using the parallel polynomial circuit shown in FIG. 1. The results of the multiply and add operations continue to accumulate over time and can use up a large amount of space. The pebbling technique may be implemented by selecting a suitable number of pebbles. As an example, in an instance where there are eight pebbles, 64 rounds of the Newton-Raphson method may be implemented and there will be an overhead of 369. Additional details for implementing the pebbling strategy as part of the parallel polynomial circuit use are described previously with respect to Table 1.

Once all of the intervals are completed, as determined in step 680, the output, including the circuit for approximation of function f, may be provided.

FIG. 7 shows an example system environment for implementing aspects of the technology and the methods described in the present disclosure. System environment includes a quantum computing system 710 and a classical computing system 730, which is coupled to a remote computer 750. Quantum computing system 710 consumes the quantum circuits generated using the classical computing system 730. Quantum computing system 710 may include a quantum processor 712 and measuring/monitoring devices 714. In one example, quantum processor 712 and measuring/monitoring devices 714 may be configured to operate in a cryogenic environment (e.g., between 4 Kelvin and 77 Kelvin) such that quantum processor 712 may perform superconducting operations. Quantum processor 712 may execute quantum circuits that are compiled by classical computing system 730. The compiled quantum circuits may be communicated to quantum processor 712 via bus 706.

With continued reference to FIG. 7, classical computing system 730 may include communication interface(s) 732, classical processor 734, and memory 736. Memory 736 may include a compiler unit 738 including libraries and other programs or code to compile a high-level description of a quantum algorithm into quantum circuits. Compiler unit 738 may further include programs or code to execute the steps described with respect to FIGS. 3-5. Thus, compiler unit 738 may include programs or code that when executed by classical processor 734 may perform the various methods described with respect to the present disclosure. In one example, the code shown in Tables 1, 2, and 3 may be stored in memory 736 either as part of compiler unit 738 or separately. The high-level description of the quantum algorithm may be stored in memory 736 or in memory 752, which is associated with remote computer 750.

For quantum chemistry or machine learning applications, non-smooth functions are also required. As an example, the inverse square root can be used in both examples, namely for the calculation of the Coulomb potential and to determine the reciprocal when employing HHL for quantum machine learning.

In classical computing, inverse square roots appear in computer graphics and the term “fast inverse square root” is often used: it labels the procedure to approximate the inverse square root using bit-operations on the floating-point representation of the input. The code ultimately performs a Newton-Raphson iteration in order to improve upon a pretty accurate initial guess, which it finds using afore-mentioned bit-operations. Loosely speaking, the bit-operations consist of a bit-shift to divide the exponent by two in order to approximate the square root, followed by a subtraction of this result from a magic number, effectively negating the exponent and correcting the mantissa, which was also shifted together with the exponent. The magic number can be chosen using an auto-tuning procedure and varies depending on the objective function being used. This provides an extremely good initial guess for the Newton iteration at very low cost.

In an example implementation, a similar procedure is used to compute the inverse square root using fixed-point arithmetic. While the floating-point representation is not used, a low-cost initial guess is found, which allows for a small number of Newton iterations to be sufficient (i.e., 2-4 iterations). This includes determining the position of the first one in the bit-representation of the input, followed by an initialization which involves a case distinction on the magic number to use. The three constants were tuned such that the error peaks near powers of two in FIG. 8 vanish. The peaks appear due to the fact that the initial guess takes into account the location of the first one but completely ignores the actual magnitude of the input. For example, all inputs in [1,2) yield the same initial guess. The error plot with tuned constants is depicted in FIG. 9. The errors were evaluated for N=2000 (equidistant) points in the interval [1/N, 5] using m ∈ {2,3,4} Newton iterations and corresponding bit sizes n ∈ {25,35,55}. The fixed-point position is p=12, in order to ensure that no overflow occurs for small inputs. One can clearly observe that an entire Newton iteration can be saved when aiming for a given L_∞-error.

Arcsine may be implemented as a combination of polynomial evaluation and (inverse) square root. Approximating the arcsine using only a polynomial allows for a good approximation in [0,0.5], but not near 1 (where it diverges). The Cephes math library may remedy this problem by adding a case distinction, employing a “double-angle identity” for x≥0.5. This may require computing the square root, which can be achieved by first calculating the inverse square root, followed by x·1/√{square root over (x)}. A resulting error plot from the implementation of the arcsine function is shown in FIG. 10, which is on [0,1] for N=2000 points of the example reversible implementation of the arcsine using m ∈ {3,4,5} Newton iterations for calculating the inverse square root. The fixed-point position is chosen to be p=2 and total bit size n was chosen to be in {35,50,55}. The oscillations stem from the minimax polynomial which is used to approximate the arcsine on [0,0.5].

Note that certain applications may allow to a trade-off between the error in the arcsine with, e.g., probability of success by rescaling the input such that the arcsine needs to be computed only for values in [0,0.5]. This would allow one to remove the case-distinction and the subsequent calculation of the square root. As a result, one could evaluate the arcsine at a cost that is similar to the implementation costs of sin/cos. Estimates for the Toffoli and qubit counts for this case are shown in Table 4 below.

TABLE 4
					Number of
		Polynomial	Number of	Number of	Toffoli
Function	L∞ error	degree	subintervals	qubits	gates
arcsin(x)	10−5	3	2	105	4872
		4	2	131	6038
		5	2	157	7204
		6	2	183	8370
	10−7	3	3	134	7784
		4	2	166	9419
		5	2	199	11250
		6	2	232	13081
	10−9	3	6	159	11264
		4	3	197	13138
		5	3	236	15672
		6	2	274	17938

In fixed-point arithmetic, one represents numbers x using n bits as

$x=\underset{p}{\underbrace{x_{n-1}\dots x_{n-p}}}·\underset{n-p}{\underbrace{x_{n-p-1}\dots x_0}},$

where x_i ∈ {0,1} is the i-th bit of the binary representation of x, and the point position p denotes the number of binary digits to the left of the binary point. We choose both the total number of bits n and the point position p to be constant over the course of a computation. As a consequence, over- and underflow errors are introduced, while keeping the required bit-size from growing with each operation.

Unless indicated otherwise, a fixed-point addition implementation is used, which keeps the bit-size constant. This amounts to allowing over- and underflow, while keeping the registers from growing with each operation. Multiplication can be performed by repeated-addition-and-shift, which can be seen from

x·y=x _n−12ⁿ⁻¹y+ . . . +x₀2⁰y,

where x=Σ_ix_i2ⁱ with x_i ∈ {0,1} denotes the binary expansion of the n-bit number x. Thus, for i ∈ {0, . . . , n−1}, 2^i−(n−p)y is added to the result register (which is initially zero) if x_i=1. This can be implemented using n controlled additions on 1, 2, . . . , n bits if one allows for pre-truncation. Instead of computing the 2n-bit result and copying out the first n bits before uncomputing the multiplication again, the additions can be executed on a subset of the qubits, ignoring all bits beyond the scope of the n-bit result. Thus, each addition introduces an error of at most

${\epsilon }_A=\frac{1}{2^{n-p}}.$

Since there are (at most) n such additions, the total error is

$\epsilon =\frac{n}{2^{n-p}},$

a factor n larger than using the costly approach mentioned above.

Negative multipliers are dealt with by substituting the controlled addition by a controlled subtraction when conditioning on the most significant bit because it has negative weight x_MSB=−2ⁿ⁻¹ in two's-complement notation. The multiplicand is assumed to be positive throughout, which removes the need for conditional inversions of input and output (for every multiplication), thus tremendously reducing the size of circuits that require many multiplications such as, e.g., polynomial evaluation.

The square of a number can be calculated using the same approach as for multiplication. Yet, one can save (almost) an entire register by only copying out the bit being conditioned on prior to performing the controlled addition. Then the bit can be reset using another CNOT gate, followed by copying out the next bit and performing the next controlled addition. The gate counts are identical to performing

|x|0|0|x|x|0|x|x|x²|x|x²|0,

while allowing to save n−1 qubits.

Next, resource estimates for polynomial evaluation are described. The evaluation of a degree d polynomial requires an initial multiplication a_d·x, an addition of a_d−1, followed by d−1 multiply-accumulate instructions. The total number of Toffoli gates is thus equal to the cost of d multiply-accumulate instructions. Furthermore, d+1 registers are required for holding intermediate and final result(s) if no in-place adder is used for the last iteration (and no non-trivial pebbling strategy is applied). Other strategies may be employed in order to reduce the number of ancilla registers, at the cost of a larger gate count, see Tables 5-8 for examples of resource estimates for a number of functions.

TABLE 5
					Number of
		Polynomial	Number of	Number of	Toffoli
Function	L∞	degree	subintervals	qubits	gates
tanh(x)	10−5	3	15	136	12428
		4	9	169	13768
		5	7	201	15492
		6	5	234	17544
	10−7	3	50	166	27724
		4	23	205	23095
		5	14	244	23570
		6	10	284	26037
	10−9	3	162	192	77992
		4	59	236	41646
		5	30	281	35460
		6	19	327	36578

TABLE 6
					Number of
		Polynomial	Number of	Number of	Toffoli
Function	L∞ error	degree	subintervals	qubits	gates
esp(−x2)	10−5	3	11	132	10884
		4	7	163	12141
		5	5	195	14038
		6	4	226	15863
	10−7	3	32	161	20504
		4	15	199	19090
		5	10	238	21180
		6	7	187	23254
	10−9	3	97	187	49032
		4	36	231	32305
		5	19	275	30234
		6	12	319	31595

TABLE 7
					Number of
		Polynomial	Number of	Number of	Toffoli
Function	L∞ error	degree	subintervals	qubits	gates
sin(x)	10−5	3	2	132	6188
		4	2	141	7679
		5	2	169	9170
		6	2	197	10661
	10−7	3	3	142	9444
		4	2	176	11480
		5	2	211	13720
		6	2	246	15960
	10−9	3	7	167	13432
		4	3	207	15567
		5	2	247	18322
		6	2	288	21321

TABLE 8
					Number of
		Polynomial	Number of	Number of	Toffoli
Function	L∞ error	degree	subintervals	qubits	gates
exp(−x)	10−5	3	11	116	8106
		4	6	143	8625
		5	5	171	10055
		6	4	198	11245
	10−7	3	31	149	17304
		4	15	184	15690
		5	9	220	16956
		6	7	255	18662
	10−9	3	97	175	45012
		4	36	216	28302
		5	19	257	25721
		6	12	298	26452

Note that all multiplications can be carried out assuming x>0, i.e. x can be conditionally inverted prior to the polynomial evaluation (and the pseudo-sign bit is copied out). The sign is then absorbed into the coefficients. Before adding a_i into the y_i−1x-register, it is inverted conditioned on the sign-bit of x being set if the coefficient corresponds to an odd power. This is done because it is cheaper to implement a fixed-point multiplier, which can only deal with y_i−1 being negative.

The Toffoli gate count of multiplying two n-bit numbers is

$\begin{array}{c}{T}_{\mathrm{mul}}\left(n,p\right)={\sum }_{i=0}^{p-1}T\mathrm{cadd}\left(n-i\right)+{\sum }_{i=1}^{n-p}T\mathrm{cadd}\left(n-i\right)\\ ={\sum }_{i=0}^{p-1}3\left(n-i\right)+{\sum }_{i=1}^{n-p}3\left(n-i\right)+3n\\ =\frac{3}{2}{n}^2+3\mathrm{np}+\frac{3}{2}n-3p^2+3p\end{array}$

if one uses the controlled addition circuit by Takahashi et al., which requires 3n+3 Toffoli gates to (conditionally) add two n-bit numbers. The subsequent addition can be implemented using the addition circuit by Takahashi et al. , featuring 2n−1 Toffoli gates. Thus, the total cost of a fused multiply-accumulate instruction is

T _{fma(n,p)=3/2n +3np +7/2n−3p 2 +3p−1.}

Therefore, the total Toffoli count for evaluating a degree d polynomial is

T _{poly(n,d,p)=3/2n 2 d+3npd+7/2nd−3p 2 d+3pd−d.}

Evaluating M polynomials in parallel for piecewise polynomial approximation requires only

$n+\underset{2}{\lceil \log }M\rceil$

additional qubits (since one n-qubit register is required to perform the addition in the last iteration, which is no longer just a constant)

$\mathrm{and}2M\underset{2}{\lceil \log }M\rceil -\mathrm{controlled}$

NOT gates, which can be performed in parallel with the multiplication. This increases the circuit size by

T _extra(M)=2M(4┌log₂ M┐−8)

Toffoli gates per multiply-accumulate instruction, since a k-controlled NOT can be achieved using 4(k−2) Toffoli gates and k−2 dirty ancilla qubits, which are readily available in this construction.

The label register l can be computed using 1 comparator per subinterval by:

I _i =[a _i , a _{i i+1} ), a₀<a₁< . . . <a_M−1.

The comparator stores its output into one extra qubit, flipping it to 1 if x≤a_i+1. The label register is then incremented from i−1 to i, conditioned on this output qubit still being 0 (indicating that x>a_i). Incrementing l can be achieved using CNOT gates applied to the qubits that correspond to ones in the bit-representation of (i−1) ⊕ i. Finally, the comparator output qubit is uncomputed again. This procedure is carried out M times for i=0, . . . , M−1 and requires 1 additional qubit. The number of extra Toffoli gates for this label initialization is

T _label(M,n)=M·2T_cmp(n)=4Mn,

where, as a comparator, we use the CARRY-circuit, which needs 2n Toffoli gates to compare a classical value to a quantum register, and another 2n to uncompute the output and intermediate changes to the n required dirty ancilla qubits.

In total, the parallel polynomial evaluation circuit thus requires

$\begin{array}{c}{T}_{pp}\left(n,d,p,M\right)=T_{poly}\left(n,d,p\right)+d·T_{extra}\left(M\right)+T_{\mathrm{label}}\left(M,n\right)\\ =\frac{3}{2}{n}^{2}d+3n\mathrm{pd}+{\frac{7}{2}}^{\mathrm{nd}}-3p^{2}d+3pd-d+3\mathrm{Md}\left(4\lceil {\log }_{2}M\rceil -8\right)+4\mathrm{Mn}\end{array}$ $\mathrm{Toffoli}\mathrm{gates}\mathrm{and}\left(d+1\right)n+\underset{2}{\lceil \log }M\rceil +1\mathrm{qubits}.$

The inverse square root, i.e.,

$f\left(x\right)=\frac{1}{\sqrt{x}}$

can be computed efficiently using Newton's method. The iteration looks as follows:

$x_{n+1}=x_n\left(1.5-\frac{a{x}_n^2}{2}\right),$

where a is the input and

$x_n{\to }^{n\to \infty }\frac{1}{\sqrt{a}}$

if the initial guess is sufficiently close to the true solution.

In some examples, finding a good initial guess

$x_0\approx \frac{1}{\sqrt{a}}$

for Newton's zero-finding routine is useful for (fast) convergence. A crude approximation which turns out to be sufficient is the following:

$\frac{1}{\sqrt{a}}={\left(2^{\underset{2}{\log a}}\right)}^{-\frac{1}{2}}=2^{-\frac{\underset{2}{\log a}}{2}}\approx 2^{\lfloor -\frac{\lfloor \underset{2}{\log a}\rfloor }{2}\rceil }={\tilde{x}}_0,$

where

$\lfloor \underset{2}{\log a}\rfloor$

can be determined by tinning the first “1” when traversing the bit-representation of a from left to right (MSB to LSB). While the space requirement for {tilde over (x)}₀ is in

$𝒪\left(\underset{2}{\log }n\right),$

such a representation would be impractical for the first Newton round. Furthermore, noting that the first iteration on {tilde over (x)}₀=2^k leads to

$\begin{array}{cc}{\tilde{x}}_1=2^k\left(1.5-\frac{a{2}^{2k}}{2}\right)\text{=:}{x}_0,& \left(4\right)\end{array}$

one can directly choose this x₀ as the initial guess.

The preparation of x₀ can be achieved using (n−1)+n+1 ancilla qubits, which must be available due to the space requirements of the subsequent Newton steps. The one ancilla qubit is used as a flag indicating whether the first “1” from the left has already been encountered. For each iteration i ∈ {n−1, . . . , 1, 0}, one determines whether the bit a_i is 1 and stores this result r_i in one of the n work qubits, conditioned on the flag being unset.

Then, conditioned on r_i=1, the flag is flipped, indicating that the first “1” has been found. If r_i=1, the x₀-register is initialized to the value in (4) as follows: using CNOTs, the x₀-register can be initialized to the value 1.5 shifted by

$k=\frac{p-2i}{2},$

where p denotes the binary point position of the input, followed by subtracting the (3k−1)-shifted input a from x₀, which may require up to n−1 ancilla qubits.

In order to improve the quality of the first guess for numbers close to 2^k for some k ∈ , one can tune the constant 1.5 in ( ), i.e., turn it into a function C(k) of the exponent k. This increases the overall cost of calculating x₀ merely by a few CNOT gates but allows to save an entire Newton iteration even when only distinguishing three cases, namely

$\begin{array}{cc}C\left(k\right):=\{\begin{array}{cc}1.613,& k<0\\ 1.5,& k=0\\ 1.62,& k>0\end{array}.& \left(5\right)\end{array}$

The Newton iteration may include computing x_n+1 from x_n by

$x_{n+1}=x_n\left(1.5-\frac{a{x}_n^2}{2}\right),$

which can be achieved as shown in Table 9 below:

TABLE 9
1. Compute the square of xn into a new register.
2. Multiply xn2 by the shifted input to obtain axn2/2.
3. Initialize another register to 1.5 and subtract axn2/2.
4. Multiply the result by xn to arrive at xn+1.
5. Uncompute the three intermediate results.

The circuit of one such Newton iteration is depicted in FIG. 11. SQR 1110 computes the square of previous iterate x_n into an empty result-register, which is then multiplied by the input a (MUL 1120), followed by subtracting (SUB 1140) this intermediate result from the value 1.5 initialized by block 1130. Finally, the next iterate, i.e.,

$x_{n+1}=x_n\left(1.5-\frac{a{x}_n^2}{2}\right)$

can be computed by multiplying this intermediate result by x_n using MUL 1150. All temporary results may then be cleared by running the appropriate operations in reversed order (using SQR 1190, MUL 1170, SET 1180, and ADD 1160).

Therefore, for m Newton iterations, this requires m+3 n-qubit registers if no pebbling is done on the Newton iterates, i.e., if all x_i are kept in memory until the last Newton iteration has been completed.

Computing the initial guess for the fast inverse square root requires n controlled additions of two n-bit numbers plus 2n Toffoli gates for checking/setting the flag (and uncomputing it again). Thus, the Toffoli count for the initial guess is

T _{init(n)=nT cadd(n)+2n=3n 2 +5n.}

Each Newton iteration features squaring, a multiplication, a subtraction, a final multiplication (yielding the next iterate), and then an uncomputation of the three intermediate results. In total, 5 multiplications and 2 additions are employed (of which 2 multiplications and 1 addition are run in reverse), which yields the Toffoli count

$\begin{array}{c}{T}_{\mathrm{iter}}\left(n,p\right)=5T_{\mathrm{mul}}\left(n,p\right)+2T_{\mathrm{add}}\left(n\right)\\ =\frac{15}{2}{n}^2+15np+\frac{23}{2}n-15p^2+15p-2.\end{array}$

The number of Toffoli gates for the entire Newton procedure (without un-computing the iterates) for m iterations thus reads

$\begin{array}{c}{T}_{\mathrm{invsqrt}}\left(n,m,p\right)=T_{\mathrm{init}}\left(n\right)+m{T}_{\mathrm{iter}}\left(n,p\right)\\ =n^2\left(\frac{15}{2}m+3\right)+15npm+n\left(\frac{23}{2}m+5\right)-15p^{2}m+15\mathrm{pm}-2m.\end{array}$

Since each Newton iteration requires 3 ancilla registers (which are cleaned up after each round) to produce the next iterate, the total number of qubits is n(m+4), where one register holds the initial guess x₀.

Note that this is an upper bound on the required number of both qubits and Toffoli gates. Since Newton converges quadratically, there is no need to perform full additions and multiplications at each iteration. Rather, the number of bits n used for the fixed point representation should be an (increasing) function of the Newton iteration.

The square root can be calculated using

$\sqrt{x}=x·\frac{1}{\sqrt{x}},$

i.e., at a cost of an additional multiplication into a new register. Note that this new register would be required anyway when copying out the result and running the entire computation in reverse, in order to clear registers holding intermediate results. Thus, the total number of logical qubits remains unchanged.

While sin(x) and cos(x) may be easily approximated using, e.g., polynomials as described, their inverses are not. The main difficulty arises near ±1, where

$\frac{d\mathrm{arc}\sin \left(x\right)}{dx}=\frac{1}{\sqrt{1-x^2}}$

diverges.

Therefore, an alternative representation of arc sin(x) for larger values of x is used, e.g.,

$\begin{array}{c}\arcsin \left(x\right)=\frac{\pi }{2}-\mathrm{arc}\cos \left(x\right)\\ =\frac{\pi }{2}-\arcsin \left(\sqrt{1-x^2}\right).\end{array}$

Applying the double-argument identity to the last expression yields

$\begin{array}{cc}\arcsin \left(x\right)=\frac{\pi }{2}-2\arcsin \left(\sqrt{\frac{1-x}{2}}\right),& \left(6\right)\end{array}$

a very useful identity which was already used in a classical math library called Cephes. The interval is again partitioned, using a minimax polynomial to approximate arc sin(x) for x ∈ [0,0.5), and the transformation in (6) for x ∈ [0.5,1].

The inverse square root implementation is used to compute √{square root over (z)} for

$z=\frac{1-x}{2},$

which satisfies z ∈ [0,0.25], for x ∈ [0.5,1]. Therefore, the fixed point position has to be chosen large, as the inverse square root diverges for small x. The multiplication by x after this computation takes care of the singularity and, since most bits of low-significance of

$\sqrt{\frac{1}{x}}$

will cause underflow for small x, a shifted version of the inverse square root can be computed. This optimization reduces the number of extra bits required during the evaluation of the inverse square root.

In many applications, evaluating arc sin(x) only on the interval [0,0.5] may be sufficient. In such cases, the cost is much lower since this can be achieved using our parallel polynomial evaluation circuit. The Toffoli counts for this case are noted in Table 4.

In one example, the Arcsine is implemented as a combination of polynomial evaluation and the inverse square root to extend the polynomial approximation on [0,0.5] to the entire domain [0,1] employing the double-argument identity above.

First, the (pseudo) sign-bit of x is copied out and x is conditionally inverted (modulo two's-complement) to ensure x≥0. Since there are plenty of registers available, this can be achieved by conditionally initializing an extra register to 1 and then using a normal adder to increment x by one, where x denotes the bit- or one's-complement of x.

Since x ∈ [0,1], one can determine whether x<0.5 using just one Toffoli gate (and 4 NOT gates). The result of this comparison is stored in an ancilla qubit denoted by a. z=(1−x)/2 can be computed using an adder (run in reverse) acting on x shifted by one and a new register, after having initialized it to 0.5 using a NOT gate.

Then, conditioned on ā (i.e., on a being 0), this result is copied into the polynomial input register p_in and, conditioned on a, x is squared into p_in. After having applied the polynomial evaluation circuit shown in FIG. 1 (which un-computes intermediate results) to this input, p_in can be uncomputed again, followed by computing the square root of z.

Then, the result of the polynomial evaluation is multiplied by either √{square root over (z)} or x, which can be achieved using 2n controlled swaps and one multiplier. The final transformation of the result consists of an initialization to π/2 followed by a subtraction, both conditioned on ā, and a copy conditioned on a. Finally, the initial conditional inversion of x can be undone after having (conditionally) inverted the output.

Following this procedure, the Toffoli count for this arcsine implementation on n-bit numbers using m Newton iterations for calculating √{square root over (z)} and a degree-d polynomial to approximate arc sin(x) on [0,0.5] can be written as

$\begin{array}{c}{T}_{\arcsin }=3T_{\mathrm{inv}}+\left(2T_{\mathrm{poly}}-T_{\mathrm{fma}}\right)+2T_{\mathrm{csquare}}+T_{\mathrm{mul}}+T_{\mathrm{cadd}}+\left(2T_{\mathrm{invsqrt}}+T_{\mathrm{mul}}\right)+5n+2+T_{\mathrm{add}}\\ =3T_{add}+2T_{\mathrm{poly}}+3T_{\mathrm{mul}}+T_{\mathrm{cadd}}+2T_{\mathrm{invsqrt}}+9n+2\\ =d\left(3n^2+n\left(6p+7\right)-6\left(p-1\right)p-2\right)+m\left(n\left(15n+30p+23\right)-30p\left(p-1\right)-4\right)+9\left(n+1\right)p+\frac{9}{2}n\left(n+7\right)+6n^2+28n-9p^2+2\end{array}$

where T_inv(n) denotes the Toffoli count for computing the two's-complement of an n-bit number and T_{csquare(n,p)=Tmul(n,p)+2n} is the number of Toffoli gates required to perform a conditional squaring operation. Furthermore, 2n Toffoli gates are needed to achieve the conditional n-bit swap operation (twice), and another 3n are used for (conditional) copies.

All circuits were implemented at the gate level and tested using a reversible simulator extension to LIQUID.

A summary of the required resource for implementing tanh(x), exp(−x²), and sin(x) can be found in Tables 5-8. For each function, one set of parameters was implemented reversibly at the level of Toffoli gates in order to verify the proposed circuits.

The convergence of the reversible fast inverse square root implementation with the number of Newton iterations are shown in FIG. 9, where the bit sizes and point positions have been chosen such that the roundoff errors do not interfere significantly with the convergence.

For all practical purposes, choosing between 3 and 5 Newton iterations should be sufficient. The effect of tuning the constants in the initial guess (see Eqn. 5) can be seen when comparing FIG. 8 to FIG. 9: The initial guess is obtained from the location of the first non-zero in the bit-representation of the input, which results in large rounding-effects for inputs close to an integer power of two.

Tuning the initial guess results in almost uniform convergence, which saves an entire Newton iteration for a given L_∞-error.

The square root converges better than the inverse square root for small values, which can be expected, since

$\sqrt{x}=x·\frac{1}{\sqrt{x}}$

has a regularizing effect for small x. The error after m Newton iterations when using n bits for the fixed point representation is depicted in FIG. 5. Additionally, the initial guess could be improved by tuning the constants in Eqn. 4 such that the error is minimal after multiplying

$x·\frac{1}{\sqrt{x}},$

instead or just optimizing ror the inverse square root itself.

The example implementation of Arcsine uses both the polynomial evaluation and square root subroutines. The oscillatory behavior which can be seen in FIG. 8 is typical for minimax approximations. For x>0.5, the resolution is lower due to the wider range of

$\frac{1}{\sqrt{x}},$

which was accounted for by calculating a shifted version of the inverse square root. FIG. 12 shows absolute error on [0,5] for N=2000 equidistant points of the example reversible implementation of the square root for m ∈ {2,3,4} Newton iterations and corresponding bit sizes n ∈ {25,35,50}. The fixed-point position is chosen to be p=5.

While this saves a few qubits (to the left of the binary point), the reduced number of qubits to the right of the binary point fail to resolve the numbers as well, which manifests itself by bit-noise for x>0.5 in FIG. 3. The degrees of the minimax approximation were chosen to be 7, 13, and 17 for m=3,4, 5, respectively. Since arc sin(x) is an odd function, this amounts to evaluating a degree 3, 6, and 8 polynomial in x², followed by a multiplication by x.

In conclusion, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include determining a polynomial interpolation for a set of sub-intervals corresponding to the input. The method may further include constructing a quantum circuit for performing, in parallel, polynomial evaluation corresponding to each of the set of sub-intervals.

The method may further comprise determining an overall achievable precision for any input value that falls into any of the set of intervals. The method may further comprise partitioning the input into the set of sub-intervals by determining whether an input value falls into any of the set of sub-intervals and labeling the input value using a label indicative of one of the set of sub-intervals, where the sub-intervals may be chosen by pre-processing the input.

The polynomial evaluation may be performed in a nested fashion. The polynomial evaluation may be performed by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficients to the output.

In another example, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include determining a polynomial interpolation of an initial degree d for an initial set of sub-intervals corresponding to an input. The method may further include determining whether a required target precision is achievable based on the initial degree d and the set of sub-intervals. The method may further include if the target precision is achievable, then constructing a quantum circuit for performing, in parallel, polynomial evaluation corresponding to each of the set of sub-intervals. The method may further include if the target precision is not achievable, then either incrementing the initial degree d or subdividing the initial set of subintervals into additional intervals until the target precision is achievable and then constructing the quantum circuit.

The method may further include partitioning the input into the initial set of sub-intervals by determining whether an input value falls into any of the initial set of sub-intervals and labeling the input value using a label indicative of one of the initial set of sub-intervals. The initial set of sub-intervals may be chosen by pre-processing the input.

The polynomial evaluation may be performed in a nested fashion. The polynomial evaluation may be performed by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficients to the output.

The input may be divided into the initial set of sub-intervals to reduce an approximation error associated with the polynomial evaluation. At least a subset of the initial set of intervals may be fused together to reduce a quantum cost associated with the polynomial evaluation.

In yet another example, the present disclosure relates to a processor-implemented method for performing an evaluation of a polynomial corresponding to an input. The method may include for each interval of a set of sub-intervals corresponding to the input, determining whether to perform the evaluation of the polynomial in a first mode or a second mode. The method may further include, in the first mode, performing the evaluation of the polynomial by: determining a polynomial interpolation of at most degree d for a first subset of the set of sub-intervals corresponding to the input, and constructing a quantum circuit for performing polynomial evaluation corresponding to each of the first subset of the set of sub-intervals. The method may further include, in the second mode, performing the evaluation of the polynomial using Newton's method for a second subset of the set of sub-intervals corresponding to the input, wherein the second subset is different from the first subset.

The method may further include determining an overall achievable precision for any input value that falls into any of the set of intervals. The method may further include partitioning the input into the set of sub-intervals by determining whether an input value falls into any of the set of sub-intervals and labeling the input value using a label indicative of one of the set of sub-intervals.

The polynomial evaluation may be performed by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficient with the output. The iterations corresponding to the Newton's method may be implemented using a result of a pebbling strategy to meet a resource requirement corresponding to the quantum circuit. The resource requirement may be a number of available qubits corresponding to the quantum circuit. The resource requirement may be a circuit depth corresponding to the quantum circuit.

It is to be understood that the methods, modules, and components depicted herein are merely exemplary. Alternatively, or in addition, the functionally described herein can be performed, at least in part, by one or more hardware logic components. For example, and without limitation, illustrative types of hardware logic components that can be used include Field-Programmable Gate Arrays (FPGAs), Application-Specific Integrated Circuits (ASICs), Application-Specific Standard Products (ASSPs), System-on-a-Chip systems (SOCs), Complex Programmable Logic Devices (CPLDs), etc. In an abstract, but still definite sense, any arrangement of components to achieve the same functionality is effectively “associated” such that the desired functionality is achieved. Hence, any two components herein combined to achieve a particular functionality can be seen as “associated with” each other such that the desired functionality is achieved, irrespective of architectures or inter-medial components. Likewise, any two components so associated can also be viewed as being “operably connected,” or “coupled,” to each other to achieve the desired functionality.

The functionality associated with some examples described in this disclosure can also include instructions stored in a non-transitory media. The term “non-transitory media” as used herein refers to any media storing data and/or instructions that cause a machine to operate in a specific manner. Exemplary non-transitory media include non-volatile media and/or volatile media. Non-volatile media include, for example, a hard disk, a solid-state drive, a magnetic disk or tape, an optical disk or tape, a flash memory, an EPROM, NVRAM, PRAM, or other such media, or networked versions of such media. Volatile media include, for example, dynamic memory, such as, DRAM, SRAM, a cache, or other such media. Non-transitory media is distinct from, but can be used in conjunction with transmission media. Transmission media is used for transferring data and/or instruction to or from a machine. Exemplary transmission media include coaxial cables, fiber-optic cables, copper wires, and wireless media, such as radio waves.

Furthermore, those skilled in the art will recognize that boundaries between the functionality of the above described operations are merely illustrative. The functionality of multiple operations may be combined into a single operation, and/or the functionality of a single operation may be distributed in additional operations. Moreover, alternative embodiments may include multiple instances of a particular operation, and the order of operations may be altered in various other embodiments.

Although the disclosure provides specific examples, various modifications and changes can be made without departing from the scope of the disclosure as set forth in the claims below. Accordingly, the specification and figures are to be regarded in an illustrative rather than a restrictive sense, and all such modifications are intended to be included within the scope of the present disclosure. Any benefits, advantages, or solutions to problems that are described herein with regard to a specific example are not intended to be construed as a critical, required, or essential feature or element of any or all the claims.

Furthermore, the terms “a” or “an,” as used herein, are defined as one or more than one. Also, the use of introductory phrases such as “at least one” and “one or more” in the claims should not be construed to imply that the introduction of another claim element by the indefinite articles “a” or “an” limits any particular claim containing such introduced claim element to inventions containing only one such element, even when the same claim includes the introductory phrases “one or more” or “at least one” and indefinite articles such as “a” or “an.” The same holds true for the use of definite articles.

Unless stated otherwise, terms such as “first” and “second” are used to arbitrarily distinguish between the elements such terms describe. Thus, these terms are not necessarily intended to indicate temporal or other prioritization of such elements.

Claims

1.-20. (canceled)

21. A processor-implemented method for performing an evaluation of a polynomial, the method comprising: receiving an input comprising a polynomial function and a set of sub-intervals;for each interval of a set of sub-intervals corresponding to the input, determining whether to perform the evaluation of the polynomial in a first mode or a second mode;in the first mode, performing the evaluation of the polynomial by: determining a polynomial interpolation of at most degree d for a first subset of the set of sub-intervals corresponding to the input, andusing the processor, generating an output including the quantum circuit for performing polynomial evaluation corresponding to each of the first subset of the set of sub-intervals by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficients with the output; andin the second mode, performing the evaluation of the polynomial using Newton's method for a second subset of the set of sub-intervals corresponding to the input, wherein the second subset is different from the first subset.

22. The processor-implemented method of claim 21, further comprising determining an overall achievable precision for any input value that falls into any of the set of intervals.

23. The processor-implemented method of claim 21, further comprising partitioning the input into the set of sub-intervals by determining whether an input value falls into any of the set of sub-intervals and labeling the input value using a label indicative of one of the set of sub-intervals.

24. The processor-implemented method of claim 21, wherein the polynomial evaluation is performed in a nested fashion.

25. The processor-implemented method of claim 21, further comprising determining a number of rounds of the Newton's method.

26. The processor-implemented method of claim 21, wherein iterations corresponding to the Newton's method are implemented using a result of a pebbling strategy to meet a resource requirement corresponding to the quantum circuit.

27. The processor-implemented method of claim 26, wherein the resource requirement is a number of available qubits corresponding to the quantum circuit.

28. The processor-implemented method of claim 26, wherein the resource requirement is a circuit depth corresponding to the quantum circuit.

29. A processor-implemented method for performing an evaluation of a polynomial, the method comprising: receiving an input comprising a polynomial function and a set of sub-intervals;for each interval of a set of sub-intervals corresponding to the input, determining whether to perform the evaluation of the polynomial in a first mode or a second mode;in the first mode, performing the evaluation of the polynomial by: determining a polynomial interpolation of at most degree d for a first subset of the set of sub-intervals corresponding to the input, andusing the processor, generating an output including the quantum circuit for performing polynomial evaluation corresponding to each of the first subset of the set of sub-intervals by determining whether an input value of the polynomial function falls into any of the set of sub-intervals and computing a label register having at least one label indicative of at least one of the set of sub-intervals; andin the second mode, performing the evaluation of the polynomial using Newton's method for a second subset of the set of sub-intervals corresponding to the input, wherein the second subset is different from the first subset.

30. The processor-implemented method of claim 29, further comprising determining an overall achievable precision for any input value that falls into any of the set of intervals.

31. The processor-implemented method of claim 29, wherein the polynomial evaluation is performed by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficients to the output.

32. The processor-implemented method of claim 29, wherein the polynomial evaluation is performed in a nested fashion.

33. The processor-implemented method of claim 29, wherein iterations corresponding to the Newton's method are implemented using a result of a pebbling strategy to meet a resource requirement corresponding to the quantum circuit.

34. The processor-implemented method of claim 33, wherein the resource requirement is a number of available qubits corresponding to the quantum circuit.

35. The processor-implemented method of claim 33, wherein the resource requirement is a circuit depth corresponding to the quantum circuit.

36. A processor-implemented method for performing an evaluation of a polynomial, the method comprising: receiving an input comprising a polynomial function and a set of sub-intervals;for each interval of a set of sub-intervals corresponding to the input, determining whether to perform the evaluation of the polynomial in a first mode or a second mode, wherein the set of sub-intervals are chosen by pre-processing the input;in the first mode, performing the evaluation of the polynomial by: determining a polynomial interpolation of at most degree d for a first subset of the set of sub-intervals corresponding to the input, andusing the processor, generating an output including the quantum circuit for performing polynomial evaluation corresponding to each of the first subset of the set of sub-intervals by determining whether an input value of the polynomial function falls into any of the set of sub-intervals and computing a label register having at least one label indicative of at least one of the set of sub-intervals; andin the second mode, performing the evaluation of the polynomial using Newton's method for a second subset of the set of sub-intervals corresponding to the input, wherein the second subset is different from the first subset.

37. The processor-implemented method of claim 36, further comprising determining an overall achievable precision for any input value that falls into any of the set of intervals.

38. The processor-implemented method of claim 36, wherein the polynomial evaluation is performed by applying a gate to load a next set of coefficients into a quantum register, while multiplying a previous set of coefficients by a variable to produce an output, and then adding the next set of coefficients to the output.

39. The processor-implemented method of claim 36, wherein the polynomial evaluation is performed in a nested fashion.

40. The processor-implemented method of claim 36, wherein iterations corresponding to the Newton's method are implemented using a result of a pebbling strategy to meet a resource requirement corresponding to the quantum circuit.