## 2 Numerical Approximation of RLFI
The GBFA interpolant of \(f\) is given by  
\[I_{n}f(t)=f_{0,n}^{\top}\hat{G}_{0,n}^{\lambda}\{t\}, \tag{7}\]where \(\mathcal{L}_{k}^{\lambda}(t)\) is defined as  
\[\mathcal{L}_{k}^{\lambda}(t)=\hat{\omega}_{k}^{\lambda}\text{trp}\Big{(}\hat{ \lambda}_{0:n}^{\lambda\ \ \ +}\Big{)}\Big{(}\hat{G}_{0:n}^{\lambda\ \ \ +}\Big{[}\hat{t}_{n,k}^{\lambda}\Big{]}\odot\hat{G}_{0:n}^{\lambda}\{t_{1}^{ \lambda}\}\Big{)},\quad\forall k\in\mathbb{J}_{n}^{+}, \tag{8}\]  
with normalization factors and Christoffel numbers:  
\[\hat{\lambda}_{j}^{\lambda}=\frac{\pi 2^{1-4\lambda}\Gamma(j+2 \lambda)}{j!\Gamma^{2}(\lambda)(j+\lambda)},\] (9) \[\hat{\omega}_{k}^{\lambda}=1/\bigg{[}\text{trp}\Big{(}\hat{\lambda }_{0:n  
We term \(\,{}^{E}\mathbf{Q}_{n}^{a}\) the \"\(\star\)th-order FSGIM\" for the RLFI and \(\,{}^{E}\mathbf{Q}_{n}^{a}\) the \"\(\star\)th-order FSGIM Generator.\" Equation (17) is preferred for computational efficiency. For a detailed dimensional analysis of the matrix expressions in Equation (15), including the explicit form of \(\mathcal{I}_{1}^{(y)}\hat{G}_{0:n}^{\lambda}[t(1-y^{1/\alpha})\mathbf{I}_{n+1}]\) and clarification of the dimensions and reshaping operation in Equation (16), see Appendix C.  
To compute \(\mathcal{I}_{1}^{(y)}\hat{G}_{j}^{\lambda}[t(1-y^{1/\alpha})]\), we can use the SGIRV \(\hat{\mathbf{P}}=\frac{1}{2}\mathbf{P}\) with SGG nodes \(\hat{\mathbf{I}}_{n_{q}}^{\lambda_{q}}\):  
\[\mathcal{I}_{1}^{(y)}\hat{G}_{j}^{\lambda}[t(1-y^{1/\alpha})]\approx\hat{ \mathbf{P}}\hat{G}_{j}^{\lambda}\left(t\left(1-\hat{\mathbf{I}}_{n_{q}}^{ \lambda_{q}^{q}^{q}/\lambda}\right)\right),\quad\forall j\in\mathbb{J}_{n}^{ \perp},\,t\in\mathbf{\Omega}_{1}, \tag{21}\]  
cf. (35) (Algorithm 6 or 7). Formula (21) represents the \((n_{q},\lambda_{q})\)-GBFA quadrature used for the numerical partial calculation of the RLFI. We denote the approximate \(\star\)th-order RLFI of a function at point \(t\), computed using Equation (21) in conjunction with either Equation (17) or Equation (20), by \(\,{}^{n,\lambda,E}_{n_{q},\lambda_{q}}\gamma_{t}^{a}\) (The parameters in the left superscripts and subscripts of \(\,{}^{n,\lambda,E}_{n_{q},\lambda_{q}}\gamma_{t}^{a}\) are the same parameters used by the numerical method: \(n\) and \(\lambda\) indicate the degree and index of the SG polynomial \(\hat{G}_{n}^{\lambda}(x)\), while \(n_{q}\) and \(\lambda_{q}\) specify the parameters of the polynomial basis used in the quadrature scheme. \(E\) is a letter distinguishing this discrete operator from other operators in the literature, referring to the first initial of the author's surname). Figure 1 provides a visual summary of the GBFA method's workflow, illustrating the seamless integration of interpolation, transformation, and quadrature steps. This schematic highlights the method's flexibility, as the tunable parameters \(\lambda\) and \(\lambda_{q}\) allow practitioners to tailor the approximation to specific problem characteristics, optimizing both accuracy and computational efficiency. The alternative path of precomputing the FSGIM, as indicated by the dashed arrow, underscores the method's suitability for applications requiring repeated evaluations. **Remark 1**. : _The selection of \(f\in L^{2}(\Omega_{1})\) for defining the RLFI is motivated by both mathematical rigor and practical relevance. The kernel \((t-\tau)^{\alpha-1}\) exhibits a singularity at \(\tau=t\), necessitating integrability of the product for the RLFI to be well-defined. The space \(L^{2}(\Omega_{1})\), characterized by \(\int_{0}^{1}|f(t)|^{2}\,dt<\infty\), ensures sufficient regularity to guarantee convergence. Since the singularity is integrable for \(\alpha>0\), the RLFI remains finite when \(f\in L^{2}(\Omega_{1})\). Moreover, \(L^{2}\)-integrability is consistent with many physical models in viscoelasticity, anomalous diffusion, and control theory, making it a natural function space for modeling.