A method for floating-point program error optimization based on expression rewriting and error-free transformation

Abstract

The application discloses a kind of based on expression rewrites and error-free conversion floating point program error optimization method, belong to floating point operation technical field, to reduce the error of floating point calculation;Including the following steps: through the fusion static and dynamic interval analysis accurate derivation variable in floating point program range of values;Based on abstract syntax tree analysis, identify error sensitive substructure that can easily lead to numerical cancellation, interval across zero or order imbalance;For error sensitive node, using error-free splitting algorithm, by integer-decimal or bit width division mode decomposition operand, construct compensation expression to retain more effective numerical bits;Combined with rewrite rule library, design three kinds of strategies such as constant decomposition, variable splitting and operation sequencing, generate mathematical semantic equivalent but lower error optimization expression, finally generate lower target program than initial program floating point calculation error.The application does not need to improve floating point number precision type, can effectively reduce operation error, improve the numerical stability of program as a whole.

Description

Technical Field

[0001] This invention relates to a floating-point program error optimization method, specifically a floating-point program error optimization method based on expression rewriting and error-free transformation, belonging to the field of floating-point arithmetic technology. Background Technology

[0002] In modern high-performance computing, scientific simulation, engineering modeling, and other fields, floating-point programs serve as the core of numerical computation, undertaking a large number of real number operations. Because computer systems use finite precision to represent real numbers, floating-point numbers inevitably introduce rounding errors during representation and calculation. Without effective control, these errors can lead to numerical instability or even invalid results, severely impacting the reliability of upper-level application systems.

[0003] Currently, various optimization and detection methods for floating-point errors have been proposed. Mainstream techniques include: error detection and patching methods, improving floating-point precision types, using correct rounding techniques, and expression rewriting methods based on mathematical equivalence transformations. Among these, expression rewriting has received widespread attention due to its advantages of not requiring changes to program semantics and having low overhead. Representative tools such as Herbie can rewrite expressions within a specified input range to reduce errors; tools like FPTaylor calculate the upper bound of errors using symbolic Taylor expansion; and Herbgrind combines dynamic analysis to identify error-sensitive regions. These tools have achieved certain results in expression-level error optimization.

[0004] However, existing technologies still have the following shortcomings in practical engineering applications: 1) Most expression optimization tools rely on fixed intervals or sampling points, making it difficult to adapt to the dynamic changes in variable intervals during the execution of real floating-point programs; 2) Mathematical equivalent transformations performed in the program may not necessarily reduce errors within local intervals, and may even amplify errors; 3) There is a lack of active identification and intervention mechanisms for error-sensitive sub-expressions in the floating-point program structure, resulting in poor optimization effects on ill-conditioned structures; 4) Expression rewriting usually fails to combine specific error propagation mechanisms and compensation theories, making it difficult to fundamentally suppress error sources.

[0005] Therefore, there is an urgent need for a floating-point program error optimization method that can systematically identify and reconstruct floating-point expressions that are prone to errors, thereby improving the overall numerical accuracy and stability of the program. Summary of the Invention

[0006] The purpose of this invention is to address the problems of unstable accuracy, insufficient error analysis, and weak structural intervention capabilities in existing floating-point program error optimization methods. It provides a floating-point program error optimization method based on expression rewriting and error-free transformation. By integrating static and dynamic interval analysis, structural pre-analysis, and error compensation mechanisms, it effectively identifies and reconstructs error-sensitive expressions, thereby improving the overall numerical stability and computational accuracy of the program.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: a floating-point program error optimization method based on expression rewriting and error-free transformation, comprising the following steps: S1, Interval Analysis Stage: For floating-point programs with interval information as input, a combination of static analysis and dynamic sampling is used to perform mixed interval analysis on the variables in the floating-point program and derive the value range of each variable in the floating-point program. Static analysis derives the dependencies between variables based on Abstract Syntax Tree (AST) and Control Flow Graph (CFG), constructs a global interval model covering all executable paths, propagates the initial intervals of variables, and obtains static interval estimates of variables. Dynamic analysis obtains the actual value range of variables during actual operation by instrumenting and sampling the target program, and uses affine arithmetic to optimize the dynamic sampling results to tighten the interval boundaries. Finally, by fusing the results of static and dynamic analysis, the precise set of intervals of variables is output. S2, Structural Pre-analysis Stage: Using the variable interval information generated after fusion in step S1 as input, for each operation node in the expression, the abstract syntax tree structure of the expression is analyzed, and high-risk nodes in the error-sensitive area of ​​the floating-point expression are identified. Key substructures that may amplify errors are screened and marked as high-error-risk structures, including difference amplification structures that cancel out errors due to small numerical differences, sign-sensitive structures where positive and negative values ​​cross the zero interval, and order-of-magnitude imbalance structures that cause underflow risk due to the interaction of constants and extremely small variables. Finally, all the identification results are used to form a set of error-sensitive nodes with high errors. S3, Error Compensation Decomposition Stage: Static analysis of the floating-point program source code was performed. For the error-sensitive nodes identified in step S2, the SplitAM algorithm without error transformation was used to split the operands into multiple sub-items, generating high-order principal values ​​and low-order compensation items, retaining more effective numerical bits, and reducing the accumulation of rounding errors. S4. Expression rewriting stage and optimization process: Based on the error compensation decomposition results in step S3, optimization is performed using the expression rewriting tool Herbie and the rewriting rule library. Three custom rewriting strategies are introduced: constant decomposition strategy, variable decomposition strategy, and operation sorting strategy. Each strategy generates multiple candidate expressions while ensuring mathematical equivalence. Interval arithmetic is used to evaluate their numerical errors within a specific variable range, and a score is calculated using an error function. Finally, the expression with the lowest error score is selected as the optimal rewriting result, completing the optimization process and outputting the error-optimized floating-point program.

[0008] In step S1, the static interval analysis process is as follows: the initial interval of the input variables is extracted from the floating-point program entry point. The initial interval is specified by the user or set by default to the legal range of floating-point representation. Then, forward data flow analysis is performed along the control flow graph, and the initial interval of the variables is propagated according to the topological order of the basic blocks or by iterating to the fixed point. For each arithmetic operation node, the interval is propagated conservatively according to its basic interval arithmetic rules. For division and nonlinear operations, the corresponding interval operation rules are combined and the case of crossing zero is considered. Finally, the static interval estimates of each intermediate variable and the output variable are obtained.

[0009] For interval operations, variables are represented as intervals [l, u] to handle uncertainties in computation. The rules for interval operations are as follows: Addition is [x] + [y] = [x] + [y ... l + y l , x u + y u Subtraction is [x] - [y] = [x] l - y u , x u - y l The multiplication is [x] * [y] = [min(S), max(S)], where S = {x} l y u , x u y u , x u y l , x l y u Interval division follows the same rules as multiplication.

[0010] In step S1, the dynamic interval analysis process is as follows: Instrumentation positions are set at key program points, including loop bodies, conditional branch exits, and before and after function calls; during dynamic sampling, a uniform sampling strategy is used to randomly generate test inputs within the static interval; simultaneously, for high-risk areas identified by static analysis, a priority sampling strategy is introduced to improve coverage of error-sensitive intervals; through multiple executions of the instrumentation program, runtime samples of variables are collected to form empirical intervals; affine arithmetic is introduced to optimize the representation of the dynamic sampling results, representing variables in an affine form containing central values ​​and error terms to characterize the linear correlation between variables, further tightening the interval boundaries.

[0011] In step S1, after calculating the static and dynamic intervals, the final variable interval is generated by interval fusion: let the static interval be ([s_l, s_u]) and the dynamic interval be ([d_l, d_u]). The final interval ([max(s_l, d_l), min(s_u, d_u)]) is obtained by taking the intersection of the two intervals. If the dynamic interval is sufficiently sampled and strictly contained within the static interval, the dynamic result is used first, or the numerical interval with a narrower range and a boundary that is closer to the actual value is selected by combining the sampling coverage and uncertainty assessment.

[0012] In step S2, the screening process for key substructures is as follows: First, it is determined whether there is a risk of numerical cancellation among the operands in each subexpression. If the center values ​​of the intervals of two operands are similar and the subtraction result is close to zero, the structure is marked as a difference amplification structure. Second, it is detected whether the intervals of the expression sub-items cross zero. If there are intervals containing positive and negative signs, the structure is identified as a sign-sensitive structure. In addition, it is also detected that the combined structure of constants and variables with unbalanced magnitudes participates in the operation, and the operation nodes that may cause underflow or rounding to mask the contribution of variables and cause errors are identified. Such structures are marked as unbalanced magnitude structures. All identified error-sensitive nodes will be used as candidate compensation positions, and their syntax tree path, variable name, interval, and magnitude information will be recorded for reference in subsequent steps.

[0013] In the screening of the key substructure, the threshold for the error evaluation value of the two operands is set in the range of 0~0.5ulp, and when the result of the floating-point operation error is greater than 0.5ulp, the node is defined as an error-sensitive node.

[0014] In step S3, the SplitAM algorithm enhances the precision of variables prone to errors in floating-point addition and multiplication operations. The SplitAM algorithm includes two splitting modes: "integer-fractional" and "bit-level segmentation." The optimal splitting method is selected based on the interval characteristics of the variables: 1) Using the "integer-fractional" splitting mode, for addition operations, the operand is split into an integer part and multiple fractional parts of different magnitudes to avoid the loss of effective bits caused by exponent alignment; for multiplication operations, the operand is split into a principal value part and an error compensation part to improve the precision of the product result; 2) Using the "bit-level segmentation" mode, the mantissa is split according to a preset bit width based on the improved SplitR algorithm to generate a set of sub-variables with more balanced numerical magnitudes; all split variables will be uniformly identified as high-bit principal values ​​and low-bit compensation terms, and participate in the construction of subsequent expressions.

[0015] In step S3, the specific implementation of the error-compensated splitAM algorithm is as follows: Function: Split floating-point data into appropriate floating-point arrays Input: variable a, operation type op_type, model, parameter r Output: the split array [] 1. if mode == "integer-frac" then 2. int_part = floor(a) 3. frac = a - int_part 4. if op_type == "add" then 5. [x2..xn] =extractRange(frac, [1, 2^(-r)],n) 8. return [int_part, x2, ..., xn] 9. else if op_type == "mul" then 10. return [int_part, frac] 11. else if mode == "split-by-bits" then 12. if op_type == "add" then 13. [x1, x2, ..., xn] = SplitR(a, r_bits) 14. return [x1, x2, ..., xn] 15. else if op_type == "mul" then 16. [x_hi, x_lo] =SplitR(a, r_bits) 17. return [x_hi, x_lo] 18. endif.

[0016] In step S4, the constant decomposition strategy is as follows: the floating-point constant in the expression is decomposed into several sub-items that are close to the magnitude of the variable, so that the magnitude of the value is closer to the adjacent variable and the exponent alignment error is reduced; the constRewrgte algorithm extracts the constant in the expression and combines the variable interval information to generate a variety of constant decomposition forms to construct candidate polynomials. The specific implementation of constant decomposition using the constRewrite algorithm is as follows: Function: Reduces floating-point expression errors by converting constant values ​​to constant form. Input: floating-point expression `expr`, variable range file `intervals` Output: The optimized expression opt_expr 1. c = ExtractConstants(expr) 2. decomposition= GenerateDecompositions(c, intervals, constants) 3. best_expr = expr 4. min_score = ∞ 5. for each d in decompositions do: 6. candidate_expr = ApplyDecomposition(expr, d) 7. candidate_ast = ParseAST(candidate_expr) 8. result_range = EvaluateInterval(candidate_ast, intervals) 9. round_error = EstimateError(candidate_ast, intervals) 10. score = Error(result_range, round_error) 11. if score <min_score then: 12. min_score = score 13. opt_expr = candidate_expr 14. return opt_expr.

[0017] The variable decomposition process is as follows: using the error compensation decomposition result in step S3 to retain the principal value and compensation term, reducing error propagation through structural adjustment, the varRewrite algorithm extracts variables from the expression and identifies sensitive operation nodes that may cause large errors by combining interval information; the SplitAM algorithm is used to decompose the variables into high-order principal values ​​and low-order compensation terms for sensitive operation nodes, and the expression is reconstructed by combining the compensation terms. The specific implementation of the variable splitting varRewrite algorithm is as follows: Function: Reduces floating-point expression errors by changing the variable form. Input: floating-point expression `expr`, variable range file `intervals` Output: The optimized expression opt_expr 1. vars = extractVariables(expr) 2. criticalNodes = IdentifySensitiveNodes(vars, intervals) 3. for each node in criticalNodes do 4. if node.op == "+" or node.op == "-" then 5. node = SplitAM(node, add, mode, r) 6. elif node.op == "*" then 7. node = SplitAM(node, mul, mode, r) 8. endif 9. endfor 10. rewrittenExpr = herbieRewrite(expr) 11. return rewrittenExpr.

[0018] The sorting process is as follows: The orderRewrite algorithm first extracts the variables in the expression and obtains their numerical range [l, u] through the interval file. It calculates the center value of the variables and sorts them to ensure that variables with similar numerical magnitudes participate in the operation first. Then, it reconstructs the abstract syntax tree (AST) of the expression based on the sorting result and evaluates the error change before and after rewriting by combining interval arithmetic. If the error decreases significantly after optimization, the sorting result is accepted; otherwise, it reverts to the original expression.

[0019] The specific implementation of the sorting algorithm, orderRewrite, is as follows: Function: Adjusts the order of operations in floating-point representation to reduce errors. Input: floating-point expression `expr`, variable range file `intervals` Output: The optimized expression opt_expr 1. variables = extractVariables(expr) 2. centers = ∅ 3. for each var in variables do 4. interval = lookupInterval(intervals, var) 5. μ = (interval.lower + interval.upper) / 2 6. centers[var] = μ 7. endfor 8. sortedVars = sortByValue(centers) 9. sortedExpr = rebuildExpression(expr, sortedVars) 10. Return sortedExpr.

[0020] The beneficial effects of this invention are: 1) The method of this invention integrates static and dynamic interval analysis, expression structure pre-analysis and error compensation mechanism, which can systematically identify error-sensitive substructures in floating-point programs. By combining expression rewriting and error-free transformation techniques, it can effectively reduce the accumulation and propagation of rounding errors in floating-point operations and improve the overall numerical stability of the program.

[0021] 2) In the method of the present invention, by analyzing the abstract syntax tree structure of the expression, three types of error-sensitive subtrees are identified, including numerical cancellation, interval crossing zero, and order of magnitude imbalance. This step is based on operand interval, structural pattern matching and sign change judgment to extract potential high-error nodes and construct a set of nodes to be optimized.

[0022] 3) The method of this invention adopts three types of expression rewriting strategies: constant decomposition, variable splitting, and operation sorting. The three strategies correspond to three rewriting algorithms: constRewrite, varRewrite, and orderRewrite. Finally, the optimal rewritten expression is selected and written back to the program through error evaluation, thereby improving the accuracy of the expression in order to optimize program error.

[0023] 4) The variable splitting strategy in this invention is based on the idea of ​​error-free transformation and integrates the compensation algorithm SplitAM. While maintaining the original mathematical semantics, it generates equivalent expressions with smaller numerical errors. This method can effectively repair most high-error calculation results, improve the error suppression capability of the program and the flexibility of expression structure reconstruction. Furthermore, the SplitAM algorithm is used to split variables into high-order principal values ​​and low-order compensation terms, which enhances the robustness of the expression in ill-conditioned structures. The operation sorting strategy rearranges the operation order in the syntax tree according to the center value of the variable interval, avoiding the amplification of rounding errors caused by excessive numerical differences.

[0024] 5) The present invention tested 10 functions in the GNU Scientific Computing Library (GSL), which fully demonstrates that the method of the present invention has a good error correction effect on floating-point programs in some practical applications, without bringing significant performance overhead, effectively reducing calculation errors in intervals, and achieving a good error defect correction effect on the function gsl_sf_legendre_P2. Attached Figure Description

[0025] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a graph showing the ill-conditioned interval error of gsl_sf_legendre_P3 obtained from an embodiment of the method of the present invention. Figure 3 This is a graph showing the error effect of gsl_sf_legendre_P3 in the [0,1] interval obtained in an embodiment of the method of the present invention; Figure 4 This is a comparison chart showing the changes in the computation speed of the test function obtained from an embodiment of the method of the present invention. Detailed Implementation

[0026] The present invention will be further explained and described below with reference to the accompanying drawings and specific embodiments.

[0027] Example:Figure 1 As shown, this invention provides a floating-point program error optimization method based on expression rewriting and error-free transformation, comprising the following steps: S1, Interval Analysis Stage: For floating-point programs with interval information as input, a hybrid approach combining static analysis and dynamic sampling is used to perform mixed interval analysis on the variables in the floating-point program, deriving the value range of each variable. Static analysis derives the dependencies between variables based on Abstract Syntax Trees (ASTs) and Control Flow Graphs (CFGs), constructs a global interval model covering all executable paths, propagates the initial intervals of variables, and obtains static interval estimates of the variables. Dynamic analysis obtains the actual value range of variables during actual operation by instrumenting and sampling the target program, and uses affine arithmetic to optimize the dynamic sampling results and tighten the interval boundaries. Finally, the precise set of intervals for the variables is output by fusing the results of static and dynamic analysis.

[0028] The interval analysis stage aims to accurately derive the value range of each variable in the floating-point program, providing a reliable input interval basis for subsequent error-sensitive structure identification, error compensation, and expression rewriting. This stage integrates two complementary methods, static analysis and dynamic analysis, to model and constrain the variable value range, achieving accurate interval estimation.

[0029] The static interval analysis process is as follows: the initial interval of the input variables is extracted from the floating-point program entry point. The initial interval is specified by the user or set by default to the legal range of floating-point representation. Then, forward data flow analysis is performed along the control flow graph, and the initial interval of the variables is propagated according to the topological order of the basic blocks or by iterating to the fixed point. For each arithmetic operation node, the interval is conservatively propagated according to its basic interval arithmetic rules. For division and nonlinear operations, the corresponding interval operation rules are combined and the case of crossing zero is considered. Finally, the static interval estimates of each intermediate variable and output variable are obtained.

[0030] For interval operations, variables are represented as intervals [l, u] to handle uncertainties in computation. The rules for interval operations are as follows: Addition is [x] + [y] = [x] + [y ... l + y l , x u + y u Subtraction is [x] - [y] = [x] l - y u , x u - y l The multiplication is [x] * [y] = [min(S), max(S)], where S = {x} l y u , x u y u , xu y l , x l y u Interval division follows the same rules as multiplication.

[0031] The dynamic interval analysis process is as follows: Instrumentation positions are set at key program points, including loop bodies, conditional branch exits, and before and after function calls; during dynamic sampling, a uniform sampling strategy is used to randomly generate test inputs within the static interval; simultaneously, for high-risk areas identified by static analysis, a priority sampling strategy is introduced to improve coverage of error-sensitive intervals; through multiple executions of the instrumented program, runtime samples of variables are collected to form empirical intervals; affine arithmetic is introduced to optimize the representation of the dynamic sampling results, representing variables in an affine form containing central values ​​and error terms to characterize the linear correlation between variables, further tightening the interval boundaries.

[0032] After calculating the static and dynamic intervals, the final variable interval is generated through interval fusion: Let the static interval be ([s_l, s_u]) and the dynamic interval be ([d_l, d_u]). A conservative strategy of taking the intersection is used to obtain the final interval ([max(s_l, d_l), min(s_u, d_u)]). If the dynamic interval is sufficiently sampled and strictly contained within the static interval, the dynamic result is prioritized, or a numerical interval with a narrower range and boundaries that better match the actual values ​​is selected based on a combination of sampling coverage and uncertainty assessment. The final set of output variable intervals serves as the input basis for the subsequent structural pre-analysis stage.

[0033] S2, Structural Pre-analysis Stage: Before transforming the expressions, it is necessary to obtain a set of expressions with high errors. Taking the variable interval information generated after fusion in step S1 as input, for each operation node in the expression, the abstract syntax tree structure of the expression is analyzed. High-risk nodes in the error-sensitive areas of the floating-point expression are identified, and key substructures that may amplify errors are screened and marked as high-error-risk structures. These include difference amplification structures where small numerical differences lead to the cancellation of errors, sign-sensitive structures where positive and negative values ​​cross the zero interval, and order-of-magnitude imbalance structures where the interaction between constants and extremely small variables causes the risk of underflow. Finally, all the identification results are used to form a set of error-sensitive nodes with high errors.

[0034] The screening process for key substructures is as follows: First, it determines whether there is a risk of numerical cancellation among operands in each subexpression. If the center values ​​of the intervals of two operands are similar and the subtraction result approaches zero, the structure is marked as a difference amplification structure. Second, it checks whether the intervals of expression terms cross zero. If there are intervals containing positive and negative signs, the structure is identified as a sign-sensitive structure. Furthermore, it detects combinations of constants and variables with unbalanced magnitudes, identifying computational nodes that may cause underflow or rounding to mask variable contributions and lead to errors; such structures are marked as unbalanced magnitude structures. All identified error-sensitive nodes are used as candidate compensation locations, and their syntax tree path, variable name, interval, and magnitude information are recorded for reference in subsequent compensation insertion and expression rewriting stages. This analysis stage can effectively identify ill-conditioned structures in advance, improving subsequent optimization results.

[0035] In the screening of critical substructures, the threshold for the error evaluation value of the two operands is set in the range of 0~0.5ulp, and when the result of the floating-point operation error is greater than 0.5ulp, the node is defined as an error-sensitive node.

[0036] S3, Error Compensation Decomposition Stage: Static analysis of the floating-point program source code is performed. For the error-sensitive nodes identified in step S2, the SplitAM algorithm with no error transformation is used to split the operands into multiple sub-items, generating high-order principal values ​​and low-order compensation terms, retaining more valid numerical bits and reducing the accumulation of rounding errors. In the error compensation splitting stage, for the error-sensitive nodes identified in the structural pre-analysis stage, compensation structures are inserted through variable splitting to reduce calculation errors and maintain the original mathematical semantics.

[0037] The multiple combinations generated in this stage can effectively reduce error propagation, significantly improve the stability and accuracy of the rewritten expression, and provide higher quality candidate structures for the expression rewriting stage.

[0038] S4. Expression rewriting stage and optimization process: The expression rewriting stage aims to reduce floating-point calculation errors through structural transformations without altering the program's mathematical semantics. Based on the error compensation decomposition results from step S3, optimization is performed using the expression rewriting tool Herbie and a rewriting rule library. Three custom rewriting strategies are introduced: constant decomposition strategy, variable decomposition strategy, and operation sorting strategy. Each strategy generates multiple candidate expressions while ensuring mathematical equivalence. Interval arithmetic is used to evaluate their numerical errors within a specific variable range, and a score is calculated using an error function. Finally, the expression with the lowest error score is selected as the optimal rewriting result, completing the optimization process and outputting the error-optimized floating-point program.

[0039] The constant decomposition strategy involves decomposing floating-point constants in an expression into several sub-terms close to the magnitude of variables, making their numerical magnitudes closer to adjacent variables and reducing exponent alignment errors. The constRewrgte algorithm extracts constants from the expression and combines them with variable interval information to generate various constant decomposition forms to construct candidate polynomials.

[0040] The variable splitting process is as follows: using the error compensation splitting result in step S3 to retain the principal value and compensation term, and reducing error propagation through structural adjustment, the varRewrite algorithm extracts variables from the expression and identifies sensitive operation nodes (such as addition, subtraction, and multiplication) that may cause large errors by combining interval information; for sensitive operation nodes, the SplitAM algorithm is used to decompose the variables into high-order principal values ​​and low-order compensation terms, thereby suppressing the propagation of rounding errors. The splitting strategy is adaptively adjusted according to the operation type (addition / subtraction or multiplication) and parameters to generate a more structurally stable expression; then, the extended Herbie rewrite rule is called to perform algebraic rewriting on the new expression containing compensation terms, and the result with smaller error is retained as the final optimized expression.

[0041] For example, the expression `gsl_sf_legendre_P3` in GSL is used to verify the optimization effect of the transformed expression error. By comparing the codes of five different methods shown in Table 1 (original code 1.a, FMA instruction optimized code 1.b, Herbie optimized ill-conditioned interval 1.c, normal interval 1.d, and the method of this invention code 1.e), using the results from the MPFR library as the accurate value, the difference in the number of bits between the result and the accurate value is calculated, as shown in the appendix. Figure 2 The figure shows the error distribution of gsl_sf_legendre_P3 under different methods in the ill-conditioned interval; as shown in the appendix. Figure 3 The figure shows the error distribution of gsl_sf_legendre_P3 within the commonly used interval [0,1]. As can be seen from the figure, the method of this invention can significantly reduce the tail proportion of the error distribution and increase the proportion of errors concentrated in a smaller range in both intervals.

[0042] Table 1 The sorting process is as follows: The orderRewrite algorithm first extracts the variables in the expression and obtains their numerical range [l, u] through the interval file. It calculates the center value of the variables and sorts them to ensure that variables with similar numerical magnitudes participate in the operation first. Then, it reconstructs the abstract syntax tree (AST) of the expression based on the sorting result and evaluates the error change before and after rewriting by combining interval arithmetic. If the error decreases significantly after optimization, the sorting result is accepted; otherwise, it reverts to the original expression.

[0043] Experimental Results and Analysis: 1. Experimental conditions and preparation The experiment of this invention was conducted on the Linux 4.4.15-deepin-aere sw_64 operating system, SW1621 64-bit 16-core processor, 16GB of memory, and the compiler was swgcc-8.3.0, with all compilation options using the default settings.

[0044] This invention selects 10 commonly used single-parameter double-type functions in GSL as test cases to evaluate the effectiveness of the optimization method in function error control. Since GSL functions use different computation paths within different variable intervals, randomly generated input data within the overall interval is insufficient to cover all error regions and lacks representativeness. Therefore, this invention, based on variable interval division, subdivides the input domain into multiple sub-intervals and randomly generates an equal number of test data within each sub-interval. Table 2 below shows the overall interval and the branch intervals of the random input in this invention.

[0045] Table 2 2. Analysis of Experimental Results The experimental subjects were derived from the GNU Scientific Computing Library (GSL). Ten commonly used SpecCn functions and basic functions were selected for testing. IEEE 754 standard double-type data was used as input, the Multi-Precision Computing Library (MPFR) was used as the benchmark, and the proportion of output errors within 1 ULP (a) was used as the evaluation criterion. Specific optimization results are shown in Tables 3 and 4 and Appendix. Figure 4 As shown.

[0046] Table 3 shows the percentage of the floating-point program's output across the four error intervals before and after optimization (unit: %). Table 4 shows the percentage of the four error intervals before and after gsl_sf_airy_Ai optimization (unit: %). From Table 3, Table 4 and Appendix Figure 4The function execution speed test results shown demonstrate that the error correction method proposed in this invention can effectively reduce the calculation error of floating-point programs while maintaining essentially the same program execution efficiency. Table 3 shows that, overall, after optimization, the proportion of errors within the [0,1] range for most functions is significantly improved, with an average improvement of 8.57%. Simultaneously, the proportion of larger errors (error > 8) is reduced in several functions, with an average reduction of 2.01%, indicating that the optimization method has a significant effect on reducing function errors and making the calculation results more accurate. Before optimization, the proportion of errors within the [0,1] range for the gsl_sf_legendre_P2 function was 53.77%, which increased to 75.92% after optimization; the proportion of errors > 8 decreased from 6.19% to 0. This indicates that the optimization has a very significant effect on improving the accuracy of this function, effectively reducing the occurrence of large errors. Since the gsl_sf_legendre_P2 function is also a polynomial function, the program incorporates error compensation in the corrected formula, correcting the large errors in specific intervals. Before optimization, the percentage of errors in the [0, 1] range of the gsl_sf_log_1plusx function was 74.78%, which was significantly improved to 84.12% after optimization. Furthermore, the optimization reduced the error to 0 in the (4, 8] range and the error > 8, demonstrating significant optimization effectiveness.

[0047] For the functions gsl_sf_airy_Ai and gsl_sf_psi, the high error ratio was not effectively improved. More detailed interval detection revealed, as shown in Table 4, that the proportion of high errors remained high in certain branches, resulting in poor repair performance. The repair procedure of this invention primarily improves the mathematical expressions in the program, aiming to enhance their accuracy and stability, thereby reducing overall error. However, even with improved average computational accuracy through expression transformation, ill-conditioned problems cannot be completely avoided. These errors are usually attributed to forward errors, whose condition number characteristics are independent of program execution and only related to the properties of the function itself. Therefore, simply improving program execution accuracy is insufficient to reduce such errors. Thus, the method of this invention is not suitable for handling this type of forward error.

[0048] This invention treats a program as a set of expressions and combines interval analysis with error-free transformation techniques, focusing on expression transformation to reduce the error of individual expressions, thereby achieving error optimization in floating-point programs. Experiments on the GNU Scientific Computing Library (GSL) have fully demonstrated that this method has a good accuracy repair effect on floating-point programs in some practical applications without incurring additional time overhead. It can effectively reduce the computational error of programs and has achieved good repair results for the function gsl_sf_legendre_P2.

[0049] The above description is only used to illustrate the technical solution of the present invention and is not intended to limit it. Any other modifications or equivalent substitutions made by those skilled in the art to the technical solution of the present invention, as long as they do not depart from the spirit and scope of the technical solution of the present invention, should be covered within the scope of the claims of the present invention.

Claims

1. A floating-point program error optimization method based on expression rewriting and error-free transformation, characterized in that: Includes the following steps: S1, Interval Analysis Stage: For floating-point programs with interval information as input, a combination of static analysis and dynamic sampling is used to perform mixed interval analysis on the variables in the floating-point program and derive the value range of each variable in the floating-point program. Static analysis derives the dependencies between variables based on Abstract Syntax Tree (AST) and Control Flow Graph (CFG), constructs a global interval model covering all executable paths, propagates the initial intervals of variables, and obtains static interval estimates of variables. Dynamic analysis obtains the actual value range of variables during actual operation by instrumenting and sampling the target program, and uses affine arithmetic to optimize the dynamic sampling results to tighten the interval boundaries. Finally, by fusing the results of static and dynamic analysis, the precise set of intervals of variables is output. S2, Structural Pre-analysis Stage: Using the variable interval information generated after fusion in step S1 as input, for each operation node in the expression, the abstract syntax tree structure of the expression is analyzed, and high-risk nodes in the error-sensitive area of ​​the floating-point expression are identified. Key substructures that may amplify errors are screened and marked as high-error-risk structures, including difference amplification structures that cancel out errors due to small numerical differences, sign-sensitive structures where positive and negative values ​​cross the zero interval, and order-of-magnitude imbalance structures that cause underflow risk due to the interaction of constants and extremely small variables. Finally, all the identification results are used to form a set of error-sensitive nodes with high errors. S3, Error Compensation Decomposition Stage: Static analysis of the floating-point program source code was performed. For the error-sensitive nodes identified in step S2, the SplitAM algorithm without error transformation was used to split the operands into multiple sub-items, generating high-order principal values ​​and low-order compensation items, retaining more effective numerical bits, and reducing the accumulation of rounding errors. S4. Expression rewriting stage and optimization process: Based on the error compensation decomposition results in step S3, optimization is performed by combining the expression rewriting tool Herbie and the rewriting rule library. At the same time, three types of custom rewriting strategies are introduced: constant decomposition strategy, variable decomposition strategy and operation sorting strategy. Each strategy generates multiple candidate expressions while ensuring mathematical equivalence. The numerical error of each expression within a specific variable range is evaluated using interval arithmetic, and a score is calculated using an error function. Finally, the expression with the lowest error score is selected as the optimal rewrite result, and the optimized floating-point program is output after the optimization process is completed.

2. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S1, the static interval analysis process is as follows: the initial interval of the input variables is extracted from the floating-point program entry point. The initial interval is specified by the user or set by default to the legal range of floating-point representation. Then, forward data flow analysis is performed along the control flow graph, and the initial interval of the variables is propagated according to the topological order of the basic blocks or by iterating to the fixed point. For each arithmetic operation node, the interval is propagated conservatively according to its basic interval arithmetic rules. For division and nonlinear operations, the corresponding interval operation rules are combined and the case of crossing zero is considered. Finally, the static interval estimates of each intermediate variable and the output variable are obtained.

3. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 2, characterized in that: For interval operations, variables are represented as intervals [l, u] to handle uncertainties in computation. The rules for interval operations are as follows: Addition is [x] + [y] = [x] + [y ... l + y l , x u + y u Subtraction is [x] - [y] = [x] l - y u ,x u - y l The multiplication is [x] * [y] = [min(S), max(S)], where S = {x} l y u , x u y u , x u y l , x l y u Interval division follows the same rules as multiplication.

4. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S1, the dynamic interval analysis process is as follows: Instrumentation positions are set at key program points, including loop bodies, conditional branch exits, and before and after function calls; during dynamic sampling, a uniform sampling strategy is used to randomly generate test inputs within the static interval; simultaneously, for high-risk areas identified by static analysis, a priority sampling strategy is introduced to improve coverage of error-sensitive intervals; through multiple executions of the instrumentation program, runtime samples of variables are collected to form empirical intervals; affine arithmetic is introduced to optimize the representation of the dynamic sampling results, representing variables in an affine form containing central values ​​and error terms to characterize the linear correlation between variables, further tightening the interval boundaries.

5. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S1, after calculating the static and dynamic intervals, the final variable interval is generated by interval fusion: let the static interval be ([s_l, s_u]) and the dynamic interval be ([d_l, d_u]). The final interval ([max(s_l, d_l), min(s_u, d_u)]) is obtained by taking the intersection of the two intervals. If the dynamic interval is sufficiently sampled and strictly contained within the static interval, the dynamic result is used first, or the numerical interval with a narrower range and a boundary that is closer to the actual value is selected by combining the sampling coverage and uncertainty assessment.

6. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S2, the screening process for key substructures is as follows: First, it is determined whether there is a risk of numerical cancellation among the operands in each subexpression. If the center values ​​of the intervals of two operands are similar and the subtraction result is close to zero, the structure is marked as a difference amplification structure. Second, it is detected whether the intervals of the expression sub-items cross zero. If there are intervals containing positive and negative signs, the structure is identified as a sign-sensitive structure. In addition, it is also detected that the combined structure of constants and variables with unbalanced magnitudes participates in the operation, and the operation nodes that may cause underflow or rounding to mask the contribution of variables and cause errors are identified. Such structures are marked as unbalanced magnitude structures. All identified error-sensitive nodes will be used as candidate compensation positions, and their syntax tree path, variable name, interval, and magnitude information will be recorded for reference in subsequent steps.

7. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 6, characterized in that: In the screening of the key substructure, the threshold for the error evaluation value of the two operands is set in the range of 0~0.5ulp, and when the result of the floating-point operation error is greater than 0.5ulp, the node is defined as an error-sensitive node.

8. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S3, the SplitAM algorithm performs precision enhancement processing on variables that are prone to errors in addition and multiplication operations in floating-point operations. The SplitAM algorithm includes two splitting modes: "integer-fractional" and "digit segmentation". The optimal splitting method is selected based on the interval characteristics of the variable: 1) The "integer-fractional" splitting mode is adopted. For addition operations, the operand is split into an integer part and multiple fractional parts of different orders of magnitude to avoid the loss of effective bits caused by exponent alignment. For multiplication operations, the operand is split into a main value part and an error compensation part to improve the precision of the product result. 2) The "bit segmentation" mode is adopted. Based on the improved SplitR algorithm, the mantissa is segmented according to the preset bit width to generate a set of sub-variables with more balanced numerical magnitudes. All the split variables will be uniformly identified as high-bit principal values ​​and low-bit compensation terms, and participate in the construction of subsequent expressions.

9. The floating-point program error optimization method based on expression rewriting and error-free transformation according to claim 1, characterized in that: In step S4, the constant decomposition strategy is as follows: the floating-point constant in the expression is decomposed into several sub-items that are close to the magnitude of the variable, so that the magnitude of the value is closer to the adjacent variable and the exponent alignment error is reduced; the constRewrgte algorithm extracts the constant in the expression and combines the variable interval information to generate a variety of constant decomposition forms to construct candidate polynomials. The variable decomposition process is as follows: using the error compensation decomposition result in step S3 to retain the principal value and compensation term, reducing error propagation through structural adjustment, the varRewrite algorithm extracts variables from the expression and identifies sensitive operation nodes that may cause large errors by combining interval information; the SplitAM algorithm is used to decompose the variables into high-order principal values ​​and low-order compensation terms for sensitive operation nodes, and the expression is reconstructed by combining the compensation terms. The sorting process is as follows: The orderRewrite algorithm first extracts the variables in the expression and obtains their numerical range [l, u] through the interval file. It calculates the center value of the variables and sorts them to ensure that variables with similar numerical magnitudes participate in the operation first. Then, it reconstructs the abstract syntax tree (AST) of the expression based on the sorting result and evaluates the error change before and after rewriting by combining interval arithmetic. If the error decreases significantly after optimization, the sorting result is accepted. Otherwise, revert to the original expression.

Images