/* POTM ENTRY: ticketysplit */ /* Your Name: Vincent Goffin */ /* Your email: vjg@research.att.com */ /* compile instructions: make ticketysplit */ #include #include /* rand48 */ #include #include #include #include typedef unsigned char tup_t; /* use short for N_MAX > 67 */ /* * Solves the "King's Lottery" problem: The King draws 1 ticket of 7 numbers * from a pool of numbers. Print the smallest collection of tickets of * 7 numbers such that at least 1 ticket will match at least three of the * King's numbers, for every possible King's draw. * * The problem is reduced to 3 ditinct covering problems of the form: * * Given N (pool/3) numbers, find the smallest collection of "blocks" * (tickets) of * B (7) numbers such that all M-tuples (triples) of these * are covered by at least one block. * * References: * * I used ideas from the article "New Constructions for Covering Designs" * by D. M. Gordon, G. Kuperberg and O. Patashnik. * at http://sdcc12.ucsd.edu/~xm3dg/cover.html * * I used some code from Art Owen's netlib "oa" package for finite fields and * orthogonal latin squares. * * vgoffin@att.com * Sat Feb 28 21:49:15 EST 1998 */ /* problem parameters */ #define N_MAX 100 /* max size of number pool */ #define B_MAX 16 /* max size of block */ #define M_MAX 3 /* max size of match dimension */ #define min(a,b) ((a) < (b)? (a): (b)) int g_user_sys_time = 0; int g_half_way = 0; /* options struct */ typedef struct opt_s { int oa_tuck; /* bool */ int gr_accell; /* accellerator heuristic */ int gr_slower; /* correction to above heuristic */ long random_seed; int cycle; /* cycle seed covering */ int pack; /* pack tickets to reduce jokers */ int improve; /* toss final result */ } Opt; /* a block covering */ typedef struct Cov_s { int N; int B; int M; int maxtuples; int tuples_per_block; int **blocks; int nblocks; int maxblocks; Opt opt; } Cov; /* a search context */ typedef struct Con_s { int* block; /* the B entry block */ int* offset; /* array of offsets for speed */ char noffset; /* number of offsets used */ char level; /* all nums expected to be < 128 */ char begin; char end; char shared_tuples; char target_tuples; char accell; char slower; char only_one; char use_triple_mark; } Con; /* * tuple cover (used tuples) * U.tuples[i] = n if tuple @ i is covered by a n blocks */ struct Tuples_s { int M; int N; int ntuples; tup_t* tuples; tup_t* triple_mark; /* heuristic speed-up */ } U; /* function declarations */ void initOpt(Opt*); void fillSolution(Cov*); void greedySolution(Cov*); int packSolution(Cov*); int tossSolution(Cov*); void genCov(Cov*, Opt*, int N, int B, int M); void initCov(Cov*, Opt*, int N, int B, int M); void embedCov(Cov* pDestC, Cov*); void stretchCov(Cov* pDestC, Cov*); void importCov(Cov* toC, Cov* froC); void cycleCov(Cov*, int nCyc); int simplifyCov(Cov*); int scrubBlock(Cov*, int *block, int nB, int* new_block); void addBlock(Cov*, int* block); int removeBlock(Cov*, int pos); void repackCov(Cov*); int greedySearch(Cov*, Con ct); void lookaheadSharedTuples(Con* pCn, int* lahead); int countSharedTuples(Cov*, int* block, int nB); void removeAllBlocks(Cov*); void freeCov(Cov*); void clearTripleMark(); void recalcTupleArray(Cov*); int countAllTuples(Cov*); int markBlockTuples(Cov*, int* block, int nB); int clearBlockTuples(Cov*, int* block, int nB); int getUniqueTuples(Cov*, int pos, int* utuples); int partialTupleMatch(Cov*, int* utuples, int nUT, int j, int* jokPos, int* jokVal); void startBlock(int* block, int nB); int nextBlock(int* block, int nB, int nN, int* pIndex); void sortBlock(int* block, int nB); void dumpBlocks(Cov*, int offset); void markFlatTuples(tup_t* usedTuples, int* flat, int nF, int nB, int nN); void incrContext(Con*, Cov*, int begin, int shared_tuples); /* timer & misc */ void initTimer(); int checkTime(); /* special constructions */ void cyclicCovering(Cov*); void finiteGeometryCovering(Cov*, char fg_type, int prime, int pow, int K); void orthogonalArrayCovering(Cov*, int B, int primepow); void combinedCovering(Cov*, int N, int ty, int B1, int M1, int B2, int M2); main(int argc, char* argv[]) { int N = 0; /* size of number pool */ int B = 7; /* block dimension */ int M = 3; /* match dimension, 2 or 3 */ int i; Opt opt; int split1, split2, split3; Cov cov1, cov2, cov3; if (argc != 2) { return 1; } N = atoi(argv[1]); if (N < 7 || N > 200) { return 1; } /* split pool into 3 smaller pools: split1 >= split2 >= split3 */ if (N <= 11) { split1 = 7; split2 = 0; split3 = 0; } else if (N <= 16) { split1 = 7; split2 = 7; split3 = 0; } else if (N <= 21) { split1 = 7; split2 = 7; split3 = 7; } else if (N < 27) { split1 = 7; split2 = (N - split1)/2; split3 = N - split1 - split2; } else { split1 = N/3; split2 = (N - split1)/2; split3 = N - split1 - split2; } if (split1 < split2) { i = split2, split2 = split1, split1 = i; } if (split2 < split3) { i = split2, split2 = split3, split3 = i; } if (split1 < split3) { i = split3, split3 = split1, split1 = i; } initTimer(); initOpt(&opt); /* big one first */ genCov(&cov1, &opt, split1, B, M); dumpBlocks(&cov1, 0); g_half_way = 1; if (split2 >= 3) { if (split2 == split1) { cov2 = cov1; } else { genCov(&cov2, &opt, split2, B, M); } dumpBlocks(&cov2, N - split3 - split2); } if (split3 >= 3) { if (split3 == split2) { cov3 = cov2; } else { genCov(&cov3, &opt, split3, B, M); } dumpBlocks(&cov3, N - split3); } return 0; } /* * genCov() * generate the covering; * potentially recursive as large coverings may combine smaller ones */ void genCov(Cov* pC, Opt* pO, int N, int B, int M) { Cov tmpC; Opt opt; opt = *pO; initCov(pC, &opt, N, B, M); tmpC.opt = opt; switch (M) { case 2: switch (B) { case 4: if (N == 15) { opt.random_seed += 1L; opt.improve += 1; } break; case 5: if (N == 15) { finiteGeometryCovering(&tmpC, 'p', 2, 4, 2); importCov(pC, &tmpC); freeCov(&tmpC); opt.cycle += 1; opt.improve += 1; } else if (N == 25) { finiteGeometryCovering(&tmpC, 'a', 5, 3, 1); importCov(pC, &tmpC); freeCov(&tmpC); } break; case 6: if (N == 31) { finiteGeometryCovering(&tmpC, 'p', 5, 3, 1); importCov(pC, &tmpC); freeCov(&tmpC); } break; } break; case 3: switch (B) { case 5: if (N == 10) { cyclicCovering(pC); opt.improve += 1; } break; case 7: if (N < 13) { ; } else if (N == 13) { cyclicCovering(pC); } else if (N < 16) { finiteGeometryCovering(&tmpC, 'p', 2, 4, 2); importCov(pC, &tmpC); freeCov(&tmpC); } else if (N == 16) { combinedCovering(pC, 'e', 15, 6, 2, 7, 3); opt.random_seed += 13; opt.improve += 1; } else if (N == 17) { combinedCovering(pC, 's', 15, 5, 2, 7, 3); } else if (N == 18) { combinedCovering(pC, 's', 15, 4, 2, 7, 3); /* 19-22 */ } else if (N < 23) { finiteGeometryCovering(&tmpC, 'a', 3, 4, 2); importCov(pC, &tmpC); freeCov(&tmpC); if (N == 19) { opt.cycle += 17; opt.improve += 4; opt.random_seed += 20; } else if (N == 20) { opt.cycle += 7; opt.improve += 2; } else if (N == 21) { opt.random_seed += 10; opt.cycle += 19; opt.improve += 1; } else if (N == 22) { /* TBD: > 5mins! */ opt.cycle += 16; opt.random_seed += 10; opt.improve += 3; } /* 23-28 */ } else if (N < 29) { /* derived steiner sys */ genCov(&tmpC, &opt, 22, 6, 3); importCov(pC, &tmpC); if (N == 23) { opt.cycle += 9; opt.improve += 4; opt.random_seed += 22; } else if (N == 24) { opt.cycle += 15; opt.improve += 2; opt.random_seed += 12; } else if (N == 25) { /* re-embed */ embedCov(pC, &tmpC); opt.cycle += 22; opt.improve += 5; opt.gr_accell += 1; } else if (N == 26) { embedCov(pC, &tmpC); opt.cycle += 11; opt.improve += 1; } else if (N == 27) { embedCov(pC, &tmpC); opt.cycle += 12; opt.improve += 1; opt.random_seed += 10; } else if (N == 28) { embedCov(pC, &tmpC); opt.cycle += 15; } freeCov(&tmpC); /* 29-31 */ } else if (N < 32) { finiteGeometryCovering(&tmpC, 'p', 2, 5, 2); importCov(pC, &tmpC); freeCov(&tmpC); if (N == 29) { opt.cycle += 18; opt.improve += 1; } } else if (N == 32) { combinedCovering(pC, 's', 31, 6, 2, 7, 3); opt.pack = 0; } else if (N == 33) { combinedCovering(pC, 's', 31, 5, 2, 7, 3); opt.pack = 0; } else if (N == 34) { combinedCovering(pC, 's', 31, 4, 2, 7, 3); opt.pack = 0; /* 35-41 */ } else if (N < 42) { orthogonalArrayCovering(&tmpC, B, 5); importCov(pC, &tmpC); /* re-embed */ embedCov(pC, &tmpC); freeCov(&tmpC); opt.gr_accell += (N != 35); opt.gr_slower += (N != 35); /* 42-52 */ } else if (N < 53) { tmpC.opt.oa_tuck = 1; orthogonalArrayCovering(&tmpC, B, 7); importCov(pC, &tmpC); /* re-embed */ embedCov(pC, &tmpC); freeCov(&tmpC); opt.gr_accell += (N >= 51); opt.gr_slower += (N >= 51); /* 53-58 */ } else if (N < 59) { tmpC.opt.oa_tuck = 1; orthogonalArrayCovering(&tmpC, B, 8); importCov(pC, &tmpC); embedCov(pC, &tmpC); freeCov(&tmpC); opt.gr_accell += (N >= 57); opt.gr_slower += (N >= 57); /* 59-62 */ } else if (N < 63) { tmpC.opt.oa_tuck = 1; orthogonalArrayCovering(&tmpC, B, 9); importCov(pC, &tmpC); embedCov(pC, &tmpC); freeCov(&tmpC); } else if (N == 63) { finiteGeometryCovering(&tmpC, 'p', 2, 6, 2); importCov(pC, &tmpC); freeCov(&tmpC); } else if (N == 64) { combinedCovering(pC, 's', 63, 6, 2, 7, 3); opt.pack = 0; } else if (N == 65) { combinedCovering(pC, 's', 63, 5, 2, 7, 3); opt.pack = 0; } else if (N == 66) { combinedCovering(pC, 's', 63, 4, 2, 7, 3); opt.pack = 0; } else if (N == 67) { combinedCovering(pC, 's', 65, 5, 2, 7, 3); opt.pack = 0; } break; } break; } recalcTupleArray(pC); pC->opt = opt; fillSolution(pC); simplifyCov(pC); } /* * combinedCovering() * linear combination of 2 smaller coverings to produce a larger one */ void combinedCovering(Cov* pC, int ty, int N, int B1, int M1, int B2, int M2) { Cov aC, bC; Opt opt; initOpt(&opt); genCov(&aC, &opt, N, B1, M1); genCov(&bC, &opt, N, B2, M2); importCov(pC, &bC); if (opt.cycle) { cycleCov(pC, opt.cycle); } if (ty == 's') { stretchCov(pC, &aC); } else { embedCov(pC, &aC); } freeCov(&aC); freeCov(&bC); } void fillSolution(Cov* pC) { Opt opt = pC->opt; int i; srand48(opt.random_seed); if (opt.cycle) { cycleCov(pC, opt.cycle); } if (U.ntuples < pC->maxtuples) { greedySolution(pC); } if (opt.pack) { while (packSolution(pC)) { } } for (i = 0; i < opt.improve; ++i) { if (!tossSolution(pC)) { break; } } } /* * greedySolution() * perform the greedy search * including the multipass speedup heuristic using accell & slower options */ void greedySolution(Cov* pC) { Con cn; int block[B_MAX]; int offset[512]; /* re-dim if needed */ memset(&cn, 0, sizeof(cn)); cn.block = block; cn.offset = offset; cn.end = pC->N - (pC->B - 1); cn.accell = pC->opt.gr_accell; cn.slower = pC->opt.gr_slower; cn.use_triple_mark = (cn.accell > 0); cn.only_one = 0; cn.target_tuples = min(pC->tuples_per_block, pC->maxtuples - U.ntuples); do { while (cn.target_tuples > 0) { greedySearch(pC, cn); if (U.ntuples >= pC->maxtuples) { break; } cn.target_tuples--; } if (cn.use_triple_mark) { clearTripleMark(); } if (U.ntuples >= pC->maxtuples) { break; } /* * += 1 does well, * += 2 can do better but at a much steeper cost */ cn.target_tuples += 1; } while (cn.accell-- != 0); } void initOpt(Opt* pO) { memset(pO, 0, sizeof(Opt)); pO->random_seed = 1L; /* several constructions depend on this value */ pO->pack = 1; } int calcTuples(int n, int m) /* (n*(n-1)*(n-2) ...) /(1*2*3 ...) */ { int i, binCoef = n; for (i = 1; i < m; ++i) { binCoef *= (n - i); binCoef /= (i + 1); } return binCoef; } void initCov(Cov* pC, Opt* pO, int N, int B, int M) { memset(pC, 0, sizeof(Cov)); pC->N = N; pC->B = B; pC->M = M; pC->maxtuples = calcTuples(N, M); pC->tuples_per_block = calcTuples(B, M); if (pO && pO != &pC->opt) { pC->opt = *pO; } } /* * greedySearch() * recursively generate ordered blocks and test each block for fit. * This routine is usually called within a loop with ever decreasing * values of target_tuples that ensures the greedy best fit. * return 1 if only one fit is needed; else continue generating * blocks to build a complete cover. */ int greedySearch(Cov* pC, Con cn) { int i; int all_shared_tuples, nblocks, shared_lahead[N_MAX]; Con new_cn; if (U.ntuples >= pC->maxtuples) { return 1; } if (cn.level == pC->B) { /* * bottom of recursion; * block fits as well as is required, add it to list */ int newtuples; int cpy_block[B_MAX]; scrubBlock(pC, cn.block, pC->B, cpy_block); newtuples = markBlockTuples(pC, cpy_block, pC->B); U.ntuples += newtuples; addBlock(pC, cpy_block); return cn.only_one? 1: 0; } lookaheadSharedTuples(&cn, shared_lahead); for (i = cn.begin; i < cn.end; ++i) { all_shared_tuples = cn.shared_tuples + shared_lahead[i]; if (pC->tuples_per_block - all_shared_tuples < cn.target_tuples) { continue; } cn.block[cn.level] = i; /* * keep current nblocks as a change indicator */ nblocks = pC->nblocks; /* * recurse to fill the block entry by entry; * but only if there's a chance of meeting the shared tuple target * (cn.slower = 1 for a better, but slower, covering) */ if (cn.use_triple_mark && cn.level + 1 == 3 + cn.accell) { int* b = cn.block; int toff = (b[0]*pC->N + b[1])*pC->N + b[2]; if (all_shared_tuples + U.triple_mark[toff] > pC->tuples_per_block - cn.target_tuples + cn.slower) { continue; } } new_cn = cn; /* binary copy: keep cn.block, cn.offset */ incrContext(&new_cn, pC, i, all_shared_tuples); if (greedySearch(pC, new_cn) || U.ntuples >= pC->maxtuples) { return 1; } if (pC->nblocks > nblocks) { if (i < cn.end - 1) { /* * a new block was generated; * review shared counts for the current block */ cn.shared_tuples = countSharedTuples(pC, cn.block, cn.level); cn.begin = i; lookaheadSharedTuples(&cn, shared_lahead); } } else if (cn.use_triple_mark && cn.level + 1 == 3) { /* * did not produce an insertion so we have a lower limit * on the availability of free tuples */ int* b = cn.block; int toff = (b[0]*pC->N + b[1])*pC->N + b[2]; U.triple_mark[toff] = pC->tuples_per_block - cn.target_tuples - all_shared_tuples + 1; } } /* cn.block[cn.level] = pC->N; */ return 0; } void addBlock(Cov* pC, int* block) /* block must have length to match what pC expects */ { int i, j; int *p; if (pC->nblocks >= pC->maxblocks) { if (pC->maxblocks == 0) { pC->maxblocks = 256; pC->blocks = (int **)malloc(pC->maxblocks*sizeof(int*)); } else { pC->maxblocks *= 2; pC->blocks = (int **)realloc(pC->blocks, pC->maxblocks*sizeof(int*)); } } p = (int *)malloc(pC->B*sizeof(int)); memcpy(p, block, pC->B*sizeof(int)); pC->blocks[pC->nblocks] = p; pC->nblocks++; } /* * lookaheadSharedTuples() * look ahead to calculate the coming shared tuples */ void lookaheadSharedTuples(Con* pCn, int* lahead) { int i, j; memset(lahead + pCn->begin, 0, sizeof(int)*(pCn->end - pCn->begin)); for (i = 0; i < pCn->noffset; ++i) { tup_t* tuples = U.tuples + pCn->offset[i]; for (j = pCn->begin; j < pCn->end; ++j) { lahead[j] += (tuples[j] != 0); } } } /* * doCountSharedTuples() * countSharedTuples() * count how many of this block's tuples are already covered. */ int doCountSharedTuples(int* block, int nB, int imin, int imax, int nN, int toff) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { count += (tuples[block[j]] > 0); } } } else { count += doCountSharedTuples(block, nB, i+1, imax+1, nN, ioff); } } } return count; } int countSharedTuples(Cov* pC, int* block, int nB) { return doCountSharedTuples(block, nB, 0, nB - pC->M + 1, pC->N, 0); } /* * doMarkBlockTuples() * markBlockTuples() * mark all the tuples for this block by incrementing the tuple counter array. */ int doMarkBlockTuples(int* block, int nB, int imin, int imax, int nN, int toff) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { count += (tuples[block[j]] == 0); tuples[block[j]]++; } } } else { count += doMarkBlockTuples(block, nB, i+1, imax+1, nN, ioff); } } } return count; } int markBlockTuples(Cov* pC, int* block, int nB) { return doMarkBlockTuples(block, nB, 0, nB - pC->M + 1, pC->N, 0); } /* * doclearBlockTuples() * clearBlockTuples() * clear all tuples associated with this block from the g_tuples array * return the number of tuple holes left. */ int doClearBlockTuples(int* block, int nB, int imin, int imax, int nN, int toff) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { tuples[block[j]]--; count += (tuples[block[j]] == 0); } } } else { count += doClearBlockTuples(block, nB, i+1, imax+1, nN, ioff); } } } return count; } int clearBlockTuples(Cov *pC, int* block, int nB) { return doClearBlockTuples(block, nB, 0, nB - pC->M + 1, pC->N, 0); } /* * countAllTuples() * docountAllTuples() * go thru the tuple array and count the non-zero entries */ int doCountAllTuples(int imin, int imax, int nN, int toff) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { int ioff = (toff + i)*nN; if (imax == nN-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nN; ++j) { count += (tuples[j] > 0); } } else { count += doCountAllTuples(i+1, imax+1, nN, ioff); } } return count; } int countAllTuples(Cov* pC) { return doCountAllTuples(0, pC->N - pC->M + 1, pC->N, 0); } /* * doNeededEntries() * scrubBlock() * identify unneeded entries in block (jokers) and set them to pC->N; * resort ticket if needed */ int doNeededEntries(int* block, int nB, int imin, int imax, int nN, int toff, int* needed, int* hist, int nhist) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; hist[nhist] = i; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { if (tuples[block[j]] == 0) { int h; needed[j]++; for (h = 0; h <= nhist; ++h) { needed[hist[h]]++; } count++; } } } } else { count += doNeededEntries(block, nB, i+1, imax+1, nN, ioff, needed, hist, nhist+1); } } } return count; } int scrubBlock(Cov* pC, int* block, int nB, int* cpy_block) { int i; int nN = pC->N; int needed[B_MAX]; int hist[B_MAX]; int newtuples = 0; int nscrubbed = 0; memset(needed, 0, sizeof(needed)); memcpy(cpy_block, block, nB*sizeof(int)); newtuples = doNeededEntries(block, nB, 0, nB-pC->M+1, nN, 0, needed, hist, 0); /* remove entries not needed */ for (i = 0; i < nB; ++i) { if (!needed[i]) { if (cpy_block[i] != pC->N) { ++nscrubbed; } cpy_block[i] = pC->N; } } if (nscrubbed) { sortBlock(cpy_block, nB); } return newtuples; } /* * doSimplifyBlock() * simplifyBlock() * simplifyCov(Cov* pC) * double check if any jokers are present */ int doUnneededEntries(int* block, int nB, int imin, int imax, int nN, int toff, int* needed, int* hist, int nhist) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; hist[nhist] = i; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { if (tuples[block[j]] <= 1) { int h; needed[j]++; for (h = 0; h <= nhist; ++h) { needed[hist[h]]++; } count++; } } } } else { count += doUnneededEntries(block, nB, i+1, imax+1, nN, ioff, needed, hist, nhist+1); } } } return count; } int simplifyBlock(Cov* pC, int pos) { int i, j; int tupindex[128]; int ntupindex; int new_needed[B_MAX]; int needed[B_MAX]; int hist[B_MAX]; int nJokers = 0; int nN = pC->N; int nB = pC->B; int* block = pC->blocks[pos]; int blockmap[N_MAX]; memset(needed, 0, sizeof(needed)); doUnneededEntries(block, nB, 0, nB-pC->M+1, nN, 0, needed, hist, 0); /* jokerize unused numbers */ for (i = 0; i < nB; ++i) { if (needed[i] == 0) { nJokers++; break; } } if (nJokers) { clearBlockTuples(pC, block, nB); for (i = 0; i < nB; ++i) { if (needed[i] == 0) { block[i] = pC->N; ++nJokers; } } /* re-sort block */ sortBlock(block, nB); markBlockTuples(pC, block, nB); } return nJokers; } /* * transform as many entries as possible into jokers */ int simplifyCov(Cov* pC) { int i, nJokers = 0; for (i = 0; i < pC->nblocks; ++i) { if (pC->blocks[i]) { nJokers += simplifyBlock(pC, i); } } return nJokers; } void freeTupleArray() { if (U.tuples) { free(U.tuples); U.tuples = 0; } if (U.triple_mark) { free(U.triple_mark); U.triple_mark = 0; } U.ntuples = 0; U.N = U.M = 0; } /* * packSolution() * remove as many blocks as possible by making use of the joker entries * A JOKER is a block entry that can be set to an arbitrary value * without affecting the covering. */ int packSolution(Cov* pC) { int i, j, k, m, n; int jokPos, jokVal; int block[B_MAX]; int nJokers = 0; int utuples[128]; int nUniqTuples; int *iB, *jB; int nRemoved = 0; if (!(pC->M == 3 || pC->M == 2)) { return 0; } /* * start from last block with the least unique-tuples. * Identify the unique elements and try to un-unique them * by setting joker entries in other blocks. */ for (i = pC->nblocks - 1; i >= 0; --i) { if (!(iB = pC->blocks[i])) { /* = */ continue; } if ((nJokers = simplifyCov(pC)) == 0) { break; } /* find the unique tuples in this block */ nUniqTuples = getUniqueTuples(pC, i, utuples); if (nUniqTuples == 0) { /* * can happen because of block rearrangements farther down */ removeBlock(pC, i); ++nRemoved; if (++nRemoved % 10 == 0) { repackCov(pC); } continue; } /* * check how many entries can be removed using jokers from other cards */ for (j = pC->nblocks - 1; j >= 0 && nUniqTuples != 0; --j) { /* keep only blocks with jokers */ if (j == i || !(jB = pC->blocks[j]) || jB[pC->B-1] != pC->N) { continue; } /* * if block j contains 2 numbers from at least one * unique tuple of block i, it can be used to reduce tiket i's * unique tuple count by at least 1 */ if (partialTupleMatch(pC, utuples, nUniqTuples, j, &jokPos, &jokVal)) { jB = pC->blocks[j]; clearBlockTuples(pC, jB, pC->B); jB[jokPos] = jokVal; sortBlock(jB, pC->B); markBlockTuples(pC, jB, pC->B); nUniqTuples = getUniqueTuples(pC, i, utuples); /* re-use the j-block if it has more than 1 joker */ ++j; /* backwards */ } } if (nUniqTuples == 0) { removeBlock(pC, i); iB = 0; if (++nRemoved % 10 == 0) { repackCov(pC); } /* } else { clearBlockTuples(pC, iB, pC->B); scrubBlock(pC, iB, pC->B, block); memcpy(iB, block, pC->B*sizeof(int)); markBlockTuples(pC, iB, pC->B); */ } } repackCov(pC); return nRemoved; } /* * partialTupleMatch() * ticket at bpos has joker entries. see if by setting the joker we can * match any tuples in utuples[]. * TBD: needs to be generalized wrt pC->M */ int partialTupleMatch(Cov* pC, int* utuples, int nUT, int bpos, int* jokPos, int* jokVal) { int u, i, j; int a, b, c, d; int match = 0; int* block = pC->blocks[bpos]; int nN = pC->N; int nJok = 0; for (i = pC->B-1; i >= 0 && block[i] == nN; --i) { ++nJok; } for (u = 0; u < nUT; ++u) { int toff = utuples[u]; switch (pC->M) { case 4: break; case 3: /* * remember that the blocks are sorted */ a = (toff/(nN*nN)) % nN; b = (toff/nN) % nN; c = toff % nN; for (i = 0; i < pC->B-nJok-1; ++i) { if (block[i] != a && block[i] != b) { continue; } for (j = i+1; j < pC->B-nJok; ++j) { if (block[j] == b || block[j] == c) { match = 1; /* by setting the joker we can make one of the unique tuples of i not-unique anymore we still have to figure out the joker setting */ if (block[i] != a) { *jokVal = a; } else if (block[j] != c) { *jokVal = c; } else { *jokVal = b; } break; } } } break; case 2: a = (toff/nN) % nN; b = toff % nN; for (i = 0; i < pC->B - 1; ++i) { if (block[i] == a || block[i] == b) { match = 1; /* by setting the joker we can make one of the unique tuples of i not-unique anymore we still have to figure out the joker setting */ if (block[i] == a) { *jokVal = b; } else { *jokVal = a; } break; } } break; } } if (match) { /* return first joker in B */ *jokPos = pC->B - nJok; return 1; } return 0; } /* * tossSolution() * remove random block(s) and try again... * return 1 at first improvement, 0 after timeout otherwise. */ int tossSolution(Cov* pC) { int i, j, n; int delta; Cov impC; int utuples[128]; int nUniqTuples; int block[B_MAX]; int cpy_block[B_MAX]; int timeLimit = (g_half_way? 590 : 300); initCov(&impC, &pC->opt, pC->N, pC->B, pC->M); importCov(&impC, pC); for (i = 0; checkTime() < timeLimit; ++i) { n = drand48()*(impC.nblocks - .5); memcpy(block, impC.blocks[n], impC.B*sizeof(int)); U.ntuples -= removeBlock(&impC, n); repackCov(&impC); scrubBlock(&impC, block, impC.B, cpy_block); U.ntuples += markBlockTuples(&impC, cpy_block, impC.B); addBlock(&impC, cpy_block); while (packSolution(&impC) > 0) { ; } /* * return at 1st improvement */ if (impC.nblocks < pC->nblocks) { importCov(pC, &impC); freeCov(&impC); return 1; } drand48(); /* compatibility... */ } freeCov(&impC); return 0; } /* * doGetUniqueTuples() * getUniqueTuples() * return the tuples that this block is the only one to cover. */ int doGetUniqueTuples(int* block, int nB, int imin, int imax, int nN, int toff, int* utuples) { int i, j; int count = 0; for (i = imin; i < imax; ++i) { if (block[i] < nN) { int ioff = (toff + block[i])*nN; if (imax == nB-1) { tup_t* tuples = U.tuples + ioff; for (j = i+1; j < nB; ++j) { if (block[j] < nN) { if (tuples[block[j]] == 1) { utuples[count++] = ioff + block[j]; } } } } else { count += doGetUniqueTuples(block, nB, i+1, imax+1, nN, ioff, utuples + count); } } } return count; } int getUniqueTuples(Cov* pC, int pos, int* utuples) { return doGetUniqueTuples(pC->blocks[pos], pC->B, 0, pC->B - pC->M + 1, pC->N, 0, utuples); } /* * remove the block from the list and clear its tuple marks. * leave a null pointer behind * the list must be subsequently re-packecd! * return the number of holes introduced in the tuple cover. */ int removeBlock(Cov* pC, int pos) { int count = 0; count = clearBlockTuples(pC, pC->blocks[pos], pC->B); free(pC->blocks[pos]); pC->blocks[pos] = 0; return count; } void repackCov(Cov* pC) { int i, j; for (i = 0, j = 0; i < pC->nblocks; ++i) { if (pC->blocks[i]) { pC->blocks[j] = pC->blocks[i]; j++; } } pC->nblocks = j; } /* * clear C * free block memory but not container */ void removeAllBlocks(Cov* pC) { int i; if (pC->blocks) { for (i = 0; i < pC->nblocks; ++i) { if (pC->blocks[i]) { free(pC->blocks[i]); pC->blocks[i] = 0; } } } /* don't set blocks or maxblocks to 0! */ pC->nblocks = 0; } void freeCov(Cov* pC) { removeAllBlocks(pC); free(pC->blocks); memset(pC, 0, sizeof(Cov)); } /* * importCov() * copy froC blocks to toC, adjusting block length/content along the way */ void importCov(Cov* pToC, Cov* pFroC) { int i, j; int block[B_MAX]; int minB = pFroC->B; if (pToC->B < minB) { minB = pToC->B; } /* clear previous blocks */ removeAllBlocks(pToC); for (i = 0; i < pFroC->nblocks; ++i) { for (j = 0; j < minB; ++j) { block[j] = pFroC->blocks[i][j]; if (block[j] > pToC->N) { block[j] = pToC->N; } } for (; j < pToC->B; ++j) { block[j] = pToC->N; } /* no need to resort block */ addBlock(pToC, block); } } /* * embedCov() * general purpose way of fitting a smaller solution into a larger space * TBD: in several cases this does less well than a stretchCov()... why? */ void embedCov(Cov* pDestC, Cov* pC) { int i, j, k, l; int block[B_MAX]; int offset[512]; /* increase as needed */ Con cn; if (pC->B >= pDestC->B) { return; /* not ready for this */ } memset(&cn, 0, sizeof(cn)); cn.block = block; cn.offset = offset; recalcTupleArray(pDestC); for (i = 0; i < pC->nblocks; ++i) { memcpy(block, pC->blocks[i], pC->B*sizeof(int)); /* force a fit by removing as many tail entries as needed */ for (j = pC->B-1; j >= 0; --j) { if (block[j] < pDestC->N - (pDestC->B - j - 1)) { break; } } cn.level = j + 1; cn.begin = block[j] + 1; cn.end = pDestC->N - (pDestC->B - cn.level - 1); cn.shared_tuples = countSharedTuples(pDestC, block, cn.level); cn.target_tuples = pDestC->tuples_per_block - cn.shared_tuples; cn.use_triple_mark = 0; cn.only_one = 1; cn.noffset = 0; /* * harder to initialize because we don't start from sratch */ switch(pDestC->M) { case 3: for (k = 0; k < cn.level-1; ++k) { for (l = k+1; l < cn.level; ++l) { cn.offset[cn.noffset++] = (cn.block[k]*pDestC->N + cn.block[l])*pDestC->N; } } break; case 2: for (k = 0; k < cn.level; ++k) { cn.offset[cn.noffset++] = cn.block[k]*pDestC->N; } break; } while (!greedySearch(pDestC, cn)) { if (--cn.target_tuples == 0) { break; } } } } /* * stretchCov() * extend each ticket of pC and add it to pDestC */ void stretchCov(Cov* pDestC, Cov* pC) { int block[B_MAX]; int i, j, k; if (pC->B >= pDestC->B) { return; /* not ready for this */ } for (i = 0; i < pC->nblocks; ++i) { memcpy(block, pC->blocks[i], pC->B*sizeof(int)); k = pC->N; for (j = 0; j < pC->B; ++j) { if (block[j] == k) { /* 1st joker */ break; } } for (; j < pDestC->B; ++j) { block[j] = k; if (k < pDestC->N) { ++k; } } addBlock(pDestC, block); } } /* * cycleCov() * just cycle the lettering of the covering -- ideally a noop but if * a greedy search is performed to fill the coverage, this can * have an effect */ void cycleCov(Cov* pC, int nCyc) { int i, j; int block[B_MAX]; for (i = 0; i < pC->nblocks; ++i) { memcpy(block, pC->blocks[i], pC->B*sizeof(int)); for (j = 0; j < pC->B; ++j) { if (block[j] != pC->N) { pC->blocks[i][j] = (block[j] + pC->N + nCyc) % pC->N; } } sortBlock(pC->blocks[i], pC->B); } recalcTupleArray(pC); } /*---------------------------------------------------------------------------- * Tuple Array Routines */ void clearTupleArray() { int tsize = ipow_sum(U.N, U.M)*sizeof(tup_t); memset(U.tuples, 0, tsize); } void newTupleArray(int N, int M) { int tsize; /* remove previous */ freeTupleArray(); U.N = N; U.M = M; tsize = ipow_sum(N, M)*sizeof(tup_t); U.tuples = (tup_t *)malloc(tsize); memset(U.tuples, 0, tsize); /* * used by all M */ tsize = ipow_sum(N, 3)*sizeof(tup_t); U.triple_mark = (tup_t *)malloc(tsize); memset(U.triple_mark, 0, tsize); } /* * recalcTupleArray() * use the provided block suite to rebuild the tuple array from scratch */ void recalcTupleArray(Cov* pC) { int i; if (U.N != pC->N || U.M != pC->M) { newTupleArray(pC->N, pC->M); } clearTupleArray(); U.ntuples = 0; for (i = 0; i < pC->nblocks; ++i) { U.ntuples += markBlockTuples(pC, pC->blocks[i], pC->B); } } /* * clearTripleMark() * suports a heuristic used in greedy algor. */ void clearTripleMark() { int tsize = ipow_sum(U.N, 3)*sizeof(tup_t); memset(U.triple_mark, 0, tsize); } /*---------------------------------------------------------------------------- * Cyclic covering */ void addCyclicBlocks(Cov* pC, int* block); void cyclicCovering(Cov* pC) /* Generate the best possible cyclic covering. */ { int block[B_MAX]; int sav_block[B_MAX]; int most_tuples = 0; startBlock(block, pC->B); while (U.ntuples < pC->maxtuples && nextBlock(block, pC->B, pC->N, 0)) { if (block[0] != 0) { break; /* translation */ } removeAllBlocks(pC); recalcTupleArray(pC); addCyclicBlocks(pC, block); if (U.ntuples > most_tuples) { most_tuples = U.ntuples; memcpy(sav_block, block, sizeof(block)); } } if (U.ntuples < pC->maxtuples) { removeAllBlocks(pC); recalcTupleArray(pC); addCyclicBlocks(pC, sav_block); } } void addCyclicBlocks(Cov* pC, int* block) { int i, j; int cyc_block[B_MAX]; int cpy_block[B_MAX]; memcpy(cyc_block, block, pC->B*sizeof(int)); for (i = 0; i < pC->N; ++i) { sortBlock(cyc_block, pC->B); scrubBlock(pC, cyc_block, pC->B, cpy_block); U.ntuples += markBlockTuples(pC, cpy_block, pC->B); addBlock(pC, cpy_block); for (j = 0; j < pC->B; ++j) { cyc_block[j] = (cyc_block[j] + 1) % pC->N; } } } /*----------------------------------------------------------------------------- * Finite Field Construction * The code to generate finite fiels comes from the oa package in netlib, * by Art Owen. An original header comment is reproduced below. */ /* These programs construct and manipulate orthogonal arrays. They were prepared by Art Owen Department of Statistics Sequoia Hall Stanford CA 94305 They may be freely used and shared. This code comes with no warranty of any kind. Use it at your own risk. I thank the Semiconductor Research Corporation and the National Science Foundation for supporting this work. */ /* * used by finite geometries and orthogonal arrays */ #define Q_MAX 125 /* max size of prime pow */ #define P_MAX 16 /* max degree of any polynomial */ /* * expression for x^n = polynomial_of_order_less_than_n * These are primitive irreducible polynomials that are * used to define finite field (gf) "extensions" on finite subfields. * The subfields are of prime order, the etensions aore of prime power order * the irred polynomial */ int p2n2[] = { 1, 1 }; /* gf(2^2) = gf(4) */ int p2n3[] = { 1, 0, 1 }; /* gf(2^3) = gf(8) */ int p2n4[] = { 1, 0, 0, 1 }; /* gf(2^4) = gf(16) */ int p2n5[] = { 1, 0, 0, 1, 0 }; /* gf(2^5) = gf(32) */ int p2n6[] = { 1, 0, 0, 0, 0, 1 }; /* gf(2^6) = gf(64) */ int p3n2[] = { 1, 2 }; /* gf(3^2) = gf(9) */ int p3n3[] = { 2, 0, 1 }; /* gf(3^3) = gf(27) */ int p3n4[] = { 1, 0, 0, 2 }; /* gf(3^4) = gf(81) */ int p5n2[] = { 3, 4 }; /* gf(5^2) = gf(25) */ int p5n3[] = { 3, 0, 4 }; /* gf(5^3) = gf(125) */ int p7n2[] = { 4, 6 }; /* gf(7^2) = gf(49) */ void int_to_poly(int q, int n, int k, int pol[]); int poly_to_int(int q, int n, int pol[]); void poly_sum(int p, int n, int* p1, int* p2, int* sum); void poly_prod(int p, int n, int* xton, int* p1, int* p2, int* prod); int *xton; int *ffe[Q_MAX], ffebuf[Q_MAX*Q_MAX]; int *add[Q_MAX], addbuf[Q_MAX*Q_MAX]; int *mul[Q_MAX], mulbuf[Q_MAX*Q_MAX]; /* * initFiniteField(): * initialize a finite field on q elements; * there are only finite field for primes and prime powers. * return 1 if there was a finite field and it was initialized, * else return 0. */ int initFiniteField(int q) { /* q: number of symbols */ int p; /* prime */ int n; /* power: q = p^n */ int i, j, k; int click; int poly[P_MAX]; /* make a design of size q*q*q over q symbols, strength 3 and ncol columns using a Bush construction based on finite field extensions */ /* note: the expressions don't require a polynomial rep for fields of prime order, and are much simpler */ switch (q) { /* fields of prime order: none used */ case 2: case 3: case 5: case 6: case 7: case 11: case 13: case 17: xton = 0; p = q, n = 1; break; /* extended fields: */ case 4: xton = p2n2, p = 2, n = 2; break; case 8: xton = p2n3, p = 2, n = 3; break; case 9: xton = p3n2, p = 3, n = 2; break; case 16: xton = p2n4, p = 2, n = 4; break; case 25: xton = p5n2, p = 5, n = 2; break; case 27: xton = p3n3, p = 3, n = 3; break; case 32: xton = p2n5, p = 2, n = 5; break; case 49: xton = p7n2, p = 7, n = 2; break; case 64: xton = p2n6, p = 2, n = 6; break; case 81: xton = p3n4, p = 3, n = 4; break; case 125: xton = p5n3, p = 5, n = 3; break; default: return 0; } for (i = 0; i < q; i++) { ffe[i] = ffebuf + i*q; add[i] = addbuf + i*q; mul[i] = mulbuf + i*q; } /* generate ffe, the polynomial representation of each element of the extended finite field */ memset(ffebuf, 0, q*q*sizeof(int)); for (i = 1; i < q; i++) { for (click = 0; ffe[i-1][click] == (p-1); click++) { ffe[i][click] = 0; } ffe[i][click] = ffe[i-1][click] + 1; for (j = click+1; j < n; j++) { ffe[i][j] = ffe[i-1][j]; } } /* * build the add and mul tables by adding and muling the polynomials, * then converting the result back to ints */ for (i = 0; i < q; i++) { for (j = 0; j < q; j++) { poly_sum(p, n, ffe[i], ffe[j], poly); add[i][j] = poly_to_int(p, n, poly); poly_prod(p, n, xton, ffe[i], ffe[j], poly); mul[i][j] = poly_to_int(p, n, poly); } } return 1; } void int_to_poly(int p, int n, int k, int pol[]) { int i; for (i = 0; i < n; ++i) { pol[i] = k % p; k /= p; } } int poly_to_int(int p, int n, int pol[]) { int i; int k = 0; for (i = 1; i <= n; ++i) { k *= p; k += pol[n-i]; } return k; } void poly_sum(int p, int n, int* p1, int* p2, int* sum) { int i; for (i = 0; i < n; i++) { sum[i] = (p1[i] + p2[i]) % p; } } /* * poly_prod(): * Set prod = p1*p2 with coefficients modulo p, * and use xton to reduce x^n and higher powers to lower ones. */ void poly_prod(int p, int n, int* xton, int* p1, int* p2, int* prod) { int i, j; int expandedProd[32]; /* expanded product, using powers > x^n */ for (i = 0; i < 2*n-1; i++) { expandedProd[i] = 0; } for (i = 0; i < n; i++) { for (j = 0; j < n; j++) { expandedProd[i+j] += p1[i]*p2[j]; } } for (i = 2*n-2; i>n-1; i--) { for (j = 0; j < n; j++) { expandedProd[i-n+j] += xton[j]*expandedProd[i]; } } for (i = 0; i < n; i++) { prod[i] = expandedProd[i] % p; } } /*---------------------------------------------------------------------------- * Finite geometry covering */ #define F_MAX B_MAX /* size of flat won't exceed size of block */ int initFiniteField(int q); void flatFromPoints(Cov* pC, int K, int p, int* gPts, int nGPts, int nFPts); int linearDepend(int testPt, int* basePts, int nBPts, int prim); int ipow(int a, int b); int ipow_sum(int a, int b); int qBinomialCoeff(int q, int n, int k); int* gf_to_gm; /* geom index to field value */ int* gm_to_gf; /* field value to geom index */ int* gf_equiv; /* field value's equivalence class representative */ tup_t *usedTuples; /* * Seed the initial block list with a covering generated from * k-flats of the (affine, projective) finite geometry on GF(n-1, p) */ void finiteGeometryCovering(Cov* pC, char fg_type, int p, int n, int K) { int i, j, k; int nGPts; /* points in geometry */ int nFPts; /* points per k-flat */ int nKF; /* number of k-flats */ int m = n-1; /* p^n numbers correspond to a GF(n-1, p) field */ int utSize; int changes; int N = ipow(p, n); int points[B_MAX]; int tIndex; initFiniteField(N); if (fg_type != 'a') { nGPts = ipow_sum(p, m); nFPts = ipow_sum(p, K); nKF = qBinomialCoeff(p, m+1, K+1); } else { nGPts = ipow(p, m); nFPts = ipow(p, K); nKF = ipow(p, m-K)*qBinomialCoeff(p, m, K); } initCov(pC, 0, nGPts, nFPts, K+1); /* find representatives of the points of the geometry; geom pts are equiv classes of non-null (m+1)-vectors over GF(p) */ gf_to_gm = (int*) malloc( N*sizeof(int) ); gm_to_gf = (int*) malloc( nGPts*sizeof(int) ); gf_equiv = (int*) malloc( N*sizeof(int) ); k = 0; gf_to_gm[0] = -1; /* skip the null vector */ gf_equiv[0] = 0; for (i = 1; i < N; ++i) { gf_equiv[i] = i; } changes = 1; while (changes) { changes = 0; for (i = 1; i < N; ++i) { for (j = 2; j < p; ++j) { /* over base field only */ if (mul[j][i] < gf_equiv[i]) { gf_equiv[i] = mul[j][i]; changes++; } } } } gf_to_gm[0] = -1; k = 0; for (i = 1; i < N; ++i) { if (fg_type=='a' && i%p == 0) { /* affine condition */ gf_to_gm[i] = -1; continue; } if (i == gf_equiv[i]) { gm_to_gf[k] = i; gf_to_gm[i] = k; ++k; } else { gf_to_gm[i] = gf_to_gm[gf_equiv[i]]; } } utSize = ipow_sum(nGPts, K+1)*sizeof(tup_t); usedTuples = (tup_t*) malloc(utSize); memset(usedTuples, 0, utSize); /* * generate K+1 points, numerically ordered * then for each tuple print out the k-flat. */ startBlock(points, K+1); while (nextBlock(points, K+1, nGPts, &tIndex)) { if (!usedTuples[tIndex]) { flatFromPoints(pC, K, p, points, nGPts, nFPts); } } free(usedTuples); free(gf_to_gm); free(gm_to_gf); } void flatFromPoints(Cov* pC, int K, int p, int* gPts, int nGPts, int nFPts) /* k+1 independent points determine a k-flat; a k-flat contains a number of points; given k+1 points, determine all the points in the k-flat; */ { int i, t = 0; int flat[F_MAX]; for (i = 0; i < nGPts; i++) { if (linearDepend(i, gPts, K+1, p)) { flat[t++] = i; } } /* save only non-degenerate k-flats */ if (t == nFPts) { markFlatTuples(usedTuples, flat, nFPts, K+1, nGPts); sortBlock(flat, pC->B); addBlock(pC, flat); } } /*---------------------------------------------------------------------------- * Linear dependence over a finite field */ int doLinearDepend(int testPt, int* basePts, int nBPts, int prim, int level, int fval) { int i; if (level == nBPts) { return (testPt == gf_to_gm[fval]); } for (i = 0; i < prim; ++i) { int gfPt = gm_to_gf[basePts[level]]; int fv = add[fval][ mul[gfPt][i] ]; if (doLinearDepend(testPt, basePts, nBPts, prim, level + 1, fv)) { return 1; } } return 0; } int linearDepend(int testPt, int* basePts, int nBPts, int prim) { return doLinearDepend(testPt, basePts, nBPts, prim, 0, 0); } int ipow(int a, int b) /* ipow() = a^b */ { int i, p = 1; for (i = 0; i < b; ++i) { p *= a; } return p; } int ipow_sum(int a, int b) /* ipow_sum() = 1 + a + a^2 + ... a^b */ { int i, s = 1; for (i = 0; i < b; ++i) { s *= a; s += 1; } return s; } int qBinomialCoeff(q, n, k) /* q binomial coefficient */ { int i, coeff = 1; for (i = n; i >= n-(k-1); --i) { coeff *= ipow(q, i) - 1; } for (i = k; i >= 1; --i) { coeff /= ipow(q, i) - 1; } return coeff; } /*---------------------------------------------------------------------------- * Orthogonal array covering * The idea is to modify lists of mutually orthogonal latin squares (mols) * to produce q^3 blocks, q a prome power. * I thought this up but probably someone has done this the right way already. * If you know something about this please let me know!. * The code to generate the mols is from Art Owen's netlib oa package. * The method seems to go back to Bose and Bush, ca 1955. */ /* * seed the initial block list with a partial covering generated * form an orthogonal array. An oa of strength 3 has no repeated triples in * matching columns when rows are compared to each other (+ columns <=> rows) * The oa used here is built using a Bush construction based on finite fields. * If I add a constant to each column I can make symbols nunique to each * column; the resulting rows will share no common triple, anywhere. * I can then use the rows as blocks. */ void orthogonalArrayCovering(Cov* pC, int bmax, int q) { int i, j, k; int strength = 3; int block[2*B_MAX]; int poly[P_MAX]; Opt opt = pC->opt; if (bmax > q+1) { bmax = q+1; } initCov(pC, 0, q*q, bmax, 3); /* * the block generated can be as long as q+1 numbers; * for q = 3, only 4 numbers of the block are determined, leaving 3 to be * determined by embedding. */ if (!initFiniteField(q)) { return; } /* * design of size q*q*q over q symbols, strength 3, and pC->B columns */ for (i = 0; i < q*q*q; i++) { int_to_poly(q, strength, i, poly); block[0] = poly[2]; for (j = 1; j < bmax; j++) { k = poly[2]; k = mul[k][j-1]; k = add[k][ poly[1] ]; k = mul[k][j-1]; k = add[k][ poly[0] ]; block[j] = j*q + k; } addBlock(pC, block); } /* * and now tuck the q^3 block into the unit tetrahedron... */ if (opt.oa_tuck) { Cov tuckC; int i, j, k, m; int map[2*B_MAX]; Opt opt; initOpt(&opt); genCov(&tuckC, &opt, 2*q, bmax, 3); /* pairs of 0...q */ for (i = 0; i < bmax-1; ++i) { for (j = i+1; j < bmax; ++j) { for (k = 0; k < q; ++k) { map[k] = q*i + k; } for (k = 0; k < q; ++k) { map[q+k] = q*j + k; } for (k = 0; k < tuckC.nblocks; ++k) { for (m = 0; m < bmax; ++m) { if (tuckC.blocks[k][m] >= tuckC.N) { /* jokers */ break; } block[m] = map[tuckC.blocks[k][m]]; } for (; m < bmax; ++m) { block[m] = pC->N; } addBlock(pC, block); } } } } } /*---------------------------------------------------------------------------- * Block Stuff */ void startBlock(int* block, int nB) /* the returned block must be passed thru nextBlock() to be true beginning */ { int i; for (i = 0; i < nB; ++i) { block[i] = i; } block[nB-1]--; } /* * return the next tuple of points in sorted lexicographic order. * also calc the index if requested; * return 0 to indicate end is reached. */ int nextBlock(int* block, int nB, int nN, int* pIndex) { int i; for (i = nB - 1; i >= 0; --i) { block[i]++; if (block[i] < nN - (nB - 1 - i)) { for (; i < nB-1; ++i) { block[i+1] = block[i] + 1; } if (pIndex) { int j, n = 0; for (j = 0; j < nB; ++j) { n *= nN; n += block[j]; } *pIndex = n; } return 1; } } return 0; } void sortBlock(int* block, int nB) { int i, j; for (i = 0; i < nB - 1; ++i) { for (j = i + 1; j < nB; ++j) { if (block[j] < block[i]) { int b = block[i]; block[i] = block[j]; block[j] = b; } } } } void markFlatTuples(tup_t* usedTuples, int* flat, int nF, int nB, int nN) { int i; int findx[F_MAX]; startBlock(findx, nB); while(nextBlock(findx, nB, nF, 0)) { int uindx = 0; for (i = 0; i < nB; ++i) { uindx *= nN; uindx += flat[findx[i]]; } usedTuples[uindx]++; } } /*---------------------------------------------------------------------------- * Timer & Misc Stuff */ void initTimer() { struct tms tmsbuf; times(&tmsbuf); g_user_sys_time = tmsbuf.tms_utime + tmsbuf.tms_stime; } int checkTime() { struct tms tmsbuf; times(&tmsbuf); return (tmsbuf.tms_utime + tmsbuf.tms_stime - g_user_sys_time)/CLK_TCK; } /*---------------------------------------------------------------------------- * Context Stuff */ /* * increment context for recursion */ void incrContext(Con* pCn, Cov* pCv, int begin, int shared_tuples) { int j; pCn->level++; pCn->begin = begin + 1; pCn->end++; pCn->shared_tuples = shared_tuples; switch(pCv->M) { case 3: for (j = 0; j < pCn->level-1; ++j) { pCn->offset[pCn->noffset++] = (pCn->block[j]*pCv->N + begin)*pCv->N; } break; case 2: pCn->offset[pCn->noffset++] = begin*pCv->N; break; } } void printBlock(int *block, int nB, int nN, int offset) { int i, k; int used[3*N_MAX]; memset(used, 0, sizeof(used)); for (i = 0; i < nB; ++i) { if (block[i] < nN) { k = block[i] + offset; } else { for (k = 0; used[k]; ++k) { ; } } printf("%s%d", (i? " " : ""), k + 1); used[k] = 1; } printf("\n"); } void dumpBlocks(Cov* pC, int offset) { int i; for (i = 0; i < pC->nblocks; ++i) { printBlock(pC->blocks[i], pC->B, pC->N, offset); } }