{'fraction_non_alphanumeric': 0.03591532355063748, 'fraction_numerical': 0.021361558816454174, 'mean_word_length': 4.598033934823593, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Benchmarking 50-Photon Gaussian Boson Sampling on the Sunway TaihuLight

Boson sampling is expected to be one of an important milestones that will demonstrate quantum supremacy. The present work establishes the benchmarking of Gaussian boson sampling (GBS) with threshold detection based on the Sunway TaihuLight supercomputer. To achieve the best performance and provide a competitive scenario for future quantum computing studies, the selected simulation algorithm is fully optimized based on a set of innovative approaches, including a parallel scheme and instruction-level optimizing method. Furthermore, data precision and instruction scheduling are handled in a sophisticated manner by an adaptive precision optimization scheme and a DAG-based heuristic search algorithm, respectively. Based on these methods, a highly efficient and parallel quantum sampling algorithm is designed. The largest run enables us to obtain one Torontonian function of a 100 x 100 submatrix from 50-photon GBS within 20 hours in 128-bit precision and 2 days in 256-bit precision.

implementation, instantaneous quantum polynomial time circuit and RCS are based on quantum bits, while boson sampling and GBS are based on bosons, such as single photons. Since the proposal of GBS greatly simplifies its quantum implementation, quantum sampling based on bosons is expected to be an important milestone in demonstrating quantum supremacy [1], [13].
This work establishes the quantum supremacy frontier for the boson-based quantum sampling problem. To achieve the best performance and provide a competitive scenario for quantum computers, the selected algorithm is fully optimized based on the Sunway TaihuLight supercomputer, which is one of the most powerful classical computers in the world. Our major contributions are as follows: • an effective parallel scheme is proposed to reduce the requirements of cache-level storage and obtain an almost ideal load balance of the selected boson sampling algorithm.
• an instruction-based and multiple-precision optimizing scheme that enables customizable data precision modes is designed to guarantee sufficient accuracy of the simulation; an adaptive framework based on upper-and lower-bound estimation is further applied to determine the best precision of the configuration.
• in terms of instruction-level optimization, a DAGbased heuristic search algorithm is adopted to achieve optimal instruction scheduling for all kernels.

RELATED WORKS
In 2019, Google first claimed the quantum supremacy by using 53 superconducting quantum bits in an RCS experiment [14]. The quantum device took 200 seconds to sample one instance of a quantum circuit a million times, while classical benchmarking was expected to need 10,000 years arXiv:2009.01177v1 [cs.DC] 2 Sep 2020 on the Summit supercomputer. To confirm and establish this quantum supremacy, there have been many developments in the classical benchmarking algorithm [15], [16], [17], [18], [19], [20] before and after the experiments. The latest result found by IBM indicates that the calculation can be run on the same supercomputer in less than two and a half daysrather than 10,000 years [21], [22]. Furthermore, by leveraging secondary storage of supercomputer, the sampling probabilities of all the 2 53 states can be parallelly computed. For large number of samples, the average classical calculation time for a single sample will be 2 53 times faster. Therefore, a more reliable state-of-the-art benchmarking on classical computers is crucial to prove the quantum supremacy.
Today, many studies have been based on the RCS problem. In general, the classical simulators of RCS problems are categorized into three classes. The first class and the second class are the direct evolutions of the quantum state [15], [23], [24], [25], [26] and the perturbation of stabilizer circuits [27], [28], [29], respectively. The tensor network contraction class [18], [19], [30] is the most suitable method for current flop-oriented architectures such as on the Fugaku [31] and Summit [32] supercomputers. In addition, several hybrid algorithms [16], [21], [33] show their great performance in achieving a simulation or benchmark with ∼ 50 quantum bits. There are also some efforts [24], [34] that attempt to build up the benchmarking of the quantum supremacy based on the RCS problem on the Sunway TaihuLight supercomputer.
The boson sampling problem, which was introduced by Aaronson and Arkhipov [35], is historically the first protocol to conclusively demonstrate the quantum supremacy. Several algorithms have been designed and experiments have been conducted referring to such approaches [36], [37], [38], [39], [40]. However, experimentally, it is actually difficult to scale the boson sampling to high photon numbers owing to its intrinsic photon cost. GBS [41] and GBS with threshold detection [12] were recently proposed to obtain higher-order photon numbers than those of boson sampling. There are several classical benchmarks for GBS. Based on the Titan supercomputer, Brajesh Gupt et al. [42] proposed a benchmark for GBS with threshold detection. However, due to the limitations of the method, it can only handle approximately 22 clicks (or photons) when using the entire Titan supercomputer.
Different from previous works, we establish the first classical benchmarking for GBS with threshold detection based on calculation of the Torontonian function, which has not been used in previous works. Our design is able to calculate a single Torontonian for 50-photon GBS in less than one day.

Selected Algorithm: Boson Sampling
In the standard boson sampling experiment, N indistinguishable single photons are sent to an M -port linear optical network, and the output scattered photons are detected by N single-photon detectors. The probabilities of these Nphoton events are related to the permanent function of an N × N sampling matrix, which is proven to be #P-hard for classical computers. However, a large-scale experiment suffers from the formidable challenge of preparing N perfect single photons [36], [38], [39], [43], [44], [45], [46], [47].  A practical improvement was proposed to change the input single photons to nonclassical Gaussian light and use single-photon detectors without photon-number resolution to register N -click events, which is called GBS with threshold detection [11], [12], [48], [49]. In this case, the probability of N -click events is the Torontonian function of a 2N × 2N sampling matrix, which is expected to still be exponentially hard for classical computers.
In a GBS experiment, before photon detection, the output Gaussian state is described by a covariance matrix . The sampling matrix A is then defined as A = I − −1 , and the probability p(S) of an N -click event S = s 1 , s 2 , ..., s M is determined by , where ρ is the output quantum state of a covariance matrix , S is the threshold click pattern with s k ∈ 0, 1 and M k=1 = N , and A {S} is the 2N × 2N submatrix of A by selecting the k-th and (M + k)-th row and column when s k = 1. The matrix function T or(·) is the Torontonian function, which is defined as , where P N is the power set of 1, 2, ..., N , and hence, there are 2 N determinant terms in the summation [12]. The computational complexity of the Torontonian function is O(2 N ), which is the same as the permanent function in standard boson sampling with N photons. In this work, our major task is to calculate the Torontonian function-which has exponential complexity.

Sunway Architecture
To achieve the best simulation performance, this work selects the Sunway TaihuLight supercomputer as the classical platform. The system consists of 40,960 many-core processors called SW26010s and is able to provide a peak performance of 125 PFlops and a sustainable performance of 93 PFlops. Each SW26010 processor consists of four core groups (CGs), which are shown in Figure 2. Each CG has a memory controller (MC), a management processing element (MPE) and 64 computing processing elements (CPEs).  As shown in Figure 2, the MPE is a fully functional 64bit RISC general-purpose core that contains a 32-KB L1 data cache and 256-KB L2 instruction/data cache. The CPE is a simplified 64-bit RISC computing-oriented core and has a private 16-KB L1 instruction cache and 64-KB scratchpad memory (SPM). Direct memory access (DMA) is the major method for the CPEs to access data from global memory. Furthermore, a thread-level communication is enabled based on register communication. Thus, inside one CG, the CPE cores within the same row or column are able to communicate with each other directly.
A distinctive feature of the Sunway architecture is its 256-bit large integer instruction sets, as shown in Table  1. The first five instructions (multiplication, addition, subtraction and shifting) belong to the 256-bit large integer instruction sets. However, the multiplication instruction is not a straightforward 256-bit version; instead, it performs an unsigned multiplication of two 128-bit unsigned numbers and obtains a full 256-bit result. The rest of the instructions in the table are regular vector instructions and are useful in optimizing the boson sampling process. The parameters of instruction vshff is a, b and mask are shown in Figure 3a. a and b are two 256bit registers that contain four 64-bit numbers. mask provides information on how to construct a new register from a (the first two digits of mask) and b (the last two digits of mask). Figure 3a shows an example: based on the four digits of mask, 0 and 2 of A (corresponding value a0 and a2), and position 0 and 1 of B (corresponding value b0 and b1) are selected to construct a new register.
The conditional selection instruction vsellt takes three parameters a, b, and c, which are all 256-bit registers composed of four 64-bit numbers. a is the conditional register that chooses between b and c. As shown in Figure 3b, if one of the 64-bit numbers in a is less than 0, then the corresponding 64bit number of the result comes from the corresponding bits of b and vice versa.

PARALLEL SCHEME DESIGN
This section introduces the parallel scheme of solving Equation 2, including partition and storage strategies for load balancing and the LDM capacity issues.

Partition Strategy
For Equation 2, the best dimension to partition is the enumeration of Z due to its embarrassing parallelism. For convenience, we map each Z to a number mask "mask z " ranging from 1 to 2 N − 1. If i ∈ Z, the i-th bit in the binary representation of mask z is "1" and vice versa.
A naive partition strategy can be directly obtained by evenly dividing the 1 ∼ 2 N − 1 mask into every CPE from different CGs. However, while two continuous numbers must have different numbers of "one" in the binary representation, the matrix size and computing load of each CPE are different. Therefore, the simple and naive strategy brings a severe thread-level load imbalance.

Algorithm 1 Calculate T orontonian(A)
Require: Matrix A Ensure: A is Hermitian positive definite 1: function T orontonian(A) 2: result ← 0 3: To avoid this situation, we design a new partition strategy with high load balancing. We first divide the workloads of matrices into N parts according to the matrix size |Z|. In each part, we divide all masks evenly and continuously into every process. Algorithm 1 describes the Torontonian function with the new partition strategy. Line 3 indicates the enumerating of N parts; in Line 4 and Line 5, the masks assigned to a process are determined; in Line 7, the calculating matrix A Z is generated by the current mask; Line 8 is the determinant calculation that is the most expensive part; Line 9 is the update of the result; and in Line 10, the next mask is derived from the current mask, and then, the next loop begins. In a certain part, each CPE processes a matrix of the same size, obtaining the best load balance.  17: cto ← COUNT TAIL ONE(mask) 18:

Algorithm 2 Get Matrix Mask
x return masknext 22: end function 23: In the naive strategy, we only need to calculate mask z = 2 N * rank/nproc to obtain the first mask of a certain process (or obtain the k-th mask) and to calculate mask z +1 to obtain the next mask. nproc denotes the total number of processes, and rank denotes the rank of a certain process.
However, the two functions for obtaining the mask, the pseudocodes of which are shown in Algorithm 2, are rather complicated in the new strategy. The get kth mask function is based on a bit-by-bit algorithm that determines the position of every "1" from high to low in binary representation. In each round of the outer loop, the position of a "1" is determined. In the inner loop, we try each possible position one by one according to the comparison results of the combinatorial enumeration "cnt" and the current "k". If cnt is greater than "k", then we set the position to "1" and begin the next round of the outer loop; otherwise, we subtract cnt from "k". The main task of the get next mask function is to move the consecutive "1"s at the tail of the binary representation to a proper position. For instance, consider that the current mask is "010011". To obtain the next mask "001110", we need to move the last two "1"s to the fourth position (beginning with 1) and move the first "1" to the next position. Therefore, a bit operation algorithm is designed. cto and ctz represent the number of consecutive ones and zeros at the tail, respectively. cto is used to count the consecutive "1"s at the tail, while ctz is used to indicate the proper position.
Such a partition scheme is also applicable to thread-level parallelisms with trivial modifications.

Storage Strategy
Neither the input nor the calculating matrix can be fit to the 64 KB LDM space if the scale N is too large. For instance, when N is 45 based on 256-bit multiple precision, the storage requirement of a matrix is (2N ) × (2N ) × sizeof(complex) = 90 * 90 * 64 = 506.25 KB, and far exceeds the LDM capacity. In addition, when reserving some LDM space to store the input matrix to avoid extra DMA copies, the insufficiency of capacity can be even worse. Therefore, we propose a storage strategy to lower the LDM requirement of both the input matrix A and calculating matrix A Z .
For A, we use lossy compression and a matrix symmetric. As scientific measurement data, there are only three to five significant digits in each element of matrix A. Therefore, a 16-bit number (with an error of ±1.525×10 −5 ) is sufficient to describe the measurement data. Additionally, in our implementation, 32-bit lossy compression is supported for higher-precision requirements.
Another reason that brings storage improvements is the usage of a matrix symmetric. Owing to the physical properties of the matrix A, if we represent matrix A in the form of a block matrix with four blocks as then the following three symmetries are satisfied: Now, only half of A 00 and A 01 need to be stored, and thus, only approximately 1/4 of the LDM space is required. For calculating the matrix A Z , we use the Hermitian symmetric and distribute the matrix into several CPEs to save the LDM space. Since A is Hermitian positive definite, it is not difficult to prove that I − A Z is also Hermitian positive definite, and in each round of row reduction, the matrix is always Hermitian. Thus, we can save half the space of the matrix when performing calculations, i.e., the upper triangular part of the matrix, which saves half of the LDM space.
It is still not enough to put the whole calculating matrix into LDM. The matrix must be scattered into different CPEs, which brings an extra cost when a CPE needs data from other CPEs. Fortunately, the register communication mechanism of the Sunway architecture provides an effective and efficient way to exchange data among different CPEs. Figure 4 illustrates our distribution strategy and the usage of register communication when the matrix is scattered to two CPEs. Each row of the matrix is alternately assigned to two CPEs. In the upper panel, the row with pivot is located at CPE 0. In this round, CPE 0 sends the row with pivot to CPE 1 via register communication. Then, each CPE multiplies this row by a proper scalar and adds the scaled row to other rows to finish reduction. In the next round of row reduction, as shown in the lower panel, the row with pivot comes to CPE 1, and thus, the data need to be transferred from CPE 1 to CPE 0. The process of reduction is the same as in the previous round.
Owing to the two strategies for the calculating matrix, the LDM is now enough to store all the required data. When N is 45 and the matrix is scattered to eight CPEs, each CPE is occupied 506.25 KB/2/8 = 31.64 KB LDM space, where 506.25 KB is mentioned in the beginning of this subsection, the division factor 2 is due to Hermitian symmetry, and 8 is due to the utilization of eight CPEs. In addition, when using 32-bit compression, the size of the input matrix is (2N ) * (2N ) * 2 * 4/4 = 16.17 KB. Therefore, the total LDM requirement is approximately 31.64 KB + 16.17 KB = 47.81 KB, which is less than the capacity limit of 64 KB.

Motivation
For real-world data, the result of the Torontonian function is a very small decimal, which results in a high-precision requirement. For instance, when N is 45, the result is 1.44 × 10 −25 . According to Equation 2, the result T or(A) is accumulated from the summation of massive reciprocal square roots of the determinant. We can prove that the value of each reciprocal square root is greater than 1. Thus, at least 25 decimal digits (roughly equivalent to 83 binary bits) is required to achieve adequate accuracy in the determinant calculation. According to the IEEE 754 standard [50], the double-precision type, which has only 53 binary bits of significant precision, is not competent to represent these numbers. Therefore, a customized multiple-precision type must be designed for the Torontonian function. To design such a precision type, we first consider the upper and lower bound analysis of the absolute value of intermediate results. The actual upper bound depends on the summation order of the Torontonian function. If all positive numbers are summed first and all negative numbers are added afterward, then the intermediate results will be very large and proportional to the times of summation (i.e., O(2 N )). However, the result of the Torontonian function is guaranteed to be a very small number. Thus, the overflow of the summation can be omitted. In other words, the part of the intermediate results greater than 1 can be truncated.
The upper bound, which cannot be omitted, is reached when calculating the reciprocal of pivot in the Gaussian elimination. The value of pivot is difficult to obtain, but fortunately, it can be proved to be smaller than its corresponding determinant, so we can replace all of them with the value of the determinant. Furthermore, it can be shown that the determinant of I − A is less than any of its principal minors, which means that the upper bound can be estimated by the reciprocal of the determinant of I − A. To make an estimation without precise computation, I − A can be calculated with this approximation with |Z| = N : , where k is a correction factor depending on matrix A that has an observed range of 0.0009 ∼ 0.0035, and a is the representative value of the diagonal element a ii . In realworld experiments, the range of a is 0.16±0.06. In addition, the value of the determinant when N = 60 will reach 10 −6 , which means that the upper bound can be estimated to be , where k and a are the same as in Equation 4. When k, a, and N are not too large (k < 0.0035, a < 0.2, N < 60), we have kN/(1 − a) 1 − a. Therefore, the approximation is close to (1/(1 − a) − 1) N , which can be guaranteed to be smaller than the actual result. Then, our approximation can be trusted to be a lower bound of the actual result.
With different k and a, there will be some variances in the value of the Torontonian function. Taking the range of these parameters into consideration, the lower bound can reach 10 −61 at N = 60.
According to the bound analysis, the range of intermediate results falls into a range of 10 −61 to 10 7 , which can be represented by a large fixed-point number. Therefore, we choose a fixed-point number as our precision type, which is much faster than a floating-point number.

Multiple-Precision Design
The large integer instruction set in the Sunway architecture is very suitable for implementing fixed-point numbers. As mentioned in Section 3.2, the Sunway processor supports 256-bit large integer operations, which implies that up to 256-bit fixed-point calculations can be implemented by hardware instructions directly or almost directly. According to the bound analysis, 256-bit fixed-point precision is sufficient for up to N = 60, which is an impossible task for a classical computer.
However, it must be noted that there is no 256-bit large integer multiplication in the architecture. The largest multiplication instruction is the 128-bit unsigned multiplication. Constructing the 256-bit signed multiplication requires numerous instructions based on the supported 128-bit multiplication. It is too expensive to perform the 256-bit signed multiplication.
To decrease the cost of multiplication, we propose an adjustable multiple-precision design. In addition to the 256bit precision, we mainly provide a precision mode of 128bit fixed-point numbers, which is much less expensive in terms of multiplication. The 128-bit precision mode is not competent for all cases but is much faster than the 256-bit precision mode when it works.
For the scaling factor, we just use a global factor for every fixed-point number. Unifying the scaling factor is beneficial to the performance efficiency.

Operation Design
In the determinant solver of the Torontonian function, O(2 N N 3 ) times of multiplications are required, which is the bottleneck of the whole calculation. Constructing signed multiplication based on the unsigned multiplication is the crucial challenge in multiple-precision design. A simple method is to convert negative numbers to positive numbers by using if statements. However, the if statements may result in a failure of branch predictions, which is a hazard to performance. Thus, avoiding if statements is considered the best way to enhance performance.
Take a 128-bit multiplication between 128-bit signed numbers a and b as an example. For convenience, we assume that a and b are both negative. The complements of a and b are expressed as (2 128 − |a|) and (2 128 − |b|). The unsigned multiplication is represented below: (2 128 −|a|) * (2 128 −|b|) = 2 256 −(|a|+|b|) * 2 128 +|a| * |b| (6) Therefore, the 256-bit result produced by unsigned multiplication is |a| * |b| − (|a| + |b|) * 2 128 , in which the term of 2 256 is truncated. The result is not equal to the correct one, which is |a| * |b|. It worth noting that if a and b are 256-bit signed numbers, then the second term of Equation 6 is (|a| + |b|) * 2 256 , which can be truncated. This indicates that the 256-bit results of the signed multiplication can be taken from the unsigned multiplication between two 256-bit signed numbers.
However, as mentioned above, there is no 256-bit multiplication instruction in the Sunway architecture. Instead, three times of 128-bit multiplication are required, i.e., a * b = a  return c 9: end function Algorithm 3 describes our 128-bit signed multiplication design. The first line describes 128-bit unsigned multiplication. The second and third lines are prepared for conditional selection instructions. By vshff instruction, the sign bit of a is extracted and put into high 128 bits of a sign . The low 128 bits of a sign is 0. In the next line, according to a sign , the vectorized conditional selection instruction vsellt is applied to make a choice between b and 0. Last, c is calculated and scaled in Line 7.
Therefore, in our method, a total of eight instructions are adequate for 128-bit signed multiplication. The 256bit signed multiplication is much more difficult but shares the same idea. In our implementation, 19 instructions are sufficient.
Division, which is used for calculating the scalar applied in the elimination, is the second expensive operation. Although only O(2 N N 2 ) times of divisions are required, the time complexity of division is much higher than that of multiplication. In our work, we use O(2 N N ) times of reciprocal calculations and O(2 N N 2 ) times of multiplications, instead of using division directly. In general, the reciprocal optimization results in a decrease in precision; however, our multiple precision design is enough to withstand the loss of accuracy. We choose trial division to implement the reciprocal. The time complexity is O(b), where b represents the bit length of the multiple-precision fixed-point number.
The last operation required is the reciprocal square root. It has the same time complexity with division but is only called O(2 N ) times, and thus, it is evidently not a bottleneck in the Torontonian function. Similar to division, a bitwise method is applied to the reciprocal square root.

Adaptive Precision Strategy
The final task of our precision design is to determine the scaling factor and bit length of the multiple-precision type (i.e., 128 bits or 256 bits). If the scaling factor is too small or the bit length is too short, then the results may have an overflow or insufficient significant digits. In contrast, if the bit length is too long, then the performance efficiency will be affected. In this subsection, we propose a strategy to select a proper scaling factor and bit length automatically.
The scaling factor depends on the upper bound of all intermediate results. As explained in Subsection 5.1, the upper bound can be represented by the reciprocal of determinant of I − A. As long as we calculate the determinant of I − A in advance, the upper bound can be obtained, hence the scaling factor. The choice of bit length is much more complicated and depends on several factors including both the upper bound and the lower bound. In our adaptive precision strategy, Equation 5 is chosen to approximate the trusted lower bound.
In addition, bits owing to the accumulated error should be taken into account. Since each element has N 2 times of elementary arithmetic in a determinant calculation, we assume that one time of determinant calculation will produce an error of αN 2 , where α is a correction factor. Meanwhile, 2 N times of determinant calculations are required. Thus, we use α2 N N 2 to estimate the accumulated error, where α is taken as 0.5 by some experiments.
Last but not least, the significant digits required (usually for physical experiments, this is three decimal digits) should be considered.
The final bit length can be represented as follows: , where B denotes the total bit length required, B lower denotes the part of the bit length derived from the upper bound (or scaling factor), B upper denotes the part derived from the upper bound, B sgn denotes the part derived from the significant digit requirement, and B accum denotes the part derived from the accumulated error. According to the relationship between B and 128, the 128-bit precision mode or 256-bit precision mode is chosen.

Framework
To obtain the optimal instruction scheduling, the corresponding hardware features must be considered carefully. For each CPE from the Sunway processor, there are two hardware pipelines with different functions. The first pipeline mainly takes charge of floating-point or vector instructions, and the second pipeline mainly takes charge of memory access and branch instructions. In addition, integer scalar instructions are processed by both pipelines. In our multiple-precision design, most instructions are large integer instructions, which belong to vector instructions and are thereby taken charge of by the first pipeline. Only about 15% of the instructions belong to the second pipeline in each kernel.
The performance (or the pipeline stall) of the first pipeline determines the overall performance. Two hazards affect the pipeline efficiency in the Sunway architecture. One is the read-after-write (RAW) data hazard; the other is the writeback stage structural hazard, which conflicts with the writeback of general-purpose registers. Given both hazards, it is hard to find an optimal instruction scheduling directly. Thus, we turn to optimization methods to find the optimal instruction scheduling.
For the instructions assigned to the second pipeline, we manually insert them into the scheduled sequence of instructions that belong to the first pipeline. To eliminate computation dependence and avoid data hazards, software pipeline technology is applied. As a result, these inserted instructions do not increase the number of total clocks in our case. In other words, all time cost comes from the first pipeline.

Heuristic Search Based on DAG Constraints
In our instruction scheduling framework, Stage 2 is the critical stage that determines the final performance. A heuristic search method is proposed to address this stage. We choose the A-Star search algorithm [51] as the framework. To search in heuristic order, the A-Star algorithm selects the next state that minimizes f (n) = g(n) + h(n), where n is the next state be selected, g(n) is the cost of the selected state set, and h(n) is a heuristic function that estimates the lowest cost from the current state set to the goal. In our work, g(n) is the execution clocks of all selected instructions, and h(n) is the number of remaining instructions.
However, only an A-Star search is insufficient with regard to the instruction scheduling with a scale of more than 30 instructions. More constraints to prune are necessary. Using a directed acyclic graph (DAG) to represent constraints is very suitable to the instruction scheduling. In general, a node in a DAG represents an instruction; an edge a → b in a DAG indicates that b has data dependence on a. With regard to the heuristic methods based on DAGs, a previous work [52] surveyed three DAG construction algorithms and twenty-six proposed heuristics methods and evaluated their performance exhaustively.
Different from existing methods, a node in a DAG represents a sequence of instructions in our method. Due to the properties of complex operations, instructions in the kernel are symmetrical to some extent. For example, computing of the real and imaginary parts is symmetrical. Moreover, if the operation is fused multiply-add (FMA), then the instructions can be divided into 4 symmetrical parts. For each symmetrical part, the corresponding instructions are the same instructions with different operands. We assign the corresponding instructions from different parts into the same DAG nodes. Then, we can set an order for these corresponding instructions distinguished by the operands, and all nodes share the same order. Setting such an order does not lead to a loss of the optimal solution because (1) swapping corresponding instructions in some DAG nodes  (58) cannot affect the structural hazard in the writeback stage since all corresponding instructions have the same clocks, and (2) data hazards cannot be reduced by swapping since the best choice is for all DAG nodes to share the same order. Figure 5 shows the main part of the instruction DAG of a 128-bit reduction. In the top left yellow panel, the order within the node is determined as "crr", "cri", "cir", and "cii", which are four different combinations of the real and imaginary parts of two complex numbers. The remaining DAG nodes all maintain the same internal order. In addition, in the middle orange panel, we can put eight instructions into this node due to the symmetry between "b, ar sign " and "a, br sign " rather than the complex symmetry. Through the optimization, the efficiency of the heuristic search is increased many times, and thus, we can obtain the optimal results. Table 2 lists the clocks of different methods in different computing parts. "128 Outer" represents the clocks of the instruction block in the outer loop, which is called O(N 2 ) times (while the inner loop (or the reduction of Gaussian elimination) is called O(N 3 ) times) when using 128-bit precision mode. For the three different methods, "Origin" represents results without any scheduling; "Manual" represents results by a redundant removal and a relatively careful manual scheduling; and "Optimal" represents the optimal results achieved by heuristic search described in this section. The number in brackets indicates the number of instructions, and the number outside of the brackets indicates the clocks.
We can see that the optimal results are much more efficient than the results of the previous two methods, showing our brilliant pipeline performance of approximately 90% efficiency. This result implies that for our multiple-precision design, we achieve 90% of the theoretical optimal floatingpoint performance.

Performance Results
This part demonstrates the performance results. Unless otherwise specified, we use 39,286 nodes (157,144 processes) of Sunway TaihuLight for experimentation, which is almost the entire capacity of the supercomputer (a total of 40,960 nodes). Figure 6 shows the execution time of both precision modes based on a log-transformed distribution. In the beginning, the execution time is nonlinear because the workload occupancy is insufficient. Then, the time is almost linear in the end due to the high scalability.

Time to Solution
To obtain a sample, it usually needs to calculate about 100 probabilities of the candidate samples using Markov  Fig. 7: CP U × T ime results of various GBS classical benchmarks. Red line indicates this work when using 256-bit precision mode. Blue, yellow, green, and pink lines indicate the works of [53], [54], [55], and [42], respectively.

CPU X Time (s) New GBS Algorithm by BJ Wu GBS with TD by N. Quesada Hafnian Computing by A. Bjorklund GBS with TD on Titan by B. Gupt Torontonian Computing on SW with 256 bits
Chain Monte Carlo (MCMC) sampling method [40]. To calculate one Torontonian probability for a 50-click photon detecting event, the execution time in our benchmark is 73,773s (about 20 hours) for 128-bit mode and 170,891s (about 2 days) for 256-bit mode, respectively. Figure 7 shows a comparison with other works that handle GBS. In the figure, we can see that except for our work (red line) and the work of [55] (green line), most works cannot process up to 50 photons. Compared with work represented by the green line, our work is hundreds of times faster. The time-to-solution results show the excellent performance of our method.

FLOPS
For FLOPS counting, we count all arithmetic operations in the Torontonian function. As a standard, solving a linear equation needs 2/3n 3 + O(n 2 ) FLOPS with a size n matrix [56]. However, we need to make adjustments to the determinant calculation of the Torontonian function. First, as mentioned above, the input matrix is guaranteed to be an Hermitian positive definite complex matrix. Thus, the cost of calculating a complex number should be considered, and half of the operations can be saved by Chomsky decomposition (or a modified Gaussian elimination). In addition, the matrix size n is usually small in a Torontonian function, so the term O(n 2 ) should not be omitted. Therefore, we count the operations of the linear solver accurately, as shown in Formula 8.
Bringing this formula into the Torontonian function and using some combinatorial math skills, we obtain the total FLOPS: The highest FLOPS did not occur at N = 50. This is because when N becomes larger, the LDM space requirement becomes larger. Thus, as mentioned in Subsection 4.2, a calculating matrix is scattered to more CPEs, which results in additional overhead (explained in Subsection 7.1.5). This issue is especially severe in 256-bit mode.

Scalability
For strong scalability, we choose the data scale N = 36 and use 4,096 to 131,072 processes. Figure 9 shows the strong scalability results. As the figure shows, both precision modes achieve an almost linear strong scalability.

Load Balance
This part reports the load balance performance with regard to the partition strategy described in Subsection 4.1. Table 3 shows the result in 128-bit precision mode when N = 40.
Regardless of the average, maximum, and minimum values, they are all 55.6G clocks. Thus, as long as the data scale is large enough, our partition strategy can achieve an almost completely balanced load. Figure 10 shows the percentage of time cost of each kernel when N = 50. As shown in Algorithm 1, GENAZ corresponds to Line 6, which is used to generate the matrix A Z based on the current mask. INV, SYNC, REGCOM, and ELIMI all correspond to Line 7, i.e., the function   Figure 10 we can see the time cost of SYNC in 256 bits is much higher than that for 128 bits since the matrix is scattered to more CPEs in 256 bits. This is the main factor that affects performance when N becomes larger.

G E N A Z IN V S Y N C R E G C O M E L IM I S Q R T IN V
ELIMI is related to the inner loop of the linear solver, which is composed of O(2 N N 3 ) times of arithmetic and thus is the bottleneck that takes up 70% ∼ 90% of the time.
Last, SQRTINV corresponds to Line 8, which computes the reciprocal square root.

Comparisons with Other Architectures
We use the FLOPS metric to compare the performance with other architectures. For the x86 architecture (Xeon E5-2680 v3) and the ARM architecture (Kunpeng 920-4826 ARM v8), we directly run a series of additions without data dependence (called kernel "ADD") based on several efficient libraries such as DD and GMP. For the GPU platform (Tesla C2050), we choose the CAMPARY library as our counterpart, which is one of the best multiple precision libraries in GPU. We directly use data supported by work [57] as the GPU results. The work [58] proposes a new multiple precision algorithm that is faster than CAMPARY when requiring greater than 106 bits of precision. However, in this work, all libraries are tested based on a single-precision version on the NVIDIA GTX architecture, which is much slower than the double-precision version.
Although the kernels used by each architecture are different, we guarantee that the Torontonian function used in our work is the most difficult one. Table 4 lists the FLOPS results among all architectures. Only the GPU architecture is comparable to our algorithm on the Sunway architecture. Since the Tesla C2050 is a relatively old version of the GPU, we also choose the Tesla V100 GPU (7.066 TFLOPS peak performance). For this purpose, we convert the performance results based on the ratio of the peak performance of two GPUs. The results of 106 bits and 212 bits are 204 GFLOPS and 13.5 GFLOPS, respectively. Although the 106-bit mode is much better now, the 212-bit mode, which must be applied when the data scale is very large, is over three times slower than our algorithm.
Consequently, our algorithm based on the Sunway architecture has an advantage compared with other architectures.

DISCUSSION
In addition to introducing a state-of-the-art benchmark based on GBS with threshold detection, our work provides some new insights from the computing perspective. Our partition strategy is applicable to other platforms. Our heuristic search for instruction scheduling is suitable for other applications with a hand-write assembly kernel. With minor modifications, our idea can also be applied to other platforms, especially for those platforms with a write-back stage structural hazard.
The multiple-precision fixed-point design is one of our major contributions. Based on a architecture-specific large integer instruction set, a set of fast operators is achieved to support the Torontonian function, which includes addition, subtraction, multiplication, division, and reciprocal square root. Our design can be easily applied to other scientific applications requiring high accuracy. In addition to the Sunway architecture, any current or future machine with such a specific instruction set can benefit. In our future work, we may increase and improve the operators to form an efficient and cross-platform fixed-point multiple-precision library based on the hardware instructions.
Last, our program has good portability by replacing our customized precision design with existing multipleprecision libraries under other platforms.

CONCLUSION
GBS with threshold detection is one of the most feasible ways to demonstrate quantum supremacy. Based on the Sunway TaihuLight supercomputer, which contains suitable architecture for this problem, we established a state-ofthe-art classical benchmark waiting for its future quantum computing counterpart.
To simultaneously achieve almost optimal performance and sufficient accuracy, a series of methods including a partition strategy with excellent load balancing, optimal instruction scheduling based on a heuristic search, multipleprecision fixed-point design based on an architecturespecific instruction set, and adaptive precision mode selection according to upper-and lower-bound estimates was proposed and applied.
As a result, sustained performance of 2.78 PFLOPS for 128-bit precision and 1.27 PFLOPS for 256-bit precision was achieved with a proper N . The largest run enabled us to obtain one Torontonian function of a 100 × 100 submatrix from 50-photon GBS within 20 hours with 128-bit precision and 2 days in 256-bit precision.