defininig periodic functions

Librairie C++ de calcul formel/ C++ symbolic computation library

Modérateur : xcasadmin

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

defininig periodic functions

Message par lukamar » dim. sept. 03, 2017 8:44 am

Hi,

I would like to suggest a way for defining non-trigonometric periodic functions. Namely, if f(x) is defined on [a,b], then F(x)=f(x-T*floor((x-a)/T)), where T=b-a, is periodic with period T where the segment of graph of f on [a,b] is repeated. This could be automated with the following code:

Code : Tout sélectionner

gen _periodic(const gen & g,GIAC_CONTEXT) {
    if (g.type==_STRNG && g.subtype==-1) return g;
    if (g.type!=_VECT || g.subtype!=_SEQ__VECT)
        return gentypeerr(contextptr);
    vecteur & gv = *g._VECTptr;
    if (gv.size()!=4 && gv.size()!=2)
        return gensizeerr(contextptr);
    gen & e=gv[0],x,a,b;
    if (e.type!=_SYMB)
        return gentypeerr(contextptr);
    vecteur vars(*_lname(e,contextptr)._VECTptr);
    if (vars.empty())
        return e;
    if (gv.size()==2) {
        if (!gv[1].is_symb_of_sommet(at_equal))
            return gentypeerr(contextptr);
        vecteur & fl=*gv[1]._SYMBptr->feuille._VECTptr;
        if ((x=fl[0]).type!=_IDNT || !fl[1].is_symb_of_sommet(at_interval))
            return gentypeerr(contextptr);
        vecteur & ab=*fl[1]._SYMBptr->feuille._VECTptr;
        a=ab[0];
        b=ab[1];
    }
    else {
        x=gv[1];
        if (x.type!=_IDNT)
            return gentypeerr(contextptr);
        if (find(vars.begin(),vars.end(),x)==vars.end())
            return e;
        a=gv[2];
        b=gv[3];
    }
    gen T(b-a);
    if (!is_strictly_positive(T,contextptr))
        return gentypeerr(contextptr);
    gen p(subst(e,x,x-T*_floor((x-a)/T,contextptr),false,contextptr));
    return _unapply(makesequence(p,x),contextptr);
}
static const char _periodic_s []="periodic";
static define_unary_function_eval (__periodic,&_periodic,_periodic_s);
define_unary_function_ptr5(at_periodic,alias_at_periodic,&__periodic,0,true);
Now the above periodic function is obtained by calling

Code : Tout sélectionner

periodic(f(x),x,a,b)
or

Code : Tout sélectionner

periodic(f(x),x=a..b)
For example:

Code : Tout sélectionner

f:=periodic(x^2,x,-1,1);
plot(f(x),x=-5..5)

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: defininig periodic functions

Message par parisse » dim. sept. 03, 2017 6:54 pm

Good idea, but I think it would be more coherent to return an expression since we feed an expression as 1st argument.

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: defininig periodic functions

Message par lukamar » dim. sept. 03, 2017 8:22 pm

You are right, most of Giac functions follow the same principle. Then the last line of _periodic should be just

Code : Tout sélectionner

return p;

lukamar
Messages : 331
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

Re: defininig periodic functions

Message par lukamar » dim. mai 12, 2019 5:44 pm

Hello Bernard,
I noticed a bug in periodic command. On the line 6195 in intg.cc there is a check

Code : Tout sélectionner

if (e.type!=_SYMB)
return gentypeerr(contextptr);
It should be removed (i.e. lines 6195-6196 deleted) because it prevents entering e.g. x as the first argument.

By the way, I also coded a better implementation of Yen's algorithm using heaps (it is much faster). The changed files are graphe.cc/h and graphtheory.cc in my Github repository. Can you please synchronize the code?

parisse
Messages : 5731
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Re: defininig periodic functions

Message par parisse » lun. mai 13, 2019 7:20 am

My source is now updated, thank you!

Répondre