IEEE 754 double-precision floating-point number rapid printing method and system based on AVX-512 instruction set

Abstract

The invention belongs to the technical field of floating-point number printing, and discloses an IEEE 754 double-precision floating-point number rapid printing method based on an AVX-512 instruction set, and a floating-point number printing algorithm d2sci based on the AVX-512 instruction set is excellent in performance and accuracy. Information is kept, the printing result can be analyzed into the original floating-point number, and no information is lost. And the printing result is close to the binary representation of the input floating-point number, and the maximum error is within the error range of the IEEE 754 double-precision floating-point number. And the performance is excellent: compared with the existing printing algorithm, the performance is averagely improved by about 4-20 times, and the time required for processing a large amount of floating-point number printing in engineering application is greatly reduced.

Description

Technical Field

[0001] The present invention belongs to the technical field of floating point number printing, and in particular relates to an IEEE 754 double-precision floating point number fast printing method and system based on the AVX-512 instruction set. Background Art

[0002] Floating-point numbers are a core component of modern computers, and most computers use the IEEE 754 standard for internal storage. In scenarios such as scientific computing, data analysis, and industrial software, it's necessary to convert floating-point numbers to text in decimal scientific notation (known as floating-point printing) for easier human readability or data storage. This requirement exists in common industrial software applications. For example, in computational fluid dynamics (CFD), grid coordinates are typically stored using double-precision floating-point numbers, with either ASCII or binary formats available. When ASCII is used to store grid coordinate data, the floating-point numbers must be converted to ASCII (a string in decimal scientific notation). This conversion is computationally intensive and is sometimes referred to as floating-point rendering or visualization. Efficiently and accurately printing floating-point numbers has long been a key topic in computer science. Striking a balance between information losslessness (the converted decimal scientific notation is equivalent to the original floating-point number), correct rounding (the printed result is as close as possible to the binary representation of the input floating-point number within the valid length), and high performance is a hot topic of research.

[0003] In 2010, Loitsch proposed the Grisu algorithm, a classic floating-point number to string conversion algorithm that achieves efficient conversion while ensuring accuracy. In the article, the author introduced three versions of the algorithm: (1) Grisu: By caching the integer power of 10 in a local array, the mantissa and exponent in decimal scientific notation are calculated and the output is fixed length, but the result may not be the shortest output form. (2) Grisu2: By adding extra bits to achieve the shortest output, about 99.8% of floating-point numbers can get the shortest output, while the remaining approximately 0.2% will degenerate to the Grisu3 algorithm. (3) Grisu3: Based on Grisu2, Grisu3 uses higher-precision operations, and about 99.5% of the inputs produce the shortest output. For the remaining 0.5%, other perfect algorithms are used to produce the shortest output. The Grisu algorithm has not yet fully utilized the vector instruction set of modern processors, so there is still potential for performance improvement. In 2016, Andrysco et al. proposed three versions of the Errol algorithm based on Grisu: (1) Errol1: Faster and more optimal, with the fastest speed and producing optimal solutions for most data; (2) Errol2: Almost optimal, producing optimal solutions for almost all floating-point numbers, but still not for some data; (3) Errol3: Always optimal, producing optimal solutions for all floating-point numbers. Although Errol3 is slightly slower than Grisu3 in the worst case, it is still about 5.2 times faster than the previous complete algorithm.

[0004] In 2019, Leonid Yuriev developed the erthink algorithm based on Grisu and open-sourced it at https: / / github.com / erthink / erthink.git. This algorithm ensures that the output can always be converted to the original value, and for more than 99.963% of IEEE 754 double-precision values, the converted string representation is the shortest, meaning that for less than 0.037% of values, one digit is added. Furthermore, for less than 0.06% of double-precision values, the last digit differs from the ideal nearest digit by ±1. Compared to Ryu's algorithm, the implementation code size is significantly smaller, and calculating each digit takes approximately 16 to 17 clock cycles. In 2020, Junekey Jeon proposed the Grisu-Exact algorithm based on Grisu2, which improves information fidelity, shortest representation, and correct rounding. For numbers with fewer digits (such as 2 or 6 digits), its average performance is faster than Ryu's algorithm.

[0005] In 2019, Raffaello Giulietti proposed the Schubfach algorithm, an efficient method for converting double-precision floating-point numbers into the shortest, correctly rounded decimal representation. This algorithm uses fixed-precision integer operations to replace large integer calculations, pre-calculate multiplication factors of powers of 10, and correctly round off boundaries to improve algorithm efficiency. In 2021, Junekey Jeon proposed the Dragonbox algorithm based on the Schubfach algorithm, and drew inspiration from Grisu, Grisu2, and Grisu-Exact. The algorithm tries to avoid the expensive operation of multiplying 128-bit integers by 64-bit integers in the Schubfach algorithm, at the cost of more branches and constant division operations. Test results for single-precision and double-precision floating-point numbers show that the algorithm outperforms the Ryu, Grisu-Exact, and Schubfach algorithms.

[0006] Through the above analysis, the problems and defects of the existing technology are as follows:

[0007] (1) Information loss: After conversion to decimal scientific notation, a small number of values ​​cannot be represented or the reverse conversion results in a value that is not equal to the original floating-point value.

[0008] (2) Rounding error: For the last few digits of the mantissa, the conversion may result in a loss of precision, resulting in only approximate equality in some cases. In most cases, this will not affect the calculation results. However, for high-precision scenarios such as scientific computing and numerical simulation, it will cause additional errors or incorrect results in the final problem solution.

[0009] (3) Long calculation time: Existing algorithms often need to make a trade-off between calculation speed and result accuracy. In order to obtain completely accurate results, a long calculation time is usually required. Summary of the Invention

[0010] In view of the problems existing in the prior art, the present invention provides an IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set.

[0011] The present invention is implemented as follows: a fast printing method for IEEE 754 double-precision floating-point numbers based on the AVX-512 instruction set, which mainly calculates the values ​​of a and b in the decimal scientific notation aEb corresponding to the IEEE 754 double-precision floating-point number, that is,

[0012] x=a×10 6 (1)

[0013] The main steps include:

[0014] Step 1, calculate b;

[0015] Calculate the decimal scientific notation aEb in b, that is, calculate a floating point number x

[0016] Step 2, calculate a;

[0017] Calculate a in decimal scientific notation aEb, where a is a decimal number with up to 17 decimal digits and the range of a is [1,10). Multiply a by 10 16 And round it off to get the value you want to print, which is in the range of [10 16 ,10 17 );

[0018] Step 3, output a and b;

[0019] According to the calculated a and b, output their corresponding ASCII codes.

[0020] Furthermore, the calculation b is:

[0021] This is converted into the following two questions.

[0022] (1) Solution

[0023] Question 1: Given a non-zero positive floating point number x, find

[0024] Solution: Assume Then there is For convenience, remember Then we have:

[0025] 2 z ≤x<2 z+1 (2)

[0026] Taking the common logarithm of both sides of the formula, we have:

[0027] zlg2≤lgx=e 10 <(z+1)lg2 (3)

[0028] Rounding the formula down, we have:

[0029]

[0030] Therefore There are only two possible values, namely or Therefore:

[0031]

[0032] (2) Solve for z;

[0033] In addition to zero, infinity, and NaN, in order to calculate the formula, it is necessary to find the exponent z of the floating-point number; at this time, it can be divided into two cases for discussion: normalized numbers and denormalized numbers.

[0034] ① Normalized number

[0035] When the exponent is not all 0 or all 1, simply subtract the bias value 1023 from the exponent exp, and you have:

[0036]

[0037] ②Non-normalized number

[0038] Convert the denormalized number into Problem 2 and calculate the exponent z based on the number of significant digits of the mantissa frac.

[0039] Question 2: For a denormalized number x with all exponents 0 and arbitrary mantissa, find

[0040] Solution: According to the IEEE 754 standard, the formula for calculating the true value x of a double-precision floating-point number is as follows:

[0041] x=2 1-1023 ×(frac×2 -52 )=frac×2 -1074 (7)

[0042] in, and Taking the base 2 logarithm of both sides of the formula, we have:

[0043]

[0044] The expression value x of a denormalized double-precision positive floating-point number with all exponents 0 and arbitrary mantissa b , its high 12 bits (1-bit sign bit and 11-bit exponent part exp) are all 0, and the low 52 bits (mantissa part frac) are uncertain; if clz is used to represent x b The number of leading 0s (excluding the sign bit) is:

[0045]

[0046] Substituting the formula into the equation, we have:

[0047]

[0048] Furthermore, the calculation a is:

[0049] According to the formula, we have:

[0050] a=x×10 -b

[0051] According to the definition of the mantissa a, its value range is a decimal in [1,10); for ease of processing, it needs to be converted into an integer when printing, that is, the actual number to be printed Multiply by 10 16 The results are:

[0052]

[0053] Calculation We need to calculate the multiplication of two numbers. To improve efficiency, we introduce a custom floating-point representation. A floating-point number consists of the following two parts:

[0054] ① Mantissa part: uses 64-bit unsigned integer F;

[0055] ②Exponent part: uses 32-bit signed integer E;

[0056] Among them, the exponent part and the mantissa part are both in binary form;

[0057] At this time, the custom expression value of the floating point number Defined as:

[0058]

[0059] The true value x of the floating-point number can be calculated using the following formula:

[0060] x=F×2 E (12)

[0061] Two multiplication formulas using custom floating-point numbers are defined as follows;

[0062] Definition: Two floating-point numbers with true values ​​x1 and x2, whose custom expression values ​​are and Then floating point multiplication is defined as:

[0063]

[0064] Since F1 and F2 are 64-bit unsigned integers, the result of F1×F2 requires 128 bits to be fully recorded. In order to uniformly adopt the new method to represent the multiplication results, the algorithm only retains the upper 64 bits. The final approximate result is:

[0065]

[0066] Right now:

[0067]

[0068] Compute the printed value by converting x and 10-b+16 to a custom floating point representation This can be converted into the following question;

[0069] (1) Calculate the custom form of the floating-point true value x

[0070] Problem 3: Convert an IEEE 754 standard double-precision floating-point number x to a custom floating-point number Format;

[0071] Solution: Note x b To express the value of the double-precision floating-point number x without the sign bit, there are two cases:

[0072] ① All order codes are 0, non-normalized number

[0073] Similarly, use clz to represent x b The number of leading 0s; at this time, for the exponent E, according to the formula:

[0074] E=-1011-clz

[0075] The last digit F is:

[0076] F=x b <<clz

[0077] ② The order code is not all 0 or all 1, normalized number

[0078] According to the formula, the exponent E is:

[0079] E=(x b >>52)-1023

[0080] According to the IEEE 754 standard, the mantissa F is:

[0081] F=(x b <<11)|2 63

[0082] In summary, custom floating point numbers The calculation formula is:

[0083]

[0084] Among them, "and" are left shift and right shift operations respectively, which are bitwise OR operations;

[0085] (2) Custom form for calculating integer powers of 10

[0086] Question 4: For any integer power of 10 Calculate its custom form

[0087] untie:

[0088] ①Calculate the order code part; but:

[0089] a) when e 10 =0, then e2=0;

[0090] b) When e 10 ≠0, then e2=e 10 ×log210; because And log210 is an irrational number, so

[0091] Then there is And satisfy:

[0092]

[0093] Rounding down the above formula, we have:

[0094]

[0095] because are all integers, then:

[0096]

[0097] For case (1), the formula still holds; in summary:

[0098]

[0099] 2. Calculate the mantissa. To improve accuracy, fix the highest bit of the 64-bit mantissa F to 1. The corresponding exponent and mantissa are:

[0100]

[0101] In order to obtain a more accurate mantissa F, first calculate the 65-bit result F t , so we have:

[0102]

[0103] a) when e 10 When <0, there are:

[0104]

[0105] b) When e 10 ≥20, Then we have:

[0106]

[0107] c) When 0≤e 10 ≤19 o'clock, Then we have:

[0108]

[0109] Finally, according to F t The last digit determines F. If it is 1, it is rounded, and if it is 0, it is not rounded. Then we have:

[0110]

[0111] Among them, & is the bitwise AND operation;

[0112] (3) Calculate and print the value

[0113] If the custom binary representation of the input floating point number x is Denoted as {Fx,Ex}, the custom binary representation of the integer power of 10 Denoted as {Fn,En}, then substituting the formula and into, we have:

[0114]

[0115] E a =(E x -63)+(E n -63)+64=E x +E n -62<0

[0116] After simplification, round the result to get the printed value integer have:

[0117]

[0118] Furthermore, the outputs a and b are:

[0119] The print value of the mantissa a is, and There are 17 decimal digits in total; they are divided into two parts: the high 9-bit part high9 and the low 8-bit part low8; the high 9-bit part high9 can correspond to three groups of consecutive 3-digit decimal digits, and the print output ASCII codes of these 3 groups of numbers can be directly obtained by looking up the table and output them to the print buffer; while looking up the table, the ASCII code of the low8 part is calculated in real time and the result is output to the print buffer; for the exponent b, b∈[-324, 308], which can be directly obtained by looking up the table.

[0120] Another object of the present invention is to provide an IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set, comprising:

[0121] The calculation b module is used to calculate the b in the decimal scientific notation aEb, that is, to calculate a floating point number x;

[0122] The calculation module a is used to calculate a in decimal scientific notation aEb. a is a decimal number with up to 17 decimal digits. The range of a is [1,10). Multiply a by 10 16 And round it off to get the value you want to print, which is in the range of [10 16 ,10 17 );

[0123] The output a and b modules are used to output the corresponding ASCII codes based on the calculated a and b.

[0124] Another object of the present invention is to provide a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set.

[0125] Another object of the present invention is to provide a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to execute the steps of the IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set.

[0126] Another object of the present invention is to provide an information data processing terminal, which is used to implement the IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set.

[0127] In combination with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present invention are as follows:

[0128] First, the floating-point number printing algorithm d2sci based on the AVX-512 instruction set proposed in the present invention performs well in both performance and accuracy.

[0129] Information preservation: The print result of this algorithm is a fixed 17-digit decimal digit, which can be parsed into the original floating-point number without information loss.

[0130] Correct rounding: This algorithm converts integer powers of 10 to custom floating-point numbers and calculates printed values. It has been proven that the maximum error range introduced is Within the IEEE 754 double-precision floating point error range.

[0131] Excellent Performance: Leveraging the vector processing capabilities of modern processors, this approach significantly reduces branch instructions and mitigates the penalties for branch prediction failures by optimizing instruction selection and loop structures. Compared to existing printing algorithms, this approach delivers an average performance improvement of approximately 4 to 20 times, significantly reducing the time required to process large floating-point number prints in engineering applications.

[0132] Second, in common industrial software applications, there is a need to convert IEEE 754 double-precision floating-point numbers into character strings and store them in files or output them to the screen. For example, in the field of computational fluid dynamics (CFD), double-precision floating-point numbers are usually used in CAE software to represent grid coordinate points. When a high-precision grid with 1 billion degrees of freedom is involved in the solution, all parameters and grid data need to be frequently output in ASCII code to the screen and intermediate result files (the size of the intermediate result file output each time is about 30G) during the solution process. At this time, it is necessary to perform a floating-point number conversion to ASCII code. This conversion is a computationally intensive task. Improving the efficiency of floating-point number conversion helps to shorten the parallel solution time overall. BRIEF DESCRIPTION OF THE DRAWINGS

[0133] Figure 1 This is a flow chart of a method for quickly printing IEEE 754 double-precision floating-point numbers based on the AVX-512 instruction set provided by an embodiment of the present invention.

[0134] Figure 2 This is a structural block diagram of an IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set provided by an embodiment of the present invention.

[0135] Figure 3 This is a diagram of the IEEE 754 standard 64-bit floating-point number format provided by an embodiment of the present invention.

[0136] Figure 4 This is a flowchart of the algorithm processing provided by an embodiment of the present invention.

[0137] Figure 5 This is a graph of the average execution time of printing random single floating-point numbers for each algorithm provided by the embodiments of the present invention.

[0138] Figure 6 This is a graph of the average execution time of each algorithm provided in the embodiments of the present invention for printing a random single floating-point number with 1 to 17 significant digits. DETAILED DESCRIPTION

[0139] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0140] An embodiment of the present invention provides an IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set, characterized in that the method is executed by a processor that supports the AVX-512 instruction set and includes the following steps:

[0141] S1. Calculate the scientific notation exponent b of floating-point numbers in batches using AVX-512 parallel instructions.

[0142] S2. Express the floating point number x as a×10 b , multiply a by 10 to the power of 16 and round it off with the help of SIMD instructions to get the printed integer value;

[0143] S3. Quickly map a and b to a print output format through table lookup and vector ASCII code conversion instructions.

[0144] The exponent b is extracted by the hardware instruction auxiliary module from the IEEE 754 floating-point number exponent and inferred based on the parallel lg estimation instruction. Finally, a SIMD vector compare and select instruction is used to complete the exponent range judgment and determine the integer b.

[0145] An exponential z extraction acceleration circuit is provided in the processor for:

[0146] a) Under normalized numbers, directly decode the floating point exponent and subtract the bias value 1023;

[0147] b) For denormalized numbers, the hardware clz (Count Leading Zeros) instruction is used to quickly count the leading zeros in the mantissa and estimate the exponent z based on this count, shortening the logic branches and judgment paths.

[0148] The calculation of the mantissa a is completed by a custom floating-point multiplication unit implemented in FPGA or ASIC, which supports truncation of the high 64 bits of the 64-bit integer multiplication result. In conjunction with the exponent register update and mantissa rounding mechanism, the high-precision integer a is generated.

[0149] The custom floating-point multiplication unit includes:

[0150] a) a mantissa multiplication logic array for calculating the upper 64 bits of the product of two 64-bit unsigned integers;

[0151] b) exponent accumulator logic for adding two 32-bit signed integers;

[0152] c) Mantissa rounding controller, used to detect the low bit of the product and determine whether to round it up or down.

[0153] The output module includes parallel lookup table logic and ASCII code generation circuit. The lookup table logic reads three groups of ASCII codes corresponding to the upper 9 digits from the three-segment compressed lookup table ROM and writes them into the character buffer at high speed; the lower 8 bits are calculated in real time by the pipeline parallel generation circuit.

[0154] Integer powers of 10 are pre-calculated and stored in a high-speed read-only buffer. The processor or chip quickly obtains the corresponding custom floating-point format {F, E} through indexing, reducing the reliance on lg and power function calculations at runtime.

[0155] The IEEE 754 double-precision floating-point fast printing system based on the AVX-512 instruction set includes the following modules, all implemented through software and hardware collaboration:

[0156] a) Exponent calculation module, which extracts floating-point exponent b based on hardware decoding and AVX-512 SIMD estimation instructions;

[0157] b) A custom multiplication module that calculates the product of floating-point numbers and powers of 10 and rounds off the mantissa based on a custom accelerator;

[0158] c) Table lookup output module, integrating three-segment high and low bit table lookup structure and ASCII code converter, supporting pipeline parallel output.

[0159] The index calculation module includes:

[0160] a) Hardware floating-point decoder for parsing the IEEE 754 format to extract the exp field;

[0161] b) a leading zero detector to handle exponential calculation of denormalized numbers;

[0162] c) A table lookup multiplication factor module is used to quickly map z to the decimal exponent b value corresponding to lg2×z.

[0163] The table lookup output module includes:

[0164] a) Compression coding ROM lookup table, used for three-segment lookup table mapping of the high-order ASCII code output of the mantissa a;

[0165] b) a pipeline parallel decimal encoder for calculating the low-order mantissa into ASCII code;

[0166] c) an index b lookup table mapping unit, used to convert the integer b into a corresponding exponential character output format.

[0167] like Figure 1 As shown, an embodiment of the present invention provides an IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set, comprising the following steps:

[0168] S101, calculate b;

[0169] Calculate the decimal scientific notation aEb in b, that is, calculate a floating point number x

[0170] S102, calculate a;

[0171] Calculate a in decimal scientific notation aEb, where a is a decimal number with up to 17 decimal digits and the range of a is [1,10). Multiply a by 10 16 And round it off to get the value you want to print, which is in the range of [10 16 ,10 17 );

[0172] S103, output a and b;

[0173] According to the calculated a and b, output their corresponding ASCII codes.

[0174] The calculation b provided in the embodiment of the present invention is:

[0175] (1) Solution

[0176] Question 1: Given a non-zero positive floating point number x, find

[0177] Solution: Assume Then there is For convenience, remember Then we have:

[0178] 2 z ≤x<2 z+1 (twenty two)

[0179] Taking the common logarithm of both sides of the formula, we have:

[0180] zlg2≤lgx=e 10 <(z+1)lg2 (23)

[0181] Rounding the formula down, we have:

[0182]

[0183] Therefore There are only two possible values, namely or Therefore:

[0184]

[0185] (2) Solve for z;

[0186] In addition to zero, infinity, and NaN, in order to calculate the formula, we need to find the exponent z of the floating-point number. In this case, we can discuss it in two cases: normalized number and denormalized number.

[0187] 1 Normalized number

[0188] When the exponent is not all 0 or all 1, simply subtract the bias value 1023 from the exponent exp, and you have:

[0189] z=exp-1023 (26)

[0190] 2Denormalized number

[0191] Convert the denormalized number to problem 2 and calculate the exponent z based on the number of significant digits of the mantissa frac;

[0192] Question 2: For a denormalized number x with all exponents 0 and arbitrary mantissa, find

[0193] Solution: According to the IEEE 754 standard, the formula for calculating the true value x of a double-precision floating-point number is as follows:

[0194] x=2 1-1023 ×(frac×2-52)=frac×2 -1074 (27)

[0195] in, And frac∈[1,2 52 -1]; taking the base 2 logarithm of both sides of the formula, we have:

[0196]

[0197] The expression value x of a denormalized double-precision positive floating-point number with all exponents 0 and arbitrary mantissa b , its high 12 bits (1-bit sign bit and 11-bit exponent part exp) are all 0, and the low 52 bits (mantissa part frac) are uncertain; if clz is used to represent x b The number of leading 0s (excluding the sign bit) is:

[0198]

[0199] Substituting the formula into the equation, we have:

[0200]

[0201] The calculation a provided in the embodiment of the present invention is:

[0202] According to the formula, we have:

[0203] a=x×10 -b

[0204] According to the definition of the mantissa a, its value range is a decimal in [1,10); for ease of processing, it needs to be converted into an integer when printing, that is, the actual number to be printed Multiply by 10 16 The results are:

[0205]

[0206] Calculation We need to calculate the multiplication of two numbers. To improve efficiency, we introduce a custom floating-point representation. A floating-point number consists of the following two parts:

[0207] ① Mantissa part: uses 64-bit unsigned integer F;

[0208] ②Exponent part: uses 32-bit signed integer E;

[0209] Among them, the exponent part and the mantissa part are both in binary form;

[0210] At this time, the custom expression value of the floating point number Defined as:

[0211]

[0212] The true value x of the floating-point number can be calculated using the following formula:

[0213] x=F×2 E (32)

[0214] Two multiplication formulas using custom floating-point numbers are defined as follows;

[0215] Definition: Two floating-point numbers with true values ​​x1 and x2, whose custom expression values ​​are and Then floating point multiplication is defined as:

[0216]

[0217] Since F1 and F2 are 64-bit unsigned integers, the result of F1×F2 requires 128 bits to be fully recorded. In order to uniformly adopt the new method to represent the multiplication results, the algorithm only retains the upper 64 bits. The final approximate result is:

[0218]

[0219] Right now:

[0220]

[0221] Compute the printed value by converting x and 10-b+16 to a custom floating point representation This can be converted into the following question;

[0222] (1) Calculate the custom form of the floating-point true value x

[0223] Problem 3: Convert an IEEE 754 standard double-precision floating-point number x to a custom floating-point number Format;

[0224] Solution: Note x b To express the value of the double-precision floating-point number x without the sign bit, there are two cases:

[0225] The first-order code is all 0, which is a non-normalized number.

[0226] Similarly, use clz to represent x b The number of leading 0s; at this time, for the exponent E, according to the formula:

[0227] E=-1011-clz

[0228] The last digit F is:

[0229] F=x b <<clz

[0230] The 2nd order code is not all 0 or all 1, normalized number

[0231] According to the formula, the exponent E is:

[0232] E=(x b >>52)-1023

[0233] According to the IEEE 754 standard, the mantissa F is:

[0234] F=(x b <<11)|2 63

[0235] In summary, custom floating point numbers The calculation formula is:

[0236]

[0237] Among them, "and" are left shift and right shift operations respectively, which are bitwise OR operations;

[0238] (2) Custom form for calculating integer powers of 10

[0239] Question 4: For any integer power of 10 Calculate its custom form untie:

[0240] ①Calculate the order code part; but:

[0241] a) when e 10 =0, then e2=0;

[0242] b) When e 10 ≠0, then e2=e 10 ×log210; because And log210 is an irrational number, so

[0243] Then there is And satisfy:

[0244]

[0245] Rounding down the above formula, we have:

[0246]

[0247] because are all integers, then:

[0248]

[0249] For case (1), the formula still holds; in summary:

[0250]

[0251] ② Calculate the mantissa. To improve accuracy, the highest bit of the 64-bit mantissa F is fixed to 1. The corresponding exponent and mantissa are:

[0252]

[0253] In order to obtain a more accurate mantissa F, first calculate the 65-bit result F t , so we have:

[0254]

[0255] a) when e 10 When <0, there are:

[0256]

[0257] b) When e 10 ≥20, Then we have:

[0258]

[0259] c) When 0≤e 10 ≤19 o'clock, Then we have:

[0260]

[0261] Finally, according to F t The last digit determines F. If it is 1, it is rounded, and if it is 0, it is not rounded. Then we have:

[0262]

[0263] Among them, & is the bitwise AND operation;

[0264] (3) Calculate and print the value

[0265] If the custom binary representation of the input floating point number x is Denoted as {Fx,Ex}, the custom binary representation of the integer power of 10 Denoted as {Fn,En}, then substituting the formula and into, we have:

[0266]

[0267] E a =(E x -63)+(E n -63)+64=E x +E n -62<0

[0268] After simplification, round the result to get the printed value integer have:

[0269]

[0270] Outputs a and b provided by this embodiment of the present invention are:

[0271] The print value of the mantissa a is, and There are 17 decimal digits in total; they are divided into two parts: the high 9-bit part high9 and the low 8-bit part low8; the high 9-bit part high9 can correspond to three groups of consecutive 3-digit decimal digits, and the print output ASCII codes of these 3 groups of numbers can be directly obtained by looking up the table and output them to the print buffer; while looking up the table, the ASCII code of the low8 part is calculated in real time and the result is output to the print buffer; for the exponent b, b∈[-324, 308], which can be directly obtained by looking up the table.

[0272] like Figure 2 As shown, an embodiment of the present invention provides an IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set, including:

[0273] The calculation b module is used to calculate the b in the decimal scientific notation aEb, that is, to calculate a floating point number x;

[0274] The calculation module a is used to calculate a in decimal scientific notation aEb. a is a decimal number with up to 17 decimal digits. The range of a is [1,10). Multiply a by 10 16 And round it off to get the value you want to print, which is in the range of [10 16 ,10 17 );

[0275] The output a and b modules are used to output the corresponding ASCII codes based on the calculated a and b.

[0276] Another object of the present invention is to provide a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set.

[0277] Another object of the present invention is to provide a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to execute the steps of the IEEE 754 double-precision floating-point number fast printing method based on the AVX-512 instruction set.

[0278] Another object of the present invention is to provide an information data processing terminal, which is used to implement the IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set.

[0279] The present invention is specifically implemented:

[0280] 1.1 Algorithm Design Ideas

[0281] For Figure 3 The printing problem of the IEEE 754 double-precision floating-point number shown can be converted into the problem of efficiently calculating the values ​​of a and b in decimal scientific notation aEb, that is,

[0282] x=a×10 b (42)

[0283] After conversion, the sign bit in decimal scientific notation is the same as the sign bit in binary. Since negative numbers only require the sign bit to be inverted, the following discussion will focus on non-zero positive floating-point numbers. Special cases such as positive and negative infinity, NaN (Not a Number), and zero require identification and handling during algorithm implementation.

[0284] First, the exponent and mantissa of the input floating-point number are extracted through bitwise operations. The mantissa is converted to a normalized binary form to obtain the significand frac; the exponent is converted to a decimal integer by the algorithm and the corresponding offset adjustment is performed to exp according to the provisions of the IEEE 754 standard. At this time, the mantissa a and exponent b in decimal scientific notation will be determined by the exponent exp and mantissa frac of the input value x. The process of this algorithm (d2sci) is as follows: Figure 4 As shown, it mainly includes the following three steps.

[0285] (1) Calculate b

[0286] Calculate the decimal scientific notation aEb in b, that is, calculate a floating point number x

[0287] (2) Calculate a

[0288] Calculate a in decimal scientific notation aEb, where a is a decimal number with up to 17 decimal digits and the range of a is [1,10). Multiply a by 10 16 And round it off to get the value you want to print, which is in the range of [10 16 ,10 17 ).

[0289] (3) Output a and b

[0290] According to the calculated a and b, output their corresponding ASCII codes.

[0291] 1.2 Calculation of b

[0292] The calculation of b in aEb can be abstracted as the solution to the following problem.

[0293] 1.2.1 Solution

[0294] Question 1: Given a non-zero positive floating point number x, find

[0295] Solution: Assume Then there is For convenience, remember Then we have:

[0296] 2 z ≤x<2 z+1 (43)

[0297] Taking the common logarithm of both sides of the formula, we have:

[0298] zlg2≤lgx=e 10 <(z+1)lg2 (44)

[0299] Rounding the formula down, we have:

[0300]

[0301] Therefore There are only two possible values, namely or Therefore:

[0302]

[0303] 1.2.2 Solving for z

[0304] In addition to zero, infinity, and NaN, in order to calculate the formula, we need to find the exponent z of the floating-point number. In this case, we can discuss two cases: normalized numbers and denormalized numbers.

[0305] (1) Normalized number

[0306] When the exponent is not all 0 or all 1, simply subtract the bias value 1023 from the exponent exp, and you have:

[0307] z=exp-1023 (47)

[0308] (2) Non-normalized numbers

[0309] For non-normalized numbers, the problem can be transformed into Problem 2, and the exponent z is calculated based on the number of significant digits of the mantissa frac.

[0310] Question 2: For a denormalized number x with all exponents 0 and arbitrary mantissa, find

[0311] Solution: According to the IEEE 754 standard, the formula for calculating the true value x of a double-precision floating-point number is as follows:

[0312] x=2 1-1023 ×(frac×2 -52 )=frac×2 -1074 (48)

[0313] in, And frac∈[1,2 52 -1]. Taking the base 2 logarithm of both sides of the formula, we have:

[0314]

[0315] The expression value x of a denormalized double-precision positive floating-point number with all exponents 0 and arbitrary mantissa b , its high 12 bits (1-bit sign and 11-bit exponent part exp) are all 0, and the low 52 bits (frac part) are uncertain. bThe number of leading 0s (excluding the sign bit) is:

[0316]

[0317] Substituting the formula into the equation, we have:

[0318]

[0319] 1.3 Calculate a

[0320] According to the formula, we have:

[0321] a=x×10 -b

[0322] According to the definition of the mantissa a, its value range is a decimal in [1,10). In order to facilitate processing, it needs to be converted into an integer when printing, that is, the actual number to be printed is Multiply by 10 16 The results are:

[0323]

[0324] Calculation We need to calculate the multiplication of two numbers. To improve efficiency, we introduce a custom floating point representation. By converting x and 10-b+16 into custom floating point representations, we calculate the printed value. It can be transformed into the following problem.

[0325] 1.3.1 Custom floating point representation

[0326] To ensure calculation accuracy, this algorithm uses a custom floating-point number representation method similar to the Grisu algorithm. In this method, a floating-point number consists of the following two parts:

[0327] Mantissa part: uses 64-bit unsigned integer F.

[0328] Exponent part: uses a 32-bit signed integer E.

[0329] The exponent part and the mantissa part are both in binary form.

[0330] At this time, the custom expression value of the floating point number Defined as:

[0331]

[0332] The true value x of the floating-point number can be calculated using the following formula:

[0333] x=F×2 E (53)

[0334] The multiplication formulas for two custom floating-point numbers are defined as follows.

[0335] Definition: Two floating-point numbers with true values ​​x1 and x2, whose custom expression values ​​are and Then floating point multiplication is defined as:

[0336]

[0337] Since F1 and F2 are 64-bit unsigned integers, the result of F1×F2 requires 128 bits to be fully recorded. In order to uniformly adopt the new method to represent the multiplication result, the algorithm only retains the upper 64 bits, and the final approximate result is:

[0338]

[0339] Right now:

[0340]

[0341] 1.3.2 Calculating the Custom Form of the Floating-Point True Value x

[0342] Problem 3: Convert an IEEE 754 standard double-precision floating-point number x to a custom floating-point number Format.

[0343] Solution: Note x b To express the value of the double-precision floating-point number x without the sign bit, there are two cases:

[0344] (1) All exponents are 0, non-normalized number

[0345] Similarly, use clz to represent x b The number of leading 0s. At this time, for the exponent E, according to the formula:

[0346] E=-1011-clz

[0347] The last digit F is:

[0348] F=x b <<clz

[0349] (2) The exponent is not all 0 or all 1, the normalized number

[0350] According to the formula, the exponent E is:

[0351] E=(x b >>52)-1023

[0352] According to the IEEE 754 standard, the mantissa F is:

[0353] F=(x b <<11)|263

[0354] In summary, custom floating point numbers The calculation formula is:

[0355]

[0356] Among them, << and >> are left shift and right shift operations respectively, which are bitwise OR operations.

[0357] 1.3.3 Custom forms for calculating integer powers of 10

[0358] Question 4: For any integer power of 10 Calculate its custom form

[0359] untie:

[0360] 1. Calculate the exponent part. but:

[0361] (1) When e 10 =0, then e2=0.

[0362] (2) When e 10 ≠0, then e2=e 10 ×log210. Because And log210 is an irrational number, so Then there is And satisfy:

[0363]

[0364] Rounding down the above formula, we have:

[0365]

[0366] because are all integers, then:

[0367]

[0368] For case (1), the formula still holds. In summary:

[0369]

[0370] Calculate the mantissa. To improve accuracy, the highest bit of the 64-bit mantissa F is fixed to 1, and the corresponding exponent and mantissa are:

[0371]

[0372] In order to obtain a more accurate mantissa F, first calculate the 65-bit result F t, so we have:

[0373]

[0374] When e 10 When <0, there are:

[0375]

[0376] When e 10 ≥20, Then we have:

[0377]

[0378] When 0≤e 10 ≤19 o'clock, Then we have:

[0379]

[0380] Finally, according to F t The last digit determines F. If it is 1, it is rounded, and if it is 0, it is not rounded. Then we have:

[0381]

[0382] Among them, & is a bitwise AND operation.

[0383] 1.3.4 Calculating Print Values

[0384] If the custom binary representation of the input floating point number x is Denoted as {Fx,Ex}, the custom binary representation of the integer power of 10 Denoted as {Fn,En}, then substituting the formula and into, we have:

[0385]

[0386] E a =(E x -63)+(E n -63)+64=E x +E n -62<0

[0387] After simplification, round the result to get the printed value integer have:

[0388]

[0389] 1.4 Output a and b

[0390] The print value of the mantissa a is, and There are 17 decimal digits in total. These are split into two parts: the upper 9 bits (high9) and the lower 8 bits (low8). The upper 9 bits (high9) correspond to three consecutive groups of three decimal digits. The ASCII codes for these three groups can be directly obtained by table lookup and output to the print buffer. Simultaneously with the table lookup, the ASCII codes for the lower 8 bits are calculated in real time and output to the print buffer. For the exponent b, b∈[-324, 308], which can be directly obtained by table lookup.

[0391] 2.1 Using fast calculation method

[0392] In order to improve the operation speed, when z∈[-1074, 1023], the formula It can be quickly calculated as:

[0393]

[0394] Similarly, when e 10 ∈[-292, 340], in the formula It can be quickly calculated as:

[0395]

[0396] To calculate the number of leading zeros clz in the formula, the Intel Intrinsic function _lzcnt_u64 of the x86-64 architecture can be used for efficient operation.

[0397] 2.2 Pre-storage Technology

[0398] 2.2.1 Pre-storage floating point numbers

[0399] In the formula, z∈[-1074, 1023], then To accelerate x and To compare, a local double array is needed to store 632 double values ​​from 1e-323 to 1e308, and the operations in Section 4.1.1 can be converted to a lookup table.

[0400] 2.2.2 Precompute custom floating-point numbers that are integer powers of 10

[0401] According to the IEEE 754 standard, the smallest absolute positive value that can be represented by a double-precision floating-point number is approximately 5×10 -324 , the maximum absolute positive value is approximately 1.7976931348623158×10 308 , then there is an index And b∈[-324,308]. To calculate the formula, 10-b+16 needs to be converted to a custom floating-point number, so the integer power range to be used is [-292,340]. Therefore, in programming implementation, an array can be used to cache the values ​​of 10 to the power of -292 to 340, a total of 633 unsigned 64-bit integers, and the operations in Section 4.2.3 can be converted to a lookup table. The exponent in the custom floating-point number can be calculated through multiplication and shift operations according to Formula (44) and does not require pre-storage.

[0402] 2.2.3 Pre-calculating ASCII codes for printout

[0403] In order to improve the printing speed, three parts are pre-stored. (1) ASCII codes from 00 to 99 are pre-calculated and stored in a 16-bit integer array of length 100; (2) ASCII codes from 000 to 999 are pre-calculated and stored in a 32-bit integer array of length 1000; (3) ASCII codes in the form of "x.xx" are pre-calculated, ranging from "0.00" to "9.99", and stored in a 32-bit integer array of length 1000. Print output When , you can find the ASCII code corresponding to the value through the array subscript. There are 17 bits to be output, and the total of 6 values ​​can be calculated by "3+3+3+2+3+3" The complete print result.

[0404] 2.2.4 Precomputed exponential ASCII code

[0405] The decimal exponent range of a double-precision floating-point number is [-324, 308]. Precompute the results from "e-324\0" to "e+308\0" and store them in a 64-bit integer array of length 633. The ASCII code representation of the exponent can be determined through the array subscript.

[0406] 2.3 Acceleration using AVX-512 instruction set

[0407] In the printout When using the SIMD (Single Instruction Multiple Data) instruction built-in functions in the Intel AVX-512 instruction set, such as _mm512_set_epi64, _mm512_set1_epi64, _mm_set_epi16, _mm512_permutexvar_epi16, etc., multiple data are initialized, filled, transformed, and other operations are performed simultaneously, thereby improving the computing speed with the help of multiple computing units in the CPU.

[0408] The present invention can be applied to all occasions where IEEE 754 double-precision floating-point numbers need to be converted into fixed-length ASCII code forms, and can replace the standard implementations of sprintf(), std::ostrstream(), and doubleconv() functions in C, C++, JavaScript engines, etc. In particular, in common industrial software applications, there is a need to convert IEEE754 double-precision floating-point numbers into string form and store them in files or output them to the screen. For example, in the field of computational fluid dynamics (CFD), double-precision floating-point numbers are usually used in CAE software to represent grid coordinate points. During the solution process, it is necessary to frequently output all parameters and grid data in ASCII code form to the screen and intermediate result files, and this conversion work is a computationally intensive task. The application of this algorithm can effectively improve the efficiency of floating-point number conversion and help to shorten the parallel solution time as a whole.

[0409] Table 1 lists six different implementations of this algorithm (d2sci). d2sci_lut represents the original implementation using only a lookup table. d2sci_sse, d2sci_avx512_lut, and d2sci_avx512 implement three optimization strategies. In these implementations, only one output buffer is used to cache the printed ASCII codes. d2sci_32 and d2sci_32v use the AVX-512 instruction set to print 32 double-precision floating-point numbers in parallel. The only difference between them is the number of output buffers, which affects the memory usage during algorithm runtime.

[0410] Table 1 Description of different implementation versions of this algorithm

[0411]

[0412] This algorithm is compared with the performance of the representative algorithms shown in Table 2. Figure 5 The table shows a comparison of the time taken by different printing algorithms for printing completely random double-precision floating-point numbers, measured in nanoseconds. For completely random floating-point numbers, traditional floating-point printing algorithms, such as ryu, take an average of approximately 25.29 nanoseconds to print a single floating-point number. In comparison, our original d2sci_lut algorithm has an average execution time of only 4.53 nanoseconds, a performance improvement of approximately 5.6 times. When processing large numbers of floating-point numbers, the d2sci_32 and d2sci_32v algorithms have an average execution time of 1.24 nanoseconds per number, equivalent to approximately 7-8 CPU clock cycles, a 20.4-fold performance improvement compared to the ryu algorithm. This significant performance improvement is primarily due to the use of the AVX-512 instruction set and algorithmic optimizations, which reduce the number of instructions and computational complexity.

[0413] Table 2 All algorithms selected for benchmarking

[0414]

[0415] Figure 6 The average execution time of different algorithms at different precision levels in decimal scientific notation is shown. All six implementations exhibit nearly constant execution time across different precision levels. In contrast, the performance of existing algorithms fluctuates depending on precision, with the general trend being that execution time increases with increasing precision. This phenomenon is due to the fact that our algorithm outputs fixed-length fractions and is therefore less sensitive to changes in precision. Overall, our approach demonstrates significant performance advantages over existing methods.

[0416] It should be noted that the embodiments of the present invention can be implemented by hardware, software, or a combination of software and hardware. The hardware portion can be implemented using dedicated logic; the software portion can be stored in a memory and executed by an appropriate instruction execution system, such as a microprocessor or dedicated design hardware. Those skilled in the art will appreciate that the above-mentioned devices and methods can be implemented using computer-executable instructions and / or contained in processor control code, for example, such as a carrier medium such as a disk, CD or DVD-ROM, a programmable memory such as a read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuits such as very large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., can also be implemented by software executed by various types of processors, or can be implemented by a combination of the above-mentioned hardware circuits and software, such as firmware.

[0417] The above description is only a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions and improvements made by any technician familiar with this technical field within the technical scope disclosed by the present invention and within the spirit and principles of the present invention should be covered by the scope of protection of the present invention.

Claims

1. A fast printing method for IEEE 754 double-precision floating-point numbers based on the AVX-512 instruction set, characterized in that: The method is executed by a processor supporting the AVX-512 instruction set and includes the following steps: S1. Calculate the scientific notation exponent b of floating-point numbers in batches using AVX-512 parallel instructions. S2. Express the floating-point number x as a×10^b, multiply a by 10 to the power of 16 using SIMD instructions, and round it off to obtain the printed integer value. S3. Quickly map a and b to a print output format through table lookup and vector ASCII code conversion instructions.

2. The method according to claim 1, wherein The exponent b is extracted by the hardware instruction auxiliary module from the exponent of the IEEE754 floating point number and is inferred based on the parallel lg estimation instruction, and finally the exponent range judgment and the determination of the integer b are completed through a SIMD vector comparison and selection instruction.

3. The method according to claim 2, wherein An exponential z extraction acceleration circuit is provided in the processor for: a) Under normalized numbers, directly decode the floating point exponent and subtract the bias value 1023; b) For denormalized numbers, the hardware clz (Count Leading Zeros) instruction is used to quickly count the leading zeros in the mantissa and estimate the exponent z based on this count, shortening the logic branches and judgment paths.

4. The method according to claim 1, wherein The calculation of the mantissa a is completed by a custom floating-point multiplication unit implemented in FPGA or ASIC, which supports truncation of the high 64 bits of the 64-bit integer multiplication result, and cooperates with the exponent register update and mantissa rounding mechanism to complete the generation of high-precision integer a.

5. The method according to claim 4, wherein The custom floating-point multiplication unit includes: a) a mantissa multiplication logic array for calculating the upper 64 bits of the product of two 64-bit unsigned integers; b) exponent accumulator logic for adding two 32-bit signed integers; c) Mantissa rounding controller, used to detect the low bit of the product and determine whether to round it up or down.

6. The method according to claim 1, wherein The output module includes parallel lookup table logic and ASCII code generation circuit. The lookup table logic reads three groups of ASCII codes corresponding to the upper 9 digits from the three-segment compressed lookup table ROM and writes them into the character buffer at high speed; the lower 8 bits are calculated in real time by the pipeline parallel generation circuit.

7. The method according to claim 1, wherein The integer powers of 10 are pre-calculated and stored in a high-speed read-only buffer. The processor or chip quickly obtains the corresponding custom floating-point format {F, E} through indexing, reducing the calculation dependence on lg and power functions during runtime.

8. An IEEE 754 double-precision floating-point number fast printing system based on the AVX-512 instruction set, characterized in that: Includes the following modules, all implemented through software and hardware collaboration: a) Exponent calculation module, which extracts floating-point exponent b based on hardware decoding and AVX-512 SIMD estimation instructions; b) A custom multiplication module that calculates the product of floating-point numbers and powers of 10 and rounds off the mantissa based on a custom accelerator; c) Table lookup output module, integrating three-segment high and low bit table lookup structure and ASCII code converter, supporting pipeline parallel output.

9. The system according to claim 8, wherein The index calculation module includes: a) Hardware floating-point decoder for parsing the IEEE 754 format to extract the exp field; b) a leading zero detector to handle exponential calculation of denormalized numbers; c) A table lookup multiplication factor module is used to quickly map z to the decimal exponent b value corresponding to lg2×z.

10. The system according to claim 8, wherein The table lookup output module includes: a) Compression coding ROM lookup table, used for three-segment lookup table mapping of the high-order ASCII code output of the mantissa a; b) a pipeline parallel decimal encoder for calculating the low-order mantissa into ASCII code; c) an index b lookup table mapping unit, used to convert the integer b into a corresponding exponential character output format.

Images