/* This (DPE) C source file was automatically generated by CS version 1.4, with cs::HOLONOMIC = FALSE and cs::LAZY = FALSE, from the following specification: {T = Prod(Z, Sequence(T))}, in the standard form: {T4 = ProdPrdSeq(T, T3), T = Prod(Z, T1), T5 = ProdSeq(T, T3), T1 = Union(T2, T5), Z = Atom(Z), T3 = Union(T2, T4), T2 = Epsilon(T2)}. You can get CS from http://dept-info.labri.u-bordeaux.fr/~dutour/CS. Please report any bug, improvement to the authors of CS: cs@labri.u-bordeaux.fr Sylvie Corteel - LRI Alain Denise - LRI Isabelle Dutour - LaBRI Paul Zimmermann - INRIA */ /* This is a C Dynamic Module for MuPAD */ #include #include #include #ifdef SUN4 #include #define TOINFP fpsetround(FP_RP) #define TOINFM fpsetround(FP_RM) #define Double double #elif IRIX64 #include #define TOINFP swapRM(ROUND_TO_PLUS_INFINITY) #define TOINFM swapRM(ROUND_TO_MINUS_INFINITY) #define Double double #elif ALPHA #include #define TOINFP write_rnd(FP_RND_RP) #define TOINFM write_rnd(FP_RND_RM) #define Double volatile double #elif AIX #include #define TOINFP fp_swap_rnd(FP_RND_RP) #define TOINFM fp_swap_rnd(FP_RND_RM) #define Double double #elif SUNOS #include char *out; #define TOINFP ieee_flags("set","direction","positive",&out) #define TOINFM ieee_flags("set","direction","negative",&out) #define Double double #elif HP700 #define TOINFP fpsetround(FP_RP) #define TOINFM fpsetround(FP_RM) #define Double double #elif LINUX #include #define TOINFP __setfpucw(0x1b7f) #define TOINFM __setfpucw(0x177f) #define Double double #endif typedef struct { Double lo,hi; int e; } DPE; #define MAXPREC 16 /* IEEE 754 mantissa length in digits */ Double eps; void DPEinit() { eps=1.0; while (eps/2.0 != 0.0) eps/=2.0; } void DPEprint(DPE *x) { printf("%1.16f*10^%d..%1.16f*10^%d\n",x->lo,x->e,x->hi,x->e); } /* x <- y */ void DPEset(DPE *x, double y) { x->lo=x->hi=y; x->e=0; } /* x <- y */ void DPEcopy(DPE *x, DPE *y) { x->lo=y->lo; x->hi=y->hi; x->e=y->e; } void DPEnormal(DPE *x) { if (fabs(x->hi)>=10.0) { while (fabs(x->hi)>=10.0) { TOINFM; x->lo /= 10; TOINFP; x->hi /= 10; x->e++; } } else while (x->lo!=0.0 && fabs(x->lo)<=1.0) { TOINFM; x->lo *= 10; TOINFP; x->hi *= 10; x->e--; } } /* x <- y*z, ONLY for y and z positive */ void DPEmul(DPE *x, DPE *y, DPE *z) { TOINFM; x->lo = y->lo * z->lo; TOINFP; x->hi = y->hi * z->hi; x->e = y->e + z->e; DPEnormal(x); } /* x <- y*r */ void DPEmulr(DPE *x, DPE *y, double r) { Double tmp; if (r<0) { TOINFM; tmp = x->lo; x->lo = r*x->hi; TOINFP; x->hi = r*tmp; } else { TOINFM; x->lo = r*y->lo; TOINFP; x->hi = r*y->hi; } x->e = y->e; DPEnormal(x); } /* x <- n*x */ void DPEmuli(DPE *x, long n) { Double tmp; if (n<0) { TOINFM; tmp = x->lo; x->lo = n*x->hi; TOINFP; x->hi = n*tmp; } else { TOINFM; x->lo *= n; TOINFP; x->hi *= n; } DPEnormal(x); } /* x <- x/n */ void DPEdivi(DPE *x, long n) { Double tmp; if (n<0) { TOINFM; tmp = x->lo; x->lo = x->hi/n; TOINFP; x->hi = tmp/n; } else { TOINFM; x->lo /= n; TOINFP; x->hi /= n; } DPEnormal(x); } /* x <- x+y, destructive on y */ void DPEadd(DPE *x, DPE *y) { int d=x->e - y->e; if (d>0) { if (d>MAXPREC) { y->lo=0.0; y->hi=eps; } /* y too small */ else while (d-->0) { TOINFM; y->lo /= 10.0; TOINFP; y->hi /= 10.0; } } else if (d<0) { if (d<-MAXPREC) /* x too small */ { x->lo=y->lo; x->hi=y->hi; y->lo=0.0; y->hi=eps;} else while (d++<0) { TOINFM; x->lo /= 10.0; TOINFP; x->hi /= 10.0; } x->e = y->e; } else {} /* d=0 */ TOINFM; x->lo += y->lo; TOINFP; x->hi += y->hi; DPEnormal(x); } int DPEcmp(DPE *x, DPE *y) /* supposes normalized and positive */ { int d=x->e - y->e; if (d>0) return(1); else if (d<0) return(-1); else /* d=0 */ if (x->lo>y->hi) return(1); else if (x->hilo) return(-1); else return(0); } extern double drand48(); #define Nan (-1.0) #define NOINIT(A) ((A)->lo==Nan) #define INIT(A) ((A)->lo!=Nan) DPE s,tmp; DPE *T4,*countT4(int); MTcell drawT4(int); DPE *T,*countT(int); MTcell drawT(int); DPE *T5,*countT5(int); MTcell drawT5(int); DPE *T1,*countT1(int); MTcell drawT1(int); DPE *Z,*countZ(int); MTcell drawZ(int); DPE *T3,*countT3(int); MTcell drawT3(int); DPE *T2,*countT2(int); MTcell drawT2(int); int Nmax = -1; #define ntT 0 #define ntT1 1 #define ntT2 2 #define ntT3 3 #define ntZ 4 #define ntT4 5 #define ntT5 6 DPE *(*count[])(int) = { countT, countT1, countT2, countT3, countZ, countT4, countT5 }; MTcell (*draw[])(int) = { drawT, drawT1, drawT2, drawT3, drawZ, drawT4, drawT5 }; MTcell DPE2list(DPE *x) { //MFeval(MFtext2expr("old_DIGITS:=DIGITS")); //MFeval(MFtext2expr("DIGITS:=16")); MTcell e = MFlong( x->e ); MTcell l = MFnewList(2); MTcell p = MFcall("_power",2,MFlong(10),MFcopy(e)); MFsetList(&l,0,MFcall("_mult",2,MFdouble( x->lo ),MFcopy(p))); MFsetList(&l,1,MFcall("_mult",2,MFdouble( x->hi ),MFcopy(p))); MFsig(l); //MFeval(MFtext2expr("DIGITS:=old_DIGITS")); //MFeval(MFtext2expr("unassign(old_DIGITS)")); return(l); } void init(int n) { int i; srand48(getpid()); DPEinit(); //MFprintf("init n=%d\n",n); if (n<1) n=1; T4=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<1;i++) DPEset(T4+i,0.0); for (;i<=n;i++) T4[i].lo=Nan; T=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<1;i++) DPEset(T+i,0.0); for (;i<=n;i++) T[i].lo=Nan; T1=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<0;i++) DPEset(T1+i,0.0); for (;i<=n;i++) T1[i].lo=Nan; Z=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<1;i++) DPEset(Z+i,0.0); for (;i<=n;i++) Z[i].lo=Nan; T3=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<0;i++) DPEset(T3+i,0.0); for (;i<=n;i++) T3[i].lo=Nan; T2=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<0;i++) DPEset(T2+i,0.0); for (;i<=n;i++) T2[i].lo=Nan; T5=T4; Nmax=n; } void reinit(int n) { int i; //MFprintf("reinit n=%d\n",n); T4=(DPE*)realloc(T4, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) T4[i].lo=Nan; T=(DPE*)realloc(T, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) T[i].lo=Nan; T1=(DPE*)realloc(T1, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) T1[i].lo=Nan; Z=(DPE*)realloc(Z, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) Z[i].lo=Nan; T3=(DPE*)realloc(T3, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) T3[i].lo=Nan; T2=(DPE*)realloc(T2, (n+1)*sizeof(DPE)); for (i=Nmax+1;i<=n;i++) T2[i].lo=Nan; T5=T4; Nmax=n; } void countint(int i, int n) { if (n<0) MFerror("Bad size\n"); //MFprintf("Nmax=%d n=%d\n",Nmax,n); if (Nmax==-1) init(n); if (Nmaxlo<=0.0) /* se bloque quand on fait B[n].lo !! */ MFerror("There is no such structure of this size\n"); MFreturn( MFeval( draw[i](n) ) ); } MFEND /* T4 = ProdPrdSeq(T, T3) */ DPE* countT4(int n) { int k; if (NOINIT(T4+n)) { DPEmul(T4+n,countT(1),countT3(n-1)); for (k=2;k<=n;k++) { DPEmul(&tmp,countT(k),countT3(n-k)); DPEadd(T4+n,&tmp); } } return(T4+n); } /* T = Prod(Z, T1) */ DPE* countT(int n) { int k; if (NOINIT(T+n)) { DPEmul(T+n,countZ(1),countT1(n-1)); } return(T+n); } /* T5 = ProdSeq(T, T3) */ DPE* countT5(int n) { int k; return(countT4(n)); } /* T1 = Union(T2, T5) */ DPE* countT1(int n) { int k; if (NOINIT(T1+n)) { DPEcopy(T1+n,countT2(n)); DPEcopy(&tmp,countT5(n)); DPEadd(T1+n,&tmp); } return(T1+n); } /* Z = Atom(Z) */ DPE* countZ(int n) { int k; if (NOINIT(Z+n)) { DPEset(Z+n,(n==1) ? 1.0 : 0.0); } return(Z+n); } /* T3 = Union(T2, T4) */ DPE* countT3(int n) { int k; if (NOINIT(T3+n)) { DPEcopy(T3+n,countT2(n)); DPEcopy(&tmp,countT4(n)); DPEadd(T3+n,&tmp); } return(T3+n); } /* T2 = Epsilon(T2) */ DPE* countT2(int n) { int k; if (NOINIT(T2+n)) { DPEset(T2+n,(n==0) ? 1.0 : 0.0); } return(T2+n); } /* T4 = ProdPrdSeq(T, T3) */ MTcell drawT4(int n) { int k,c; MTcell res; DPEmulr(&s,countT4(n),drand48()); for (k=1;;k++) { DPEmul(&tmp,countT(k),countT3(n-k)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T4+n))==0) {MFerror("FAIL\n");} else if (c==1) break; k=n+1-k; DPEmul(&tmp,countT(k),countT3(n-k)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T4+n))==0) {MFerror("FAIL\n");} else if (c==1) break; k=n+1-k; } if ((n-k)>0) res=MFnewExprSeq(2, drawT(k), drawT3(n-k) ); else res=drawT(k); MFsig( res ); return( res ); } /* T = Prod(Z, T1) */ MTcell drawT(int n) { int k,c; MTcell res; res=MFnewExpr(3, MFident("Prod"), drawZ(1), drawT1(n-1) ); MFsig( res ); return( res ); } /* T5 = ProdSeq(T, T3) */ MTcell drawT5(int n) { int k,c; MTcell res; DPEmulr(&s,countT5(n),drand48()); for (k=1;;k++) { DPEmul(&tmp,countT(k),countT3(n-k)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T5+n))==0) {MFerror("FAIL\n");} else if (c==1) break; k=n+1-k; DPEmul(&tmp,countT(k),countT3(n-k)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T5+n))==0) {MFerror("FAIL\n");} else if (c==1) break; k=n+1-k; } if ((n-k)>0) res=MFnewExpr(3, MFident("Sequence"), drawT(k), drawT3(n-k) ); else res=MFnewExpr(2, MFident("Sequence"), drawT(k) ); MFsig( res ); return( res ); } /* T1 = Union(T2, T5) */ MTcell drawT1(int n) { int k,c; DPEmulr(&s,countT1(n),drand48()); DPEcopy(&tmp,countT2(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T1+n))==0) {MFerror("FAIL\n");} else if (c==1) { return( drawT2(n) ); } else { DPEcopy(&tmp,countT5(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T1+n))==0) {MFerror("FAIL\n");} else if (c==1) { return( drawT5(n) ); } else {MFerror("FAIL\n");} } } /* Z = Atom(Z) */ MTcell drawZ(int n) { int k,c; return( MFident("Z") ); } /* T3 = Union(T2, T4) */ MTcell drawT3(int n) { int k,c; DPEmulr(&s,countT3(n),drand48()); DPEcopy(&tmp,countT2(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T3+n))==0) {MFerror("FAIL\n");} else if (c==1) { return( drawT2(n) ); } else { DPEcopy(&tmp,countT4(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T3+n))==0) {MFerror("FAIL\n");} else if (c==1) { return( drawT4(n) ); } else {MFerror("FAIL\n");} } } /* T2 = Epsilon(T2) */ MTcell drawT2(int n) { int k,c; return( MFident("Epsilon") ); }