Skip to content
Find file
Fetching contributors…
Cannot retrieve contributors at this time
286 lines (235 sloc) 7.74 KB
/* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
/* License in license.txt. */
/* */
/* Verbs: Complex-Valued Scalar Functions */
#include "j.h"
static Z zj={0,1};
static Z z1={1,0};
static D hypoth(D u,D v){D p,q,t; MMM(u,v); R INF(p)?inf:p?(t=q/p,p*sqrt(1+t*t)):0;}
static ZF1(jtzjx){Z z; z.re=-v.im; z.im= v.re; R z;}
static ZF1(jtzmj){Z z; z.re= v.im; z.im=-v.re; R z;}
Z zrj0(D a){Z z; z.re=a; z.im=0.0; R z;}
ZS1(jtzconjug, zr=a; zi=-b;)
ZS2(jtzplus, zr=a+c; zi=b+d;)
ZS2(jtzminus, zr=a-c; zi=b-d;)
ZF1(jtztrend){D a,b,t;Z z;
a=v.re; b=v.im;
if(ZOV(v)){a/=2; b/=2;}
t=hypoth(a,b);
if(t<inf){if(!t)++t; z.re=a/t; z.im=b/t;}
else switch((INF(a)?2:0)+INF(b)){
case 1: z.re=0.0; z.im=SGN(b); break;
case 2: z.re=SGN(a); z.im=0.0; break;
case 3: ZASSERT(0,EVNAN);
}
R z;
}
ZF2(jtztymes){D a,b,c,d;Z z;
a=u.re; b=u.im; c=v.re; d=v.im;
z.re=TYMES(a,c)-TYMES(b,d);
z.im=TYMES(a,d)+TYMES(b,c);
R z;
}
ZF2(jtzdiv){ZF2DECL;D t;
if(ZNZ(v)){
if(ABS(c)<ABS(d)){t=a; a=-b; b=t; t=c; c=-d; d=t;}
a/=c; b/=c; d/=c; t=1+d*d; zr=(a+(b&&d?b*d:0.0))/t; zi=(b-(a&&d?a*d:0))/t;
}else if(ZNZ(u))switch(2*(0>a)+(0>b)){
case 0: if(a> b)zr= inf; else zi= inf; break;
case 1: if(a>-b)zr= inf; else zi=-inf; break;
case 2: if(a<-b)zr=-inf; else zi= inf; break;
case 3: if(a< b)zr=-inf; else zi=-inf;
}
ZEPILOG;
}
ZF1(jtznegate){R zminus(zeroZ,v);}
D zmag(Z v){R hypoth(v.re,v.im);}
B jtzeq(J jt,Z u,Z v){D a=u.re,b=u.im,c=v.re,d=v.im,p,q;
if(a==c&&b==d)R 1;
if(ZEZ(u)||ZEZ(v)||!jt->ct||(0>a!=0>c&&0>b!=0>d))R 0;
if(ZOV(u)||ZOV(v)){a/=2; b/=2; c/=2; d/=2;}
if(ZUN(u)||ZUN(v)){a*=2; b*=2; c*=2; d*=2;}
p=hypoth(a,b); q=hypoth(c,d);
R p!=inf && q!=inf && hypoth(a-c,b-d)<=jt->ct*MAX(p,q);
}
ZF1(jtzfloor){D p,q;
ZF1DECL;
zr=jfloor(a); p=a-zr;
zi=jfloor(b); q=b-zi;
if(1<=p+q+jt->ct)if(p>=q)++zr; else ++zi;
ZEPILOG;
}
ZF1(jtzceil){R znegate(zfloor(znegate(v)));}
ZF2(jtzrem){D a,b,d;Z q;
if(ZEZ(u))R v;
ZASSERT(!ZINF(v),EVNAN);
if(INF(u.re)&&!u.im&&!v.im){
if(u.re==inf )R 0<=v.re?v:u;
if(u.re==infm)R 0>=v.re?v:u;
}
ZASSERT(!ZINF(u),EVNONCE);
d=u.re*u.re+u.im*u.im;
a=u.re*v.re+u.im*v.im; q.re=tfloor(a/d);
b=u.re*v.im-u.im*v.re; q.im=tfloor(b/d);
R zminus(v,ztymes(u,q));
}
ZF2(jtzgcd){D a,b;Z t,z;
ZASSERT(!(ZINF(u)||ZINF(v)),EVNAN);
while(ZNZ(u)){
t=zrem(u,v); if(t.re>v.re/2)t.re-=v.re; if(t.im>v.im/2)t.im-=v.im;
v.re=u.re; v.im=u.im; u.re=t.re; u.im=t.im;
}
z.re=a=v.re; z.im=b=v.im;
switch(2*(0>a)+(0>b)){
case 0: if(!a){z.re= b; z.im=0;} break;
case 1: z.re=-b; z.im= a; break;
case 2: if(!b){z.re=-a; z.im=0;}else{z.re= b; z.im=-a;} break;
case 3: z.re=-a; z.im=-b;
}
R z;
}
ZF2(jtzlcm){ZASSERT(!(ZINF(u)||ZINF(v)),EVNAN); R ZEZ(u)||ZEZ(v) ? zeroZ : ztymes(u,zdiv(v,zgcd(u,v)));}
ZF1(jtzexp){D a,b,c,s,t;Z z;
a=v.re; b=v.im;
if(a<EMIN)z.re=z.im=0.0;
else{
ZASSERT(-THMAX<b&&b<THMAX,EVLIMIT);
c=cos(b); s=sin(b);
if(a<=EMAX){t=exp(a); z.re=t*c; z.im=t*s;}
else{
if(!c)z.re=0; else{t=a+log(ABS(c)); t=EMAX<t?inf:exp(t); z.re=0>c?-t:t;}
if(!s)z.im=0; else{t=a+log(ABS(s)); t=EMAX<t?inf:exp(t); z.im=0>s?-t:t;}
}}
R z;
}
ZF1(jtzlog){ZF1DECL;
zr=b?log(hypoth(a,b)):INF(a)?inf:a?log(hypoth(a,b)):-inf;
zi=a||b?atan2(b,a):0;
#if SY_WINCE_MIPS && !defined(WIN32_PLATFORM_PSPC)
if(!b) zi=a<0?PI : 0; /* atan2(0,v) fails in mips handheld wince - _9^2.5-1.5*/
#endif
ZEPILOG;
}
ZF2(jtzpow){ZF2DECL;D m;I n;
if(!a&&!b){z.re=d?0:0>c?inf:!c; z.im=0; R z;}
if(!d&&IMIN<c&&c<=IMAX&&(n=(I)jfloor(c),c==n)){
if(0>n){u=zdiv(z1,u); n=-n;}
z=z1;
while(n){if(1&n)z=ztymes(z,u); u=ztymes(u,u); n>>=1;}
R z;
}
z=zexp(ztymes(v,zlog(u)));
if(!b&&!d){
m=jfloor(c);
if(0>a&&c>m&&c==0.5+m)z.re=0;
if(c==m)z.im=0;
}
R z;
}
ZF1(jtzsqrt){D p,q,t;
ZF1DECL;
MMM(a,b);
if(p){
t=0.5*q/p; zr=sqrt(ABS(a/2)+p*sqrt(0.25+t*t)); zi=b/(zr+zr);
if(0>a){t=ABS(zi); zi=0>b?-zr:zr; zr=t;}
}
ZEPILOG;
}
/* See Abramowitz & Stegun, Handbook of Mathematical Functions, */
/* National Bureau of Standards, 1964 6. */
static ZF1(jtzsin){D a,b,c,s;Z z;
a=v.re; b=v.im;
ZASSERT(-THMAX<a&&a<THMAX,EVLIMIT);
s=sin(a); z.re=s?s*(b<-EMAX2|| EMAX2<b?inf:cosh(b)):0.0;
c=cos(a); z.im=c?c*(b<-EMAX2?infm:EMAX2<b?inf:sinh(b)):0.0;
R z;
} /* 4.3.55 */
static ZF1(jtzcos){D a,b,c,s;Z z;
a=v.re; b=v.im;
ZASSERT(-THMAX<a&&a<THMAX,EVLIMIT);
c=cos(a); z.re=c? c*(b<-EMAX2|| EMAX2<b?inf:cosh(b)):0.0;
s=sin(a); z.im=s?-s*(b<-EMAX2?infm:EMAX2<b?inf:sinh(b)):0.0;
R z;
} /* 4.3.56 */
static ZF1(jtztan){R zdiv(zsin(v),zcos(v));}
static ZF1(jtzp4){R zsqrt(zplus(z1,ztymes(v,v)));}
static ZF1(jtzm4){R 1e16<hypoth(v.re,v.im)?v:ztymes(zplus(v,z1),zsqrt(zdiv(zminus(v,z1),zplus(v,z1))));}
static ZF1(jtzsinh){R zmj(zsin(zjx(v)));} /* 4.5.7 */
static ZF1(jtzcosh){R zcos(zjx(v));} /* 4.5.8 */
static ZF1(jtztanh){R v.re<-TMAX?zrj0(-1.0):TMAX<v.re?z1:zdiv(zsinh(v),zcosh(v));}
static ZF1(jtzp8){R zsqrt(ztymes(zplus(zj,v),zminus(zj,v)));}
static ZF1(jtzasinh){R 0>v.re ? znegate(zasinh(znegate(v))) : zlog(zplus(v,zp4(v)));}
static ZF1(jtzacosh){Z z;
z=zlog(zplus(v,zm4(v)));
if(0>=z.re){z.re=0; z.im=ABS(z.im);}
R z;
}
static ZF1(jtzatanh){R ztymes(zrj0((D)0.5),zlog(zdiv(zplus(z1,v),zminus(z1,v))));}
static ZF1(jtzatan){ZF1DECL;
if(!b&&(a<-1e13||1e13<a))R zrj0(0<a?PI/2.0:-PI/2.0);
z=zmj(zatanh(zjx(v)));
if(!b)z.im=0;
R z;
} /* 4.4.22 */
static ZF1(jtzasin){R !v.im&&-1<=v.re&&v.re<=1?zrj0(asin(v.re)):zmj(zasinh(zjx(v)));} /* 4.4.20 */
static ZF1(jtzacos){R zminus(zrj0(PI/2.0),zasin(v));}
static ZF1(jtzarc){D x,y;Z t,z;
z.re=z.im=0;
t=ztrend(v); x=t.re; y=t.im;
if(0!=x||0!=y)z.re=atan2(y,x);
#if SY_WINCE_MIPS && !defined(WIN32_PLATFORM_PSPC)
if(!y) z.re=x<0?PI : 0; /* atan2(0,v) fails in mips handheld wince - 12 o. _3 */
#endif
R z;
}
ZF2(jtzcir){D r;I x;Z z;
z=zeroZ;
r=u.re;
x=(I)jfloor(0.5+r);
ZASSERT(-12<=r&&r<=12&&FEQ(x,r)&&!u.im,EVDOMAIN);
switch(x){
default: ZASSERT(0,EVDOMAIN);
case 0: R zsqrt(ztymes(zplus(z1,v),zminus(z1,v)));
case 1: R zsin(v); case -1: R zasin(v);
case 2: R zcos(v); case -2: R zacos(v);
case 3: R ztan(v); case -3: R zatan(v);
case 4: R zp4(v); case -4: R zm4(v);
case 5: R zsinh(v); case -5: R zasinh(v);
case 6: R zcosh(v); case -6: R zacosh(v);
case 7: R ztanh(v); case -7: R zatanh(v);
case 8: R zp8(v); case -8: R znegate(zp8(v));
case 9: z.re=v.re; R z; case -9: R v;
case 10: z.re=zmag(v); R z; case -10: R zconjug(v);
case 11: z.re=v.im; R z; case -11: R zjx(v);
case 12: R zarc(v); case -12: R zexp(zjx(v));
}}
B jtztridiag(J jt,I n,A a,A x){I i,j,n1=n-1;Z*av,d,p,*xv;
av=ZAV(a); xv=ZAV(x); d=xv[0];
for(i=j=0;i<n1;++i){
ASSERT(ZNZ(d),EVDOMAIN);
p=zdiv(xv[j+2],d);
xv[j+3]=d=zminus(xv[j+3],ztymes(p,xv[j+1]));
av[i+1]= zminus(av[i+1],ztymes(p,av[i] ));
j+=3;
}
ASSERT(ZNZ(d),EVDOMAIN);
i=n-1; j=AN(x)-1; av[i]=zdiv(av[i],d);
for(i=n-2;i>=0;--i){j-=3; av[i]=zdiv(zminus(av[i],ztymes(xv[j+1],av[i+1])),xv[j]);}
R 1;
}
DF1(jtexppi){A z;B b;D r,th,y;I k;Z*v,t;
F1RANK(0,jtexppi,0);
if(!(CMPX&AT(w)))R expn1(pix(w));
v=ZAV(w); r=exp(PI*v->re); y=v->im; if(b=0>y)y=-y;
th=y-2*(I)(y/2); k=(I)(2*th); if(k!=2*th)k=-1; else if(b&&k)k=4-k;
if(!(0<=k&&k<=3))R expn1(pix(w));
switch(k){
case 0: t.re= r; t.im= 0; break;
case 1: t.re= 0; t.im= r; break;
case 2: t.re=-r; t.im= 0; break;
case 3: t.re= 0; t.im=-r; break;
}
GA(z,CMPX,1,0,0); *ZAV(z)=t; R z;
} /* special code for ^@o. */
ZF1(jtznonce1){ZASSERT(0,EVNONCE);}
ZF2(jtznonce2){ZASSERT(0,EVNONCE);}
Something went wrong with that request. Please try again.