<commit_before>/*
 * Copyright (c) 2007, 2008, 2009, 2010 <NAME>
 * Copyright (c) 2007, 2008, 2009, 2010 Virginia Polytechnic
 *                                      Institute & State University
 *
 * This file is part of liquid.
 *
 * liquid is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * liquid 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 liquid.  If not, see <http://www.gnu.org/licenses/>.
 */

//
// Elliptic filter design
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "liquid.internal.h"

#define LIQUID_DEBUG_ELLIP_PRINT   0

// forward declarations

// Landen transformation (one iteration)
float landenf(float _k);

// compute elliptic integral K(k) for _n recursions
void ellipkf(float _k,
             unsigned int _n,
             float * _K,
             float * _Kp);

// elliptic degree
float ellipdegf(float _N,
                float _k1,
                unsigned int _n);

// Elliptic cd() function (_n recursions)
float cdef(float _u,
           float _k,
           unsigned int _n);

// ***************************************************************

// Landen transformation (one iteration)
float landenf(float _k)
{
    float kp = sqrtf(1.0f - _k*_k);
    return (1.0f - kp) / (1.0f + kp);
}

// compute elliptic integral K(k) for _n recursions
void ellipkf(float _k,
             unsigned int _n,
             float * _K,
             float * _Kp)
{
    float kn = _k;
    float knp = sqrtf(1.0f - _k*_k);
    float K  = 0.5f*M_PI;
    float Kp = 0.5f*M_PI;
    unsigned int i;
    for (i=0; i<_n; i++) {
        kn  = landenf(kn);
        knp = landenf(knp);

        K  *= (1.0f + kn);
        Kp *= (1.0f + knp);
    }
    *_K  = K;
    *_Kp = Kp;
}

// elliptic degree
float ellipdegf(float _N,
                float _k1,
                unsigned int _n)
{
    // compute K1, K1p from _k1
    float K1, K1p;
    ellipkf(_k1,_n,&K1,&K1p);

    // compute q1 from K1, K1p
    float q1 = expf(-M_PI*K1p/K1);

    // compute q from q1
    float q = powf(q1,1.0f/_N);

    // expand numerator, denominator
    unsigned int m;
    float b = 0.0f;
    for (m=0; m<_n; m++)
        b += powf(q,(float)(m*(m+1)));
    float a = 0.0f;
    for (m=1; m<_n; m++)
        a += powf(q,(float)(m*m));

    float g = b / (1.0f + 2.0f*a);
    float k = 4.0f*sqrtf(q)*g*g;
    return k;
}

// Elliptic cd() function (_n recursions)
float cdef(float _u,
           float _k,
           unsigned int _n)
{
    float wn = cosf(_u*M_PI*0.5f);
    float wn_inv = 1.0f / wn;
    float kn;
    unsigned int i;
    for (i=_n; i>0; i--) {
        unsigned int j;
        kn = _k;
        for (j=0; j<i; j++)
            kn = landenf(kn);

        wn = 1.0f / wn_inv;
        wn_inv = 1.0f / (1.0f + kn) * (wn_inv + kn*wn);
    }

    return 1.0f / wn_inv;
}<commit_msg>ellip: simplifying cd(), sn() methods and inverse; needs testing<commit_after>/*
 * Copyright (c) 2007, 2008, 2009, 2010 <NAME>
 * Copyright (c) 2007, 2008, 2009, 2010 Virginia Polytechnic
 *                                      Institute & State University
 *
 * This file is part of liquid.
 *
 * liquid is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * liquid 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 liquid.  If not, see <http://www.gnu.org/licenses/>.
 */

//
// Elliptic filter design
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "liquid.internal.h"

#define LIQUID_DEBUG_ELLIP_PRINT   0

// forward declarations

// Landen transformation (one iteration)
float landenf(float _k);

// compute elliptic integral K(k) for _n recursions
void ellipkf(float _k,
             unsigned int _n,
             float * _K,
             float * _Kp);

// elliptic degree
float ellipdegf(float _N,
                float _k1,
                unsigned int _n);

// 
float complex ellip_pqf(float complex _u,
                       float complex _wM,
                       float _k,
                       unsigned int _M);

// 
float complex ellip_apqf(float complex _u,
                        float complex _w0,
                        float _k,
                        unsigned int _M);

// elliptic cd() function (_n recursions)
float complex ellip_cdf(float complex _u,
                        float complex _k,
                        unsigned int _n);

// elliptic inverse cd() function (_n recursions)
float complex ellip_acdf(float complex _u,
                         float complex _k,
                         unsigned int _n);


// elliptic sn() function (_n recursions)
float complex ellip_snf(float complex _u,
                        float complex _k,
                        unsigned int _n);

// elliptic inverse sn() function (_n recursions)
float complex ellip_asnf(float complex _u,
                         float complex _k,
                         unsigned int _n);

// ***************************************************************

// Landen transformation (one iteration)
float landenf(float _k)
{
    float kp = sqrtf(1.0f - _k*_k);
    return (1.0f - kp) / (1.0f + kp);
}

// compute elliptic integral K(k) for _n recursions
void ellipkf(float _k,
             unsigned int _n,
             float * _K,
             float * _Kp)
{
    float kn = _k;
    float knp = sqrtf(1.0f - _k*_k);
    float K  = 0.5f*M_PI;
    float Kp = 0.5f*M_PI;
    unsigned int i;
    for (i=0; i<_n; i++) {
        kn  = landenf(kn);
        knp = landenf(knp);

        K  *= (1.0f + kn);
        Kp *= (1.0f + knp);
    }
    *_K  = K;
    *_Kp = Kp;
}

// elliptic degree
float ellipdegf(float _N,
                float _k1,
                unsigned int _n)
{
    // compute K1, K1p from _k1
    float K1, K1p;
    ellipkf(_k1,_n,&K1,&K1p);

    // compute q1 from K1, K1p
    float q1 = expf(-M_PI*K1p/K1);

    // compute q from q1
    float q = powf(q1,1.0f/_N);

    // expand numerator, denominator
    unsigned int m;
    float b = 0.0f;
    for (m=0; m<_n; m++)
        b += powf(q,(float)(m*(m+1)));
    float a = 0.0f;
    for (m=1; m<_n; m++)
        a += powf(q,(float)(m*m));

    float g = b / (1.0f + 2.0f*a);
    float k = 4.0f*sqrtf(q)*g*g;
    return k;
}

// 
float complex ellip_pqf(float complex _u,
                       float complex _wM,
                       float _k,
                       unsigned int _M)
{
    float complex wn = _wM;
    float complex wn_inv = 1.0f / wn;
    float kn;
    unsigned int i;
    for (i=_M; i>0; i--) {
        unsigned int j;
        kn = _k;
        for (j=0; j<i; j++)
            kn = landenf(kn);

        wn = 1.0f / wn_inv;
        wn_inv = 1.0f / (1.0f + kn) * (wn_inv + kn*wn);
    }

    return 1.0f / wn_inv;
}

// 
float complex ellip_apqf(float complex _u,
                         float complex _w0,
                         float _k,
                         unsigned int _M)
{
    float complex wn = _w0;
    float complex kn;
    unsigned int i;
    for (i=1; i<=_M; i++) {
        unsigned int j;
        kn = _k;
        for (j=0; j<i; j++)
            kn = landenf(kn);

        wn = 2*wn / ((1.0f+kn)*(1.0f + sqrtf(1.0f - kn*kn*wn*wn)));
    }

    return wn;
}

// elliptic cd() function (_n recursions)
float complex ellip_cdf(float complex _u,
                        float complex _k,
                        unsigned int _n)
{
    float complex wM = cosf(_u*M_PI*0.5f);
    return ellip_pqf(_u,wM,_k,_n);
}

// elliptic inverse cd() function (_n recursions)
float complex ellip_acdf(float complex _u,
                         float complex _k,
                         unsigned int _n)
{
    float complex w0 = cosf(_u*M_PI*0.5f);
    return ellip_apqf(_u,w0,_k,_n);
}

// elliptic sn() function (_n recursions)
float complex ellip_snf(float complex _u,
                        float complex _k,
                        unsigned int _n)
{
    float complex wM = sinf(_u*M_PI*0.5f);
    return ellip_pqf(_u,wM,_k,_n);
}

// elliptic inverse sn() function (_n recursions)
float complex ellip_asnf(float complex _u,
                         float complex _k,
                         unsigned int _n)
{
    float complex w0 = sinf(_u*M_PI*0.5f);
    return ellip_apqf(_u,w0,_k,_n);
}

{'fraction_non_alphanumeric': 0.12591443138993572, 'fraction_numerical': 0.024163156727998226, 'mean_word_length': 1.859904912836767, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5769850', 'n_tokens_mistral': 3253, 'n_tokens_neox': 3024, 'n_words': 966}