Mercurial > hg > jgplsrc

comparison vbang.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: ! */
5
6 #include "j.h"
7
8
9 static Z z1={1, 0};
10
11 static D coeff[]={
12 0.0, 1.0, 0.5772156649015329, -0.6558780715202538,
13 -0.0420026350340952, 0.1665386113822915, -0.0421977345555443,
14 -0.009621971527877, 0.007218943246663, -0.0011651675918591,
15 -0.0002152416741149, 0.0001280502823882, -0.0000201348547807,
16 -0.0000012504934821, 0.000001133027232, -0.0000002056338417,
17 6.116095e-9, 5.0020075e-9, -1.1812746e-9,
18 1.043427e-10, 7.7823e-12, -3.6968e-12,
19 5.1e-13, -2.06e-14, -5.4e-15,
20 1.4e-15, 1.0e-16
21 };
22
23 static I terms=sizeof(coeff)/sizeof(D);
24
25 static Z jtzhorner(J jt,I n,D*c,Z v){Z s;D*d=n+c;
26 s=zeroZ;
27 DO(n, s=zplus(zrj0(*--d),ztymes(v,s)););
28 R s;
29 }
30
31 static D dgps(D v){D*d=terms+coeff,s=0.0; DO(terms, s=*--d+v*s;); R 1/s;}
32 /* Abramowitz & Stegun, 6.1.34 */
33
34 static Z jtzgps(J jt,Z z){R zdiv(z1,zhorner(terms,coeff,z));}
35
36 static D jtdgamma(J jt,D x){B b;D t;
37 t=1.0; b=x==jfloor(x);
38 if(b&&0>=x){ASSERT(x>x-1,EVLIMIT); R x==2*jfloor(x/2)?inf:infm;}
39 if(0<=x) while(1<x){t*=--x; if(t==inf)R inf;}
40 else {while(0>x){t*=x++; if(t==inf)R 0.0;} t=1.0/t;}
41 R b?t:t*dgps(x);
42 } /* gamma(x) using recurrence formula */
43
44 static Z jtzgrecur(J jt,Z z){Z t;
45 t=z1;
46 if(0<=z.re) while( 0.5<z.re){--z.re; t=ztymes(t,z); if(t.re==inf)R t; }
47 else {while(-0.5>z.re){t=ztymes(t,z); ++z.re; if(t.re==inf)R zeroZ;} t=zdiv(z1,t);}
48 R ztymes(t,zgps(z));
49 } /* gamma(z) using recurrence formula */
50
51 static Z jtzgauss(J jt,D n,Z z){D d=1/n;Z p,t;
52 if(1>=n)R zgrecur(z);
53 p=ztymes(zpow(zrj0(2*PI),zrj0((1-n)/2)),zpow(zrj0(n),zminus(z,zrj0(0.5))));
54 t=zdiv(z,zrj0(n));
55 DO((I)n, p=ztymes(p,zgrecur(t)); t.re+=d;);
56 R p;
57 } /* Abramowitz & Stegun, 6.1.20 */
58
59 static Z jtzstirling(J jt,Z z){Z p,q;
60 static D c[]={1.0, 1.0/12, 1.0/288, -139.0/51840, -571.0/2488320},
61 e=2.718281828459045235360287;
62 p=ztymes(zsqrt(zdiv(zrj0(2*PI),z)),zpow(zdiv(z,zrj0(e)),z));
63 q=zhorner(5L,c,zdiv(z1,z));
64 R ztymes(p,q);
65 } /* Abramowitz & Stegun, 6.1.37 */
66
67 static Z jtzgamma(J jt,Z z){D y=ABS(z.im);
68 R !y?zrj0(dgamma(z.re)):20<y?zstirling(z):zgauss(ceil(y/0.8660254),z);
69 }
70
71 AMON(factI, D,I, *z=dgamma(1.0+(D)*x);)
72 AMON(factD, D,D, *z=_isnan(*x)?*x:dgamma(1.0+*x);)
73 AMON(factZ, Z,Z, *z=zgamma(zplus(z1,*x));)
74
75
76 #define PQLOOP(expr) while(n&&h&&h!=inf&&h!=infm){h*=expr; --n;}
77
78 static D pq(D h,D m,D*c,D*d){D x=*c,y=*d;I n=(I)MIN(m,IMAX);
79 if(0>=m)R h;
80 switch(2*(0>x)+(0>y)){
81 case 0: if(x!= y)PQLOOP(x--/y--); break;
82 case 1: if(x!=-y)PQLOOP(x--/y++)else if(m>2*jfloor(0.5*m))h=-h; break;
83 case 2: if(x!=-y)PQLOOP(x++/y--)else if(m>2*jfloor(0.5*m))h=-h; break;
84 case 3: if(x!= y)PQLOOP(x++/y++); break;
85 }
86 if(0>=*c)*c+=m; else *c-=m;
87 if(0>=*d)*d+=m; else *d-=m;
88 R h;
89 }
90
91 static I signf(D x){R 0<=x||1<=x-2*jfloor(0.5*x)?1:-1;}
92 /* sign of !x */
93
94 static D jtdbin(J jt,D x,D y){D c,d,e,h=1.0,p,q,r;I k=0;
95 c=y; if(0<=c)p=jfloor(c); else{k+=4; ++c; p=jfloor(-c);}
96 d=y-x; if(0<=d)q=jfloor(d); else{k+=2; ++d; q=jfloor(-d);}
97 e=x; if(0<=e)r=jfloor(e); else{k+=1; ++e; r=jfloor(-e);}
98 switch(k){
99 case 0: h=pq(h,q,&c,&d); h=pq(h,r,&c,&e); break;
100 case 1: h=pq(h,p,&c,&d); h=pq(h,r,&e,&d); --e; break;
101 case 2: h=pq(h,p,&c,&e); h=pq(h,q,&d,&e); --d; break;
102 case 5: h=pq(h,p,&e,&c); h=pq(h,q,&e,&d); --c; --e; break;
103 case 6: h=pq(h,p,&d,&c); h=pq(h,r,&d,&e); --c; --d; break;
104 case 7: h=pq(h,q,&d,&c); h=pq(h,r,&e,&c); --c; --d; --e; break;
105 }
106 if(!h)R 0;
107 if(h==inf||h==infm)R inf*signf(x)*signf(y)*signf(y-x);
108 R h*dgamma(1+c)/(dgamma(1+d)*dgamma(1+e));
109 } /* x and y-x are not negative integers */
110
111 static D ibin(D x,D y){D d=MIN(x,y-x),p=1;
112 DO((I)d, p*=y--/d--; if(p==inf)R p;);
113 R jfloor(0.5+p);
114 } /* x and y are non-negative integers; x<=y */
115
116 static Z jtzbin(J jt,Z x,Z y){Z a,b,c;
117 a=zgamma(zplus(z1,y));
118 b=zgamma(zplus(z1,x));
119 c=zgamma(zplus(z1,zminus(y,x)));
120 R zdiv(a,ztymes(b,c));
121 }
122
123 #define MOD2(x) ((x)-2*jfloor(0.5*(x)))
124
125 static D jtbindd(J jt,D x,D y){B id,ix,iy;D d;
126 if(_isnan(x))R x; else if(_isnan(y))R y;
127 d=y-x;
128 id=d==jfloor(d);
129 ix=x==jfloor(x);
130 iy=y==jfloor(y);
131 switch(4*(ix&&0>x)+2*(iy&&0>y)+(id&&0>d)){
132 default: ASSERTSYS(0,"bindd");
133 case 5: /* 1 0 1 */ /* Impossible */
134 case 0: /* 0 0 0 */
135 case 2: /* 0 1 0 */ R ix&&iy?ibin(x,y):dbin(x,y);
136 case 3: /* 0 1 1 */ R (MOD2(x)?-1:1)*ibin(x,x-y-1);
137 case 6: /* 1 1 0 */ R (MOD2(d)?-1:1)*ibin(-1-y,-1-x);
138 case 1: /* 0 0 1 */
139 case 4: /* 1 0 0 */
140 case 7: /* 1 1 1 */ R 0;
141 }} /* P.C. Berry, Sharp APL Reference Manual, 1979, p. 132 */
142
143 static Z jtbinzz(J jt,Z x,Z y){B id,ix,iy;D rd,rx,ry;Z d;
144 if(!x.im&&!y.im)R zrj0(bindd(x.re,y.re));
145 d=zminus(y,x);
146 rd=d.re; id=rd==jfloor(rd)&&0==d.im;
147 rx=x.re; ix=rx==jfloor(rx)&&0==x.im;
148 ry=y.re; iy=ry==jfloor(ry)&&0==y.im;
149 switch(4*(ix&&0>rx)+2*(iy&&0>ry)+(id&&0>rd)){
150 default: ZASSERT(0,EVSYSTEM);
151 case 5: /* 1 0 1 */ /* Impossible */
152 case 0: /* 0 0 0 */
153 case 2: /* 0 1 0 */ R zbin(x,y);
154 case 3: /* 0 1 1 */ R zrj0((MOD2(rx)?-1:1)*ibin(rx,rx-ry-1));
155 case 6: /* 1 1 0 */ R zrj0((MOD2(rd)?-1:1)*ibin(-1-ry,-1-rx));
156 case 1: /* 0 0 1 */
157 case 4: /* 1 0 0 */
158 case 7: /* 1 1 1 */ R zeroZ;
159 }}
160
161
162 ANAN(binDD, D,D,D, bindd)
163 ANAN(binZZ, Z,Z,Z, binzz)