Mercurial > hg > jgplsrc

comparison ve.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: Elementary Functions (Arithmetic, etc.) */
5
6 #include "j.h"
7 #include "vasm.h"
8
9 #define DIVI(u,v) (u||v ? u/(D)v : 0.0)
10 #define DIVBB(u,v) (v?u:u?inf:0.0)
11
12 #define TYMESBX(u,v) (u?v:0)
13 #define TYMESXB(u,v) (v?u:0)
14 #define TYMESID(u,v) (u ?u*v:0)
15 #define TYMESDI(u,v) ( v?u*v:0)
16 #define TYMESDD(u,v) (u&&v?u*v:0)
17
18 AOVF( plusII, I,I,I, PLUSVV, PLUS1V, PLUSV1)
19 AOVF(minusII, I,I,I, MINUSVV,MINUS1V,MINUSV1)
20 AOVF(tymesII, I,I,I, TYMESVV,TYMES1V,TYMESV1)
21
22 APFX( plusIO, D,I,I, PLUSO)
23 APFX(minusIO, D,I,I, MINUSO)
24 APFX(tymesIO, D,I,I, TYMESO)
25
26 AIFX( plusBB, I,B,B, + ) /* plusII */ AIFX( plusBD, D,B,D, + )
27 /* plusII */ /* plusII */ AIFX( plusID, D,I,D, + )
28 AIFX( plusDB, D,D,B, + ) AIFX( plusDI, D,D,I, +) ANAN( plusDD, D,D,D, PLUS)
29 ANAN( plusZZ, Z,Z,Z, zplus )
30
31 AIFX(minusBB, I,B,B, - ) /* minusII */ AIFX(minusBD, D,B,D, - )
32 /* minusII */ /* minusII */ AIFX(minusID, D,I,D, - )
33 AIFX(minusDB, D,D,B, - ) AIFX(minusDI, D,D,I, -) ANAN(minusDD, D,D,D, MINUS)
34 ANAN(minusZZ, Z,Z,Z, zminus)
35
36 /* andBB */ APFX(tymesBI, I,B,I, TYMESBX) APFX(tymesBD, D,B,D, TYMESBX)
37 APFX(tymesIB, I,I,B, TYMESXB) /* tymesII */ APFX(tymesID, D,I,D, TYMESID)
38 APFX(tymesDB, D,D,B, TYMESXB) APFX(tymesDI, D,D,I, TYMESDI) APFX(tymesDD, D,D,D, TYMESDD)
39 ANAN(tymesZZ, Z,Z,Z, ztymes )
40
41 APFX( divBB, D,B,B, DIVBB) APFX( divBI, D,B,I, DIVI) APFX( divBD, D,B,D, DIV)
42 APFX( divIB, D,I,B, DIVI ) APFX( divII, D,I,I, DIVI) APFX( divID, D,I,D, DIV)
43 APFX( divDB, D,D,B, DIVI ) APFX( divDI, D,D,I, DIVI) ANAN( divDD, D,D,D, DIV)
44 ANAN( divZZ, Z,Z,Z, zdiv )
45
46 /* orBB */ APFX( minBI, I,B,I, MIN) APFX( minBD, D,B,D, MIN)
47 APFX( minIB, I,I,B, MIN) APFX( minII, I,I,I, MIN) APFX( minID, D,I,D, MIN)
48 APFX( minDB, D,D,B, MIN) APFX( minDI, D,D,I, MIN) APFX( minDD, D,D,D, MIN)
49 APFX( minSS, SB,SB,SB, SBMIN)
50
51 /* andBB */ APFX( maxBI, I,B,I, MAX) APFX( maxBD, D,B,D, MAX)
52 APFX( maxIB, I,I,B, MAX) APFX( maxII, I,I,I, MAX) APFX( maxID, D,I,D, MAX)
53 APFX( maxDB, D,D,B, MAX) APFX( maxDI, D,D,I, MAX) APFX( maxDD, D,D,D, MAX)
54 APFX( maxSS, SB,SB,SB, SBMAX)
55
56
57 static D jtremdd(J jt,D a,D b){D q,x,y;
58 if(!a)R b;
59 ASSERT(!INF(b),EVNAN);
60 if(a==inf )R 0<=b?b:a;
61 if(a==infm)R 0>=b?b:a;
62 q=b/a; x=tfloor(q); y=tceil(q); R teq(x,y)?0:b-a*x;
63 }
64
65 ANAN(remDD, D,D,D, remdd)
66 ANAN(remZZ, Z,Z,Z, zrem )
67
68 static I jtremid(J jt,I a,D b){D r;I k;
69 ASSERT(a&&-9e15<b&&b<9e15,EWOV);
70 r=b-a*floor(b/a); k=(I)r;
71 ASSERT(k==r,EWOV);
72 R k;
73 }
74
75 APFX(remID, I,I,D, remid)
76
77 static I remii(I a,I b){I r; R!a?b:(r=b%a,0<a?r+a*(0>r):r+a*(0<r));}
78
79 AHDR2(remII,I,I,I){I u,v;
80 if(1==n) DO(m, *z++=remii(*x,*y); x++; y++; )
81 else if(b)DO(m, u=*x++; if(0<=u&&!(u&(u-1))){--u; DO(n, *z++=u&*y++;);}
82 else DO(n, *z++=remii( u,*y); y++;))
83 else DO(m, v=*y++; DO(n, *z++=remii(*x, v); x++; ));
84 }
85
86
87 static I igcd1(I a,I b){R a?igcd1(b%a,a):b;}
88
89 static I jtigcd(J jt,I a,I b){
90 if(a>IMIN&&b>IMIN){a=ABS(a); b=ABS(b);}
91 else{
92 if(a==b||!a||!b){jt->jerr=EWOV; R 0;}
93 if(a==IMIN){b=ABS(b); a=-(a+b);}else{a=ABS(a); b=-(a+b);}
94 }
95 R a?igcd1(b%a,a):b;
96 }
97
98 D jtdgcd(J jt,D a,D b){D a1,b1,t;
99 a=ABS(a); b=ABS(b); if(a>b){t=a; a=b; b=t;}
100 ASSERT(inf!=b,EVNAN);
101 if(!a)R b;
102 a1=a; b1=b;
103 while(remdd(a1/jfloor(0.5+a1/a),b1)){t=a; a=remdd(a,b); b=t;}
104 R a;
105 } /* D.L. Forkes 1984; E.E. McDonnell 1992 */
106 #if SY_64
107 #if SY_WIN32
108 static I jtilcm(J jt,I a,I b){C er=0;I b1,d,z;
109 if(a&&b){RZ(d=igcd(a,b)); b1=b/d; TYMESVV(1L,&z,&a,&b1); if(er)jt->jerr=EWOV; R z;}else R 0;
110 }
111 #else
112 static I jtilcm(J jt,I a,I b){LD z;I b1,d;
113 if(a&&b){RZ(d=igcd(a,b)); b1=b/d; z=a*(LD)b1; if(z<IMIN||IMAX<z)jt->jerr=EWOV; R (I)z;}else R 0;
114 }
115 #endif
116 #else
117 static I jtilcm(J jt,I a,I b){D z;I b1,d;
118 if(a&&b){RZ(d=igcd(a,b)); b1=b/d; z=a*(D)b1; if(z<IMIN||IMAX<z)jt->jerr=EWOV; R (I)z;}else R 0;
119 }
120 #endif
121
122 #define GCDIO(u,v) (dgcd((D)u,(D)v))
123 #define LCMIO(u,v) (dlcm((D)u,(D)v))
124
125 static D jtdlcm(J jt,D a,D b){ASSERT(!(INF(a)||INF(b)),EVNAN); R a&&b?a*(b/dgcd(a,b)):0;}
126
127 APFX(gcdIO, D,I,I, GCDIO)
128 APFX(gcdII, I,I,I, igcd )
129 APFX(gcdDD, D,D,D, dgcd )
130 APFX(gcdZZ, Z,Z,Z, zgcd )
131
132 APFX(lcmII, I,I,I, ilcm )
133 APFX(lcmIO, D,I,I, LCMIO)
134 APFX(lcmDD, D,D,D, dlcm )
135 APFX(lcmZZ, Z,Z,Z, zlcm )
136
137
138 #define GETD {d=*wv++; if(!d){z=0; break;}}
139 #define INTDIVF(c,d) (0>c==0>d?c/d:c%d?c/d-1:c/d)
140 #define INTDIVC(c,d) (0>c!=0>d?c/d:c%d?c/d+1:c/d)
141
142 F2(jtintdiv){A z;B b,flr;I an,ar,*as,*av,c,d,j,k,m,n,p,p1,r,*s,wn,wr,*ws,*wv,*zv;
143 RZ(a&&w);
144 an=AN(a); ar=AR(a); as=AS(a); av=AV(a);
145 wn=AN(w); wr=AR(w); ws=AS(w); wv=AV(w); b=ar>=wr; r=b?wr:ar; s=b?as:ws;
146 ASSERT(!ICMP(as,ws,r),EVLENGTH);
147 if(an&&wn){m=prod(r,s); n=prod(b?ar-r:wr-r,r+s);}else m=n=0;
148 GA(z,INT,b?an:wn,b?ar:wr,s); zv=AV(z);
149 d=wn?*wv:0; p=0<d?d:-d; p1=d==IMIN?p:p-1; flr=XMFLR==jt->xmode;
150 if(!wr&&p&&!(p&p1)){
151 k=0; j=1; while(p>j){++k; j<<=1;}
152 switch((0<d?0:2)+(flr?0:1)){
153 case 0: DO(n, *zv++=*av++>>k;); break;
154 case 1: DO(n, c=*av++; *zv++=0< c?1+((c-1)>>k):(c+p1)>>k;); break;
155 case 2: DO(n, c=*av++; *zv++=c>IMIN?-c>>k:-(-c>>k);); break;
156 case 3: DO(n, c=*av++; *zv++=0<=c?-(c>>k):1+(-(1+c)>>k););
157 }}else if(flr){
158 if(1==n) DO(m, c=*av++; GETD; *zv++=INTDIVF(c,d); )
159 else if(b)DO(m, GETD; DO(n, c=*av++; *zv++=INTDIVF(c,d););)
160 else DO(m, c=*av++; DO(n, GETD; *zv++=INTDIVF(c,d););)
161 }else{
162 if(1==n) DO(m, c=*av++; GETD; *zv++=INTDIVC(c,d); )
163 else if(b)DO(m, GETD; DO(n, c=*av++; *zv++=INTDIVC(c,d););)
164 else DO(m, c=*av++; DO(n, GETD; *zv++=INTDIVC(c,d););)
165 }
166 R z?z:flr?floor1(divide(a,w)):ceil1(divide(a,w));
167 } /* <.@% or >.@% on integers */
168
169
170 static F2(jtweight){RZ(a&&w); R df1(behead(over(AR(w)?w:reshape(a,w),one)),bsdot(slash(ds(CSTAR))));}
171
172 F1(jtbase1){A z;B*v;I c,d,m,n,p,r,*s,t,*x;
173 RZ(w);
174 n=AN(w); t=AT(w); r=AR(w); s=AS(w); c=r?*(s+r-1):1;
175 ASSERT(t&DENSE,EVNONCE);
176 if(c>(SY_64?63:31)||!(t&B01))R pdt(w,weight(sc(c),t&RAT+XNUM?cvt(XNUM,num[2]):num[2]));
177 if(c)m=n/c; else RE(m=prod(r-1,s));
178 GA(z,INT,m,r?r-1:0,s); x=m+AV(z); v=n+BAV(w);
179 if(c)DO(m, p=0; d=1; DO(c, if(*--v)p+=d; d+=d;); *--x=p;)
180 else memset(x-m,C0,m*SZI);
181 R z;
182 }
183
184 F2(jtbase2){I ar,*as,at,c,t,wr,*ws,wt;
185 RZ(a&&w);
186 at=AT(a); ar=AR(a); as=AS(a);
187 wt=AT(w); wr=AR(w); ws=AS(w); c=wr?*(ws+wr-1):1; t=maxtype(at,wt);
188 ASSERT(at&DENSE&&wt&DENSE,EVNONCE);
189 if(!(t&at))RZ(a=cvt(t,a));
190 if(!(t&wt))RZ(w=cvt(t,w));
191 R 1>=ar?pdt(w,weight(sc(c),a)):rank2ex(w,rank2ex(sc(c),a,0L,0L,1L,jtweight),0L,1L,1L,jtpdt);
192 }
193
194 F1(jtabase1){A d,z;B*zv;I c,n,p,r,t,*v,x;
195 RZ(w);
196 n=AN(w); r=AR(w); t=AT(w);
197 ASSERT(t&DENSE,EVNONCE);
198 if(!n||t&B01)R reshape(over(shape(w),n?one:zero),w);
199 if(!(t&INT)){
200 d=df2(num[2],maximum(one,aslash(CMAX,mag(ravel(w)))),atop(ds(CFLOOR),ds(CLOG)));
201 RZ(z=abase2(reshape(increm(d),num[2]),w));
202 R t&FL&&equ(irs1(z,0L,1L,jthead),lt(w,zero))?irs1(z,0L,1L,jtbehead):z;
203 }
204 c=x=0; v=AV(w);
205 DO(n, p=*v++; if(p==IMIN){c=SY_64?64:32; break;} x=x<p?p:x<-p?-p:x;);
206 if(!c)while(x){x>>=1; ++c;}
207 c=MAX(1,c);
208 GA(z,B01,n*c,1+r,AS(w)); *(r+AS(z))=c;
209 v=n+AV(w); zv=AN(z)+BAV(z);
210 DO(n, x=*--v; DO(c, *--zv=(B)(x&1); x>>=1;));
211 R z;
212 }
213
214 F2(jtabase2){A z;I an,ar,at,wn,wr,wt,zn;
215 RZ(a&&w);
216 an=AN(a); ar=AR(a); at=AT(a);
217 wn=AN(w); wr=AR(w); wt=AT(w);
218 ASSERT(at&DENSE&&wt&DENSE,EVNONCE);
219 if(!ar)R residue(a,w);
220 if(1==ar&&at&B01+INT&&wt&B01+INT){I*av,d,r,*u,*wv,x,*zv;
221 RZ(coerce2(&a,&w,INT));
222 RE(zn=mult(an,wn)); GA(z,INT,zn,1+wr,AS(w)); *(wr+AS(z))=an;
223 av=an+AV(a); wv=wn+AV(w); zv=zn+AV(z);
224 if(2==an&&!av[-2]&&0<(d=av[-1])){I d1,j,k;
225 k=0; j=1; while(d>j){++k; j<<=1;} d1=d-1;
226 if(d==j)DO(wn, x=*--wv; *--zv=x&d1; *--zv=x>>k;)
227 else DO(wn, x=*--wv; if(0<=x){*--zv=x%d; *--zv=x/d;}else{*--zv=d+x%d; *--zv=-1+x/d;})
228 }else DO(wn, x=*--wv; u=av; DO(an, d=*--u; *--zv=r=remii(d,x); x=d?(x-r)/d:0;););
229 R z;
230 }else{PROLOG;A y,*zv;C*u,*yv;I k;
231 F2RANK(1,0,jtabase2,0);
232 k=bp(at); u=an*k+CAV(a);
233 GA(y,at, 1, 0,0); yv=CAV(y);
234 GA(z,BOX,an,1,0); zv=an+AAV(z);
235 DO(an, MC(yv,u-=k,k); RZ(w=divide(minus(w,*--zv=residue(y,w)),y)););
236 EPILOG(ope(z));
237 }}
238
239 F1(jtintmod2){A z;B*b,*v;I k=SZI,mask,m,n,q,r,*u,*wi;
240 RZ(w);
241 n=AN(w); q=n/k; r=n%k; v=BAV(w)+!liln*(k-1);
242 GA(z,B01,n,AR(w),AS(w)); u=AV(z);
243 b=(B*)&mask; DO(k, b[i]=1;);
244 b=(B*)&m; DO(q, DO(k, b[i]=*v; v+=k;); *u++=mask&m;)
245 b=(B*)u; wi=AV(w)+q*k; DO(r, *b++=1&*wi++?1:0;);
246 R z;
247 }