fnft__misc.h

Go to the documentation of this file.

/*
* This file is part of FNFT.
*
* FNFT is free software; you can redistribute it and/or
* modify it under the terms of the version 2 of the GNU General
* Public License as published by the Free Software Foundation.
*
* FNFT is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Contributors:
* Sander Wahls (TU Delft) 2017-2018, 2023.
* Shrinivas Chimmalgi (TU Delft) 2019-2020.
* Peter J. Prins (TU Delft) 2021.
*/
 
#ifndef FNFT__MISC_H
#define FNFT__MISC_H
 
#include "fnft.h"
#include <string.h>
 
void fnft__misc_print_buf(const FNFT_UINT len, FNFT_COMPLEX const * const buf,
                          char const * const varname);
 
FNFT_REAL fnft__misc_rel_err(const FNFT_UINT len,
    FNFT_COMPLEX const * const vec_numer, FNFT_COMPLEX const * const vec_exact);
 
FNFT_REAL fnft__misc_hausdorff_dist(const FNFT_UINT lenA,
    FNFT_COMPLEX const * const vecA, const FNFT_UINT lenB,
    FNFT_COMPLEX const * const vecB);
 
FNFT_COMPLEX fnft__misc_sech(FNFT_COMPLEX Z);
 
FNFT_REAL fnft__misc_l2norm2(const FNFT_UINT N, FNFT_COMPLEX const * const Z,
    const FNFT_REAL a, const FNFT_REAL b);
 
FNFT_INT fnft__misc_filter(FNFT_UINT * const N, FNFT_COMPLEX * const vals,
    FNFT_COMPLEX * const rearrange_as_well,
    FNFT_REAL const * const bounding_box);
 
 
FNFT_INT fnft__misc_filter_nonreal(FNFT_UINT *N_ptr, FNFT_COMPLEX * const vals,
    const FNFT_REAL tol_im);
 
FNFT_INT fnft__misc_filter_inv(FNFT_UINT * const N_ptr, FNFT_COMPLEX * const vals,
    FNFT_COMPLEX * const rearrange_as_well,
    FNFT_REAL const * const bounding_box);
 
FNFT_INT fnft__misc_merge(FNFT_UINT *N_ptr, FNFT_COMPLEX * const vals,
    FNFT_REAL tol);
 
FNFT_INT fnft__misc_downsample(const FNFT_UINT D, FNFT_COMPLEX const * const q,
    FNFT_UINT * const Dsub_ptr, FNFT_COMPLEX ** qsub_ptr,
    FNFT_UINT * const first_last_index);
 
static inline FNFT_COMPLEX fnft__misc_CSINC(FNFT_COMPLEX x)
{
    const FNFT_REAL sinc_th=1.0E-8;
 
    if (FNFT_CABS(x)>=sinc_th)
        return FNFT_CSIN(x)/x;
    else
        return FNFT_CCOS(x/FNFT_CSQRT(3));
}
 
static inline FNFT_COMPLEX fnft__misc_CSINC_derivative(FNFT_COMPLEX x)
{
    const FNFT_REAL sinc_derivative_th=1.0;
    FNFT_COMPLEX returnvalue;
 
    if (FNFT_CABS(x)>=sinc_derivative_th)
        returnvalue = (FNFT_CCOS(x)-fnft__misc_CSINC(x))/x;
    else {
        returnvalue = 0.0;
        FNFT_COMPLEX x2 = x * x;
        for (FNFT_UINT n=9; n-->0; )
            returnvalue = (( n%2 ? 1.0 : -1.0) + x2/(2*n+2) * returnvalue ) / (2*n+3);
        returnvalue *= x;
    }
    return returnvalue;
}
 
FNFT_UINT fnft__misc_nextpowerof2(const FNFT_UINT number);
 
FNFT_INT fnft__misc_resample(const FNFT_UINT D, const FNFT_REAL eps_t, FNFT_COMPLEX const * const q,
    const FNFT_REAL delta, FNFT_COMPLEX *const q_new);
 
static inline void fnft__misc_matrix_mult(const FNFT_UINT n,
                                          const FNFT_UINT m,
                                          const FNFT_UINT p,
                                          FNFT_COMPLEX const * const A,
                                          FNFT_COMPLEX const * const B,
                                          FNFT_COMPLEX * const C){
    for (FNFT_UINT cn = 0; cn < n; cn++) {
        for (FNFT_UINT cp = 0; cp < p; cp++) {
            C[ cn*p + cp] = 0.0;
            for (FNFT_UINT cm = 0; cm < m; cm++)
                C[ cn*p + cp ] += A[ cn*m + cm ] * B[ cm*p + cp ];
        }
    }
    
    return;
}
 
static inline void fnft__misc_mat_mult_2x2(FNFT_COMPLEX * const U,
        FNFT_COMPLEX *const T){
 
    FNFT_COMPLEX TM[4] = { 0 };
    fnft__misc_matrix_mult(2, 2, 2, U, T, &TM[0]);
    memcpy(T,TM,4*sizeof(FNFT_COMPLEX));
    
    return;
}
 
static inline void fnft__misc_mat_mult_4x4(FNFT_COMPLEX * const U,
        FNFT_COMPLEX *const T){
 
    FNFT_COMPLEX TM[16] = { 0 };
    fnft__misc_matrix_mult(4, 4, 4, U, T, &TM[0]);
    memcpy(T,TM,16*sizeof(FNFT_COMPLEX));
    
    return;
}
 
static inline FNFT_REAL fnft__misc_legendre_poly(const FNFT_UINT n, const FNFT_REAL x){
    FNFT_UINT  i;
    FNFT_REAL P, P_1, P_2;
    if (n == 0)
        P = 1;
    else if (n == 1)
        P = x;
    else{
        P_1 = x;
        P_2 = 1;
        for (i = 2; i <= n; i++) {
            P = (2.0*i-1)*x*P_1/i -(i-1.0)*P_2/i;
            P_2 = P_1;
            P_1 = P;
        }
    }
    return P;
}
 
static inline FNFT_INT fnft__misc_normalize_vector(const FNFT_UINT len, FNFT_COMPLEX * const v)
{
    // Find max of absolute values of coefficients
    FNFT_REAL max_abs = 0.0;
    for (FNFT_UINT i=0; i<len; i++) {
        const FNFT_REAL cur_abs = FNFT_CABS( v[i] );
        if (cur_abs > max_abs)
            max_abs = cur_abs;
    }
 
    // Return if polynomial is identical to zero
    if (max_abs == 0.0)
        return 0;
 
    // Otherwise, rescale
    const FNFT_REAL a = FNFT_FLOOR( FNFT_LOG2(max_abs) );
    const FNFT_REAL scl = FNFT_POW( 2, -a );
    for (FNFT_UINT i=0; i<len; i++)
        v[i] *= scl;
 
    return (FNFT_INT) a;
}
 
#ifdef FNFT_ENABLE_SHORT_NAMES
#define misc_print_buf(...) fnft__misc_print_buf(__VA_ARGS__)
#define misc_rel_err(...) fnft__misc_rel_err(__VA_ARGS__)
#define misc_hausdorff_dist(...) fnft__misc_hausdorff_dist(__VA_ARGS__)
#define misc_sech(...) fnft__misc_sech(__VA_ARGS__)
#define misc_l2norm2(...) fnft__misc_l2norm2(__VA_ARGS__)
#define misc_filter(...) fnft__misc_filter(__VA_ARGS__)
#define misc_filter_inv(...) fnft__misc_filter_inv(__VA_ARGS__)
#define misc_filter_nonreal(...) fnft__misc_filter_nonreal(__VA_ARGS__)
#define misc_merge(...) fnft__misc_merge(__VA_ARGS__)
#define misc_downsample(...) fnft__misc_downsample(__VA_ARGS__)
#define misc_CSINC(...) fnft__misc_CSINC(__VA_ARGS__)
#define misc_CSINC_derivative(...) fnft__misc_CSINC_derivative(__VA_ARGS__)
#define misc_nextpowerof2(...) fnft__misc_nextpowerof2(__VA_ARGS__)
#define misc_resample(...) fnft__misc_resample(__VA_ARGS__)
#define misc_matrix_mult(...) fnft__misc_matrix_mult(__VA_ARGS__)
#define misc_mat_mult_2x2(...) fnft__misc_mat_mult_2x2(__VA_ARGS__)
#define misc_mat_mult_4x4(...) fnft__misc_mat_mult_4x4(__VA_ARGS__)
#define misc_legendre_poly(...) fnft__misc_legendre_poly(__VA_ARGS__)
#define misc_normalize_vector(...) fnft__misc_normalize_vector(__VA_ARGS__)
#endif
 
#endif