Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 52

PARI/GP code for Quaternions

{
my(
isreal=((x)->my(t=type(x));t=="t_INT"||t=="t_FRAC"||t=="t_REAL"),
isscalar=((x)->my(t=type(x));t=="t_INT"||t=="t_FRAC"||t=="t_REAL"||t=="t_COMPLEX"),
isquaternion=((x)->my(t=type(x));if(t=="t_COL"&&#x==4,for(i=1,4,if(!isreal(x[i]),return(0)));1,0)),
isvalid=((x)->isscalar(x)||isquaternion(x)),
toquaternion=((x)->if(isquaternion(x),x,isscalar(x),[real(x),imag(x),0,0],[])),
Qapply_=((f,x)->if(isscalar(x),f(x),isquaternion(x),my(v=x[^1],r=norml2(v)^(1/2),res=f(x[1]+I*r));concat(real(res),imag(res)*if(r,v/r,[1,0,0]~)),error("invalid type in Qapply")))
);
Qabs=((x)->if(isscalar(x),abs(x),isquaternion(x),norml2(x)^(1/2),error("invalid type in Qabs")));
Qadd=((x,y)->if(isvalid(x)&&isvalid(y),algadd(,x,y),error("invalid type in Qadd")));
Qapply=Qapply_;
Qarg=((x)->if(isscalar(x),arg(x),isquaternion(x),arg(x[1]+I*norml2(x[^1])^(1/2)),error("invalid type in Qarg")));
Qconj=((x)->if(isscalar(x),conj(x),isquaternion(x),concat(x[1],-x[^1]),error("invalid type in Qconj")));
Qdivl=((x,y)->if(isvalid(x)&&isvalid(y),algdivl(,x,y),error("invalid type in Qdivl")));
Qdivr=((x,y)->if(isvalid(x)&&isvalid(y),algdivr(,x,y),error("invalid type in Qdivr")));
Qinv=((x)->if(isvalid(x),alginv(,x),error("invalid type in Qinv")));
Qmul=((x,y)->if(isvalid(x)&&isvalid(y),algmul(,x,y),error("invalid type in Qmul")));
Qneg=((x)->if(isvalid(x),algneg(,x),error("invalid type in Qneg")));
Qpoleval=((T,x)->if(isvalid(x),algpoleval(,T,x),error("invalid type in Qpoleval")));
Qpow=((x,y)->if(type(y)=="t_INT",if(isvalid(x),return(algpow(,x,y))),isscalar(x)&&isscalar(y),return(x^y),isreal(y),if(isvalid(x),return(Qapply_((x)->x^y,x))),isreal(x),if(isvalid(y),return(Qapply_((y)->x^y,y))));error("invalid type in Qpow"));
Qreal=((x)->if(isscalar(x),real(x),isquaternion(x),x[1],error("invalid type in Qreal")));
Qsign=((x)->if(isreal(x),sign(x),isvalid(x),if(x,x/norml2(x)^(1/2),0*x),error("invalid type in Qsign")));
Qslerp=((x,y,t)->if(isvalid(x)&&isvalid(y)&&isreal(t),if(t===0,x,t===1,y,if(norml2(x)<norml2(y),[x,y,t]=[y,x,1-t]);algmul(,x,Qapply_((x)->x^t,algdivl(,x,y)))),error("invalid type in Qslerp")));
Qslerpinv=((x,y,z)->x=toquaternion(x);y=toquaternion(y);z=toquaternion(z);if(x==[]||y==[]||z==[],error("invalid type in Qslerpinv"));my(a,b,c,d,na,nb,nc,nx,nx2);nx2=norml2(x);nx=sqrt(nx2);na=x*y~;nb=x*z~;a=x*(na/nx2);b=x*(nb/nx2);c=y-a;d=z-b;na/=nx;nb/=nx;nc=sqrt(norml2(c));arg(nb+I*if(nc,(c*d~)/nc,sqrt(norml2(d))))/arg(na+nc*I));
Qsqr=((x)->if(isvalid(x),algneg(,x),error("invalid type in Qsqr")));
Qsub=((x,y)->if(isvalid(x)&&isvalid(y),algsub(,x,y),error("invalid type in Qsub")));
}

Last edited by lanxiyu (2024-04-11 00:02:54)