/* * Released under GNU General Public License December 2007 * http://www.gnu.org/copyleft/gpl.html */ // // ppsearch.c - primitive polynomial search and primitivity test program // // This program tests a binary polynomial for primitivity by taking a // primitive element and raising it to the power 2^n, modulo the // polynomial. If the original primitive element is not returned, // the polynomial fails the primitivity test. An alternate test is to // raise the candidate primitive to 2^n-1. Any primitive element raised // to 2^n-1 returns One. This test uses the 2^n form of the test because // it is slightly faster. If 2^n-1 is not prime, additional tests are required // because the sequence of polynomials could have repeated 2^n-1 / (prime) // times. For each prime factor, the primitive element is raised to the // power 2^n-1 / (prime). If one is returned for any of these additional tests, // the polynomial is not primitive. // // The polynomial exponentiation algorithm is standard square and multiply. // A small optimization is made by choosing the primitive element as a // "right shift" of the primitive polynomial. This makes the multiply part // of multiply and square relatively insignificant. // // Note: Portability has been improved. The code attempts to follow the C99 // language standard, except for the use of macros to enable MMX and XMM // code generation. Recent Microsoft, Intel, and gnu compilers support // these macros. // // email sduplichan at this domain // // Modified 2003, Jean-Luc Cooke // jlcooke -a- certainkey -point- com // // 1) Now includes LaTeX output mode. // 2) Now compiles in LINUX #include #include #include #include #include #include #include #if defined __GNUC__ #include #include #undef min #define min(X,Y) ((X) < (Y) ? (X) : (Y)) // Some gcc versions do not align XMM data properly, causing gp fault. This is a work-around #define GCC_ALIGNMENT_HELP __attribute__((aligned(16))) #elif defined (WIN32) || defined (WIN64) #include #include #define GCC_ALIGNMENT_HELP typedef unsigned __int8 uint8_t; typedef __int32 int32_t; typedef unsigned __int64 uint64_t; typedef unsigned __int32 uint32_t; #else #error === unknown compiler, add support === #endif // Set MAXBITS to the size of the largest extended integer. This needs to be at // least one bigger than the highest degree primitive polynomial search or test. // There is now little speed dependency on this setting. #ifndef MAXBITS #define MAXBITS 6400 // set to a multiple of BIG_INT_BITS #endif #define DIMENSION(array) (sizeof (array) / sizeof (array [0])) //---------------------------------------------------------------------------- // Integer type selection // // uintn_t: Several functions operate using UINTN integers. These // functions are not among the most performance sensitive. // // BIG_INT_BITS: This size should match the biggest integer type // used (MMX=64, XMM=128). The extended integer math functions // operate on chunks of this size. // #define BIG_INT_BITS 128 #if !defined (OPTIMIZE32) typedef uint64_t uintn_t; #define UINTN_BITS 64 #else // default to 64-bit optimization typedef uint32_t uintn_t; #define UINTN_BITS 32 #endif //---------------------------------------------------------------------------- #if MAXBITS & (BIG_INT_BITS - 1) # error MAXBITS must be a multiple of BIG_INT_BITS (128) #endif #define BIGINT_COUNT (MAXBITS / BIG_INT_BITS) #define UINT64_COUNT (MAXBITS / 64) #define UINT128_COUNT (MAXBITS / 128) #define UINT32_COUNT (MAXBITS / 32) #define UINT8_COUNT (MAXBITS / 8) #define UINTN_COUNT (MAXBITS / UINTN_BITS) //---------------------------------------------------------------------------- // large integer structure // typedef union { uintn_t uintn [UINTN_COUNT]; uint8_t uint8 [UINT8_COUNT]; uint32_t uint32 [UINT32_COUNT]; uint64_t uint64 [UINT64_COUNT]; __m64 m64 [UINT64_COUNT]; // MMX data, x86 only __m128i m128i [UINT128_COUNT]; // XMM data, x86 only } INTEGER GCC_ALIGNMENT_HELP; //---------------------------------------------------------------------------- // structure to hold list of numbers: 2^n-1 / (a single prime) // typedef struct { INTEGER *divisor; int count; } DIVISOR_LIST; // structure for accessing small list of irreducable polynomials typedef struct { uint32_t list [4000]; int count; } IRREDUCABLEINFO; //---------------------------------------------------------------------------- // extended integer constants static INTEGER IntegerZero; static INTEGER IntegerOne; static INTEGER IntegerTwo; //---------------------------------------------------------------------------- // // functions for sequentially returning the binomial coefficients (x choose m) // quantity of different ways to choose m of x items. Caller calls firstCombination // with m value and an array of size m. For all the additional combinations, // nextCombination is called. Array is returned with elements selected from [0, m-1]. // If an array element other than array [m-1] is changed, the global multipleBitUpdate // is incremented to flag the event. when nextCombination returns < 0, all x choose m // combinations have been generated. // static int firstCombination (int m, int *array) { int index; for (index = 0; index < m; index++) array [index] = index; return m - 1; } //---------------------------------------------------------------------------- static int nextCombination (int x, int m, int walkingIndex, int *array) { int index; if (++array [walkingIndex] > x - m + walkingIndex) { for (;;) { if (--walkingIndex < 0) return walkingIndex; if (++array [walkingIndex] <= x - m + walkingIndex) break; } for (index = walkingIndex + 1; index < m; index++) array [index] = array [walkingIndex] + index - walkingIndex; walkingIndex = m - 1; } return walkingIndex; } //---------------------------------------------------------------------------- // // copyInteger - copy the active portion of an extended integer // static void copyInteger (INTEGER *dest, INTEGER *source, int activeBits) { int index; for (index = 0; index < activeBits / 64; index++) dest->uint64 [index] = source->uint64 [index]; } //---------------------------------------------------------------------------- // // compareInteger - compare the active portions of two extended integers // static int compareInteger (INTEGER *dest, INTEGER *source, int activeBits) { int index, diff = 0; for (index = 0; index < activeBits / 64; index++) diff += dest->uint64 [index] != source->uint64 [index]; return diff; } //---------------------------------------------------------------------------- // // populationCount - return number of non-zero bits // static int populationCount (INTEGER *arg, int activeBits) { int index, bit, weight = 0; for (index = 0; index < activeBits / UINTN_BITS; index++) for (bit = 0; bit < UINTN_BITS; bit++) weight += ((arg->uintn [index] >> bit) & 1); return weight; } //---------------------------------------------------------------------------- // // logError - printf message to stdout and exit // static void logError (char *message,...) { va_list Marker; char buffer [100]; va_start (Marker, message); vsprintf(buffer, message, Marker); va_end(Marker); fprintf (stderr, "\n%s\n", buffer); exit (1); } //---------------------------------------------------------------------------- // // setbit - set a bit in an extended integer // static void setbit (INTEGER *data, int bitnumber) { data->uintn [bitnumber / UINTN_BITS] |= (uintn_t) 1 << (bitnumber % UINTN_BITS); } //---------------------------------------------------------------------------- // // roundUp - round an integer value up to a multiple of n // static int roundUp (int value, int n) { value += n - 1; value -= value % n; return value; } //---------------------------------------------------------------------------- // // highestSetBit_uint32 - finds highest set bit in 32-bit integer // static int highestSetBit_uint32 (uint32_t data) { data |= (data >> 1); data |= (data >> 2); data |= (data >> 4); data |= (data >> 8); data |= (data >> 16); data -= (data >> 1) & 0x55555555; data = (data & 0x33333333) + ((data >> 2) & 0x33333333); data = (data + (data >> 4)) & 0x0f0f0f0f; return ((data * 0x01010101) >> 24) - 1; } //---------------------------------------------------------------------------- // // highestSetBit_uint64 - finds highest set bit in 64-bit integer // static int highestSetBit_uint64 (uint64_t data) { data |= (data >> 1); data |= (data >> 2); data |= (data >> 4); data |= (data >> 8); data |= (data >> 16); data |= (data >> 32); data -= (data >> 1) & 0x5555555555555555ULL; data = (data & 0x3333333333333333ULL) + ((data >> 2) & 0x3333333333333333ULL); data = (data + (data >> 4)) & 0x0f0f0f0f0f0f0f0fULL; return ((data * 0x0101010101010101ULL) >> 56) - 1; } //---------------------------------------------------------------------------- // // highestSetBit32 - finds highest set bit in extended integer, 32-bit optimized // static int highestSetBit32 (INTEGER *data, int activeBits) { int bitno = activeBits - 32; int index = activeBits / 32 - 1; while (index > 0) { if (data->uint32 [index]) break; index--; bitno -= 32; } return highestSetBit_uint32 (data->uint32 [index]) + bitno; } //---------------------------------------------------------------------------- // // highestSetBit64 - finds highest set bit in extended integer, 64-bit optimized // static int highestSetBit64 (INTEGER *data, int activeBits) { int bitno = activeBits - 64; int index = activeBits / 64 - 1; while (index > 0) { if (data->uint64 [index]) break; index--; bitno -= 64; } return highestSetBit_uint64 (data->uint64 [index]) + bitno; } //---------------------------------------------------------------------------- // // highestSetBit - finds highest set bit in extended integer. x86 optimized // static int highestSetBit (INTEGER *data, int activeBits) { #if defined (OPTIMIZE32) return highestSetBit32 (data, activeBits); #else return highestSetBit64 (data, activeBits); #endif } //---------------------------------------------------------------------------- // // extractBit - return the value of a selected bit from a extended integer // static int extractBit (INTEGER *data, int bitNumber) { return (data->uintn [bitNumber / UINTN_BITS] >> (bitNumber % UINTN_BITS)) & 1; } //---------------------------------------------------------------------------- // // shiftLeft - shift left extended integer // static void shiftLeft (INTEGER *source, INTEGER *dest, int shiftCount, int activeBits) { int sourceIndex, destIndex, leftShift, rightShift, uintnCount; uintnCount = activeBits / UINTN_BITS; sourceIndex = uintnCount - 1 - shiftCount / UINTN_BITS; destIndex = uintnCount - 1; leftShift = shiftCount % UINTN_BITS; rightShift = UINTN_BITS - leftShift; if (shiftCount == 0) logError ("error, shift count is zero\n"); if (!leftShift) // special case: just move integers, no shifting required while (sourceIndex > 0) { dest->uintn [destIndex] = source->uintn [sourceIndex]; sourceIndex--; destIndex--; } else while (sourceIndex > 0) { dest->uintn [destIndex] = (source->uintn [sourceIndex - 0] << leftShift) | (source->uintn [sourceIndex - 1] >> rightShift); sourceIndex--; destIndex--; } dest->uintn [destIndex] = source->uintn [0] << leftShift; while (destIndex > 0) dest->uintn [--destIndex] = 0; } //---------------------------------------------------------------------------- // // shiftLeftOnce - shift left extended integer once using SSE registers // static void shiftLeftOnceSse (INTEGER *source, INTEGER *dest, int activeBits) { int sourceIndex, uint128Count; __m128i current, next, shiftedBits, recoveredBits, carryOut; uint128Count = activeBits / 128; sourceIndex = uint128Count - 1; next = source->m128i [sourceIndex]; while (sourceIndex) { current = source->m128i [sourceIndex]; next = source->m128i [sourceIndex - 1]; carryOut = _mm_srli_epi64 (next, 63); // bit 64 is bit 127 from next part carryOut = _mm_shuffle_epi32 (carryOut, 0xFE); // bit 0 is bit 127 from next part, other bits are zero recoveredBits = _mm_srli_epi64 (current, 63); // bit 0 is bit 63 of current part recoveredBits = _mm_shuffle_epi32 (recoveredBits, 0xCF); // bit 64 is bit 63 of current part, other bits are zero recoveredBits = _mm_or_si128 (recoveredBits, carryOut); // bits 0 and 64 are bits lost by the SSE shift shiftedBits = _mm_slli_epi64 (current, 1); dest->m128i [sourceIndex] = _mm_or_si128 (shiftedBits, recoveredBits); current = next; sourceIndex--; } // the final shift zero fills dest->m128i [0] = _mm_or_si128 (_mm_slli_epi64 (next, 1), _mm_shuffle_epi32 (_mm_srli_epi64 (next, 63), 0xCF)); } //---------------------------------------------------------------------------- // // shiftLeftOnce - shift left extended integer once using MMX registers // static void shiftLeftOnceMmx (INTEGER *source, INTEGER *dest, int activeBits) { int sourceIndex, uint64Count; __m64 current, next; uint64Count = activeBits / 64; sourceIndex = uint64Count - 1; while (sourceIndex) { current = source->m64 [sourceIndex]; next = source->m64 [sourceIndex - 1]; dest->m64 [sourceIndex] = _mm_or_si64 (_mm_slli_si64 (current, 1), _mm_srli_si64 (next, 63)); current = next; sourceIndex--; } // the final shift zero fills dest->m64 [0] = _mm_slli_si64 (next, 1); } //---------------------------------------------------------------------------- // // shiftLeftOnce - shift left extended integer once using 64-bit registers // static void shiftLeftOnce64 (INTEGER *source, INTEGER *dest, int activeBits) { int sourceIndex, uint64Count; uint64_t current, next; uint64Count = activeBits / 64; sourceIndex = uint64Count - 1; while (sourceIndex) { current = source->uint64 [sourceIndex]; next = source->uint64 [sourceIndex - 1]; dest->uint64 [sourceIndex] = (current << 1) | (next >> 63); current = next; sourceIndex--; } // the final shift zero fills dest->uint64 [0] = next << 1; } //---------------------------------------------------------------------------- // // shiftLeftOnce - shift left extended integer once using 32-bit registers // static void shiftLeftOnce32 (INTEGER *source, INTEGER *dest, int activeBits) { int sourceIndex, uint32Count; uint32_t current, next; uint32Count = activeBits / 32; sourceIndex = uint32Count - 1; while (sourceIndex) { current = source->uint32 [sourceIndex]; next = source->uint32 [sourceIndex - 1]; dest->uint32 [sourceIndex] = (current << 1) | (next >> 31); current = next; sourceIndex--; } // the final shift zero fills dest->uint32 [0] = next << 1; } //---------------------------------------------------------------------------- // // shiftLeftOnce - shift left extended integer once. Optimized for single shift // static void shiftLeftOnce (INTEGER *source, INTEGER *dest, int activeBits) { #if defined (OPTIMIZE32) shiftLeftOnce32 (source, dest, activeBits); #else shiftLeftOnceSse (source, dest, activeBits); #endif } //---------------------------------------------------------------------------- // // addInteger - add extended integers // static void addInteger (INTEGER *value1, INTEGER *value2, INTEGER *result, int activeBits) { int index, uint32Count, carry = 0; INTEGER total; uint32Count = activeBits / 32; for (index = 0; index < uint32Count; index++) { total.uint64 [0] = (uint64_t) value1->uint32 [index] + value2->uint32 [index] + carry; result->uint32 [index] = total.uint32 [0]; carry = total.uint32 [1]; } } //---------------------------------------------------------------------------- // // multiplyInteger - multiply extended integers // static void multiplyInteger (INTEGER *value1, INTEGER *value2, INTEGER *result, int activeBits) { int index; INTEGER temp, total; copyInteger (&total, &IntegerZero, activeBits); if (extractBit (value1, 0)) addInteger (&total, value2, &total, activeBits); for (index = 1; index < activeBits; index++) { if (!extractBit (value1, index)) continue; shiftLeft (value2, &temp, index, activeBits); addInteger (&total, &temp, &total, activeBits); } copyInteger (result, &total, activeBits); } //---------------------------------------------------------------------------- // // mul10 - multiply extended integer by 10 // static void mul10 (INTEGER *value, int activeBits) { INTEGER x8, x2; shiftLeft (value, &x8, 3, activeBits); shiftLeft (value, &x2, 1, activeBits); addInteger (&x8, &x2, value, activeBits); } //---------------------------------------------------------------------------- // // allOnes - returns extended integer with ls bits set to all ones // static INTEGER allOnes (int bits) { static INTEGER result; int index; result = IntegerZero; for (index = 0; index < bits; index++) setbit (&result, index); return result; } //---------------------------------------------------------------------------- // // skipWhiteSpace - skip spaces and tabs in a character buffer // static char *skipWhiteSpace (char *position) { while (*position == ' ' || *position == '\t') position++; return position; } //---------------------------------------------------------------------------- // // findBase - looks at ascii buffer and determines the base. A leadinf 0x forces hex. // A trailing b, o, d forces binary, octal, or decimal, respectively. // If no prefix or suffix is present, a best guess is made by scanning the digits. // static int findBase (char *position) { int base = 0, maxCharacter = '0'; char *endOfNumber; if (position [0] == '0' && tolower (position [1]) == 'x') return 16; endOfNumber = position + strspn (position, "0123456789"); if (tolower (*endOfNumber) == 'b') base = 2; if (tolower (*endOfNumber) == 'o') base = 8; if (tolower (*endOfNumber) == 'd') base = 10; // no suffix, look at digits while (position != endOfNumber) { if (maxCharacter < *position) maxCharacter = *position; position++; } if (base) { if (maxCharacter >= '0' + base) logError ("digit %c conflicts with base suffix %c", maxCharacter, *endOfNumber); } else { if (maxCharacter <= '1') base = 2; else if (maxCharacter <= '7') base = 8; else if (maxCharacter <= '9') base = 10; } return base; } //---------------------------------------------------------------------------- // // scanDecimalDigits - read a extended integer from an ascii buffer of decimal digits. // static void scanDecimalDigits (char *buffer, INTEGER *value) { char *position = buffer; INTEGER dvalue = IntegerZero; copyInteger (value, &IntegerZero, MAXBITS); for (;;) { int digit = *position++; if (!isdigit (digit)) break; mul10 (value, MAXBITS); dvalue.uintn [0] = digit - '0'; addInteger (value, &dvalue, value, MAXBITS); } } //---------------------------------------------------------------------------- // // scanDigits - read a extended integer from an ascii buffer on digits of a selectable base. // static void scanDigits (char *buffer, INTEGER *integer, int base) { static int digitBits [] = {0, 0, 1, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 4}; char *position, *endOfNumber, *startOfNumber = skipWhiteSpace (buffer); int bitno, index, bitsPerDigit; if (base == 0) base = findBase (startOfNumber); if (base == 10) { scanDecimalDigits (startOfNumber, integer); return; } if (startOfNumber [0] == '0' && tolower (startOfNumber [1]) == 'x') startOfNumber += 2; if (base == 16) endOfNumber = startOfNumber + strspn (startOfNumber, "0123456789abcdefABCDEF"); else if (base == 8) endOfNumber = startOfNumber + strspn (startOfNumber, "01234567"); else endOfNumber = startOfNumber + strspn (startOfNumber, "01"); bitsPerDigit = digitBits [base]; position = endOfNumber - 1; bitno = 0; *integer = IntegerZero; while (position >= startOfNumber) { int value; value = toupper (*position); if (value >= 'A') value = 10 + (value - 'A'); else value -= '0'; if (value >= base) logError ("invalid base %u digit (%c)", base, *position); for (index = 0; index < bitsPerDigit; index++) { if (value & 1) setbit (integer, bitno); bitno++; value >>= 1; } position--; } } //---------------------------------------------------------------------------- // // computeFactors - find factors using trial division method // static INTEGER *computeFactors (int *numberOfFactors, int order) { uint64_t n, p, end, f; int factorCount = 0; INTEGER *factorList = NULL; p = allOnes (order).uint64 [0]; n = p; end = sqrt (n); for (f = 3; f <= end; ) { if (n % f == 0) { n /= f; end = sqrt(n); factorList = realloc (factorList, sizeof (INTEGER) * (factorCount + 1)); factorList [factorCount] = IntegerZero; factorList [factorCount].uint64 [0] = f; factorCount++; } else f += 2; } factorList = realloc (factorList, sizeof (INTEGER) * (factorCount + 1)); factorList [factorCount] = IntegerZero; factorList [factorCount].uint64 [0] = n; factorCount++; *numberOfFactors = factorCount; return factorList; } //----------------------------------------------------------------------------- // // Find all the factors of 2^polynomialDegree-1 by reading them from a disk file // Return a list of 2^polynomialDegree-1 / (each single prime) // static int findFactors (int polynomialDegree, DIVISOR_LIST *divisorList, int verbose) { FILE *fp; INTEGER factor, *factorList = NULL, result, total = IntegerOne, org = allOnes (polynomialDegree); int numberOfFactors = 0, index1, index2, activeBits; char filename [100]; char *buffer = calloc (1, 10000); activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); divisorList->count = 0; divisorList->divisor = 0; sprintf (filename, "factor2n-1/%u.txt", polynomialDegree); fp = fopen (filename, "r"); if (!fp) { if (verbose) printf ("\nfailed to open %s\n", filename); if (polynomialDegree > 64) { static int warnedAlready; if (!warnedAlready) { printf ("\nfactor file missing for 2^%u - 1", polynomialDegree); if (!verbose) { warnedAlready++; printf (", further warnings suppressed"); } } return 0; } factorList = computeFactors (&numberOfFactors, polynomialDegree); total = org; // suppress sanity check for computed factorization } else { while (!feof (fp)) { if (!fgets (buffer, 10000, fp)) break; scanDecimalDigits (buffer, &factor); multiplyInteger (&total, &factor, &result, activeBits); total = result; factorList = realloc (factorList, sizeof (INTEGER) * (numberOfFactors + 1)); factorList [numberOfFactors++] = factor; } fclose (fp); } free (buffer); // do a sanity check on the factorization if (compareInteger (&org, &total, activeBits) != 0) { printf ("\n======incorrect factor file for 2^%u - 1 removed======", polynomialDegree); free (factorList); unlink (filename); return 0; } for (index1 = 0; index1 < numberOfFactors; index1++) { INTEGER result, total = IntegerOne; for (index2 = 0; index2 < numberOfFactors; index2++) { if (index2 == index1) continue; multiplyInteger (&total, &factorList [index2], &result, activeBits); total = result; } if (divisorList->count) if (memcmp (&total, &divisorList->divisor [divisorList->count - 1], sizeof (INTEGER)) == 0) continue; divisorList->divisor = realloc (divisorList->divisor, sizeof (INTEGER) * (divisorList->count + 1)); divisorList->divisor [divisorList->count++] = total; } free (factorList); return 1; } //---------------------------------------------------------------------------- // // xorInteger - exclusive-or two extended integers // uses SIMD 128-bit xor, unrolled and inline // static INTEGER *xorInteger (INTEGER *arg1, INTEGER *arg2, INTEGER *result, int activeBits) { int index; for (index = 0; index < activeBits / 128; index++) result->m128i [index] = _mm_xor_si128 (arg1->m128i [index], arg2->m128i [index]); return result; } //---------------------------------------------------------------------------- // // extract a bit field a bit from a extended integer // static int extractBits (INTEGER *data, int lsb, int msb) { int index, total = 0; for (index = lsb; index <= msb; index++) if (extractBit (data, index)) total |= 1 << (index - lsb); return total; } //---------------------------------------------------------------------------- // // binaryAscii - return in binary ascii representation of extended integer // in a static buffer. Four calls can be made before overwriting. // static char *binaryAscii (INTEGER *data, int bits) { static char buffer [4] [MAXBITS + 2]; static int cycle; char *position = buffer [cycle]; char *result = position; if (bits == 0) bits = MAXBITS; bits = highestSetBit (data, bits) + 1; if (bits == 0) bits++; while (bits) *position++ = (char) ('0' + extractBit (data, --bits)); *position++ = '\0'; if (++cycle == 4) cycle = 0; return result; } //---------------------------------------------------------------------------- // // hexAscii - return in hexadecimal ascii representation of extended integer // in a static buffer. Four calls can be made before overwriting. // static char *hexAscii (INTEGER *data, int bits) { static char buffer [4] [MAXBITS * 4 + 2]; static int cycle; int index; char *position = buffer [cycle]; char *result = position; if (bits == 0) bits = MAXBITS; bits = highestSetBit (data, bits) + 1; if (bits == 0) bits++; index = (bits + 7) / 8; while (index--) position += sprintf (position, "%02X", data->uint8 [index]); *position++ = '\0'; if ((((bits - 1) / 4) & 1) == 0) result++; if (++cycle == 4) cycle = 0; return result; } //---------------------------------------------------------------------------- // // octalAscii - return in octal ascii representation of extended integer // in a static buffer. Four calls can be made before overwriting. // static char *octalAscii (INTEGER *data, int activeBits) { static char buffer [4] [MAXBITS / 3 + 2]; static int cycle; char *position = buffer [cycle]; char *result = position; int bits; if (activeBits == 0) activeBits = MAXBITS; bits = highestSetBit (data, activeBits) + 1; if (bits == 0) bits++; bits = (bits + 2) / 3 * 3; while (bits) { *position = (char) ('0' + extractBits (data, bits - 3, bits - 1)); position++; bits -= 3; } *position++ = '\0'; if (++cycle == 4) cycle = 0; return result; } //---------------------------------------------------------------------------- // // polynomialText - return text string containing ascii description of // the polynomial represented by the extended integer // static char *polynomialText (void *data, int displayMode, int activeBits) { int index; INTEGER *polynomial = data; static char buffer [4] [MAXBITS * 12]; static int cycle; char *position = buffer [cycle]; char *result = position; if (displayMode == 'b') return binaryAscii (polynomial, activeBits); else if (displayMode == 'o') return octalAscii (polynomial, activeBits); else if (displayMode == 'h') return hexAscii (polynomial, activeBits); else if (displayMode == 'l') { position += sprintf (position, "$ "); for (index = activeBits - 1; index > 1; index--) if (extractBit (polynomial, index)) position += sprintf (position, "x^{%u} + ", index); if (extractBit (polynomial, 1)) position += sprintf (position, "x + "); position += sprintf (position, "1 $"); if (++cycle == 4) cycle = 0; return result; } else { for (index = activeBits - 1; index > 1; index--) if (extractBit (polynomial, index)) position += sprintf (position, "x^%u + ", index); if (extractBit (polynomial, 1)) position += sprintf (position, "x + "); position += sprintf (position, "1"); if (++cycle == 4) cycle = 0; return result; } } //---------------------------------------------------------------------------- // // shiftRight - shift right extended integer // static void shiftRight (INTEGER *integer, int shiftCount, int activeBits) { int destIndex, sourceIndex, rightShift, leftShift, uintnCount; uintnCount = activeBits / UINTN_BITS; destIndex = 0; sourceIndex = shiftCount / UINTN_BITS; rightShift = shiftCount % UINTN_BITS; leftShift = UINTN_BITS - rightShift; if (!rightShift) // special case: just move integers, no shifting required while (sourceIndex < uintnCount - 1) { integer->uintn [destIndex] = integer->uintn [sourceIndex] >> rightShift; sourceIndex++; destIndex++; } while (sourceIndex < uintnCount - 1) { integer->uintn [destIndex] = (integer->uintn [sourceIndex + 0] >> rightShift) | (integer->uintn [sourceIndex + 1] << leftShift); sourceIndex++; destIndex++; } integer->uintn [destIndex] = integer->uintn [uintnCount - 1] >> rightShift; while (destIndex < uintnCount - 1) integer->uintn [++destIndex] = 0; } //---------------------------------------------------------------------------- // // dividePolynomial - divide binary polynomial, coefficients are kept in extended integers // input: numerator - polynomial to divide // denominator - divisor polynomial // output: numerator - remainder polynomial // denominator - destroyed // quotient - result // static void dividePolynomial (INTEGER *numerator, INTEGER *denominator, INTEGER *quotient, int activeBits) { int numeratorPower = highestSetBit (numerator, activeBits); int denominatorPower = highestSetBit (denominator, activeBits); int denominatorShift = numeratorPower - denominatorPower; int bitsRetired; if ((signed) denominatorShift < 0) return; if (denominatorPower == -1) logError ("polynomial division by zero\n"); else if (denominatorPower == 0) { copyInteger (quotient, numerator, activeBits); copyInteger (numerator, &IntegerZero, activeBits); return; } copyInteger (quotient, &IntegerZero, activeBits); shiftLeft (denominator, denominator, denominatorShift, activeBits); for (;;) { setbit (quotient, denominatorShift); xorInteger (numerator, denominator, numerator, activeBits); bitsRetired = numeratorPower - highestSetBit (numerator, activeBits); denominatorShift -= bitsRetired; numeratorPower -= bitsRetired; if (numeratorPower < denominatorPower) break; shiftRight (denominator, bitsRetired, activeBits); } } //---------------------------------------------------------------------------- // // divideBigPolynomialBySmall - divide binary polynomial, numerator coefficient is kept in extended integer // // input: numerator - polynomial to divide // denominator - divisor polynomial, 33 bits maximum // return: - remainder polynomial // static uint32_t divideBigPolynomialBySmall (INTEGER *numerator, uint64_t denominator, int activeBits) { int numeratorPower, denominatorPower, uint32count; uint32_t remainder, mask; numeratorPower = highestSetBit (numerator, activeBits); denominatorPower = highestSetBit_uint64 (denominator); mask = denominator << (32 - denominatorPower); uint32count = numeratorPower / 32 + 1; remainder = 0; if (denominatorPower > 33) logError ("divideBigPolynomialBySmall overflow"); for (;;) { int bit, bitCount; remainder ^= numerator->uint32 [--uint32count]; bitCount = 32; if (uint32count == 0) bitCount -= denominatorPower; for (bit = 0; bit < bitCount; bit++) { int out; out = remainder >> 31; remainder <<= 1; remainder ^= ((-out) & mask); } if (uint32count == 0) break; } return remainder >> (32 - denominatorPower); } //---------------------------------------------------------------------------- // // Build a list of the first few irreducable polynomials. These are used to // quickly screen out reducable polynomials before running a more lengthy test // static void findIrreducablePolynomials (IRREDUCABLEINFO *irreducableInfo, int displayMode, int verbose) { INTEGER candidate; INTEGER numerator, quotient, denominator; int index, numeratorPower, denominatorPower, activeBits, fail; activeBits = BIG_INT_BITS; copyInteger (&candidate, &IntegerZero, activeBits); // Start with x^2 + x + 1, because candidate polynomials // are selected to avoid multiples of x + 1. candidate.uint32 [0] = 7; for (;;) { numeratorPower = highestSetBit (&candidate, activeBits); fail = 0; for (index = 0; index < irreducableInfo->count; index++) { numerator = candidate; denominator = IntegerZero; denominator.uint32 [0] = irreducableInfo->list [index]; denominatorPower = highestSetBit (&denominator, activeBits); if (denominatorPower * 2 > numeratorPower) break; dividePolynomial (&numerator, &denominator, "ient, activeBits); if (numerator.uint64 [0] == 0) fail++; if (fail) break; } if (!fail) { irreducableInfo->list [irreducableInfo->count] = candidate.uint32 [0]; irreducableInfo->count++; if (irreducableInfo->count == DIMENSION (irreducableInfo->list)) break; } candidate.uint32 [0]++; } if (verbose) { uint64_t small, big; small = irreducableInfo->list [0]; big = irreducableInfo->list [irreducableInfo->count - 1]; printf ("findIrreducablePolynomials: %s through %s\n", polynomialText (&small, displayMode, 64), polynomialText (&big, displayMode, 64)); } } //---------------------------------------------------------------------------- int findSmallPolynomialFactor (INTEGER *value, IRREDUCABLEINFO *irreducableInfo, int polynomialDegree) { uint32_t *irreducables, remainder, denominator; int index, denominatorPower, activeBits, limit; if (polynomialDegree < 96) return 0; irreducables = irreducableInfo->list; limit = highestSetBit_uint32 (polynomialDegree) + 1; activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); for (index = 0; index < irreducableInfo->count; index++) { denominator = irreducables [index]; denominatorPower = highestSetBit_uint32 (denominator); if (denominatorPower > limit) break; remainder = divideBigPolynomialBySmall (value, denominator, activeBits); if (remainder != 0) continue; return denominator; } return 0; } //---------------------------------------------------------------------------- // // when searching for primitive polynomials, this function returns the next candidate // static void nextPolynomial (INTEGER *polynomial, int polynomialWeight, int polynomialDegree) { int weight, activeBits; for (;;) { if (polynomialWeight) { int index; static int walkingIndex; static int inProgress, rowPointer [MAXBITS]; static int previousPolynomialDegree; // force reinitialization if called with different degree before previous completed if (previousPolynomialDegree != polynomialDegree) { previousPolynomialDegree = polynomialDegree; inProgress = 0; } if (!inProgress) { inProgress++; walkingIndex = firstCombination (polynomialWeight - 2, rowPointer); } else { walkingIndex = nextCombination (polynomialDegree - 1, polynomialWeight - 2, walkingIndex, rowPointer); if (walkingIndex < 0) { *polynomial = IntegerZero; inProgress = 0; return; } } *polynomial = IntegerOne; setbit (polynomial, polynomialDegree); for (index = 0; index < polynomialWeight - 2; index++) setbit (polynomial, rowPointer [index] + 1); return; } activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); addInteger (polynomial, &IntegerTwo, polynomial, activeBits); // We know that a binary primitive polynomial has an odd number of // non-zero coefficients (otherwise x+1 would divide it). So skip the evens. weight = populationCount (polynomial, activeBits); if (weight % 2 == 0) continue; break; } } //---------------------------------------------------------------------------- // // modularMultiplyPolynomial - modular multiply binary polynomial, coefficients // are kept in extended integers // static void modularMultiplyPolynomial (INTEGER *factor1, INTEGER *factor2, INTEGER *modulo, INTEGER *product, int polynomialDegree) { int index, activeBits, factor1Power; INTEGER result, factor2copy; activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); copyInteger (&factor2copy, factor2, activeBits); copyInteger (&result, &IntegerZero, activeBits); factor1Power = highestSetBit (factor1, activeBits); if (factor1Power < 0) logError ("unexpected multiply by zero\n"); activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); for (index = 0; index <= factor1Power; index++) { if (extractBit (factor1, index)) xorInteger (&result, &factor2copy, &result, activeBits); shiftLeftOnce (&factor2copy, &factor2copy, activeBits); if (extractBit (&factor2copy, polynomialDegree)) xorInteger (&factor2copy, modulo, &factor2copy, activeBits); } copyInteger (product, &result, activeBits); } //---------------------------------------------------------------------------- // // modularPowerPolynomial - modular exponentiation for binary polynomial, // coefficients are kept in extended integers // static void modularPowerPolynomial (INTEGER *primitiveElement, INTEGER *power, INTEGER *modulo, INTEGER *result, int polynomialDegree) { int bitno, totalBits, activeBits; activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); totalBits = highestSetBit (power, activeBits) - 1; copyInteger (result, primitiveElement, activeBits); for(bitno = totalBits; bitno >= 0; bitno--) { modularMultiplyPolynomial (result, result, modulo, result, polynomialDegree); if (extractBit (power, bitno)) { // a general algorithm must execute: // modularMultiplyPolynomial (result, primitiveElement, modulo, result, polynomialDegree); // But because the first argument is the right shift of the modulo, we emulate a shift register multiplier // this is not a major optimization, but it helps int lsb = extractBit (result, 0); shiftRight (result, 1, activeBits); if (lsb) xorInteger (result, primitiveElement, result, activeBits); } } } //---------------------------------------------------------------------------- // // primitivityTest - return non-zero if the polynomial is primitive // static int primitivityTest (INTEGER *polynomial, int polynomialDegree, DIVISOR_LIST *divisorList, int displayMode, int verbose) { // "quick" test to determine if the period might be 2^polynomialDegree-1 // only guaranteed for prime 2**polynomialDegree-1, but quickly screens many others INTEGER primitiveElement, power, result; int activeBits; activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); copyInteger (&power, &IntegerZero, activeBits); setbit (&power, polynomialDegree); copyInteger (&primitiveElement, polynomial, activeBits); shiftRight (&primitiveElement, 1, activeBits); modularPowerPolynomial (&primitiveElement, &power, polynomial, &result, polynomialDegree); if (compareInteger (&result, &primitiveElement, activeBits) == 0) { if (divisorList->count == 1) { if (!verbose) printf ("\n%s", polynomialText (polynomial, displayMode, activeBits)); else printf ("is primitive"); return 1; } else { int k; for (k = 0; k < divisorList->count; k++) { copyInteger (&power, &divisorList->divisor [k], activeBits); modularPowerPolynomial (&primitiveElement, &power, polynomial, &result, polynomialDegree); if (compareInteger (&result, &IntegerOne, activeBits) == 0) break; } if (k == divisorList->count) { if (!verbose) printf ("\n%s", polynomialText (polynomial, displayMode, activeBits)); else printf ("is primitive"); return 1; } else if (verbose) printf ("fails on order test %u of %u", k + 1, divisorList->count); } } else if (verbose) printf ("fails initial test"); return 0; } //---------------------------------------------------------------------------- // // findPrimitivePolynomials - search for primitive polymonials of the specified degree // static void findPrimitivePolynomials (INTEGER *polynomial, IRREDUCABLEINFO *irreducableInfo, int polynomialDegree, int targetCount, int polynomialWeight, int displayMode, int verbose) { DIVISOR_LIST divisorList; int factorizationAvailable = findFactors (polynomialDegree, &divisorList, verbose); int count = 0, activeBits; // if no factorization is available, nothing can be done if (!factorizationAvailable) return; activeBits = roundUp (polynomialDegree + 1, BIG_INT_BITS); if (verbose) printf ("\n%s", polynomialText (polynomial, displayMode, activeBits)); for (;;) { int isPrimitive; int hsb; uint32_t smallDivisor; INTEGER temp; copyInteger (&temp, polynomial, activeBits); smallDivisor = findSmallPolynomialFactor (&temp, irreducableInfo, polynomialDegree); if (smallDivisor) { signed int hsb; if (verbose) { uint64_t temp = smallDivisor; printf (" is reducable, divisor %s", polynomialText (&temp, displayMode, 64)); } nextPolynomial (polynomial, polynomialWeight, polynomialDegree); hsb = highestSetBit (polynomial, activeBits); if ((unsigned int) hsb == polynomialDegree + 1) break; if (hsb == -1) break; if (verbose) printf ("\n%s", polynomialText (polynomial, displayMode, activeBits)); continue; } else if (verbose) printf (" "); // display end of irreducability check isPrimitive = primitivityTest (polynomial, polynomialDegree, &divisorList, displayMode, verbose); count += isPrimitive; if (count >= targetCount) break; nextPolynomial (polynomial, polynomialWeight, polynomialDegree); hsb = highestSetBit (polynomial, activeBits); if (hsb == -1) break; if ((unsigned int) hsb == polynomialDegree + 1) break; if (verbose) printf ("\n%s", polynomialText (polynomial, displayMode, activeBits)); } free (divisorList.divisor); } //---------------------------------------------------------------------------- // // command line help goes here // static void helpScreen (void) { printf("use ppsearch options\n\n"); printf ("options: bits=d set polynomial size to d (decimal) bits\n"); printf (" poly=0xnnnn set initial polynomial (hex, ms bit not required)\n"); printf (" poly=nnnnb set initial polynomial (binary, ms bit not required)\n"); printf (" poly=nnnno set initial polynomial (octal, ms bit not required)\n"); printf (" poly=\"x^d+...\" set initial polynomial by non-zero coefficients & powers\n"); printf (" search=d search for next d primitive polynomials\n"); printf (" weight=d search for polynomials with d non-zero coefficients only\n"); printf (" binary results in binary\n"); printf (" octal results in octal\n"); printf (" hex results in hex\n"); printf (" loop when search completes, restart with bits+1\n"); printf (" maxbits exit loop mode when bits reaches maxbits\n"); printf (" verbose show more details\n"); exit (1); } //---------------------------------------------------------------------------- int main (int argc, char *argv []) { INTEGER polynomial = IntegerZero; int verbose = 0, loopMode = 0, polynomialWeight = 0; int displayMode = 0; int polynomialDegree = 0, targetCount = 1; int maxbits = MAXBITS; IRREDUCABLEINFO irreducableInfo = {0}; time_t startTime; memset(IntegerOne.uintn, 0, sizeof(IntegerOne)); memset(IntegerTwo.uintn, 0, sizeof(IntegerTwo)); IntegerOne.uintn[0] = 1; IntegerTwo.uintn[0] = 2; if (argc == 1) helpScreen (); while (--argc) { char *position = argv [argc]; if (strnicmp (position, "search=", 7) == 0) targetCount = atoi (position + 7); else if (strnicmp (position, "poly=x", 6) == 0) { int degree; polynomial = IntegerZero; position += 5; while (*position) { int power = 0; if (*position == '1') { power = 0; position++; } else if (tolower (*position) == 'x') { if (position [1] == '^') { position += 2; power = atoi (position); position += strspn (position, "0123456789"); } else { power = 1; position++; } } else logError ("polynomial syntax error (%s)\n", position); if (extractBit (&polynomial, power)) logError ("duplicate polynomial power\n"); setbit (&polynomial, power); position = skipWhiteSpace (position); if (*position == '+') { position = skipWhiteSpace (position + 1); continue; } if (*position == '\0') break; logError ("polynomial syntax error (%s)\n", position); } degree = highestSetBit (&polynomial, MAXBITS); if (!polynomialDegree) polynomialDegree = degree; else if (polynomialDegree != degree) logError ("bits= conflicts with poly=x^...\n"); } // this polynomial entry mode handles only polynomials else if (strnicmp (position, "poly=", 5) == 0) scanDigits (position + 5, &polynomial, 0); else if (stricmp (position, "verbose") == 0) verbose++; else if (stricmp (position, "loop") == 0) loopMode = 1; else if (strnicmp (position, "maxbits=", 8) == 0) maxbits = strtoul (position + 8, NULL, 10); else if (stricmp (position, "binary") == 0) displayMode = 'b'; else if (stricmp (position, "octal") == 0) displayMode = 'o'; else if (stricmp (position, "hex") == 0) displayMode = 'h'; else if (stricmp (position, "latex") == 0) displayMode = 'l'; else if (strnicmp (position, "weight=", 7) == 0) { polynomialWeight = atol (position + 7); if (polynomialWeight < 3) logError ("minumum weight is 3"); if (polynomialWeight % 2 != 1) logError ("weight must be odd, otherwise a factor of x+1 will be present"); } else if (strnicmp (position, "bits=", 5) == 0) polynomialDegree = atol (position + 5); else logError ("ERROR: invalid command line argument %s\n", position); } if (polynomialDegree >= MAXBITS) logError ("this build is limited to %u bits\n", MAXBITS - 1); if (polynomialDegree == 0) logError ("specify either bits= or poly=x^...\n"); findIrreducablePolynomials (&irreducableInfo, displayMode, verbose); startTime = time (NULL); // if no polynomial given, or just low bits, set ms bit setbit (&polynomial, 0); setbit (&polynomial, polynomialDegree); for (;;) { findPrimitivePolynomials (&polynomial, &irreducableInfo, polynomialDegree, targetCount, polynomialWeight, displayMode, verbose); if (!loopMode) break; polynomialDegree++; if (polynomialDegree == MAXBITS) break; if (polynomialDegree == maxbits) break; polynomial = IntegerZero; setbit (&polynomial, 0); setbit (&polynomial, 1); setbit (&polynomial, polynomialDegree); } printf("\n"); printf ("elapsed time: %lu\n", (long) (time (NULL) - startTime)); return 0; } //-------------------------end of file-----------------------------