Mercurial > hg > jgplsrc

comparison vm.c @ 0:e0bbaa717f41 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
lol J
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Mon, 25 Nov 2013 11:56:30 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e0bbaa717f41
1 /* Copyright 1990-2011, Jsoftware Inc. All rights reserved. */
2 /* License in license.txt. */
3 /* */
4 /* Verbs: "Mathematical" Functions (Irrational, Transcendental, etc.) */
5
6 #include "j.h"
7 #include "ve.h"
8
9 static D jtintpow(J jt,D x,I n){D r=1;
10 if(0>n){x=1/x; if(n==IMIN){r=x; n=IMAX;} else n=-n;}
11 while(n){if(1&n)r*=x; x*=x; n>>=1;}
12 R r;
13 } /* x^n where x is real and n is integral */
14
15 static D jtpospow(J jt,D x,D y){
16 if(0==y)R 1.0;
17 if(0==x)R 0<y?0.0:inf;
18 if(0<x){
19 if(y== inf)R 1<x?inf:1>x?0.0:1.0;
20 if(y==-inf)R 1<x?0.0:1>x?inf:1.0;
21 R exp(y*log(x));
22 }
23 if(y==-inf){ASSERT(-1>x,EVDOMAIN); R 0.0;}
24 if(y== inf){ASSERT(-1<x,EVDOMAIN); R 0.0;}
25 jt->jerr=EWIMAG;
26 R 0;
27 } /* x^y where x and y are real and x is non-negative */
28
29 #define POWXB(u,v) (v?u:1)
30 #define POWBX(u,v) (u?1.0:v<0?inf:!v)
31 #define POWII(u,v) intpow((D)u,v)
32 #define POWID(u,v) pospow((D)u,v)
33
34 APFX(powBI, D,B,I, POWBX )
35 APFX(powBD, D,B,D, POWBX )
36 APFX(powIB, I,I,B, POWXB )
37 APFX(powII, D,I,I, POWII )
38 APFX(powID, D,I,D, POWID )
39 APFX(powDB, D,D,B, POWXB )
40 APFX(powDI, D,D,I, intpow)
41 APFX(powDD, D,D,D, pospow)
42 APFX(powZZ, Z,Z,Z, zpow )
43
44 ANAN(cirZZ, Z,Z,Z, zcir )
45
46 static void jtcirx(J jt,I n,I k,D*z,D*y){D p,t;
47 NAN0;
48 switch(k){
49 default: ASSERTW(0,EWIMAG);
50 case 0: DO(n, t=*y++; ASSERTW( -1.0<=t&&t<=1.0, EWIMAG ); *z++=sqrt(1.0-t*t);); break;
51 case 1: DO(n, t=*y++; ASSERTW(-THMAX<t&&t<THMAX,EVLIMIT); *z++=sin(t);); break;
52 case 2: DO(n, t=*y++; ASSERTW(-THMAX<t&&t<THMAX,EVLIMIT); *z++=cos(t);); break;
53 case 3: DO(n, t=*y++; ASSERTW(-THMAX<t&&t<THMAX,EVLIMIT); *z++=tan(t);); break;
54 case 4: DO(n, t=*y++; *z++=t<-1e8?-t:1e8<t?t:sqrt(t*t+1.0);); break;
55 case 5: DO(n, t=*y++; *z++=t<-EMAX2?infm:EMAX2<t?inf:sinh(t);); break;
56 case 6: DO(n, t=*y++; *z++=t<-EMAX2|| EMAX2<t?inf:cosh(t);); break;
57 case 7: DO(n, t=*y++; *z++=t<-TMAX?-1:TMAX<t?1:tanh(t);); break;
58 case -1: DO(n, t=*y++; ASSERTW( -1.0<=t&&t<=1.0, EWIMAG ); *z++=asin(t);); break;
59 case -2: DO(n, t=*y++; ASSERTW( -1.0<=t&&t<=1.0, EWIMAG ); *z++=acos(t);); break;
60 case -3: DO(n, *z++=atan(*y++);); break;
61 case -4: DO(n, t=*y++; ASSERTW(t<=-1.0||1.0<=t, EWIMAG ); *z++=t<-1e8||1e8<t?t:t==-1?0:(t+1)*sqrt((t-1)/(t+1));); break;
62 case -5: p=log(2.0);
63 DO(n, t=*y++; *z++=1.0e8<t?p+log(t):-7.8e3>t?-(p+log(-t)):log(t+sqrt(t*t+1.0));); break;
64 case -6: p=log(2.0);
65 DO(n, t=*y++; ASSERTW( 1.0<=t, EWIMAG ); *z++=1.0e8<t?p+log(t):log(t+sqrt(t*t-1.0));); break;
66 case -7: DO(n, t=*y++; ASSERTW( -1.0<=t&&t<=1.0, EWIMAG ); *z++=0.5*log((1.0+t)/(1.0-t));); break;
67 case 9: DO(n, *z++=*y++;); break;
68 case 10: DO(n, t=*y++; *z++=ABS(t);); break;
69 case 11: DO(n, *z++=0.0;); break;
70 case 12: DO(n, *z++=0<=*y++?0.0:PI;); break;
71 }
72 NAN1V;
73 }
74
75 AHDR2(cirBD,D,B,D){ASSERTW(b&&1==m,EWIMAG); cirx(n,(I)*x,z,y);}
76 AHDR2(cirID,D,I,D){ASSERTW(b&&1==m,EWIMAG); cirx(n, *x,z,y);}
77
78 AHDR2(cirDD,D,D,D){I k=(I)jfloor(0.5+*x);
79 ASSERTW(k==*x,EVDOMAIN);
80 ASSERTW(b&&1==m,EWIMAG);
81 cirx(n,k,z,y);
82 }
83
84
85 F2(jtlogar2){A z;I t;
86 RZ(a&&w);
87 t=maxtype(AT(a),AT(w));
88 if(!(t&XNUM)||jt->xmode==XMEXACT){jt->xmode=XMEXACT; R divide(logar1(w),logar1(a));}
89 z=rank2ex(cvt(XNUM,a),cvt(XNUM,w),0L,0L,0L,jtxlog2a);
90 if(z)R z;
91 if(jt->jerr==EWIMAG||jt->jerr==EWIRR){RESETERR; jt->xmode=XMEXACT; R divide(logar1(w),logar1(a));}
92 R 0;
93 }
94
95 F2(jtroot){A z;I t;
96 RZ(a&&w);
97 t=maxtype(AT(a),AT(w));
98 if(!(t&XNUM))R expn2(cvt(t,w),recip(cvt(t,a)));
99 z=rank2ex(cvt(XNUM,a),cvt(XNUM,w),0L,0L,0L,jtxroota);
100 switch(jt->jerr){
101 case EWIMAG: RESETERR; R expn2(cvt(CMPX,w),recip(cvt(CMPX,a)));
102 case EWIRR: RESETERR; R expn2(cvt(FL, w),recip(cvt(FL, a)));
103 default: R z;
104 }}
105
106 F1(jtjdot1){R tymes(a0j1,w);}
107 F2(jtjdot2){R plus(a,tymes(a0j1,w));}
108 F1(jtrdot1){R expn1(jdot1(w));}
109 F2(jtrdot2){R tymes(a,rdot1(w));}
110
111
112 F1(jtpolar){RZ(w); R cvt(SPARSE&AT(w)?SFL:FL,df2(v2(10L,12L),w,qq(ds(CCIRCLE),v2(1L,0L))));}
113
114 F1(jtrect){A e,z;B b;I r,t;P*wp,*zp;Z c;
115 RZ(w);
116 t=AT(w); r=AR(w); jt->rank=0;
117 ASSERT(!AN(w)||t&NUMERIC,EVDOMAIN);
118 if(t&CMPX){GA(z,FL,2*AN(w),1+r,AS(w)); *(AS(z)+r)=2; MC(AV(z),AV(w),AN(z)*sizeof(D)); R z;}
119 else if(t&SPARSE){
120 b=1&&t&SCMPX;
121 GA(z,b?SFL:t,1,1+r,AS(w)); *(AS(z)+r)=2;
122 wp=PAV(w); zp=PAV(z);
123 if(b){e=SPA(wp,e); c=*ZAV(e); ASSERT(FEQ(c.re,c.im),EVSPARSE); SPB(zp,e,scf(c.re));}
124 else SPB(zp,e,ca(SPA(wp,e)));
125 SPB(zp,a,ca(SPA(wp,a)));
126 SPB(zp,i,ca(SPA(wp,i)));
127 SPB(zp,x,rect(SPA(wp,x)));
128 R z;
129 }else R df2(w,zero,qq(ds(CCOMMA),zero));
130 }