/* File: count.c A program for counting specific spanning subgraphs of G x P_n. Copyright by F.J. Faase */ #define VERSION "960906_gcc" /* Version history: 960906_gcc: implemented complete 3rd method 951122_gcc: take solve out of calc_req 951031_gcc: make count longword, replace 5*RADIX by RADIX2 same for orderd_matrix 951027_gcc: Option /st did not work proper any more. Modified the limit for `Error: overflow in counts' to 5*RADIX 951004_gcc: Changed name giving of types and variables accoring paper. 950930_gcc: More general solution for PROP_FT_PATH implemented. 950929_gcc: allowing /D=0,.. with /C PROP_FT_PATH toegevoegd 950906_gcc: replace s < 0 by i < 0 in ASSERT of num() added init_counting(0) in analyze_all_freq 950501_gcc: After sweeping, scan all solutions from 2 to max 950428_gcc: new algorithm, using matrix 950427_gcc: put $ around C_..^.. construct in output. bug: 2735 wrong index for M crash for /2f P_6.grf 950423_gcc: Allocatie van m[] meer dan 64K. 950422_gcc: resize did not correctly initialize the newly allocated memory. malloc probably initializes it to zero in most of the cases. 950418a_gcc: replaced problem_kind != PROP_2_FAC by constr_con in pre_next_step. bugs: strange test in resize */ #define SUN #define LONG_INT_FOR_SCANF #define PROTO(X) X #define CHECK #include #include #include #include #ifndef SUN #include #include #endif /*** DEBUGING OPTIONS ***/ #define NO_DEBUG #define NO_CHECK #define NO_NOISE 255 FILE *f_error = NULL; #ifdef DEBUG #define DEBUG_(S) S #define DEBUG_P(S) printf S #define DEBUG_N(S) PRINT_NUM_M(stdout, S) #define DEBUG_NP(S) PRINT_NUM_MP(stdout, S) #else #define DEBUG_(S) #define DEBUG_P(S) #define DEBUG_N(S) #define DEBUG_NP(S) #endif #define _DEBUG_(S) S #define _DEBUG_P(S) printf S #define _DEBUG_N(S) PRINT_NUM_M(stdout, S) #define _DEBUG_NP(S) PRINT_NUM_MP(stdout, S) #define NO_DEBUG_(S) #define NO_DEBUG_P(S) #define NO_DEBUG_N(S) #define NO_DEBUG_NP(S) #define ASSERT(C,M) if (C) { printf(M); fprintf(f_error, "%%");\ fprintf(f_error, M); fflush(f_error); { int *v = 0; *v = 0; }} #define ASSERT1(C,M,A1) if (C) { printf(M, A1); fprintf(f_error, "%%");\ fprintf(f_error, M, A1); fflush(f_error); { int *v = 0; *v = 0; }} #define ASSERT2(C,M,A1,A2) if (C) { printf(M, A1, A2); fprintf(f_error, "%%");\ fprintf(f_error, M, A1, A2); fflush(f_error); { int *v = 0; *v = 0; }} /*** ELEMENTARY TYPES ***/ typedef short bool; #define FALSE ((bool)0) #define TRUE ((bool)1) typedef unsigned char byte; typedef unsigned short int word; typedef unsigned long int longword; /*** MALLOC ***/ #define MALLOC(t) (t *)safe_malloc(sizeof(t)) #define CHAR_MALLOC(size) (char *)safe_malloc(size) #define ARRAY_MALLOC(t,size) (t *)safe_malloc(sizeof(t)*(size)) PROTO(void *safe_malloc(longword);) void *safe_malloc(size) longword size; { void *h = malloc(size); DEBUG_P(("safe_malloc(%ld)\n", (longword)size)); if (h == NULL) { printf("FATAL: Cannot allocate %ld bytes of dynamic memory\n", (longword)size); printf(" Program aborted.\n"); exit(1); } #ifdef NOISE { int i; char *s = h; for (i = 0; i < size; i++) *(s++) = NOISE; } #endif return h; } /*** INPUT GRAPH ***/ /* VERTICES INPUT GRAPH */ typedef short vertex_t; #define MAX_NR_VERTICES 10 static vertex_t nr_vertices; /* EDGES OF INPUT GRAPH */ #define MAX_NR_EDGES 40 static short nr_edges; vertex_t edges[MAX_NR_EDGES][2]; PROTO(bool read_graph(char *, FILE *);) bool read_graph(fn, out) char *fn; FILE *out; { FILE *f; vertex_t p1, p2; #ifdef LONG_INT_FOR_SCANF long int vnr_vertices, vp1, vp2; #endif nr_vertices = 0; nr_edges = 0; if ((f = fopen(fn, "r")) == NULL) return FALSE; if (out) fprintf(out, "Reading graph from file '%s'\n", fn); #ifdef LONG_INT_FOR_SCANF if (fscanf(f, " %ld", &vnr_vertices) != 1) #else if (fscanf(f, " %d", &nr_vertices) != 1) #endif { fprintf(stderr, "ERROR: Number of vertices not given"); if (out && out != stdout) fprintf(stderr, "ERROR: Number of vertices not given"); return FALSE; } #ifdef LONG_INT_FOR_SCANF nr_vertices = (vertex_t)vnr_vertices; #endif if (nr_vertices <= 0 || nr_vertices > MAX_NR_VERTICES) { fprintf(stderr, "ERROR: illegale number of vertices (%d)\n", nr_vertices); if (out && out != stdout) fprintf(out, "ERROR: illegale number of vertices (%d)\n", nr_vertices); nr_vertices = 0; return FALSE; } #ifdef LONG_INT_FOR_SCANF while (fscanf(f, " %ld %ld", &vp1, &vp2) == 2) { p1 = (vertex_t)vp1; p2 = (vertex_t)vp2; #else while (fscanf(f, " %d %d", &p1, &p2) == 2) { #endif NO_DEBUG_P(("%d %d\n",p1, p2)); if ( p1 < 0 || p1 >= nr_vertices || p2 < 0 || p2 >= nr_vertices || p1 == p2) { fprintf(stderr, "ERROR: edge (%d,%d) is ignored.\n", p1, p2); if (out && out != stdout) fprintf(out, "ERROR: edge (%d,%d) is ignored.\n", p1, p2); } else if (nr_edges >= MAX_NR_EDGES) { fprintf(stderr, "ERROR: too many edges.\n"); if (out && out != stdout) fprintf(out, "ERROR: too many edges.\n"); } else { edges[nr_edges][0] = p1; edges[nr_edges][1] = p2; nr_edges++; } } if (out) { fprintf(out, "A graph on %d vertices with %d edges has been read.\n", nr_vertices, nr_edges); fprintf(out, "The edges are:"); { short i; for (i = 0; i < nr_edges; i++) fprintf(out, " (%d,%d)", edges[i][0], edges[i][1]); } fprintf(out, "\n\n"); } fclose(f); return TRUE; } PROTO(void print_graph(FILE *, char *, char *);) void print_graph(f, name, filename) FILE *f; char *name; char *filename; { #ifdef ECHO_LINE fprintf(f, "%d", nr_vertices); { short i; for (i = 0; i < nr_edges; i++) fprintf(f, " %d %d", edges[i][0], edges[i][1]); } #endif fprintf(f, "'%s' %s", filename, VERSION); fprintf(f, "\nWe count the number of %ss, for a graph on %d vertices,", name, nr_vertices); fprintf(f, "\nand with the edges:"); { short i; for (i = 0; i < nr_edges; i++) fprintf(f, " (%d,%d)", edges[i][0] + 1, edges[i][1] + 1); } } PROTO(bool exist_edge(vertex_t, vertex_t);) bool exist_edge(a, b) vertex_t a; vertex_t b; { short i; for (i = 0; i < nr_edges; i++) if ( (edges[i][0] == a && edges[i][1] == b) || (edges[i][0] == b && edges[i][1] == a)) return TRUE; return FALSE; } /*** PROPERTY SPECIFICATION ***/ #define PROP_UNDEFINED 0 #define PROP_HAMIL_CYCLE 1 #define PROP_DOMINO_T 2 #define PROP_2_FAC 3 #define PROP_HAMIL_PATH 4 #define PROP_SPAN_TREES 5 #define PROP_SPAN_TREES13 6 #define PROP_SPEC 7 #define MAX_PROP 7 struct { char *name; char *ext; char *code; } kinds[] = { { "", "", "" } , { "Hamilton cycle", "hc", "HC" } , { "domino tilling", "dt", "DT" } , { "2 factor", "2f", "2F" } , { "Hamilton path", "hp", "HP" } , { "spanning tree", "st", "ST" } , { "spanning tree with degree 1 and 3", "st13", "ST_{1,3}" } , { "specified", "sp", "" } }; word problem_kind = PROP_UNDEFINED; /* CONSTRAINTS ON DEGREE */ bool constr_d = FALSE; #define MAX_NR_D MAX_NR_VERTICES bool allowed_d[MAX_NR_D]; byte maximum_d; /* CONSTRAINTS ON CONNECTED AND ACYCLIC */ bool constr_con = FALSE; bool constr_acycl = FALSE; bool constr_con_acycl = FALSE; /* CONSTRAINTS ON FROM-TO PATHS */ bool constr_from = FALSE, constr_to = FALSE; int from_id_path = 0, to_id_path = 0; bool constr_to_path = FALSE; /*** AUTOMORPHIC PERMUTATIONS OF THE INPUT GRAPH ***/ typedef vertex_t perm_t[MAX_NR_VERTICES]; typedef struct perm_list_T { struct perm_list_T *next; perm_t perm; } perm_list_t; perm_list_t *automorphic_perms = NULL; PROTO(bool automorphic_perm(perm_t);) bool automorphic_perm(perm) perm_t perm; { vertex_t p1, p2; if (constr_from && from_id_path != perm[from_id_path]) return FALSE; if (constr_to && to_id_path != perm[to_id_path]) return FALSE; for (p1 = 0; p1 < nr_vertices - 1; p1++) for (p2 = p1 + 1; p2 < nr_vertices; p2++) if (exist_edge(p1, p2) != exist_edge(perm[p1], perm[p2])) return FALSE; return TRUE; } PROTO(void find_automorphic_perms(word);) void find_automorphic_perms(nr_vertices) word nr_vertices; { short times[MAX_NR_VERTICES-1]; perm_t perm; short i; /* Init: */ for (i = 0; i < nr_vertices - 1; i++) times[i] = nr_vertices - i - 1; for (i = 0; i < (short)nr_vertices; i++) perm[i] = i; for(;;) { /* Check permutation is an automorphic mapping perm: */ if (automorphic_perm(perm)) { perm_list_t *new = MALLOC(perm_list_t); new->next = automorphic_perms; memcpy(new->perm, perm, sizeof(perm_t)); automorphic_perms = new; } /* Next: */ for (i = nr_vertices - 2; i >= 0 && times[i] == 0; i--) { NO_DEBUG_P(("times[%d] = %d\n", i, times[i])); times[i] = nr_vertices - i - 1; } if (i < 0) break; times[i]--; NO_DEBUG_P(("times[%d] = %d\n", i, times[i])); { short j = nr_vertices - 1; do { word h = perm[i]; perm[i] = perm[j]; perm[j] = h; i++; j--; } while (i < j); } } } PROTO(void print_permutations(FILE *, perm_list_t *);) void print_permutations(f, perm_list) FILE *f; perm_list_t *perm_list; { fprintf(f, "\nAutomorphic permutations:\n"); while(perm_list) { short i; for (i = 0; i < nr_vertices; i++) fprintf(f, "%2d", perm_list->perm[i]); fprintf(f, "\n"); perm_list = perm_list->next; } fprintf(f, "\n"); } /*** PARTIAL CONNECTION CODING (PCC) ***/ typedef byte id_t; #define UNUSED_ID 255 #define USED_ID 254 typedef id_t pcc_t[MAX_NR_VERTICES]; #define CPY_STATE(s,d) memcpy(s, d, nr_vertices*sizeof(id_t)) #define CMP_STATE(s,d) memcmp(s, d, nr_vertices*sizeof(id_t)) PROTO(void print_pcc(FILE *, pcc_t);) void print_pcc(f, pcc) FILE *f; pcc_t pcc; { vertex_t i; for (i = 0; i < nr_vertices; i++) if (pcc[i] == USED_ID) fprintf(f, "#"); else if (pcc[i] == UNUSED_ID) fprintf(f, "_"); else fprintf(f, "%1d", (short)pcc[i]); } short print_at_pos = 0; FILE *trace_f = NULL; PROTO(void print_pcc_at(FILE *, pcc_t, short);) void print_pcc_at(f, pcc, pos) FILE *f; pcc_t pcc; short pos; { if (pos < print_at_pos) { fprintf(f, "\n"); print_at_pos = 0; } for (; pos > print_at_pos; print_at_pos++) { short i; for (i = 0; i <= nr_vertices; i++) fprintf(f, " "); } print_pcc(f, pcc); fprintf(f, " "); print_at_pos++; } PROTO(void empty_pcc(pcc_t);) void empty_pcc(pcc) pcc_t pcc; { vertex_t i; for (i = 0; i < nr_vertices; i++) pcc[i] = UNUSED_ID; } PROTO(void renumber_pcc(pcc_t);) void renumber_pcc(pcc) pcc_t pcc; { id_t old_id[MAX_NR_VERTICES]; short nr_old_id = 0; vertex_t i; NO_DEBUG_P(("renumber_pcc( ")); NO_DEBUG_(print_pcc(stdout, pcc)); NO_DEBUG_P((" ) into: ")); for (i = 0; i < nr_vertices; i++) { id_t id = pcc[i]; if (id != UNUSED_ID) { bool found = FALSE; short j; for (j = 0; !found && j < nr_old_id; j++) found = old_id[j] == id; if (found) j--; else old_id[nr_old_id++] = id; pcc[i] = (id_t)j; } } NO_DEBUG_(print_pcc(stdout, pcc)); NO_DEBUG_P(("\n")); } PROTO(void normalize_pcc(pcc_t, pcc_t);) void normalize_pcc(in, norm) pcc_t in; pcc_t norm; { perm_list_t *am_p = automorphic_perms; bool first = TRUE; NO_DEBUG_P(("normalize_pcc( ")); NO_DEBUG_(print_pcc(stdout, in)); NO_DEBUG_P((" ,)\n")); while(am_p) { pcc_t try; { vertex_t i; for (i = 0; i < nr_vertices; i++) try[i] = in[am_p->perm[i]]; } renumber_pcc(try); if (first || CMP_STATE(try, norm) < 0) CPY_STATE(norm, try); first = FALSE; am_p = am_p->next; } DEBUG_P(("normalize_pcc( ")); DEBUG_(print_pcc(stdout, in)); DEBUG_P((" , ")); DEBUG_(print_pcc(stdout, norm)); DEBUG_P((" )\n")); if (trace_f != NULL) print_pcc_at(trace_f, norm, nr_edges + nr_vertices + 4); } /*** DEFINITION BIG INTEGERS ***/ #define RDIGITS 4 /* = nr digits per ndigit */ #define RADIX 10000 #define RADIX2 100000000 #define NR_NDIGITS 16 typedef struct { byte size; /* * NR_DIGITD = nr ndigits allocated */ word len; /* = nr ndigits used */ char sign; word *ndigit; /* array of (size * NR_NDIGITS) ndigits */ } big_int; /*** DIRECTED MULTI-GRAPH ***/ struct arc_list_T; typedef struct m_vertex_list_T /* vertices of directed multi-graph */ { struct m_vertex_list_T *next; pcc_t pcc; big_int old_c, new_c; byte nr_ends, first, freq; struct arc_list_T *arc_list; } m_vertex_list_t, *m_vertex_list_p; typedef struct arc_list_T /* count of arcs between two vertices */ { struct arc_list_T *next; m_vertex_list_p m_vertex; longword count; } arc_list_t, *arc_list_p; m_vertex_list_p m_vertex_list = NULL; PROTO(m_vertex_list_p find_p_to_pcc(pcc_t);) m_vertex_list_p find_p_to_m_vertex(pcc) pcc_t pcc; { m_vertex_list_p *ref_m_vertex_list = &m_vertex_list; while ( *ref_m_vertex_list && CMP_STATE((*ref_m_vertex_list)->pcc, pcc)) ref_m_vertex_list = &(*ref_m_vertex_list)->next; if (*ref_m_vertex_list == NULL) { *ref_m_vertex_list = MALLOC(m_vertex_list_t); (*ref_m_vertex_list)->next = NULL; CPY_STATE((*ref_m_vertex_list)->pcc, pcc); (*ref_m_vertex_list)->arc_list = NULL; } return *ref_m_vertex_list; } PROTO(void add_arc(pcc_t, pcc_t);) void add_arc(from_pcc, to_pcc) pcc_t from_pcc; pcc_t to_pcc; { m_vertex_list_p from_m_vertex_p = find_p_to_m_vertex(from_pcc), to_m_vertex_p = find_p_to_m_vertex(to_pcc); arc_list_p *ref_arc_list = &to_m_vertex_p->arc_list; if (trace_f != NULL) fprintf(trace_f, " *"); DEBUG_P(("add_arc( ")); DEBUG_(print_pcc(stdout, from_pcc)); DEBUG_P((" , ")); DEBUG_(print_pcc(stdout, to_pcc)); DEBUG_P((" )")); while ( *ref_arc_list && (*ref_arc_list)->m_vertex != from_m_vertex_p) ref_arc_list = &(*ref_arc_list)->next; if (*ref_arc_list == NULL) { *ref_arc_list = MALLOC(arc_list_t); (*ref_arc_list)->next = NULL; (*ref_arc_list)->m_vertex = from_m_vertex_p; (*ref_arc_list)->count = 0; } if (++(*ref_arc_list)->count >= RADIX2) { printf("Error: overflow in counts\n"); exit(1); } DEBUG_P((" okay\n")); } PROTO(longword count_of_arc(m_vertex_list_p, m_vertex_list_p);) longword count_of_arc(from_m_vertex_p, to_m_vertex_p) m_vertex_list_p from_m_vertex_p; m_vertex_list_p to_m_vertex_p; { arc_list_p arc_list = to_m_vertex_p->arc_list; while ( arc_list && arc_list->m_vertex != from_m_vertex_p) arc_list = arc_list->next; return arc_list == NULL ? 0 : arc_list->count; } /*** FINDING ARCS OF THE DIRECTED MULTI-GRAPH ***/ PROTO(id_t unique_id_in_pcc(pcc_t);) id_t unique_id_in_pcc(pcc) pcc_t pcc; { vertex_t i; bool found = FALSE; id_t max; NO_DEBUG_P(("unique_id_in_pcc( ")); NO_DEBUG_(print_pcc(stdout, pcc)); for (i = 0; i < nr_vertices; i++) if ( pcc[i] != UNUSED_ID && pcc[i] != USED_ID && (!found || max < pcc[i])) { max = pcc[i]; found = TRUE; } NO_DEBUG_P((" ) : %d\n", found ? max + 1 : 0)); return (id_t)(found ? max + 1 : 0); } PROTO(bool is_ending(pcc_t, vertex_t);) bool is_ending(pcc, a) pcc_t pcc; vertex_t a; { vertex_t i; id_t id = pcc[a]; for (i = 0; i < nr_vertices; i++) if ( i != a && pcc[i] == id) return FALSE; return TRUE; } PROTO(byte nr_ends_in_pcc(pcc_t pcc);) byte nr_ends_in_pcc(pcc) pcc_t pcc; { short i; byte nr_ends = 0; for (i = 0; i < nr_vertices; i++) if (pcc[i] != UNUSED_ID && is_ending(pcc, i)) nr_ends++; return nr_ends; } /***************\ word nr_free(pcc_t pcc, id_t id, ds_t ds, byte maximal_d) { vertex_t i; word result = 0; for (i = 0; i < nr_vertices; i++) if (pcc[i] == id && ds[i] < maximal_d) result++; return result; } \***************/ PROTO(void change_id_into(pcc_t, id_t, id_t);) void change_id_into(pcc, from_id, to_id) pcc_t pcc; id_t from_id; id_t to_id; { vertex_t i; for (i = 0; i < nr_vertices; i++) if (pcc[i] == from_id) pcc[i] = to_id; } PROTO(bool all_others_used(pcc_t, vertex_t, vertex_t);) bool all_others_used(pcc, a, b) pcc_t pcc; vertex_t a; vertex_t b; { vertex_t i; for (i = 0; i < nr_vertices; i++) if ( i != a && i != b && pcc[i] != USED_ID) return FALSE; return TRUE; } /*** DEGREES (DS) FOR VERTICES ***/ typedef byte d_t; typedef d_t ds_t[MAX_NR_VERTICES]; #define CPY_DS(s,d) memcpy(s, d, nr_vertices*sizeof(id_t)) #define CMP_DS(s,d) memcmp(s, d, nr_vertices*sizeof(id_t)) PROTO(void print_ds(FILE *, ds_t);) void print_ds(f, ds) FILE *f; ds_t ds; { vertex_t i; for (i = 0; i < nr_vertices; i++) fprintf(f, "%1d", (short)ds[i]); } /*** GENERATING VERTICES AND ARCS OF DIRECTED MULTI-GRAPHS ***/ static pcc_t cur_start_pcc; static bool first_m_vertex; /* FINDING VERTICES WITH EDGE */ void fill_in_next_pcc(new_pcc, ds, p, nr_ends) pcc_t new_pcc; ds_t ds; vertex_t p; short nr_ends; { DEBUG_P(("fill_in_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); DEBUG_P((",")); DEBUG_(print_ds(stdout, ds)); DEBUG_P((" , %d, %d)\n", (short)p, (short)nr_ends)); if (trace_f != NULL) print_pcc_at(trace_f, new_pcc, nr_edges + p + 2); while ( p < nr_vertices && !( new_pcc[p] != UNUSED_ID && allowed_d[ds[p] - 1] /*! 950929 modified: */ /*! 940129 ! toegevoegd: */ && ( !constr_con || ds[p] == 1 || !is_ending(new_pcc, p)) && (problem_kind != PROP_HAMIL_PATH || nr_ends < 2))) { DEBUG_P(("%d: %d && %d && (%d || %d || %d) && (%d || %d)\n", p, new_pcc[p] != UNUSED_ID, allowed_d[ds[p] - 1], /*! 940129_1 ! toegevoegd: */ !constr_con, ds[p] == 1, !is_ending(new_pcc, p), problem_kind != PROP_HAMIL_PATH, nr_ends < 2)); p++; } if (p >= nr_vertices) { pcc_t norm; normalize_pcc(new_pcc, norm); DEBUG_P((" %d => %d %d, %d => %d %d, %d => %d\n", constr_from, first_pcc, ds[from_id_path] == 1, constr_to, CMP_STATE(new_pcc, m_vertex_list->pcc), ds[to_id_path] == 1, constr_to_path, nr_ends < 2)); if ( ( !constr_from || !first_m_vertex || ds[from_id_path] == 1) && ( !constr_to || CMP_STATE(new_pcc, m_vertex_list->pcc) || ds[to_id_path] == 1) && (!constr_to_path || nr_ends < 2)) add_arc(cur_start_pcc, norm); } else { fill_in_next_pcc(new_pcc, ds, p + 1, nr_ends); { id_t save_id = new_pcc[p]; new_pcc[p] = UNUSED_ID; ds[p]--; fill_in_next_pcc(new_pcc, ds, p + 1, (ds[p] > 0 ? nr_ends + 1 : nr_ends - 1)); ds[p]++; new_pcc[p] = save_id; } } DEBUG_P(("fill_in_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); DEBUG_P((",")); DEBUG_(print_ds(stdout, ds)); DEBUG_P((" , %d, %d) okay\n", (short)p, (short)nr_ends)); } void construct_next_pccs(pcc, ds) pcc_t pcc; ds_t ds; { vertex_t i; pcc_t new_pcc; ds_t new_ds; id_t nr_ends = 0; id_t new_id = constr_con_acycl ? unique_id_in_pcc(pcc) : (id_t)0; CPY_STATE(new_pcc, pcc); CPY_DS(new_ds, ds); DEBUG_P(("construct_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, new_ds)); DEBUG_P((" )\n")); if (constr_con) { bool correct = TRUE; /* 950929 special case /C /D=0,.. */ if (!allowed_d[0]) for (i = 0; correct && i < nr_vertices; i++) correct = allowed_d[new_ds[i]] && new_pcc[0] == new_pcc[i]; else { bool found = FALSE; id_t val; for (i = 0; correct && i < nr_vertices; i++) if (new_ds[i] != 0) if (!allowed_d[new_ds[i]]) correct = FALSE; else if (found) correct = new_pcc[i] == val; else { found = TRUE; val = new_pcc[i]; } correct = correct && found; } if (correct) { pcc_t end_pcc; empty_pcc(end_pcc); if (trace_f != NULL) print_pcc_at(trace_f, end_pcc, nr_edges + nr_vertices + 4); if ( ( !constr_from || !first_m_vertex || ds[from_id_path] == 1) && ( !constr_to || ds[to_id_path] == 1)) add_arc(cur_start_pcc, end_pcc); } } for (i = 0; i < nr_vertices; i++) { if (allowed_d[new_ds[i] + 1]) { new_ds[i]++; if (new_ds[i] == 1) new_pcc[i] = constr_con_acycl ? new_id++ : (id_t)0; } else if ( allowed_d[new_ds[i]] && (!constr_con || !is_ending(new_pcc, i))) new_pcc[i] = UNUSED_ID; else { return; DEBUG_P(("construct_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, new_ds)); DEBUG_P((" ) okay\n")); } } if (problem_kind == PROP_HAMIL_PATH) { nr_ends = nr_ends_in_pcc(new_pcc); if (nr_ends > 2) { return; DEBUG_P(("construct_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, new_ds)); DEBUG_P((" ) okay\n")); } } else if (constr_to_path) nr_ends = nr_ends_in_pcc(new_pcc); fill_in_next_pcc(new_pcc, new_ds, (vertex_t)0, nr_ends); DEBUG_P(("construct_next_pcc( ")); DEBUG_(print_pcc(stdout, new_pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, new_ds)); DEBUG_P((" ) okay\n")); } /* FINDING ARCS WITH A VERTEX */ void construct_arc(pcc, ds, edge_nr) pcc_t pcc; ds_t ds; short edge_nr; { NO_DEBUG_P(("construct_arc( ")); NO_DEBUG_(print_pcc(stdout, pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, ds)); NO_DEBUG_P((" , %d)\n", edge_nr)); if (trace_f != NULL) print_pcc_at(trace_f, pcc, edge_nr); if (edge_nr == nr_edges) construct_next_pccs(pcc, ds); else { vertex_t a = edges[edge_nr][0], b = edges[edge_nr][1]; edge_nr++; if ( ( !constr_d || (ds[a] < maximum_d && ds[b] < maximum_d)) && ( !constr_acycl || ds[a] == 0 || ds[b] == 0 || pcc[a] != pcc[b])) { pcc_t new_pcc; CPY_STATE(new_pcc, pcc); if (ds[a] == 0) if (ds[b] == 0) { id_t new_id = constr_con_acycl ? unique_id_in_pcc(pcc) : (id_t)0; new_pcc[a] = new_id; new_pcc[b] = new_id; } else new_pcc[a] = pcc[b]; else if (ds[b] == 0) new_pcc[b] = pcc[a]; else change_id_into(new_pcc, pcc[a], pcc[b]); ds[a]++; ds[b]++; /* Optimalisatie: if ( constr_con && nr_free(new_pcc, new_pcc[a], ds, maximum_d) == 0) construct_next_pcc(new_pcc, ds); else */ construct_arc(new_pcc, ds, edge_nr); ds[a]--; ds[b]--; } construct_arc(pcc, ds, edge_nr); } NO_DEBUG_P(("construct_arc( ")); NO_DEBUG_(print_pcc(stdout, pcc)); NO_DEBUG_P((",")); NO_DEBUG_(print_ds(stdout, ds)); NO_DEBUG_P((" , %d) okay\n", edge_nr)); } void find_m_vertices_from_m_vertex(pcc) pcc_t pcc; { ds_t ds; DEBUG_P(("find_m_vertices_from_m_vertex( ")); DEBUG_(print_pcc(stdout, pcc)); DEBUG_P((" )\n")); CPY_STATE(cur_start_pcc, pcc); { vertex_t i; for (i = 0; i < nr_vertices; i++) ds[i] = pcc[i] != UNUSED_ID ? (d_t)1 : (d_t)0; } construct_arc(pcc, ds, 0); DEBUG_P(("m_vertices_from_m_vertex( ")); DEBUG_(print_pcc(stdout, pcc)); DEBUG_P((" ) okay\n")); } /*** PRINT M VERTEX LIST ***/ void print_m_vertex_list(f) FILE *f; { m_vertex_list_p sl = m_vertex_list; DEBUG_P(("print_m_vertex_list()\n")); fprintf(f, "\nGenerated transition graph:\n"); while (sl) { print_pcc(f, sl->pcc); fprintf(f, " :"); DEBUG_P(("for: ")); DEBUG_(print_pcc(stdout, sl->pcc)); DEBUG_P(("\n")); { arc_list_p bl = sl->arc_list; while (bl) { fprintf(f, " "); print_pcc(f, bl->m_vertex->pcc); fprintf(f, " (%ld)", bl->count); bl = bl->next; } } fprintf(f, "\n"); sl = sl->next; } DEBUG_P(("print_m_vertex_list() okay\n")); } /*** BIG INT PROCEDURES ***/ #define SIZE_OF_LEN(S) ((S) <= 0L ? 0L : (((S) - 1L) / (long)NR_NDIGITS + 1)) /* CONSISTENCY CHECKING */ #ifdef CHECK #define ASS_NDIGIT(S,I,V) ass_ndigit(S, I, V) #define ASS_RADIX(S,V) ass_radix(&S, V) #define NDIGIT(S,I) ndigit(S,I) #define CHECK_BIG_INT(S) check_big_int(S) PROTO(void check_big_int(big_int *);) void check_big_int(s) big_int *s; { int i; ASSERT2(s->len != s->size * NR_NDIGITS, "ERROR: inconsist %d %d\n", s->len, s->size); ASSERT1(s->len > 10000, "ERROR: s->len > 10000: %d\n", s->len); for (i = 0; i < s->len; i++) ASSERT2(s->ndigit[i] > RADIX, "ERROR: s->ndigit[%d] = %d\n", i, s->ndigit[i]); } PROTO(void ass_ndigit(big_int *, short, long);) void ass_ndigit(s, i, v) big_int *s; short i; long v; { CHECK_BIG_INT(s); ASSERT1(i >= (short)s->len, "ERROR: (ass_ndigit) i too large , %d\n", i); ASSERT1(i < 0, "ERROR: (ass_ndigit) i negative, %d\n", i); ASSERT1(v >= (long)RADIX, "ERROR: (ass_ndigit) value too large, %ld\n", v); ASSERT1(v < 0L, "ERROR: (ass_ndigit) value negative, %ld\n", v); s->ndigit[i] = (word)v; } PROTO(void ass_radix(word *, long);) void ass_radix(d, v) word *d; long v; { ASSERT1(v >= (long)RADIX, "ERROR: (ass_radix) value too large, %ld\n", v); ASSERT1(v < 0L, "ERROR: (ass_radix) value negative, %ld\n", v); *d = (word)v; } PROTO(word ndigit(big_int *, short);) word ndigit(s, i) big_int *s; short i; { CHECK_BIG_INT(s); ASSERT1(i >= (short)s->len || i < 0, "ERROR: i out of range, %d\n", i); return s->ndigit[i]; } #else #define ASS_NDIGIT(S,I,V) S->ndigit[(word)(I)] = (word)(V) #define ASS_RADIX(D,V) D = (word)(V) #define NDIGIT(S,I) S->ndigit[(word)(I)] #define CHECK_BIG_INT(S) #endif /* AUXILARY BIG INT PROCEDURES */ PROTO(word len(big_int *);) word len(a) big_int *a; { word i; CHECK_BIG_INT(a); if (a->len == 0) return 0; for (i = a->len; i > 0 && a->ndigit[i - 1] == 0; i--); return i; } #define RESIZE(D,S) { if ((word)D->size < (S)) resize(D, S); } PROTO(void resize(big_int *, long);) void resize(d, size) big_int *d; long size; { word i, *old; word old_len; DEBUG_P(("resize(%d, %ld)\n", (short)d->size, size)); ASSERT(size > 250L || size * (long)NR_NDIGITS > RADIX, "ERROR: Number to large\n"); ASSERT(size < 0L, "ERROR: negative size\n"); CHECK_BIG_INT(d); if (size <= (long)d->size) { printf("ERROR: resize - reducing size\n"); return; } old_len = d->len; old = d->ndigit; d->size = (word)size; d->len = NR_NDIGITS * (word)size; d->ndigit = ARRAY_MALLOC(word, d->len); for (i = 0; i < d->len; i++) d->ndigit[i] = 0; /* memcpy(new, d->ndigit, sizeof(word) * len);*/ for (i = 0; i < old_len; i++) ASS_NDIGIT(d, i, old[i]); for (i = old_len; i < d->len; i++) ASS_NDIGIT(d, i, 0); CHECK_BIG_INT(d); free(old); } /* ELEMENTARY BIG INT PROCEDURES */ PROTO(void free_(big_int *);) void free_(d) big_int *d; { d->sign = 1; d->size = 0; d->len = 0; free(d->ndigit); d->ndigit = NULL; } PROTO(void init(big_int *);) void init(d) big_int *d; { DEBUG_P(("init\n")); d->sign = 1; d->size = 0; d->len = 0; d->ndigit = NULL; } void print_big_int(FILE *, big_int *, bool, bool); PROTO(void assign(big_int *, big_int *);) void assign(d, s) big_int *d, *s; { word i; CHECK_BIG_INT(d); CHECK_BIG_INT(s); RESIZE(d, s->size); /* memcpy(d->ndigit, s->ndigit, sizeof(word) * s->len); */ for (i = 0; i < s->len; i++) ASS_NDIGIT(d, i, NDIGIT(s, i)); for (i = s->len; i < d->len; i++) ASS_NDIGIT(d, i, 0); CHECK_BIG_INT(d); } PROTO(bool big_int_is_zero(big_int *);) bool big_int_is_zero(bi) big_int *bi; { word i; for (i = 0; i < bi->len; i++) if (NDIGIT(bi, i) != 0) return FALSE; return TRUE; } PROTO(void assign_word(big_int *, longword);) void assign_word(bi, w) big_int *bi; longword w; { word i; CHECK_BIG_INT(bi); DEBUG_P(("assign_word(bi, %d)\n", (short)w)); RESIZE(bi, 1); ASS_NDIGIT(bi, 0, w % RADIX); ASS_NDIGIT(bi, 1, w / RADIX); for (i = 2; i < bi->len; i++) ASS_NDIGIT(bi, i, 0); } void print_big_int(f, c, right_justified, abs) FILE *f; big_int *c; bool right_justified, abs; { short i; bool lead = TRUE; if (!abs) fprintf(f, "%c", c->sign == 1 ? '+' : '-'); for (i = c->len - 1; i >= 0; i--) if (lead) if (NDIGIT(c, i) == 0 && i != 0) { if (right_justified) fprintf(f, " "); } else { if (right_justified) fprintf(f, "%*d", RDIGITS , NDIGIT(c, i)); else fprintf(f, "%d", NDIGIT(c, i)); lead = FALSE; } else fprintf(f, "%0*d", RDIGITS, NDIGIT(c, i)); } /* ADD ABSOLUTE */ PROTO(void add_abs_times_to(big_int *, big_int *, word, word);) void add_abs_times_to(s, a, sh, t) big_int *s, *a; word sh, t; { long i; word carry = 0; long len_a, len_s, len_max; DEBUG_P(("add_abs_times_to(): ")); DEBUG_(print_big_int(stdout, s, FALSE, FALSE)); DEBUG_P((" + ")); DEBUG_(print_big_int(stdout, a, FALSE, FALSE)); DEBUG_P((" * %d * RADIX ^ %d\n", t, sh)); len_a = (long)len(a) + sh; len_s = (long)len(s); RESIZE(s, SIZE_OF_LEN(len_a)); len_max = len_a > len_s ? len_a : len_s; for (i = sh; i < len_max; i++) { longword temp; temp = (long)NDIGIT(s, i) + (long)carry + ((i - sh < len_a) ? (long)NDIGIT(a, i - sh) * (long)t : 0L); DEBUG_P((" %ld temp= %d, %ld", i, (short)NDIGIT(s, i), (long)temp)); ASSERT1(temp >= RADIX2, "ERROR: temp error1 %ld\n", temp); /* { printf("s->ndigit[%ld] = %ld\n", i, (long)NDIGIT(s, i)); printf("carry: %d\n", carry); printf("a->ndigit[%ld] = %ld\n", i - sh, (i - sh < len_a) ? (long)NDIGIT(a, i - sh) : 0L); printf("t = %ld\n", (long)t); printf("* = %ld\n", ((i - sh < len_a) ? (long)NDIGIT(a, i - sh) * (long)t : 0L)); exit(1); } */ ASS_NDIGIT(s, i, temp % (long)RADIX); ASS_RADIX(carry, temp / (long)RADIX); } ASSERT1(carry > RADIX, "ERROR: carry error2 %d\n", carry); if (carry > 0) { RESIZE(s, SIZE_OF_LEN(len_max + 1)); ASS_NDIGIT(s, i, carry); } DEBUG_P(("= ")); DEBUG_(print_big_int(stdout, s, FALSE, FALSE)); DEBUG_P(("\n")); } /* COMPARE ABSOLUTE */ PROTO(int comp_abs(big_int *, big_int *);) int comp_abs(a, b) big_int *a, *b; { long sa = a->len, sb = b->len, s = sb, i; if (sa < sb) { for (i = sb - 1; i >= sa; i--) if (NDIGIT(b, i) > 0) return -1; s = sa; } else if (sb < sa) { for (i = sa - 1; i >= sb; i--) if (NDIGIT(a, i) > 0) return 1; s = sb; } for (i = s - 1; i >= 0; i--) if (NDIGIT(a, i) == NDIGIT(b, i)) /* continue */; else if (NDIGIT(a, i) > NDIGIT(b, i)) return 1; else return -1; return 0; } /* COMPARE */ PROTO(int comp(big_int *, big_int *);) int comp(a, b) big_int *a, *b; { if (a->sign < b->sign) return -1; if (b->sign < a->sign) return 1; return (a->sign == -1) ? -comp_abs(a,b) : comp_abs(a,b); } /*** CYCLE FREQUENCY OF M VERTICES ***/ PROTO(word gcd(word, word);) word gcd(a, b) word a, b; { NO_DEBUG_P(("gcd(%d,%d) = ", a, b)); if (a == 0 || b == 0) return 0; while (a != b) if (a < b) b -= a; else a -= b; NO_DEBUG_P(("%d\n", a)); return a; } PROTO(bool analyze_freq(word);) bool analyze_freq(step) word step; { m_vertex_list_p sl = m_vertex_list; bool zero_freq_found = FALSE; while (sl != NULL) { if (!big_int_is_zero(&sl->new_c)) if (sl->first == 0) sl->first = (byte)step; else if (sl->freq == 0) sl->freq = (byte)step - sl->first; else sl->freq = (byte)gcd(sl->freq, step - sl->first); if (sl->freq == 0) zero_freq_found = TRUE; sl = sl->next; } return zero_freq_found; } bool freq_analyzed = FALSE; PROTO(void calculate_nr_ends(void);) void calculate_nr_ends() { m_vertex_list_p sl = m_vertex_list; NO_DEBUG_P(("start calculate_nr_ends\n")); while (sl != NULL) { sl->nr_ends = nr_ends_in_pcc(sl->pcc); sl = sl->next; } NO_DEBUG_P(("end calculate_nr_ends\n")); } void print_freq(f) FILE *f; { m_vertex_list_p sl = m_vertex_list; NO_DEBUG_P(("start print_freq\n")); fprintf(f, "\n\nFrequences\n"); { vertex_t i; for (i = 0; i < nr_vertices; i++) fprintf(f, " "); } fprintf(f," #ends first freq\n"); while (sl != NULL) { print_pcc(f, sl->pcc); sl->nr_ends = nr_ends_in_pcc(sl->pcc); fprintf(f, " %5d %5d %4d\n", sl->nr_ends, sl->first, sl->freq); sl = sl->next; } NO_DEBUG_P(("end print_freq\n")); } /*** CALCULATING VALUES FOR C(n) ***/ void check_all_counts() { m_vertex_list_p sl = m_vertex_list; while (sl != NULL) { CHECK_BIG_INT(&sl->old_c); CHECK_BIG_INT(&sl->new_c); sl = sl->next; } } PROTO(void init_counting(int);) void init_counting(p) int p; { m_vertex_list_p sl = m_vertex_list; int i = 0; while (sl != NULL) { init(&sl->old_c); assign_word(&sl->old_c, i == p ? 1 : 0); init(&sl->new_c); assign_word(&sl->new_c, 0); sl->first = 0; sl->freq = 0; sl = sl->next; i++; } } PROTO(void pre_next_step(void);) void pre_next_step() { m_vertex_list_p sl = m_vertex_list; if (constr_con) assign_word(&sl->old_c, 0); else assign(&sl->old_c, &sl->new_c); assign_word(&sl->new_c, 0); while ((sl = sl->next) != NULL) { assign(&sl->old_c, &sl->new_c); assign_word(&sl->new_c, 0); } } PROTO(void next_step_counting(void);) void next_step_counting() { m_vertex_list_p sl = m_vertex_list; while (sl != NULL) { arc_list_p bl = sl->arc_list; while (bl != NULL) { add_abs_times_to(&sl->new_c, &bl->m_vertex->old_c, 0, bl->count); #ifdef EXT_CHECK check_all_counts(); #endif bl = bl->next; } sl = sl->next; } } #define USE_BIG_INT #ifdef USE_DOUBLE typedef double num_t; #define GLOBAL(V) double V #define TO_V(S,I,D) S[i] = convert_large_nr_to_double(D) #define INIT(S) #define FREE(S) #define ASS(D,S) (D = S) #define ASS_WORD(D,S) (D = (double)S) #define ASS_COUNT(D,S) (D = convert_large_nr_to_double(S)) #define PRINT_NUM_S(F,C) fprintf(F,"%5.0f ", C) #define PRINT_NUM_M(F,C) fprintf(F,"%5.0f ", C) #define PRINT_NUM_MP(F,C) fprintf(F,"%5.0f ", C) #define IS_ZERO(D) (D == 0.0) #define IS_ONE(D) (D == 1.0) #define ABS_LT(A,B) (fabs(A) < fabs(B)) #define POS_SIGN(A) A > 0.0 #define EXCHANGE(A,B) { double h = A; A = B; B = h; } #define ASS_ZERO(A) A = 0.0 #define ASS_ONE(A) A = 1.0 #define ASS_ADD(S,A) S += A #define ASS_SUB(S,A) S -= A #define NEG(A) A = -A #define PROD(R,A,B) R = A * B #define DIV(R,A,B) (R = A / B, TRUE) #define GCD(R,A,B) R = A double convert_large_nr_to_double(v) big_int v; { double r = 0.0; short i; for (i = v.len - 1; i >= 0; i--) r = r * RADIX + (double)v.ndigit[i]; return r; } #endif #ifdef USE_LONG typedef long num_t; #define GLOBAL(V) long V #define TO_V(S,I,D) S[I] = (long)D.ndigit[1] * RADIX + (long)D.ndigit[0] #define INIT(S) #define FREE(S) #define ASS(D,S) (D = S) #define ASS_WORD(D,S) (D = (t_v)S) #define ASS_COUNT(S) (D = (long)S.ndigit[1] * RADIX + (long)S.ndigit[0]) #define PRINT_NUM_S(F,C) fprintf(F,"%5ld ", C) #define PRINT_NUM_M(F,C) fprintf(F,"%5ld ", C) #define PRINT_NUM_MP(F,C) fprintf(F,"%5ld ", C) #define IS_ZERO(D) (D == 0L) #define IS_ONE(D) (D == 1L) #define ABS_LT(A,B) (abs(A) < abs(B)) #define POS_SIGN(A) A > 0L #define EXCHANGE(A,B) { long h = A; A = B; B = h; } #define ASS_ZERO(A) A = 0L #define ASS_ONE(A) A = 1L #define ASS_ADD(S,A) S += A #define ASS_SUB(S,A) S -= A #define NEG(A) A = -A #define PROD(R,A,B) R = A * B #define DIV(R,A,B) (R = A / B, A % B == 0L) #define GCD(R,A,B) R = gcd_long(A, B) long gcd_long(a, b) long a, b; { DEBUG_P(("gcd_long %ld %ld \n", a, b)); if (a == 0L || b == 0L) return 0L; if (a < 0L) a = -a; if (b < 0L) b = -b; DEBUG_P(("gcd %ld %ld\n", a, b)); while (a != b) if (a < b) { b = b % a; if (b == 0L) return a; } else { a = a % b; if (a == 0L) return b; } return a; } #endif #ifdef USE_BIG_INT #define num_t big_int #define GLOBAL(V) big_int V = { (byte)0, (word)0, (char)1, (word *)NULL } #define TO_V(D,I,S) assign(D + I, &S) #define INIT(S) init(&S) #define FREE(S) free_(&S) #define ASS(D,S) assign(&D, &S) #define ASS_WORD(D,S) assign_word(&D, S) #define ASS_COUNT(D,S) assign(&D, &S) #define PRINT_NUM_S(F,C) { print_big_int(F, &C, FALSE, TRUE);fprintf(F, " "); } #define PRINT_NUM_M(F,C) { print_big_int(F, &C, FALSE, FALSE);fprintf(F, " "); } #define PRINT_NUM_MP(F,C) { print_big_int(F, C, FALSE, FALSE);fprintf(F, " "); } #define IS_ZERO(D) is_word(&D, 0) #define IS_ONE(D) is_word(&D, 1) #define ABS_LT(A,B) comp_abs(&A, &B) < 0 #define POS_SIGN(A) A.sign == 1 #define EXCHANGE(A,B) exchange(&A, &B) #define ASS_ZERO(A) assign_word(&A, 0L) #define ASS_ONE(A) assign_word(&A, 1L) #define ASS_ADD(S,A) ass_add(&S, &A, FALSE) #define ASS_SUB(S,A) ass_add(&S, &A, TRUE) #define NEG(A) A.sign = -A.sign #define PROD(R,A,B) prod(&R, &A, &B) #define DIV(R,A,B) div_large(&R, NULL, &A, &B) #define GCD(R,A,B) gcd_large(&R, &A, &B) PROTO(bool is_word(num_t *, word);) bool is_word(a, w) num_t *a; word w; { int i; if (a->size == 0) return TRUE; if (NDIGIT(a, 0) != w) return FALSE; for (i = 1; i < a->len; i++) if (NDIGIT(a, i) != 0) return FALSE; return TRUE; } PROTO(void exchange(num_t *, num_t *);) void exchange(a, b) num_t *a, *b; { byte size; word len; char sign; word *ndigit; size = a->size; a->size = b->size; b->size = size; len = a->len; a->len = b->len; b->len = len; sign = a->sign; a->sign = b->sign; b->sign = sign; ndigit = a->ndigit; a->ndigit = b->ndigit; b->ndigit = ndigit; } PROTO(void sub_abs_times_to(num_t *, num_t *, word, word);) void sub_abs_times_to(s, a, sh, t) num_t *s, *a; word sh, t; { long i; word borrow = 0; long len_s, len_a; DEBUG_P(("sub_abs_times_to(")); DEBUG_NP(s); DEBUG_NP(a); DEBUG_P((",%d ,%d)\n", sh, t)); len_s = len(s); len_a = (long)a->len; DEBUG_P(("len_s %ld, len_a %ld\n", len_s, len_a)); for (i = sh; i < len_s; i++) { longword sum = (longword)borrow + ( (i - sh < len_a) ? (longword)NDIGIT(a, i - sh) * t : 0L); DEBUG_P((" s->ndigit[%ld] = %ld, sum = %ld ", i, (longword)NDIGIT(s, i), sum)); if ((longword)NDIGIT(s, i) < sum) { longword temp = sum - (longword)NDIGIT(s, i); DEBUG_P(("temp = %ld\n", temp)); ASS_RADIX(borrow, (temp - 1)/RADIX + 1); ASS_NDIGIT(s, i, (long)borrow * (long)RADIX - temp); DEBUG_P((" s->ndigit[%ld] = %ld\n", i, (longword)NDIGIT(s, i))); } else { ASS_NDIGIT(s, i, NDIGIT(s, i) - sum); DEBUG_P((" s->ndigit[%ld] = %ld\n", i, (longword)NDIGIT(s, i))); borrow = 0; } } ASSERT(borrow > 0, "ERROR: borrow error\n"); DEBUG_P(("result: ")); DEBUG_NP(s); DEBUG_P(("\n")); } PROTO(void ass_add(num_t *, num_t *, bool);) void ass_add(a, b, change_sign) num_t *a, *b; bool change_sign; /* affects both large numbers! */ { int sign = change_sign ? -b->sign : b->sign; DEBUG_P(("ass_add(); %c", a->sign == 1 ? '+' : '-')); DEBUG_NP(a); DEBUG_P((" %c %c", change_sign ? '-' : '+', b->sign == 1 ? '+' : '-')); DEBUG_NP(b); if (a->sign == sign) add_abs_times_to(a, b, 0, 1); else if (comp_abs(a, b) > 0) sub_abs_times_to(a, b, 0, 1); else { sign = -a->sign; exchange(a, b); sub_abs_times_to(a, b, 0, 1); a->sign = sign; } DEBUG_P((" = %c", a->sign == 1 ? '+' : '-')); DEBUG_NP(a); DEBUG_P(("\n")); } PROTO(void prod(num_t *, num_t *, num_t *);) void prod(r, a, b) num_t *r, *a, *b; { long i, len_a = len(a), len_b = len(b); DEBUG_P(("prod ")); DEBUG_NP(a); DEBUG_NP(b); RESIZE(r, SIZE_OF_LEN(len_a + len_b)); assign_word(r, 0); r->sign = (a->sign == b->sign) ? 1 : -1; if (len_a < len_b) { for (i = 0; i < len_a; i++) if (NDIGIT(a, i) != 0) add_abs_times_to(r, b, i, NDIGIT(a, i)); } else { for (i = 0; i < len_b; i++) if (NDIGIT(b, i) != 0) add_abs_times_to(r, a, i, NDIGIT(b, i)); } DEBUG_P(("prod ")); DEBUG_NP(a); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(r); } PROTO(int comp_abs_shift(num_t *, num_t *, int);) int comp_abs_shift(a, b, sh) num_t *a, *b; int sh; { long sa = a->len, sb = b->len + sh, s = sb, i; if (sa < sb) { for (i = sb - 1; i >= sa; i--) if ((i < sh ? 0 : NDIGIT(b, i - sh)) > 0) return -1; s = sa; } else if (sb < sa) { for (i = sa - 1; i >= sb; i--) if (NDIGIT(a, i) > 0) return 1; s = sb; } for (i = s - 1; i >= 0; i--) if (NDIGIT(a, i) == (i < sh ? 0 : NDIGIT(b, i - sh))) /* continue */; else if (NDIGIT(a, i) > (i < sh ? 0 : NDIGIT(b, i - sh))) return 1; else return -1; return 0; } GLOBAL(div_a); PROTO(bool div_large(num_t *, num_t *, num_t *, num_t *);) bool div_large(d, r, a, b) num_t *d, *r, *a, *b; { long len_a, len_b; long i; longword d1, d2; assign(&div_a, a); DEBUG_P(("div_large ")); DEBUG_NP(a); DEBUG_NP(b); len_a = len(a); len_b = len(b); DEBUG_P(("len_a %ld, len_b %ld\n", len_a, len_b)); if (d != NULL) { RESIZE(d, SIZE_OF_LEN(len_a - len_b + 2)); assign_word(d, 0); d->sign = a->sign * b->sign; } d1 = NDIGIT(b, len_b - 1) + (len_b == 1 ? 0L : 1L); d2 = NDIGIT(b, len_b - 1) * RADIX + (len_b == 1 ? 0L : (NDIGIT(b, len_b - 2) + (len_b == 2 ? 0L : 1L))); i = len_a - len_b; while (i >= 0) { DEBUG_P((" %ld: ", i)); if (len_b + i < a->len && NDIGIT(a, len_b + i) != 0) { longword n = NDIGIT((&div_a), len_b + i) * RADIX + NDIGIT((&div_a), len_b + i - 1); word q; ASS_RADIX(q, n / d1); DEBUG_P(("%4ld / %4ld = %4d ", n, d1, q)); if (d != NULL) if (NDIGIT(d, i) + q >= RADIX) { word carry = q, j = i; DEBUG_P(("-- d->ndigit[%ld] = %d, q = %d, %d\n", i, NDIGIT(d, i), q, NDIGIT(d, i) + q)); while (carry > 0 && j < d->len) { longword sum = NDIGIT(d, j) + carry; DEBUG_P((" - d->ndigit[%d] = %d, c = %d, %ld\n", j, NDIGIT(d, i), carry, sum)); ASS_NDIGIT(d, j, sum % RADIX); ASS_RADIX(carry, sum / RADIX); j++; } ASSERT1(carry > 0, "ERROR: carry in div1 %d\n", carry); } else ASS_NDIGIT(d, i, NDIGIT(d, i) + q); sub_abs_times_to(&div_a, b, i, q); } DEBUG_(else printf("skip ");) { longword n = div_a.ndigit[len_b + i - 1]*RADIX + (len_b > 1 ? div_a.ndigit[len_b + i - 2] : 0L); word q; ASS_RADIX(q, n / d2); DEBUG_P(("%4ld / %4ld = %4d ", n, d2, q)); if (d != NULL) if (NDIGIT(d, i) + q >= RADIX) { word carry = q, j = i; DEBUG_P(("-- d->ndigit[%ld] = %d, q = %d, %d\n", i, NDIGIT(d, i), q, NDIGIT(d, i) + q)); while (carry > 0 && j < d->len) { longword sum = NDIGIT(d, j) + carry; DEBUG_P((" - d->ndigit[%d] = %d, c = %d, %ld\n", j, NDIGIT(d, i), carry, sum)); ASS_NDIGIT(d, j, sum % RADIX); ASS_RADIX(carry, sum / RADIX); j++; } ASSERT1(carry > 0, "ERROR: carry in div2 %d\n", carry); } else ASS_NDIGIT(d, i, NDIGIT(d, i) + q); sub_abs_times_to(&div_a, b, i, q); } while(comp_abs_shift(&div_a, b, i) >= 0) { DEBUG_P(("%4d ", 1)); if (d != NULL) if (NDIGIT(d, i) + 1 >= RADIX) { word carry = 1, j = i; DEBUG_P(("-- d->ndigit[%ld] = %d, q = 1, %d\n", i, NDIGIT(d, i), NDIGIT(d, i) + 1)); while (carry > 0 && j < d->len) { longword sum = NDIGIT(d, j) + carry; DEBUG_P((" - d->ndigit[%d] = %d, c = %d, %ld\n", j, NDIGIT(d, i), carry, sum)); ASS_NDIGIT(d, j, sum % RADIX); ASS_RADIX(carry, sum / RADIX); j++; } ASSERT1(carry > 0, "ERROR: carry in div3 %d\n", carry); } else ASS_NDIGIT(d, i, NDIGIT(d, i) + 1); sub_abs_times_to(&div_a, b, i, 1); } i--; DEBUG_P(("\n")); } if (r != NULL) { assign(r, &div_a); DEBUG_P(("div: ")); DEBUG_NP(a); DEBUG_P((" %% ")); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(r); DEBUG_P(("\n")); } else if (!IS_ZERO(div_a)) { printf("Possible div error"); print_big_int(stdout, a, FALSE, FALSE); printf(" %% "); print_big_int(stdout, b, FALSE, FALSE); printf(" = "); print_big_int(stdout, &div_a, FALSE, FALSE); printf("\n"); } if (d != NULL) { DEBUG_P(("div: ")); DEBUG_NP(a); DEBUG_P((" / ")); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(d); DEBUG_P(("\n")); } return IS_ZERO(div_a); } GLOBAL(gcd_a); GLOBAL(gcd_b); GLOBAL(gcd_c); PROTO(void gcd_large(num_t *, num_t *, num_t *);) void gcd_large(r, a, b) num_t *r, *a, *b; { int c; DEBUG_P(("gcd_large ")); DEBUG_NP(a); DEBUG_NP(b); DEBUG_P(("\n")); ASSERT(is_word(a, 0) || is_word(b, 0), "ERROR: zero argument gcd\n"); assign(&gcd_a, a); assign(&gcd_b, b); while ((c = comp_abs(&gcd_a, &gcd_b)) != 0) { DEBUG_P(("gcd_large - ")); DEBUG_N(gcd_a); DEBUG_N(gcd_b); DEBUG_P(("\n")); if (c < 0) { if (div_large(NULL, &gcd_c, &gcd_b, &gcd_a)) { assign(r, &gcd_a); DEBUG_P(("gcd_large ")); DEBUG_NP(a); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(r); DEBUG_P(("\n")); return; } assign(&gcd_b, &gcd_c); } else { if (div_large(NULL, &gcd_c, &gcd_a, &gcd_b)) { assign(r, &gcd_b); DEBUG_P(("gcd_large ")); DEBUG_NP(a); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(r); DEBUG_P(("\n")); return; } assign(&gcd_a, &gcd_c); } } assign(r, &gcd_a); DEBUG_P(("gcd_large ")); DEBUG_NP(a); DEBUG_NP(b); DEBUG_P((" = ")); DEBUG_NP(r); DEBUG_P(("\n")); } #endif PROTO(short print_counts(FILE *, num_t *, word, bool);) short print_counts(f, array, nr_steps, as_equation) FILE *f; num_t *array; word nr_steps; bool as_equation; { short i; DEBUG_P(("print_counts()\n")); if (f != NULL) fprintf(f, "\n"); if (!as_equation && f != NULL) fprintf(f, " n H(n)\n"); init_counting(0); for (i = 1; i < (short)nr_steps; i++) { next_step_counting(); if (f != NULL) { if (as_equation) { fprintf(f, "$C(%d) = ", i); print_big_int(f, &m_vertex_list->new_c, FALSE, TRUE); fprintf(f, "$%s", i == nr_steps - 1 ? " and" : ","); } else { fprintf(f, "%3d ", i); print_big_int(f, &m_vertex_list->new_c, TRUE, TRUE); } fprintf(f, "\n"); } if (array != NULL) TO_V(array, i, m_vertex_list->new_c); freq_analyzed = !analyze_freq(i); pre_next_step(); } if (f != NULL) fprintf(f, "\n"); DEBUG_P(("print_counts() okay %d\n", i)); return i; } PROTO(void print_sequence(FILE *, char *graph_name, int problem_kind);) void print_sequence(f, graph_name, problem_kind) FILE *f; char *graph_name; int problem_kind; { short i; short pos = 10, lines = 0; DEBUG_P(("print_sequence()\n")); fprintf(f, "%%I A0000\n"); init_counting(0); fprintf(f, "%%S A0000 "); for (i = 1; lines < 3; i++) { next_step_counting(); { short l = len(&m_vertex_list->new_c); short v = l > 0 ? m_vertex_list->new_c.ndigit[l-1] : 1; short nr_digits = (l > 0 ? l : 1) * 4 - 3 + (v >= 10 ? 1 : 0) + (v >= 100 ? 1 : 0) + (v >= 1000 ? 1 : 0); fprintf(stderr, "nr_digits %d\n", nr_digits); if (nr_digits + pos > 75) { lines++; pos = 10; if (lines < 3) { if (i > 1) fprintf(f, ","); fprintf(f, "\n%%%c A0000 ", lines == 1 ? 'T' : 'U'); } } else if (i > 1) fprintf(f, ","); if (lines < 3) print_big_int(f, &m_vertex_list->new_c, FALSE, TRUE); pos += nr_digits+1; } freq_analyzed = !analyze_freq(i); pre_next_step(); } fprintf(f, "\n"); fprintf(f, "%%N A0000 Number of %s in $%s x P_n$.\n", kinds[problem_kind].name, graph_name); fprintf(f, "%%R A0000 \n"); fprintf(f, "%%K A0000 nonn\n"); fprintf(f, "%%O A0000 1,2\n"); fprintf(f, "%%A A0000 http://www.cs.utwente.nl/~faase/C/%s.%s.txt\n\n", graph_name, kinds[problem_kind].ext); DEBUG_P(("print_sequence() okay %d\n", i)); } /*** ANALYZE ALL FREQUENCIES ****/ PROTO(void analyze_all_freq(void);) void analyze_all_freq() { short i; DEBUG_P(("analyze_all_freq()\n")); if (freq_analyzed) return; /*! 950906 added: */ init_counting(0); /*! 940129_2 added: i < 100 */ for (i = 1; !freq_analyzed && i < 100; i++) { next_step_counting(); freq_analyzed = !analyze_freq(i); pre_next_step(); } DEBUG_P(("analyze_all_freq() okay\n")); } /*** BSORT ***/ static void swap(m1, m2, size) char *m1, *m2; size_t size; { while(size-- > 0) { char h; h = m1[size]; m1[size] = m2[size]; m2[size] = h; } } void bsort(arr, nr, size, comp) void *arr; size_t nr, size; #ifdef PROTO int (*comp)(const void *, const void *); #else int (*comp)(); #endif { int last, first; #ifdef QUICK int d = nr / 2 + 1; #endif DEBUG_P(("bsort()\n")); #ifdef QUICK while (d > 1) { int nr_swapped = 0, i; char *arr1 = (char *)arr, *arr2 = (char *)arr + (long)d * size; last = nr - d; for (i = 0; i < last; i++, arr1 += size, arr2 += size) if (comp(arr1, arr2) < 0) { swap(arr1, arr2, size); nr_swapped++; } if (nr_swapped == 0) d = d / 3; else { double nd = (double)nr_swapped/(double)last; nd = (1 + sqrt(nd))/2; d = (int)( d * nd ); } if (d < 1) d = 1; } #endif /* Normale left/right bubble-sort algorithme: */ first = 0; last = nr - 1; while (first < last) { int n_first, n_last, i; char *arr1 = (char *)arr + (long)size * first, *arr2 = arr1 + size; NO_DEBUG_P(("range %3d %3d forward\n", first, last)); n_last = first; for (i = first; i < last; i++, arr1 = arr2, arr2 += size) { NO_DEBUG_P((" comp: %d %d, %d %d\n", i, i+1, *((int *)arr1), *((int *)arr2))); if (comp(arr1, arr2) > 0) { swap(arr1, arr2, size); NO_DEBUG_P((" swap: %d %d\n", i, i+1)); n_last = i; } } last = n_last; if (first < last) { arr2 = (char *)arr + (long)size * last; arr1 = arr2 - size; NO_DEBUG_P(("range %3d %3d backward\n", first, last)); n_first = last; for (i = last; i > first; i--, arr2 = arr1, arr1 -= size) { NO_DEBUG_P((" comp: %d %d, %d %d\n", i-1, i, ((int *)arr)[i-1], ((int *)arr)[i])); if (comp(arr1, arr2) > 0) { swap(arr1, arr2, size); NO_DEBUG_P((" swap: %d %d\n", i-1, i)); n_first = i; } } first = n_first; } } DEBUG_P(("bsort()\n")); } /*** ORDER STATES ***/ static int by_order(e1, e2) void *e1, *e2; { m_vertex_list_p sl1 = *(m_vertex_list_p *)e1, sl2 = *(m_vertex_list_p *)e2; NO_DEBUG_P(("by order %ld %ld ", *(long int *)e1, *(long int *)e2)); NO_DEBUG_(print_pcc(stdout, sl1->pcc)); NO_DEBUG_P((" ")); NO_DEBUG_(print_pcc(stdout, sl2->pcc)); NO_DEBUG_P(("\n")); if (sl1 == m_vertex_list) return -1; if (sl2 == m_vertex_list) return 1; if (sl1->nr_ends < sl2->nr_ends) return -1; if (sl1->nr_ends > sl2->nr_ends) return 1; if (sl1->freq < sl2->freq) return -1; if (sl1->freq > sl2->freq) return 1; if (sl1->first < sl2->first) return -1; if (sl1->first > sl2->first) return 1; if (*(long int *)sl1 < *(long int *)sl2) return -1; if (*(long int *)sl1 > *(long int *)sl2) return 1; NO_DEBUG_P(("0!!!!\n")); return 1; } m_vertex_list_p *orderd_m_vertices = NULL; word nr_m_vertices = 0; PROTO(void order_m_vertices(void);) void order_m_vertices() { nr_m_vertices = 0; NO_DEBUG_P(("start order_m_vertices (count)\n")); { m_vertex_list_p sl = m_vertex_list; while (sl != NULL) { nr_m_vertices++; sl = sl->next; } } if (orderd_m_vertices) free(orderd_m_vertices); NO_DEBUG_P(("start order_m_vertices (alloc)\n")); orderd_m_vertices = ARRAY_MALLOC(m_vertex_list_p, nr_m_vertices); NO_DEBUG_P(("start order_m_vertices (transfer)\n")); { m_vertex_list_p sl = m_vertex_list; word i = 0; while (sl != NULL && i < nr_m_vertices) { orderd_m_vertices[i] = sl; i++; sl = sl->next; } } /*** NO_DEBUG_P(("start order_m_vertices (bsort) %d\n", nr_m_vertices)); if (freq_analyzed) bsort(orderd_m_vertices, nr_m_vertices, sizeof(m_vertex_list_p), by_order); NO_DEBUG_P(("end order_m_vertices \n")); ***/ } /*** PRINT ORDERD MATRIX ***/ bool logically_zero_arc(from_m_vertex_p, to_m_vertex_p) m_vertex_list_p from_m_vertex_p, to_m_vertex_p; { if ( from_m_vertex_p == m_vertex_list || to_m_vertex_p == m_vertex_list) return FALSE; if ( problem_kind == PROP_HAMIL_PATH && from_m_vertex_p->nr_ends > to_m_vertex_p->nr_ends) return TRUE; if (to_m_vertex_p->freq != (byte)gcd(to_m_vertex_p->freq, from_m_vertex_p->freq)) return TRUE; if ( to_m_vertex_p->freq > 1 && to_m_vertex_p->freq == from_m_vertex_p->freq && (to_m_vertex_p->first % to_m_vertex_p->freq) == (from_m_vertex_p->first % to_m_vertex_p->freq)) return TRUE; return FALSE; } void print_orderd_matrix(f) FILE *f; { word i; NO_DEBUG_P(("start print_orderd_matrix\n")); fprintf(f, "\n\nOrderd matrix:\n"); for (i = 0; i < nr_m_vertices; i++) { m_vertex_list_p m_vertex_list_to = orderd_m_vertices[i]; print_pcc(f, m_vertex_list_to->pcc); fprintf(f, " %1d %1d %1d : ", m_vertex_list_to->nr_ends, m_vertex_list_to->freq, m_vertex_list_to->first); { word j; for (j = 0; j < nr_m_vertices; j++) { longword nr_br = count_of_arc(orderd_m_vertices[j], m_vertex_list_to); if (logically_zero_arc(orderd_m_vertices[j], m_vertex_list_to)) if (nr_br > 0) fprintf(f, "!%2ld", nr_br); else fprintf(f, " "); else fprintf(f, " %2ld", nr_br); } } fprintf(f, "\n"); } NO_DEBUG_P(("end print_orderd_matrix\n")); } longword *orderd_matrix; #define POS(i,j) (i) + (j)*nr_m_vertices #define ORDERD_MATRIX(i,j) orderd_matrix[POS(i,j)] PROTO(static bool create_orderd_matrix(void);) static bool create_orderd_matrix() { word i; NO_DEBUG_P(("start of create_orderd_matrix\n")); orderd_matrix = ARRAY_MALLOC(longword, nr_m_vertices * nr_m_vertices); if (orderd_matrix == NULL) return FALSE; for (i = 0; i < nr_m_vertices; i++) { m_vertex_list_p m_vertex_list_to = orderd_m_vertices[i]; { word j; for (j = 0; j < nr_m_vertices; j++) ORDERD_MATRIX(i,j) = count_of_arc(orderd_m_vertices[j], m_vertex_list_to); } } NO_DEBUG_P(("end of create_orderd_matrix\n")); return TRUE; } word *perm_x, *perm_y; #define DET_VALUE(i,j) ORDERD_MATRIX(perm_x[i], perm_y[j]) #define DET_LAMBDA(i,j) (perm_x[i] == perm_y[j] ? -1 : 0) #define DET_FILLED(i,j) (perm_x[i] == perm_y[j] || DET_VALUE(i,j) != 0) #define SPACE(d) { word sp; for (sp = nr_m_vertices; sp > (d); sp--) \ printf(" "); } typedef long int result_t; #ifdef TGT double calc_det_done; time_t begin_time; #endif #ifdef TGT void calc_det_for(depth, result, step) word depth; result_t *result; double step; #else void calc_det_for(depth, result) word depth; result_t *result; #endif { word depth_1; { word i; for (i = 0; i <= depth; i++) result[i] = (result_t)0; } #ifdef DEBUG SPACE(depth); printf("calc_det_for(%d, %ld)\n", depth, (long int)result); { word x, y; SPACE(depth); printf(" "); for (x = 0; x < depth; x++) printf("%2d ", perm_x[x]); printf("\n"); for (y = 0; y < depth; y++) { SPACE(depth); printf("%d ", perm_y[y]); for (x = 0; x < depth; x++) printf("%2d%s ", DET_VALUE(x, y), DET_LAMBDA(x, y) ? "-x" : " "); printf("\n"); } } #endif if (depth == 1) { result[0] = (result_t)DET_VALUE(0, 0); result[1] = (result_t)DET_LAMBDA(0, 0); #ifdef TGT calc_det_done += step; if (++calc_det_done_s > 1000) printf("calc_det: %lf\n", calc_det_done); #endif } else { word min_x, min_x_v, min_y, min_y_v, x, y; min_x_v = depth + 1; NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("x: ")); for (x = 0; x < depth && min_x_v > 1; x++) { word c = 0; for (y = 0; y < depth; y++) if (DET_FILLED(x, y)) { c++; if (c == min_x_v) break; } if (c < min_x_v) { min_x_v = c; min_x = x; } NO_DEBUG_P(("%d ", c)); } NO_DEBUG_P(("\n")); NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("y: ")); min_y_v = depth + 1; for (y = 0; y < depth && min_y_v > 1; y++) { word c = 0; for (x = 0; x < depth; x++) if (DET_FILLED(x, y)) { c++; if (c == min_y_v) break; } if (c < min_y_v) { min_y_v = c; min_y = y; } NO_DEBUG_P(("%d ", c)); } NO_DEBUG_P(("\n")); depth_1 = depth - 1; if (min_x_v <= min_y_v) { NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("scan x = %d\n", min_x)); { word h,i; h = perm_x[min_x]; for (i = min_x; i < depth_1; i++) perm_x[i] = perm_x[i + 1]; perm_x[depth_1] = h; } for (y = 0; y < depth; y++) if (DET_FILLED(depth_1, y)) { short sign = (min_x + y) % 2 == 1 ? -1 : 1; result_t *new_result = result + depth + 1; NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("sub-scan x,y = %d, %d %d\n", depth_1, y, DET_VALUE(depth_1,y))); { word h,i; h = perm_y[y]; for (i = y; i < depth_1; i++) perm_y[i] = perm_y[i + 1]; perm_y[depth_1] = h; } calc_det_for(depth_1, new_result); { word h,i; h = perm_y[depth_1]; for (i = depth_1; i > y; i--) perm_y[i] = perm_y[i - 1]; perm_y[y] = h; } #ifdef DEBUG { word i; SPACE(depth); printf("prev-result (%d) : ", depth); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } SPACE(depth); printf("sign=%2d 1: %4d x: %4d\n", sign, DET_VALUE(depth_1, y), DET_LAMBDA(depth_1, y)); #endif if (DET_VALUE(depth_1, y) != 0) { word i; short v = sign * DET_VALUE(depth_1, y); for (i = 0; i < depth; i++) result[i] += v * new_result[i]; } #ifdef DEBUG { word i; SPACE(depth); printf("result after det : "); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } #endif if (DET_LAMBDA(depth_1, y) != 0) { word i; for (i = 0; i < depth; i++) result[i + 1] -= sign * new_result[i]; } #ifdef DEBUG { word i; SPACE(depth); printf("curr-result (%d) : ", depth); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } #endif } { word h,i; h = perm_x[depth_1]; for (i = depth_1; i > min_x; i--) perm_x[i] = perm_x[i - 1]; perm_x[min_x] = h; } } else { NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("scan y = %d\n", min_y)); { word h,i; h = perm_y[min_y]; for (i = min_y; i < depth_1; i++) perm_y[i] = perm_y[i + 1]; perm_y[depth_1] = h; } for (x = 0; x < depth; x++) if (DET_FILLED(x, depth_1)) { short sign = (x + min_y) % 2 == 1 ? -1 : 1; result_t *new_result = result + depth + 1; NO_DEBUG_(SPACE(depth)); NO_DEBUG_P(("sub-scan x,y = %d, %d %d\n", x, depth_1, DET_VALUE(x, depth_1))); { word h,i; h = perm_x[x]; for (i = x; i < depth_1; i++) perm_x[i] = perm_x[i + 1]; perm_x[depth_1] = h; } calc_det_for(depth_1, new_result); { word h,i; h = perm_x[depth_1]; for (i = depth_1; i > x; i--) perm_x[i] = perm_x[i - 1]; perm_x[x] = h; } #ifdef DEBUG { word i; SPACE(depth); printf("prev-result (%d) : ", depth); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } SPACE(depth); printf("sign=%2d 1: %4d x: %4d\n", sign, DET_VALUE(x, depth_1), DET_LAMBDA(x, depth_1)); #endif if (DET_VALUE(x, depth_1) != 0) { word i; short v = sign * DET_VALUE(x, depth_1); for (i = 0; i < depth; i++) result[i] += v * new_result[i]; } #ifdef DEBUG { word i; SPACE(depth); printf("result after det : "); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } #endif if (DET_LAMBDA(x, depth_1) != 0) { word i; for (i = 0; i < depth; i++) result[i + 1] -= sign * new_result[i]; } #ifdef DEBUG { word i; SPACE(depth); printf("curr-result (%d) : ", depth); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } #endif } { word h,i; h = perm_y[depth_1]; for (i = depth_1; i > min_y; i--) perm_y[i] = perm_y[i - 1]; perm_y[min_y] = h; } } } #ifdef DEBUG { word i; SPACE(depth); printf("result (%d) : ", depth); for (i = 0; i <= depth; i++) printf("%ld ", result[i]); printf("\n"); } #endif } void calc_det(f, inc_zero) FILE *f; bool inc_zero; { word nr_n; result_t *result; word depth = inc_zero ? nr_m_vertices : (word)(nr_m_vertices - 1); NO_DEBUG_P(("depth: %d %d \n", depth, depth * (depth + 1) / 2 + depth)); perm_x = ARRAY_MALLOC(word, depth); perm_y = ARRAY_MALLOC(word, depth); result = ARRAY_MALLOC(result_t, depth * (depth + 1) / 2 + depth); { word i; for (i = 0; i < depth; i++) { perm_x[i] = inc_zero ? i : i + 1; perm_y[i] = perm_x[i]; } } if (depth > 0) calc_det_for(depth, result); fprintf(f,"\nSolution\n"); { word i; nr_n = 0; for (i = 1; i <= depth; i++) if (result[depth - i]!= (result_t)0) nr_n = i + 1; } print_counts(f, NULL, nr_n + 1, TRUE); fprintf(f, "$C(n) = "); { word i, n; result_t sign; n = 0; sign = - result[depth]; for (i = 1; i <= depth; i++) { result_t v = result[depth - i] * sign; if (v != (result_t)0) { if (v < (result_t)0) { fprintf(f, " - "); v = -v; } else if (n > 0) fprintf(f, " + "); if (v != 1) fprintf(f ,"%ld", v); fprintf(f, "C(n-%d)", i); n++; if (n % 5 == 0) fprintf(f, "\n "); } } if (n == 0) fprintf(f, "0"); } fprintf(f, "$."); } #define M(i,j) m[(long)(i) + (long)(j) * (long)nr_s] PROTO(void fill_matrix(num_t *, word, int, int, int);) void fill_matrix(m, nr_s, start_with, nr_n, start_off) num_t *m; word nr_s; int start_with, nr_n, start_off; { short i; DEBUG_P(("fill_matrix(,%d,%d)\n", nr_s, start_with)); init_counting(start_with); for (i = 0; i < (short)nr_n + start_off; i++) { next_step_counting(); if (i >= start_off) { m_vertex_list_p sl; short i_o = i - start_off, j; for (sl = m_vertex_list, j = 0; sl != NULL; sl = sl->next, j++) { DEBUG_P(("m[%d,%d] = ", i, j)); DEBUG_N(sl->new_c); ASS_COUNT(M(i_o, j), sl->new_c); } DEBUG_P(("\n")); } pre_next_step(); } } #define NO_DIVIDE GLOBAL(rec_d); GLOBAL(rec_s); GLOBAL(rec_g); GLOBAL(rec_m1); GLOBAL(rec_m2); #define DEBUG_P_M { int i1, i2; \ for (i1 = 0; i1 < nr_l; i1++) \ { for (i2 = 0; i2 < nr_n; i2++) \ DEBUG_N(M(i2, i1)); \ DEBUG_P(("\n")); } } #define _DEBUG_P_M { int i1, i2; \ for (i1 = 0; i1 < nr_l; i1++) \ { for (i2 = 0; i2 < nr_n; i2++) \ _DEBUG_N(M(i2, i1)); \ _DEBUG_P(("\n")); } } void solve(num_t *m, int nr_n, int nr_l, int nr_s, int *nr_f, int nr_fl, int *np_t); void solve(num_t *m, int nr_n, int nr_l, int nr_s, int *nr_f, int nr_fl, int *np_t) /* m - filled matrix nr_n - number filled colums nr_l - number filled rows nr_s - number colums (used in M(i,j)) */ { int i; for (i = 0; i < nr_n && *nr_f > nr_fl; i++) { int j, k; int min_j = i; DEBUG_P(("Start with pass %d\n", i)); DEBUG_P_M /* find row >= i, with minimal non-zero value at i-th column */ for (j = i + 1; j < nr_l; j++) if ( IS_ZERO(M(i, min_j)) || ( !IS_ZERO(M(i,j)) && ABS_LT(M(i,j), M(i, min_j)))) min_j = j; /* If all the values in i-th column are zero, exit loop */ /* (Is this a correct assumption?) */ if (IS_ZERO(M(i, min_j))) break; DEBUG_P((" min %d with ", min_j)); DEBUG_N(M(i, min_j)); DEBUG_P(("\n")); /* Exchange rows, if needed */ if (i != min_j) for (k = 0; k < nr_n; k++) EXCHANGE(M(k, i), M(k, min_j)); /* For the very first row: */ if (i == 0) { bool all_zero = TRUE; /* determine gcd for the values in this row: */ for (j = 0; j < nr_n && (all_zero || !IS_ONE(rec_g)); j++) if (!IS_ZERO(M(j, 0))) if (all_zero) { ASS(rec_g, M(j, 0)); all_zero = FALSE; } else if(!IS_ONE(rec_g)) GCD(rec_g, rec_g, M(j, 0)); /* divide all numbers by gcd: */ /* (all_zero not possible!?) */ if (!all_zero && !IS_ONE(rec_g)) { for (j = 0; j < nr_n; j++) { ASSERT(!DIV(rec_d, M(j, 0), rec_g), "ERROR: div\n"); ASS(M(j, 0), rec_d); } } if (all_zero) (*nr_f)--; } DEBUG_P_M DEBUG_P(("sweeping\n")); /* Sweep all rows > i with i-th row: */ for (j = i + 1; j < nr_l; j++) if (!IS_ZERO(M(i, j))) { bool all_zero = TRUE; DEBUG_P(("sweep %d ", j)); /* determine gcd and calculate multiplier m1 and m2: */ GCD(rec_g, M(i,i), M(i, j)); ASSERT( !DIV(rec_m1, M(i, i), rec_g) || !DIV(rec_m2, M(i, j), rec_g), "ERROR: div2\n"); DEBUG_N(rec_g); DEBUG_N(rec_m1); DEBUG_N(rec_m2); DEBUG_P(("\n")); /* for all colums of the row: */ for (k = i + 1; k < nr_n; k++) #ifdef NO_DIVIDE { /* M(k, j) = M(k, j)*M(i, i) - M(k, i)*M(i, j); */ PROD(rec_s, M(k, j), rec_m1); PROD(rec_d, M(k, i), rec_m2); ASS_SUB(rec_s, rec_d); ASS(M(k, j), rec_s); #define GCD_ON_ROW #ifdef GCD_ON_ROW /* calculate gcd of the values in the row: */ if (!IS_ZERO(rec_s)) if (all_zero) { ASS(rec_g, rec_s); all_zero = FALSE; } else if(!IS_ONE(rec_g)) GCD(rec_g, rec_g, rec_s); #endif } #else M(k, j) -= M(k, i) * M(i, j) / M(i, i); #endif ASS_ZERO(M(i, j)); DEBUG_P_M #ifdef GCD_ON_ROW /* divide all values by gcd */ if (!all_zero && !IS_ONE(rec_g)) { DEBUG_P(("divide by: ")); DEBUG_N(rec_g); DEBUG_P(("\n")); for (k = i + 1; k < nr_n; k++) { ASSERT(!DIV(rec_d, M(k, j), rec_g), "ERROR: div3\n"); DEBUG_P(("div ")); DEBUG_N(M(k, j)); DEBUG_N(rec_g); DEBUG_P((" = ")); DEBUG_N(rec_d); DEBUG_P(("\n")); ASS(M(k, j), rec_d); } DEBUG_P_M } #endif if (all_zero) (*nr_f)--; } _DEBUG_P(("After sweeping\n")); DEBUG_P_M } _DEBUG_P_M *np_t = i; } GLOBAL(rec_p); bool calc_rec(FILE *, int, int, int, int, bool, bool); bool calc_rec(f, begin_at, step, end_at, start_off, skip, alt) FILE *f; int begin_at, step, end_at, start_off; bool skip, alt; { num_t *a, *m, *p; int i, j, np, np_new, np_t; int nr_s = nr_m_vertices + 1, nr_mul = 1, shift = 0, nr_a = nr_s * 2 + 5 + start_off, nr_n, nr_s_ = alt ? 2 * nr_s : nr_s, nr_c; bool found_rec = FALSE, no_div_error, only_2_folds, only_odd, only_even, previous; a = ARRAY_MALLOC(num_t, nr_a + 1); for (i = 0; i <= nr_a; i++) { INIT(a[i]); ASS_ZERO(a[i]); } if (nr_a != print_counts(stdout, a, nr_a, FALSE)) { fprintf(f, "Calculation not possible\n"); return FALSE; } { word i; only_odd = TRUE; only_even = TRUE; for (i = 1; i < nr_a; i++) if (i % 2 == 0) { if (only_odd && !IS_ZERO(a[i])) only_odd = FALSE; } else { if (only_even && !IS_ZERO(a[i])) only_even = FALSE; } } _DEBUG_P(("only_odd %d, only_even %d\n", only_odd, only_even)); if (only_odd && only_even) { fprintf(f, "$C(n) = 0$.\n"); return TRUE; } if (!alt) { if (only_odd || only_even || skip) { nr_s = (nr_s + 1) / 2; nr_mul = 2; } if (only_odd) shift = -1; } m = ARRAY_MALLOC(num_t, (long)nr_s_ * (long)nr_s); for (i = 0; i < (long)nr_s_ * (long)nr_s; i++) { INIT(m[i]); ASS_ZERO(m[i]); } _DEBUG_P(("nr_s : %d\n", nr_s)); p = ARRAY_MALLOC(num_t, nr_s + 1); for (i = 0; i < nr_s + 1; i++) { INIT(p[i]); ASS_ZERO(p[i]); } if (begin_at == -1) begin_at = nr_s; if (end_at == -1) end_at = nr_s; for (nr_n = begin_at; nr_n <= end_at && !found_rec; nr_n += step) { int nr_l = 0; DEBUG_P(("start\n")); for (i = 0; i < nr_s; i++) DEBUG_N(a[i]); if (alt) { for (i = 0; i < nr_s; i++) { int j1, j2; fill_matrix(m + nr_l * nr_s, nr_s, i, nr_n, start_off); nr_l += nr_s; solve(m, nr_n, nr_l, nr_s, &nr_l, nr_s, &np_t); /* move all non-empty rows to the top: */ for (j1 = 0, j2 = 0; j1 < nr_l; j1++) { int l; bool all_zero = TRUE; for (l = 0; l < nr_n && all_zero; l++) all_zero = IS_ZERO(M(j1, l)); if (!all_zero) { if (j1 != j2) for (l = 0; l < nr_n; l++) EXCHANGE(M(j1,l), M(j2,l)); j2++; } } nr_l = j2; } } else { for (i = 0; i < nr_n; i++) for (j = 0; j < nr_n; j++) ASS(M(i, j), a[(i + j + 1 + start_off) * nr_mul + shift]); nr_l = nr_n; } solve(m, nr_n, nr_l, nr_s, &nr_l, 0, &np_t); { int i1, i2; for (i1 = 0; i1 < np_t; i1++) { for (i2 = 0; i2 < nr_n; i2++) PRINT_NUM_S(f, M(i2, i1)); fprintf(f, "\n"); } } np_new = np_t; _DEBUG_P(("np_t %d\n", np_t)); for (np = 1; np <= np_t && !found_rec; np++) { bool all_other_zero = TRUE; no_div_error = TRUE; ASS_ONE(p[0]); for (i = 1; i <= np && no_div_error; i++) { ASS_ZERO(rec_s); _DEBUG_P(("%d: ", i)); for (j = 0; j < i; j++) { DEBUG_P(("s += p[%d] * M(%d, %d) = ", j, np - j, np - i)); _DEBUG_N(p[j]); _DEBUG_P((" * ")); _DEBUG_N(M(np - j, np - i)); _DEBUG_P(("\n")); PROD(rec_p, p[j], M(np - j, np - i)); ASS_ADD(rec_s, rec_p); _DEBUG_P(("s = ")); _DEBUG_N(rec_s); _DEBUG_P(("\n")); } NEG(rec_s); no_div_error = DIV(p[i], rec_s, M(np - i, np - i)); /* if (!no_div_error) fprintf(f, "%%div error %d %d %d\n", np_t, np, np - i);*/ if (!IS_ZERO(p[i])) { np_new = i; all_other_zero = FALSE; } _DEBUG_P(("p[%d] = ", i)); _DEBUG_N(p[i]); _DEBUG_P(("\n")); } fflush(f); if (np_new < nr_n && no_div_error && !all_other_zero) { word i; found_rec = FALSE; /* for (i = np_new * nr_mul + 2 + start_off; i < nr_a && found_rec; i++) */ for (i = nr_a - 1; i > np_new * nr_mul; i--) { word j; ASS_ZERO(rec_s); for (j = 0; j <= np_new; j++) { PROD(rec_p, p[j], a[i - j * nr_mul]); ASS_ADD(rec_s, rec_p); _DEBUG_N(p[j]); _DEBUG_P((" * ")); _DEBUG_N(a[i - j * nr_mul]); _DEBUG_P((" = ")); _DEBUG_N(rec_p); _DEBUG_P((", ")); _DEBUG_N(rec_s); _DEBUG_P(("\n")); } if (IS_ZERO(rec_s)) { found_rec = TRUE; nr_c = i; } else break; } if (found_rec) fprintf(f, "%%checked for %d to %d\n", nr_c, nr_a); else fprintf(f, "%%rec for %d %d not found\n", np_t, np); } else fprintf(f, "%%failed for %d %d\n", np_t, np); } np = np_new; } DEBUG_P(("ready\n")); /* printing result: */ if (!found_rec) return FALSE; previous = FALSE; for (i = 1; i < nr_c; i++) { if (previous) fprintf(f, ",\n"); previous = TRUE; fprintf(f, "$C(%d) = ", i); PRINT_NUM_S(f, a[i]); fprintf(f, "$"); } if (previous) fprintf(f, " and\n"); fprintf(f, "$C(n) = "); { word i, n; n = 0; for (i = 1; i <= np; i++) if (!IS_ZERO(p[i])) { if (POS_SIGN(p[i])) fprintf(f, " - "); else { if (n > 0) fprintf(f, " + "); NEG(p[i]); } if (!IS_ONE(p[i])) PRINT_NUM_S(f, p[i]); fprintf(f, "C(n-%d)", i * nr_mul); n++; if (n % 5 == 0) fprintf(f, "\n "); } if (n == 0) fprintf(f, "0"); } fprintf(f, "$."); return TRUE; } /*** PROGRAM OPTIONS ***/ #define OPT_P_INP 0 #define OPT_P_PERM 1 #define OPT_D_GEN 2 #define OPT_P_VERTICES 3 #define OPT_P_COUNTS 4 #define OPT_P_SEQ 5 #define OPT_P_MATRIX 6 #define OPT_P_SOL 7 #define OPT_A_OUT 8 #define OPT_T_TGT 9 #define MAX_OPT 9 struct { char letter; char *description; word value; } options[] = { { 'I', "print the input graph", OPT_P_INP } , { 'P', "print the permutations", OPT_P_PERM } , { 'G', "display M-vertex generation", OPT_D_GEN } , { 'V', "print M-vertex list", OPT_P_VERTICES } , { 'C', "print values of H(n)", OPT_P_COUNTS } , { 'S', "print sequence", OPT_P_SEQ } , { 'M', "print matrix of M", OPT_P_MATRIX } , { 'R', "print recurrent equation", OPT_P_SOL } , { 'A', "Append to output file", OPT_A_OUT } , { 'T', "Show time to go", OPT_T_TGT } }; #ifndef MSC static short stricmp(str1, str2) char *str1, *str2; { while (*str1 != '\0' && toupper(*str1) == toupper(*str2)) { str1++; str2++; } return (toupper(*str1) < toupper(*str2)) ? -1 : (toupper(*str1) > toupper(*str2)) ? 1 : 0; } #endif int main(argc, argv) int argc; char *argv[]; { char *in_filename = NULL, *out_filename = NULL, *graph_name = NULL; FILE *f; word sel_options = 0; int begin_at = -1, step = 1, end_at = -1, start_off = 0, only_even = 0; bool alt = FALSE, params = FALSE, do_calc_deter = FALSE; #define OPTION_SEL(o) (sel_options & (1 << o)) { bool error = FALSE; word i; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = FALSE; for (i = 1; i < argc; i++) if (argv[i][0] == '/') { word k; bool found = FALSE; for (k = 1; k <= MAX_PROP && !found; k++) if (!stricmp(argv[i] + 1, kinds[k].ext)) { problem_kind = k; found = TRUE; } NO_DEBUG_P(("/:%s %s", argv[i] + 1, found ? "found" : "not found")); if (!found) if (!stricmp(argv[i] + 1,"C")) constr_con = TRUE; else if (!stricmp(argv[i] + 1,"AC")) constr_acycl = TRUE; else if (toupper(argv[i][1]) == 'D' && argv[i][2] == '=') { char *s = argv[i]+3; NO_DEBUG_P(("scan: '%s'\n", s)); constr_d = TRUE; while (isdigit(*s)) { int d = *s++ - '0'; NO_DEBUG_P(("scan2: '%s'\n", s)); if (isdigit(*s)) d = 10 * d + *s++ - '0'; if (d < MAX_NR_D) allowed_d[d] = TRUE; NO_DEBUG_P(("scan3: '%s'\n", s)); if (*s == ',') s++; NO_DEBUG_P(("scan4: '%s'\n", s)); } } else if (toupper(argv[i][1]) == 'F' && argv[i][2] == '=') { if (sscanf(argv[i] + 3, "%d", &from_id_path) == 1) { constr_from = TRUE; from_id_path--; } else { printf("Error: no number after `/F='\n"); error = TRUE; } } else if (toupper(argv[i][1]) == 'T' && argv[i][2] == '=') { if (sscanf(argv[i] + 3, "%d", &to_id_path) == 1) { constr_to = TRUE; to_id_path--; } else { printf("Error: no number after `/T='\n"); error = TRUE; } } else error = TRUE; } else if (argv[i][0] == '-') { char *opt_chars = argv[i] + 1; while (*opt_chars != '\0') { word i; bool found = FALSE; for (i = 0; i < MAX_OPT && !found; i++) if (toupper(*opt_chars) == options[i].letter) { sel_options |= 1 << options[i].value; found = TRUE; } error = !found; opt_chars++; } } else if (argv[i][0] == '=') { switch(argv[i][1]) { case 'o': sscanf(argv[i] + 2, "%d", &start_off); params = TRUE; break; case 'A': alt = TRUE; params = TRUE; break; case 'E': only_even = 1; params = TRUE; break; case 'b': sscanf(argv[i] + 2, "%d", &begin_at); params = TRUE; break; case 's': sscanf(argv[i] + 2, "%d", &step); params = TRUE; break; case 'e': sscanf(argv[i] + 2, "%d", &end_at); params = TRUE; break; case 'd': do_calc_deter = TRUE; break; } } else if (in_filename == NULL) in_filename = argv[i]; else if (out_filename == NULL) out_filename = argv[i]; else error = TRUE; if (error || in_filename == NULL) { word i; printf("\nUsage: count [/kind] [-options] []"); printf("\n\nWhere kind is one of:"); for (i = 1; i <= MAX_PROP; i++) printf("\n %-4s - %s", kinds[i].ext, kinds[i].name); printf("\n C - connected"); printf("\n AC - acyclic"); printf("\n D=d - degree d"); printf("\n F=f - (1,f) have degree 1"); printf("\n T=t - (n,t) have degree 1"); printf("\n"); printf("\nWhere options combination of:"); for (i = 0; i < MAX_OPT; i++) printf("\n %c - %s", options[i].letter, options[i].description); printf("\n"); return 1; } } NO_DEBUG_P(("problem_kind %d\n", problem_kind)); /* Processing special problem kinds: */ switch (problem_kind) { case PROP_HAMIL_CYCLE: constr_con = TRUE; constr_acycl = FALSE; case PROP_2_FAC: { int i; constr_d = TRUE; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = i == 2; } break; case PROP_DOMINO_T: { int i; constr_d = TRUE; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = i == 1; } break; case PROP_HAMIL_PATH: constr_con = TRUE; constr_acycl = TRUE; constr_d = TRUE; { int i; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = i == 1 || i == 2; } break; case PROP_SPAN_TREES: constr_con = TRUE; constr_acycl = TRUE; constr_d = TRUE; allowed_d[1] = TRUE; { int i; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = i > 0; } break; case PROP_SPAN_TREES13: constr_con = TRUE; constr_acycl = TRUE; constr_d = TRUE; { int i; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = i == 1 || i == 3; } break; } if (!constr_d) { int i; for (i = 0; i < MAX_NR_D; i++) allowed_d[i] = TRUE; } /* Corrections on conflicting properties: */ /* 950929 removed correction: if (constr_con) allowed_d[0] = FALSE; */ if (!allowed_d[1]) { constr_acycl = FALSE; constr_from = FALSE; constr_to = FALSE; } /* Calculating derived values: */ { byte i; for (i = 0; i < MAX_NR_D; i++) if (allowed_d[i]) maximum_d = i; } constr_con_acycl = constr_con || constr_acycl; if (constr_to && constr_con && constr_acycl) { byte i; bool path = allowed_d[2]; for (i = 3; i < MAX_NR_D && path; i++) path = !allowed_d[i]; constr_to_path = path; printf("constr_t_p : %s", constr_to_path ? "Yes" : "No"); } printf("Problem: "); if (constr_con) printf("/C "); if (constr_acycl) printf("/AC "); if (constr_from) printf("/F=%d ", from_id_path + 1); if (constr_to) printf("/T=%d ", to_id_path + 1); if (constr_d) { int i; printf("/D="); for (i = 0; i < MAX_NR_D; i++) if (allowed_d[i]) printf("%d,",i); printf("\n"); } /* Detecting special problem kinds */ { bool only2 = TRUE, only1 = TRUE; int i; problem_kind = PROP_UNDEFINED; for (i = 0; (only1 || only2) && i < MAX_NR_D; i++) { if (allowed_d[i] != (i == 1)) only1 = FALSE; if (allowed_d[i] != (i == 2)) only2 = FALSE; } if (only2) problem_kind = constr_con ? PROP_HAMIL_CYCLE : PROP_2_FAC; else if (only1) problem_kind = PROP_DOMINO_T; else { bool only12 = TRUE, only13 = TRUE; for (i = 0; (only12 || only13) && i < MAX_NR_D; i++) { if (allowed_d[i] != (i == 1 || i == 2)) only12 = FALSE; if (allowed_d[i] != (i == 1 || i == 3)) only13 = FALSE; } if (only12 && constr_acycl && constr_con) problem_kind = PROP_HAMIL_PATH; else if (constr_acycl && constr_con && !allowed_d[0]) problem_kind = only13 ? PROP_SPAN_TREES13 : PROP_SPAN_TREES; } } { word i; graph_name = CHAR_MALLOC(strlen(in_filename) + 1); strcpy(graph_name, in_filename); for (i = strlen(graph_name); i > 0 && graph_name[i] != '.'; i--); if (i > 0) graph_name[i] = '\0'; } if (out_filename == NULL) { word i; out_filename = CHAR_MALLOC(strlen(in_filename) + 6); strcpy(out_filename, in_filename); for (i = strlen(out_filename); i > 0 && out_filename[i] != '.'; i--); if (i == 0) i = strlen(out_filename); out_filename[i] = '.'; strcpy(out_filename + i + 1, kinds[problem_kind].ext); } if ( !strcmp(out_filename, "stdout") || (f = fopen(out_filename, OPTION_SEL(OPT_A_OUT) ? "a" : "w")) == NULL) f = stdout; f_error = f; if (!read_graph(in_filename, stdout)) { printf("\nERROR: Cannot open file '%s'\n", in_filename); return 1; } if (OPTION_SEL(OPT_P_INP)) print_graph(f, kinds[problem_kind].name, in_filename); find_automorphic_perms(nr_vertices); if (OPTION_SEL(OPT_P_PERM)) print_permutations(f, automorphic_perms); if (OPTION_SEL(OPT_D_GEN)) { trace_f = f; fprintf(f, "\nTrace of M-vertex generation\n"); } { pcc_t begin_pcc; empty_pcc(begin_pcc); find_p_to_m_vertex(begin_pcc); } { m_vertex_list_p cur_m_vertex = m_vertex_list; first_m_vertex = TRUE; while (cur_m_vertex != NULL) { find_m_vertices_from_m_vertex(cur_m_vertex->pcc); cur_m_vertex = cur_m_vertex->next; first_m_vertex = FALSE; } } if (OPTION_SEL(OPT_D_GEN)) { fprintf(f, "\n"); } if (OPTION_SEL(OPT_P_VERTICES)) print_m_vertex_list(f); if (OPTION_SEL(OPT_P_COUNTS)) print_counts(f, NULL, 100, FALSE); if (OPTION_SEL(OPT_P_SEQ)) print_sequence(f, graph_name, problem_kind); if ( (OPTION_SEL(OPT_P_MATRIX) || OPTION_SEL(OPT_P_SOL)) && m_vertex_list->next != NULL) { analyze_all_freq(); calculate_nr_ends(); /* print_freq(f); */ order_m_vertices(); } if (OPTION_SEL(OPT_P_MATRIX)) print_orderd_matrix(f); if (!OPTION_SEL(OPT_P_SEQ)) { fprintf(f, "%% with version %s with", VERSION); { int i; for (i = 1; i < argc; i++) fprintf(f, " %s", argv[i]); } fprintf(f, "\n"); } if (do_calc_deter) { if (OPTION_SEL(OPT_P_SOL) && create_orderd_matrix()) calc_det(f, !constr_con); } else if (OPTION_SEL(OPT_P_SOL)) { fprintf(f, "For $C_{%s}^{%s}$ %s", graph_name, kinds[problem_kind].code, nr_m_vertices < 6 ? "it is known that\n" : "we found\n"); if (params) { if (!calc_rec(f, begin_at, step, end_at, start_off, (bool)only_even, alt)) fprintf(f, "%%A recurence equation could not be found\n"); } else if (!calc_rec(f, -1, 1, -1, 0, FALSE, FALSE)) fprintf(f, "%%A recurrence equation could not be found\n"); } if (!OPTION_SEL(OPT_P_SEQ)) fprintf(f, "\n\n"); if (f != stdout) fclose(f); return 0; }