{'fraction_non_alphanumeric': 0.0388651379712398, 'fraction_numerical': 0.018752429071123202, 'mean_word_length': 4.587585957292798, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast and stable computation of higher-order Hahn polynomials and Hahn moment invariants for signal and image analysis

This article presents, on the one hand, new algorithms for the fast and stable computation of discrete orthogonal Hahn polynomials of high order (HPs) based on the elimination of all gamma and factorial functions that cause the numerical fluctuations of HPs, and based on the use of appropriate stability conditions. On the other hand, a new method for the fast and numerically stable computation of Hahn moment invariants (HMIs) is also proposed. This method is mainly based on the use of new recursive relations of HPs and of matrix multiplications when calculating HMIs. To validate the efficiency of the algorithms proposed for the calculation of HPs, several signals and large images (≥4000 × 4000) are reconstructed by Hahn moments (HMs) up to the last order with a reconstruction error tending towards zero (MSE ≃ 10−10). The efficiency of the proposed method for calculating HMIs is demonstrated on large medical images (2048 × 2048) with a very low relative error (RE ≃ 10−10). Finally, comparisons with some recent work in the literature are provided.

In general, the computation of DOPs is performed recursively either with respect to the polynomial order n or with respect to the discrete variable x. Mukundan et al. discuss in [27] that the recursive computation of high-order Tchebichef polynomials with respect to the order n is unstable because of the propagation of numerical errors. To minimize this problem, Mukundan et al. use the recursive computation of Tchebichef polynomials with respect to the variable x rather than the order n. Similarly, Zhu et al. [50] propose a general form for the recursive computation of Tchebichef, Krawtchouk, Hahn, Meixner and Charlier polynomials with respect to the variable x for minimizing the propagation of numerical errors. Although the recurrence approach with respect to the variable x minimizes the propagation of numerical errors, this approach remains limited because it is unable to generate the polynomial values of high-order polynomials [1,9]. To overcome this problem for Tchebichef polynomials, Abdulhussain et al. [1] propose a method based on the sequential integration of the traditional recurrence relations with respect to order n and with respect to variable x that were introduced by Mukundan et al. [27]. However, this proposed method is not stable for higher order Tchebichef polynomials due to the problem of numerical errors propagation [5]. To overcome this problem, Camacho-bello et al. [5] propose the computation of higher-order Tchebichef polynomials via the modified Gram-Schmidt orthonormalization process (GSOP). This method succeeded to ensure the numerical stability of higher-order Tchebichef polynomials. More recently, GSOP method is used by Daoui et al. to introduce a numerically stable computation of Hahn [9], Charlier [8] and Meixner [12] polynomials.
Even if the problems of recursive computation of certain high-order polynomials are solved, the computation of high-order Hahn polynomials (HPs) is still limited by the following two major problems: (i) the numerical instability of the initial values, which are defined from gamma and factorial functions. (ii) The propagation of round-off error when computing high-order HPs, which leads to the loss of the orthogonality property of these polynomials. To solve the problem (i), new algorithms for HPs computation are proposed in this work. These algorithms are free from any gamma and factorial functions. To overcome the problem (ii), a new numerical stability method is proposed. This method that is based on appropriate stability conditions, is very low in terms of computational cost compared to the GSOP method that is presented in [9]. As result, the computation of high-order Hahn moments (HMs) becomes stable and fast when reconstructing large-size signals and images. It is worth mentioning that the high-order moments are required for better describe the content of the images [6,38].
The present work also discusses some computational aspects related to Hahn's moment invariants (HMIs) computation, with a view to overcoming the encountered problems.
Generally, the derivation process of the orthogonal moment invariants relies on the algebraic relation between the orthogonal moments and the geometric ones, to express the orthogonal moment invariants as a linear combination of geometric moment invariants. [3,14,18,20,30]. However, this process is limited by two major drawbacks: (i) the numerical instability due to gamma and factorial functions involved in the computation of the orthogonal moment invariants, and (ii) the high computational cost [3]. To overcome these problems when computing Racah moment invariants (RMIs), a fast and accurate computation method is proposed in [3]. This method is based on the use of recursive formulas of certain factors involved in RMIs computation. However, this method and others presented in [3,17,20,39] are not suitable for computing the invariant moments of large-size images (≥1024 × 1024), since these methods involve the use of certain terms depending on gamma and factorial functions when computing the image moment invariants. Motivated by overcoming these problems, a fast and numerically stable computation of large-size image moment invariants is proposed in this paper. For this purpose, significant improvements are made on the conventional computation of HMIs presented in [36]. On one hand, the proposed improvements allow to compute HMIs independently of gamma and factorial functions, which guarantees the numerical stability of the computed moment invariants. On the other hand, HMIs execution time is considerably reduced by using matrix multiplications and recursive formulas. It is important to note that the fast computation of moment invariants is relevant for real-time applications [3].
The rest of this work is organized as follows: in the second section, the conventional and the proposed methods for HPs computation are presented. The third section contains the theoretical background of Hahn moments (HMs) in the 1D and 2D domains. The fourth section presents the conventional and the proposed method for HMIs computation. In the fifth section, the results of simulations and comparisons are presented to validate the efficiency of the proposed algorithms, and the final section includes some conclusions and future works.

Conventional and proposed computation of HPs
This section presents the conventional computation methods of HPs and introduces the proposed algorithms for fast and stable HPs computation.

Conventional computation of HPs
The HPs of order n and discrete variable x are defined as follows [46]: where n, x = 0, 1, …, N − 1, α > − 1, β > − 1, and 3 F 2 (.) is a hyper-geometric function defined as: Using Eqs. (1) and (2), HPs can be explicit as follows: According to [7], 〈x〉 k is defined as follows: withS(k, i) are the Stirling numbers of the first kind, obtained by the following recurrence relation [40]: HPs satisfies the following orthogonality relation [23]: where δ nm denotes the delta of Kronecker, ω(x) is the weight function, and d 2 n represents the square norm function. These functions are defined as follows [29]: When the order n of HPs is high, the computation of these polynomials becomes unstable and expansive in terms of computation time because of gamma and factorial functions that define HPs (Eq. (1)) [9,27]. To improve the computation speed and the numerical stability of high-order HPs, researchers use the recurrence relations with respect to the order n or with respect to the variable x [46,50]. In order to derive the recurrence relations of the normalized HPs, the following form is used [46]: The following three-term recurrence relation allows the computation of the normalized with respect to the order n [50]: Where the coefficients A, B, C, D and E are defined in Table 1.
The initial values of HPs with respect to the order n are computed as follows [50]: The recursive computation of HPs is also possible with respect to the variable x by using the following three-term recurrence formula [3]: The initial values of HPs with respect to the variable x are given by [50]: The recursive computation method of HPs provided in [46,50] is limited by the problem of numerical instability of HPs initial values, either with respect to the order n or with respect to the variable x. This problem comes from ω(x) and ρ(n) functions, which are defined from gamma and factorial functions. Noting that Γ(172) = Inf under Matlab, Python, C++, etc. This Overflow problem leads to numerical instabilities of HPs when the values of the order n, the variable x and the parameters α and β become high. Therefore, the analysis of large-size signals and images by HPs is not possible.
To overcome the instability problem of high-order HPs, new algorithms are introduced in the next section.

HPs computation with respect to the order n
To avoid the numerical instability of the initial values e h α;β ð Þ 0 x ð Þ, Eq. (13) is used to develop the following formula: To eliminate Γ(x) function from the term e h α;β ð Þ 0 0 ð Þ (Eq. (18)), the following recursive formula is used: Based on Eqs. (13)- (14), the expression of e h α;β ð Þ 1 x ð Þ can be expressed as follows: The method thus developed to compute HPs initial values via Eqs. (18)-(20) is numerically stable since it is independent from gamma and factorial functions. Algorithm 1 summarizes the main steps of the proposed method for HPs computation with respect to the order n.
Algorithm 1 Proposed algorithm for HPs computation with respect to the order n.
The following subsection includes the proposed computation of HPs with respect to the variable x.

HPs computation with respect to the variable x
To ensure the numerical stability of HPs initial values e h α;β ð Þ n 0 ð Þ and e h α;β ð Þ n 1 ð Þ, Eqs. (16)- (17) are extended to obtain the formulas below.
For x > 1, the numerical stability of e h α;β ð Þ n x ð Þ is achieved by eliminating ω(x)function from Eq. (15) as follows: The fundamental steps used for HPs computation with respect to the variable x are presented in algorithm 2. Algorithm 2 Proposed algorithm for HPs computation with respect to the variable x.
It is important to mention that the computation of HPs initial values by Eqs. (20)-(21) presents the Underflow problem for high values of the order n. Indeed, Table 2 illustrates an example of this problem.
From Table 2 it can be seen that the values of e h α;β ð Þ n 0 ð Þ become zero (Underflow) when n > 2558. This problem limits the recursive computation of high-order HPs. To overcome this problem, new method is proposed (Fig. 1) for stable computation of high-order HPs with respect to the variable x. The details of this method are described in algorithm 3. Noting the proposed method is inspired from the method presented in [1], which is used for stable computation of high-order Tchebichef polynomials. Algorithm 3 Proposed Algorithm for HPs computation with respect to the variable x.
Since the computation of HPs using the proposed algorithms is recursive, this leads to the appearance and propagation of rounding error. This problem is discussed by Camacho-bello et al. in [5] when computing high-order Tchebichef polynomials. To overcome the encountered problem for Tchebichef, Meixner and Hahn polynomials, the modified Gram-Schmidt orthonormalization process (GSOP) is used [5,8,9,12]. However, the correction of numerical instability by GSOP is very time-consuming because GSOP involves a complex calculation process [5,12]. This problem limits the fast reconstruction of large-size signals and images by high-order HPs. Indeed, Table 3 shows a comparison in terms of execution time in second (s) by using algorithm 1, and GSOP method presented in [9]. The computation of HPs by these methods is performed 100 times, and then the average execution time is computed and illustrated in Table 3. It is worth mentioning that all paper tests are carried out under Matlab R2019a installed on PC of 2.4 GHz processor and 4 GB of RAM.
The results presented in Table 3 and 4 clearly indicate that the GSOP method significantly increases the computation time of HPs compared to the proposed methods (algorithms 1 and 3). This is why a new numerical stability method is introduced in the following section.

Proposed method for ensuring the numerical stability of higher-order HPs
The proposed method is based on a numerical analysis of the polynomial values in HPs matrix. Indeed, Fig. 2 present 2D plot of HPs matrix computed by algorithm 1 for n, x = 0 : 512, α = 5, β = 180. This matrix is displayed in two parts: the first one that is shown in Fig. 2 Fig. 2 (c) shows HPs matrix computed by GSOP method [9]. From Fig.1 (b), it is noticeable that the unstable values of HPs appear when n → N (high polynomial order). These values are positioned in the blue colored regions called "instability regions". After using GSOP method (Fig. 2 (c)), one can observe that the previously observed instable values become extremely low and tend towards zero. Indeed, Fig.3 shows the intervals of HPs values computed by GSOP method when n → N. Fig. 1 Steps of the proposed method for HPs computation with respect to the variable x . According to the conducted numerical analysis, appropriate "stability conditions" are introduced during the recursive computation of HPs to ensure the numerical stability of these polynomials. Indeed, the following stability condition is inserted after the three-term recurrence relation of HPs with respect to the order n.
b h Note that this condition is activated from the order n = N max /8 (lower order) where N max is the maximum order of HPs, this ensures that all unstable values are detected and eliminated (replaced by zero).
The following stability condition is implemented after the computation of each HPs value using the three-term recurrence relation with respect to the variable x: To evaluate the performance of the proposed Algorithms that including the stability conditions, HPs are computed by these Algorithms up to the order n = 8000 with x = 0 : 8000 and α = β = 50. Fig. 4 shows 2D plot of HPs matrices computed by Algorithms 1 and 3. The obtained results show that the computed HPs are stable since the polynomial values belong to the orthonormality domain [−1, 1], which justifies the efficiency of the proposed Algorithms for the stable computation of high-order HPs.
The following test presents a comparison in terms of execution time using the proposed method (algorithm 3) and the GSOP method presented in [9]. The execution time improvement ratio (ETIR (%)) defined by ETIR (%) = (1 − T 1/T 2) × 100 [23] is used as a comparison criterion, where T 1 is the execution time consumed by the proposed method, and T 2 represents the execution time consumed by GSOP method [9]. The obtained results which are displayed in Table 4 demonstrate that.
the proposed method (algorithm 3 with the stability condition) is very fast compared to the GSOP method since ETIR ≥ 99%. This is especially true when the polynomial order n is high. It is worth mentioning that the computation time of HPs is significantly reduced when using the following HPs symmetry relation defined for α = β [29]: It is also important to note that Krawtchouk (KPs) and Tchebichef (TPs) polynomials are special cases of HPs. Indeed, KPs ( e k n p ð Þ x ð Þ) are obtained from HPs when α = pt and β = (1 − p) t with t → ∞ [29]. In the case α = β = 0, HPs are reduced to TPs ( e t n x ð Þ) [22]. In the following sections, HPs will be used as kernel polynomials to define both Hahn moments (HMs) and Hahn moment invariants (HMIs) designed for large-size signal and image analysis.

Hahn moments
This section briefly presents the theoretical background of Hahn moments (HMs) in the 1D and 2D domains.
The 1D HMs of order n are defined for a 1D-signal function f(x) by the following relation: The reconstructed signal b f x ð Þ is computed from the 1D HMs as follows: The reconstruction error can be measured by the Mean Square Error (MSE) criterion, which is defined as follows: where f(x) is the original signal, and b f x ð Þ represents the reconstructed one. For a digital image f(x, y) of size N × M, the 2D HMs are defined by [50]: with n, m = 0, 1, 2, …, N − 1, x ∈ [0…N − 1], and y ∈ [0…M − 1],.
The computation time of HMs is considerably reduced by using the following matrix formula [3]: where The matrix formulation of the inverse 2D HMs is given by: is the reconstructed image of size N × M.
In order to quantify the reconstruction error by HMs, the MSE criterion is used as follows: The Peak Signal-to-Noise Ratio (PSNR) criterion in decibels (dB) is also used to quantify the reconstruction error. This criterion is defined from MSE as: where peak is the maximum value of the 2D image function.

Hahn moment invariants
This section presents both the conventional and the proposed method for deriving Hahn moment invariants (HMIs) from geometric moment invariants (GMIs).

The conventional computation of HMIs
The following relation is used to compute the geometric moments (GMs) of f (x, y) image of size N × M [14]: x n y m f x; y ð Þ ð36Þ The central geometric moments denoted μ nm are derived from GMs by the following formula [35]: The set of geometric moment invariants (GMIs) of order (n, m), denoted V nm , which are independent of translation, scaling and rotation are defined by the following relation [35]: where γ ¼ nþm 2 þ 1 and θ ¼ 1 2 tan −1 2μ 11 μ 20 −μ 02 .
In accordance with Eq. (39), the values of the angle θ are limited to −π/4 ≤ θ ≤ π/4. For angle values in different ranges, the article [40] is recommended.
HMs given by Eq. (31) can be expressed in terms of HPs (Eq. (3)) as follows [36]: (30) and (39), HMs can be expressed in terms of GMs as follows [36]: The substitution of GM nm by V nm allows to derive HMIs as follows [36]: Algorithm 4 summarizes the conventional steps of the HMIs derivation process. Algorithm 4 The conventional method for HMIs computation [36].

The conventional computation of HMIs via algorithm 4 is numerically instable because of
n functions as they are defined in terms of gamma and factorial functions. In addition, this conventional computation is time-consuming due to the use of several discrete summations. These problems limit the computation of HMIs especially for large-size images. There is therefore a clear need for a considerable improvement in HMIs computation in terms of rapidity and numerical stability.

Proposed computation of HMIs
In this section, the proposed method for fast and stable computation of HMIs is presented. This method is mainly based on the use of matrix and recursive formulas.
In order to significantly accelerate the computation of GMs, the following matrix formulation can be used: The computation time of μ nm can also be reduced through the use of the following formula: In a similar way, the computation of V nm is speed up by using the following formulation: y),x = 1 : N, y = 1 : M, x ¼ GM 10 GM 00 , y ¼ GM 01 GM 00 and the symbol ⊙ denotes the Hadamard product [4].
The following recursive relations are developed to accelerate and stabilize the computation of HMIs. Indeed, ω(x) is computed as follows: where logΓ(x) is the logarithmic gamma function [15]. This function is numerically stable for high values of the argument x. For example, logΓ(10 5 ) = 1.0513e + 06 under Matlab. The stable recursive computation of d 2 n function is given by: The fast and stable computation of A β ð Þ n and B α;β ð Þ nk functions is defined as follows The summary of the proposed method for fast and stable computation of HMIs is presented in algorithm 5.
Algorithm 5 The proposed method for fast and stable HMIs computation.

Experimental results and discussion
This section presents the simulation outcomes and the performed comparisons that validate the efficiency of the proposed algorithms in the field of large-size signal and image analysis. Initially, the numerical stability of the proposed algorithms for higher-order HPs computation is demonstrated through the reconstruction of large-size signals and images up to the last moment order. Then, the speed and numerical stability of the proposed method for HMIs computation is proved via the extraction of HMIs from large-size medical images.

Large-size signal reconstruction by higher-order HMs
In the first test, the efficiency of the proposed computation of high-order HPs with respect to the order n is demonstrated. For this purpose, HPs are computed by algorithm 1 for n, x = 0 : 8000 with α = β = 100. Noting that the parameter values are selected empirically following the study carried out in [9]. This study shows that the parameters selection case α = β allows obtaining a minimal reconstruction error (MSE) in comparison to other cases (α ≠ β). Then, an "ECG" signal of size N = 8000 is selected from the Fig. 6 Set of bio-signals selected from the database [25] that are used in the test database [25], and reconstructed by HMs up to the last moment order (n = 8000). Figure 5 shows a set of reconstructed "ECG" signals with the corresponding MSE values. On one hand, the obtained results indicate that the MSE values decrease proportionally with the increase in HMs order until the last moment order. On the other hand, it is noted that the lower MSE value is reached at the last HMs order. These results prove the numerical stability of the proposed algorithm 1, and confirm the importance of high-order HMs in terms of precise signal reconstruction.
To further validate the efficiency of algorithm 1, five bio-signals of size N = 8000 (Fig. 6) are reconstructed by HMs (with HPs are computed by algorithm 1 for α = β = 100). For comparison purposes, the test signal are also reconstructed by Tchebichef moments (TMs) [5]. Next, the average PSNR values corresponding to the reconstructed signals by both moment types are displayed in Fig. 7.
The PSNR values shown in Fig. 7 increase up to the last moment order, reflecting the efficiency and numerical stability of both TMs and HMs for the reconstruction of large-size signals. One should mention that the numerical stability of TMs is ensured via GSOP [5]. Whereas HMs numerical stability is provided by using the proposed method (stability condition). This makes HMs significantly faster to compute compared to TMs.

Large-size image reconstruction by higher-order HMs
In this section, the performance of the proposed method for the computation of high-order HPs with respect to the variable x (algorithm 3) is illustrated. In the first test, the "PIA21906" test image of size 4000 × 4000 is selected from the database [31], and then reconstructed by HMs up to the last order (4000, 4000). Noting that the used HPs parameters values in all reconstruction tests of this section are α = β = 50. Figure 8 shows a set of the reconstructed "PIA21906" images with the corresponding MSE and PSNR values. From this figure, it is noticeable that the visual similarity between the reconstructed image and the original one gradually improves until the last moment order (4000, 4000). It can also be noted that the best quality of the reconstructed image (lower MSE) is obtained at the last moment order. This is a clear indication on the efficiency of the proposed algorithm 3 for the stable reconstruction of large-sized images.
The second test presents a comparison in terms of reconstruction errors (MSE and PSNR) by using the moments of Tchebichef (TMs) [5], Krawtchouk (KMs) [47], Meixner (MMs) Fig. 9 Set of reconstructed "Lena" image by different moment types of order up to (1024, 1024) [12], Charlier (CMs) [12] and HMs (based on algorithm 3). Indeed, Fig. 9 shows the reconstructed "Lena" image by the aforementioned moments, and Fig. 11 shows the average MSE and PSNR values corresponding to the reconstructed test images shown in Fig. 10. The achieved results show on one hand that all the used moment types are stable because they allow the reconstruction of the test images until the last moment order (1024, 1024). On the other hand, it is apparent that the quality of the reconstructed images using high-order KMs is better (high PSNR) compared to other moments type, including HMs. However, KMs are unstable when the moment orders exceed (1075, 1075) [8]. For this reason, only images of size 1024 × 1024 are used in this test. Moreover, HMs are more general than TMs and KMs because the latter are special cases of HMs [22].
In the following test, the robustness of HMs against different types of noise is evaluated. For this purpose, a real medical image named "Boxer-fracture" of size 2048 × 2048 is selected from the database [32]. This image is contaminated by two noise types namely "Salt & pepper" and "Speckle" noise with density of 5% (Fig. 12). Then, these images are reconstructed by HMs, and the reconstruction error (PSNR) is displayed in Fig. 13. The achieved results show a slight decrease in the quality of the noisy reconstructed images compared to the quality of the noise-free reconstructed one, which validates the robustness of HMs in noisy conditions.

Fast and numerically stable computation of HMIs using large-size images
This section is preserved to demonstrate the effectiveness of the proposed method for fast and stable HMIs computation (algorithm 5).
The computation of HMIs via the proposed method depends on HPs local parameters (α and β). Therefore, the optimal selection of these parameters is crucial when computing HMIs. "Lena" "Pepper" "Fruits" "Baboon" "sailboat" In order to carry out this selection, one can use an empirical method, which is similar to the one described in [3]. Indeed, HMIs of a test image are computed for specific values of α and βparameters, then the test image is geometrically transformed by translation, rotation and scaling. Next, the relative error (RE) (Eq. (49)) is computed for each parameters case. Finally, the values of α and β that allow to obtain the lower RE values are therefore considered as the optimal parameters selection.
The RE criterion is defined as: where ‖.‖ denotes the Euclidean norm, Vdenotes the invariant moments of the original image, and V * represents the invariant moments of the geometrically transformed image.
To test the effectiveness of the empirical parameters selection method, one use a real X-ray medical image named "Covid 19" of size 2048 × 2048 (Fig. 14) that is selected from [32]. The test image is then translated by a vector varies from (−20, −20) to (+20, +20) with a step of (+4, +4), scaled by a factor scaling ranging from 0.5 to 2 with a step of 0.1, and rotated by a rotation angles ranging from 0°to 360°with a step of 10°. Then, HMIs of the original image and the geometrically transformed one are computed up to the order (5, 5) for four cases of α and βparameters: (A) α = 15 and β = 10, (B) α = β = 10, (C) and α = β = 5, and (D) α = β = 0.
The results of the performed test are presented in Fig. 15. From this figure, it can be noted that the cse α = β = 0 leads to lower RE values, therefore this case will be considered as the optimal parameters selection in the rest of this work. However, the use of the empirical method for the optimal parameters selection is limited due to the few number of the studied cases. It is important to note that the optimal selection of local parameters can be made via population-based algorithms such as the Sine Cosine Algorithm (SCA) [10] or by similar algorithms (PSO, ABC, TLBO…). In the following test, four grayscale medical images of size 2048 × 2048 are used (Fig. 16). These images that are selected from the database [32], are geometrically transformed, and then the corresponding average RE is displayed in Fig. 17. The outcome results indicate on one hand that the proposed computation of HMIs is stable for large-size images. On the other hand, the low RE values (RE ≃ 10 −10 ) indicate that the computed HMIs are highly accurate and robust against different geometrical transformations.
The following test is conducted to prove the fastness of the proposed HMIs computation method compared to the conventional one that is presented in [36]. Initially, the test image  Table 4.
The obtained results in Table 4 clearly demonstrate the considerable rapidity of the proposed method compared to the conventional one, especially for higher HMIs order.
Due to the limited paper length, only the advantages of the proposed method for HMIs computation have been illustrated. In future works, the proposed method (algorithm 5) can be successfully employed in various applications related images moments invariants such as pattern recognition, image classification and image zerowatermarking of large-size images.
It is important to point out that some recent methods for computing the invariant moments [3,17,20,37,39] are numerically unstable for large-size images (>1024 × 1024).For this reason, no comparison with these works is carried out. In addition, accelerating and stabilizing the computation of the previously mentioned methods might be a potentially interesting direction for further research.

Conclusion
In this paper, new algorithms for the computation of high order HMs and HMIs are presented in the field of large-sized signal and image analysis. First, special attention is paid to ensuring the numerical stability and to speed up the computation of HPs with respect to the order n and with respect to the variable x. The proposed algorithms are fully independent from gamma and factorial functions that are the source of numerical instability (Underflow and Overflow). Secondly, a new numerical stability method is proposed to detect and eliminate unstable coefficients during recursive computation of HPs. This method is based on the distribution analysis of HPs values in the polynomial matrix. As a result, the computation of HPs becomes stable and fast for any order n, any value of the variable x and for any value of the local parameters. Also new algorithm is proposed for fast and stable HMIs computation. This algorithm is based on developed recursive relations and matrix formulations. The simulation and comparison result clearly validated the stability and efficiency of the proposed methods for the large-sized signal and image analysis. However, the empirical method used for the selection of HPs parameters is limited due to the few studied cases and therefore needs to be further improved. In perspective and future work, we will focus on other applications of the proposed algorithms such as compression, classification and zero -watermarking of large-size signals and 2D/3D images. In addition, the selection of HPs parameters will be conducted based on metaheuristic algorithms.