1. add conversion from Taylor(polynomial) expansing to Pade approxim…

…ation
OstapHEP · Dec 13, 2024 · 75c3c22 · 75c3c22
1 parent 118c238
commit 75c3c22
Show file tree

Hide file tree

Showing 5 changed files with 269 additions and 49 deletions.
diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md
@@ -2,7 +2,8 @@
 
     1. add several functions, like `sec` , `csc` , `cot` , `cas`
     1. add Bring radical
-    1. add Pade intrpolation     
+    1. add Pade interpolation     
+    1. add conversion from Taylor(polynomial) expansing to Pade approximation
 
 ## Backward incompatible
 

diff --git a/source/include/Ostap/Rational.h b/source/include/Ostap/Rational.h
@@ -457,6 +457,9 @@ namespace Ostap
       // ======================================================================
     } ;
     // ========================================================================
+    /// forward declaration 
+    class Polynomial ; // forward declatation 
+    // ========================================================================
     /** @class Pade 
      *  Pade-like rational function with the optional setting 
      *  of constituent shape-fixing zeroes and poles 
@@ -587,6 +590,36 @@ namespace Ostap
 	const std::vector<std::complex<double> >& czeroes = std::vector<std::complex<double> > () ,
 	const std::vector<std::complex<double> >& cpoles  = std::vector<std::complex<double> > () ) ;      
       // =====================================================================
+      /** constructor from the polynomial/Taylor expansion 
+       *  @param p polynomial expansion 
+       *  @param n degree of P(x) 
+       *  @param m degree of Q(x) 
+       */
+      Pade
+      ( const Ostap::Math::Polynomial& p ,
+	const unsigned short           n , 
+	const unsigned short           m ) ;
+      // =====================================================================
+      /** constructor from the polynomial/Taylor expansion 
+       *  @param p polynomial expansion 
+       *  @param n degree of P(x) 
+       */
+      Pade
+      ( const Ostap::Math::Polynomial& p ,
+	const unsigned short           n ) ; 
+      // =====================================================================
+      /** constructor from the polynomial/Taylor expansion 
+       *  @param n degree of P(x) 
+       *  @param m degree of Q(x) 
+       *  @param p polynomial expansion 
+       */
+      Pade
+      ( const unsigned short       n         , 
+	const unsigned short       m         , 
+	const std::vector<double>& p         ,
+	const double               xmin = -1 ,
+	const double               xmax =  1 ) ;	
+      // =====================================================================
    public:
       // ======================================================================
       /// get x-min 
@@ -734,14 +767,10 @@ namespace Ostap
       // ======================================================================      
     }
     // ========================================================================
-
 
-    //
-    void pade_intepolation
-    ( const std::vector<double>& x , 
-      const std::vector<double>& y ,
-      const unsigned short       n ) ;
-    //
+    // Ostap::Math::Pade
+    void pade_ququ ( const std::vector<double>& f ,
+		     const unsigned short       n ) ;
 
 
   } //                                         The end of namespace Ostap::Math 

diff --git a/source/src/GSL_helpers.cpp b/source/src/GSL_helpers.cpp
@@ -1,6 +1,10 @@
 // ============================================================================
 // Include files 
 // ============================================================================
+// Ostap
+// ============================================================================
+#include "Ostap/GSL_utils.h"
+// ============================================================================
 // GSL
 // ============================================================================
 #include "gsl/gsl_linalg.h"
@@ -161,9 +165,22 @@ Ostap::GSL_Permutation::~GSL_Permutation ()
     m_permutation = nullptr ; 
   }
 }
-
-
-
+// ============================================================================
+/// print operator 
+std::ostream& operator<<( std::ostream&            s ,
+			  const Ostap::GSL_Matrix& m )
+{
+  if ( m.matrix() == nullptr ) { s << "GSL_Matrix{nullptr}" ; return s ; }
+  return Ostap::Utils::toStream ( *m.matrix() , s ) ;
+}
+// ============================================================================
+/// print operator 
+std::ostream& operator<<( std::ostream&            s ,
+			  const Ostap::GSL_Vector& v ) 
+{
+  if ( v.vector() == nullptr ) { s << "GSL_Vector{nullptr}" ; return s ; }
+  return Ostap::Utils::toStream ( *v.vector() , s ) ;
+}
 // ============================================================================
 //                                                                      The END 
 // ============================================================================
diff --git a/source/src/GSL_helpers.h b/source/src/GSL_helpers.h
@@ -4,6 +4,10 @@
 // ============================================================================
 // Include files
 // ============================================================================
+// STD/STD
+// ============================================================================
+#include <ostream>
+// ============================================================================
 // GSL 
 // ============================================================================
 #include "gsl/gsl_linalg.h"
@@ -68,13 +72,15 @@ namespace Ostap
     inline const gsl_matrix* matrix  () const { return m_matrix ; }
     // ========================================================================
     /// get the mattix element 
-    double      get    ( const unsigned short n1 , 
-                         const unsigned short n2 ) const 
+    double      get
+    ( const unsigned short n1 , 
+      const unsigned short n2 ) const 
     { return gsl_matrix_get ( m_matrix , n1 , n2 ) ; }
     /// set the matrix element 
-    void        set    ( const unsigned short n1 , 
-                         const unsigned short n2 , 
-                         const double   value    )
+    void        set
+    ( const unsigned short n1 , 
+      const unsigned short n2 , 
+      const double   value    )
     {        gsl_matrix_set ( m_matrix , n1 , n2 , value ) ; }
     // ========================================================================
   private:
@@ -125,11 +131,13 @@ namespace Ostap
     inline const gsl_vector* vector () const { return m_vector ; }
     // ========================================================================
     /// get the vector element 
-    double      get    ( const unsigned short n ) const 
+    double      get
+    ( const unsigned short n ) const 
     { return gsl_vector_get ( m_vector , n ) ; }
     /// set the vector element 
-    void        set    ( const unsigned short n , 
-                         const double   value    )
+    void        set
+    ( const unsigned short n , 
+      const double   value    )
     {        gsl_vector_set ( m_vector , n , value ) ; }
     // ========================================================================
   private:
@@ -173,7 +181,13 @@ namespace Ostap
     // ========================================================================    
   };
   // ==========================================================================
-} //                                                 The end of namespace Ostap 
+} //                                                 The end of namespace Ostap
+// ============================================================================
+/// print operator 
+std::ostream& operator<<( std::ostream& s , const Ostap::GSL_Matrix& m ) ;
+// ============================================================================
+/// print operator 
+std::ostream& operator<<( std::ostream& s , const Ostap::GSL_Vector& v ) ;
 // ============================================================================
 //                                                                      The END 
 // ============================================================================