/* 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 */ #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) #ifndef WORDS #define PROD_OPEN printf("Prod(") #define PROD_SEP printf(",") #define PROD_CLOSE printf(")") #define EPSILON printf("Epsilon") #else #define PROD_OPEN {} #define PROD_SEP {} #define PROD_CLOSE {} #define EPSILON {} #endif DPE s,tmp; DPE *T4,*countT4(int); void drawT4(int); DPE *T,*countT(int); void drawT(int); DPE *T5,*countT5(int); void drawT5(int); DPE *T1,*countT1(int); void drawT1(int); DPE *Z,*countZ(int); void drawZ(int); DPE *T3,*countT3(int); void drawT3(int); DPE *T2,*countT2(int); void drawT2(int); main(int argc,char *argv[]) { int i,n,k; if (argc<=1) { printf("Usage: a.out n [k seed] to generate k objects of size n\n"); exit(1);} n=atoi(argv[1]); k=(argc>=3) ? atoi(argv[2]) : 1; srand48((argc>=4) ? atoi(argv[3]) : getpid()); DPEinit(); T4=(DPE*)malloc((n+1)*sizeof(DPE)); for (i=0;i<1 && i<=n;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<=n;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<=n;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<=n;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<=n;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<=n;i++) DPEset(T2+i,0.0); for (;i<=n;i++) T2[i].lo=Nan; T5=T4; countT(n); if (T[n].lo<=0.0) {printf("there is no such structure of this size\n"); exit(1);} if (k==0) DPEprint(T+n); for (i=1;i<=k;i++) { drawT(n); putchar('\n'); } } /* 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) */ void drawT4(int n) { int k,c; 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) {fprintf(stderr,"FAIL\n"); exit(1);} 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) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) break; k=n+1-k; } drawT(k); if ((n-k)>0) {putchar(','); drawT3(n-k);} } /* T = Prod(Z, T1) */ void drawT(int n) { int k,c; PROD_OPEN; drawZ(1); PROD_SEP; drawT1(n-1); PROD_CLOSE; } /* T5 = ProdSeq(T, T3) */ void drawT5(int n) { int k,c; 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) {fprintf(stderr,"FAIL\n"); exit(1);} 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) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) break; k=n+1-k; } printf("Sequence("); drawT(k); if ((n-k)>0) {putchar(','); drawT3(n-k);} putchar(')'); } /* T1 = Union(T2, T5) */ void 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) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) { drawT2(n); return; } else { DPEcopy(&tmp,countT5(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T1+n))==0) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) { drawT5(n); return; } else {fprintf(stderr,"FAIL\n"); exit(1);} } } /* Z = Atom(Z) */ void drawZ(int n) { int k,c; printf("Z"); } /* T3 = Union(T2, T4) */ void 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) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) { drawT2(n); return; } else { DPEcopy(&tmp,countT4(n)); DPEadd(&s,&tmp); if ((c=DPEcmp(&s,T3+n))==0) {fprintf(stderr,"FAIL\n"); exit(1);} else if (c==1) { drawT4(n); return; } else {fprintf(stderr,"FAIL\n"); exit(1);} } } /* T2 = Epsilon(T2) */ void drawT2(int n) { int k,c; EPSILON; }