nspdh.cpp

#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/miller_rabin.hpp>
#include <boost/math/special_functions/prime.hpp>
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

#include <fstream>
#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#if __cplusplus <= 199711L && !defined(Force2011)
#include <pthread.h>
#else
#include <mutex>
#endif
#include <omp.h>

// Used to link in a C function from an object file.
#ifdef LINKASN1C
extern "C"
{
    int processBuffer(const char *textBuf, int len, char* out, int* size);
}
#endif

// Written by Jesse Daniel Mitchell (2017). 
// To-do: 
// I'll also come up with a reasonable standard format for exporting some of these values in a neutral way, either JSON or XML. [H]
// * This will support both Tri-Prime Tuples and regular DH Parameters. I could also modify it to allow really condensed signatures. (~60 Bytes to encode any set of parameters).   
// I need to consider whether I want to add support for decoding OpenSSL Files or not, (unber), for verification. [UL]
// * If I do that though, I'll need to add new code to the unber code, to support doing stuff from buffers.
// ** I think I've settled on not converting backwards using NSPDH. The high priority target will be coming up with a neutral format to convert from, with all relevant information.
// Add verification for DSA Parameters [M]
// Consider dropping the current exported XML Format, switch to your own format that will be easy to convert from. [H]
// Consider OpenSSH Moduli Files [H]

// Eventually I'll advocate for Discrete Log Cryptography over RSA, using Tuples and NSPs. 
// It's much cheaper, and can do Public Key Crypto. (I don't think Ephemeral is important). 
// Also parameters don't necessarily require entropy, allowing efficient encoding using CSPRNGs.
// While RSA would have an advantage (size of encrypted file), the differences could be mitigated 
// using CSPRNG Encoding. example format:
// [32 Byte Entropy Seed][1 Byte CSPRNG Configuration][1 Byte Hash Algorithm Configuration][2 Byte Tuple Size][2 Byte Pohlig Size][1 Byte "Moment" Count (Tuple)][1 Byte "Moment" Count (Pohlig)][4 Byte Offset (Pohlig)][4 Byte Offset (Modulus)]
// 32 + 1 + 1 + 1 + 1 + 4 + 4 + 2 + 2 = 48 Bytes to encode parameters completely.

// Also, it's quite easy to implement Ephemeral Public Key Cryptography using Discrete Logs, 
// though I don't see how it'd be necessary.
// Step 1   - s1 = Compute (pub_a) ^ priv_b mod p
// Step 2*  - s2 = Compute a standard Ephemeral Diffie-Hellman Key Exchange (could even take place on smaller params, not just the same parameters). 
// Step 3** - Combine s1 and s2. 

// *  This might theoretically be useful for semi-ephemeral key exchanges. Just rotate out a secondary public key.
// ** For step 3, I originally would've recommended multiplying the shared secrets together and moduloing (if the same parameters were used), but it's highly unnecessary.

// Alternate Scheme :
// The point of Ephemeral Public Key Crypto in this situation is to provide a layer of PFS on top of Authentication.
// It would probably be better to just perform an ephemeral key exchange and create a secure hash of the shared key 
// combined with a time stamp, and sign it. This is what we already do with DHE-DSA. 

// This is included as a workaround for older versions of g++.
#ifndef nullptr 
#define nullptr NULL
#endif

// This exists so that the code is compatible with multiple compilers.
// Mainly to help me with server access, I am supporting compilers that lack full C++11 support.
// I don't care about supporting Old Windows compilers, which is why I'm supporting only pthreads.
class Mutex
{
    private:
    #if __cplusplus <= 199711L && !defined(Force2011)
    pthread_mutex_t mut;
    #else
    std::mutex mut; 
    #endif

    public:
    Mutex()
    {
            #if __cplusplus <= 199711L && !defined(Force2011)
            pthread_mutex_init(&mut, NULL);         
            #endif
    }

    ~Mutex()
    {
        #if __cplusplus <= 199711L && !defined(Force2011)
        pthread_mutex_destroy(&mut);
        #endif
    }

    // Locks the Mutex.
    void Lock()
    {
        #if __cplusplus <= 199711L && !defined(Force2011)
        pthread_mutex_lock(&mut);
        #else
        mut.lock();
        #endif
    }

    // Unlocks the Mutex.
    void Unlock()
    {
        #if __cplusplus <= 199711L && !defined(Force2011)
        pthread_mutex_unlock(&mut);
        #else
        mut.unlock();
        #endif
    }
};

// Included as a workaround for Unix/Windows compatibility.
#if defined(_WIN32) || defined(_WIN64)
#include <windows.h>
#else
#include <unistd.h>
void Sleep(int x)
{
    // Honestly, 1024 and 1000 are close enough, so that's why I'm using that.
    // Worst case scenario is that it takes an extra 24ms for a thread to wake up. 
    // This will hardly be a noticeable factor in the scheme of things. 
    usleep(x << 10); 
}
#endif

using namespace std;

using namespace boost::random;
using namespace boost::math; 
using namespace boost::multiprecision;

using boost::format;
using boost::lexical_cast;

// Finds the bitsize. 
int logit(cpp_int val)
{   
    return msb(val)+1;
}

// Define the types of generators. //
typedef independent_bits_engine<mt19937, 8192*2, cpp_int> generator_type;
typedef independent_bits_engine<mt19937, 256, cpp_int> generator_type2; // Using a smaller generator here in hopes of speeding up the Miller-Rabin Primality test.

// Defines the available generators
generator_type gen(time(0));
generator_type2 gen2(time(0));

#define generateIt(size) gen() >> ((8192*2) - size);
#define NSPDH_SEARCH 0
#define NSPDH_POHLIG_FOUND  1
#define NSPDH_MODULUS_FOUND 2
#define NSPDH_DHPARAM (1<<2)
#define NSPDH_DSAPARAM (2<<2)

// Tests whether the given cpp_int is prime or not.
// It print out a "." when a potential prime is found. This is to
// give it more ssl-like behavior.
char fastPrimeC(const cpp_int& v)
{   
    // Iterate over a small list of primes to give a tiny boost to the primality checking
    // This is mainly used so we can output a '.' 
    for(int i = 1; i <= 1024 + 512; i++)
    {
        if(v == prime(i)) return 1;
        if(!(v % prime(i))) return 0; 
    }

    // Print out a dot to indicate that it may have potentially found a prime value.
    cout << "." << flush;

    // Perform a Miller-Rabin Test.
    return miller_rabin_test(v, 9, gen2); 
}


// Finds all the [unique] prime factors of a small integer. 
vector<int> factor(int val)
{
    int c = 0;
    vector<int> res;
    
    // While the current prime is less than or equal to the value . . . 
    while(prime(c)*prime(c) <= val)
    {
        // If the number is divisible by the prime
        if(!(val % prime(c)))
        {
            // Push the prime into the collection
            res.push_back(prime(c));

            // Completely remove the factor from the value.
            while(!(val % prime(c)))
            {
                val /= prime(c);
            }
        }

        // Go to the next prime.
        c++;
    }

    // If the value isn't 1, then the current value is prime.
    if(val != 1)
    {
        res.push_back(val);
    }

    
    // return the collection
    return res;
}

// Checks whether a given generator is valid primitive root or not.
char checkGenerator(const cpp_int& proposed, const cpp_int& modPrime, const cpp_int& phPrime, int smallVal)
{
    vector<int> facts = factor(smallVal);
    cpp_int tot = modPrime - 1;

    // Computes g^(totient/2) 
    if(powm(proposed,tot/2,modPrime) == 1) return 0;
    
    // Computes g^(totient/pohlig-hellman prime factor)
    if(powm(proposed, tot/phPrime, modPrime) == 1) return 0;

    // Computes g^(totient/p_i), which is the generator raised to the power of the totient divided by each unique prime factor.
    for(int i = 0; i < facts.size(); i++)
    {
        if(powm(proposed, tot/facts[i], modPrime) == 1) return 0; 
    }

    // This was a primitive root.
    return 1;
}


// Checks the generator, even if it's a quadratic residue. This is for the verification function.
char checkGeneratorInclusive(const cpp_int& proposed, const cpp_int& modPrime, const cpp_int& phPrime, int smallVal)
{
    // Checks if it is a Quadratic Residue.
    if(powm(proposed, (modPrime-1)/2, modPrime) == 1)
    {   
        // has second bit flag that it was a quadratic.
        return 3; // binary = 11
    }
    
    return checkGenerator(proposed, modPrime, phPrime, smallVal);
}

// Finds the number of generators of a cyclic group (totient of totient). 
// Totient(modPrime) * product of each factor (1 - 1/p_i),
// Totient(modPrime) * product of each factor ((p_i-1)/p_i)
cpp_int numberOfGenerators(const cpp_int& modPrime, const cpp_int& phPrime, int smallVal)
{
    cpp_int val = modPrime - 1;

    // Unnecessary instruction left commented for consistency. 
    // val *= 1; 
    val /= 2; 

    val *= phPrime-1; 
    val /= phPrime; 

    // Gets all the factors of the offset :)
    vector<int> facts = factor(smallVal);

    // Iterates over all the factors.
    for(int i = 0; i < facts.size(); i++)
    {
        // Skips 2 as a factor, since we already did it.
        if(facts[i] == 2) continue;
        val *= facts[i]-1;
        val /= facts[i]; 
    }
    
    return val; 
}

// Used for the XML generation. Produces something in the form tag="internal"
string quotes(const string& tag, const string& internal)
{
    return tag + "=\"" + internal + "\"";
}

// Used for XML generation, just casts the integer to a string. I could use C++11 libraries for this,
// but since I already have a Boost Dependency, I'll use the boost library for backwards compatibility. 
string quotes(const string& tag, int x)
{
    return quotes(tag, lexical_cast<string>(x));
}

// Gets the bytes for a given prime.
vector<char> getByteArray(cpp_int v)
{
    vector<char> res;
    while(v != 0)
    {
        res.push_back((v & 255).convert_to<char>());
        v >>= 8; 
    }
    return res;
}

// Prints the byte array for prime values.
void printByteArray(vector<char>& r, ostream& file)
{
    for(int i = r.size() -1; i >= 0; i--)
    {
        file << "&#x" << format("%02x")  % (int)((unsigned char)r[i]) << ";";  
    }
}

// Creates the XML File to be converted by a tool like "enber". 
void createXML(vector<cpp_int>& vec, ostream& file)
{
    // This code could be really, really buggy.
    // Thus far, it has worked quite well though. :)

    // Gets the bytes for the modulus and the generator. 
    int tl1 = 2; 

    // Computes the length of all that the sequence contains.
    int length = tl1;

    // Used to extract the data from each element.
    vector< vector<char> > data;
    vector<int> tls; 

    // Iterates through each element.
    for(int i = 0; i < vec.size(); i++)
    {
        // Gets the bytes from each.
        vector<char> x = getByteArray(vec[i]);
        data.push_back(x);

        // "Length" Container
        int tl =2;
        
        // Computes the size of the "length" container
        if((int)(data[i].size()) +1 > 127) tl++;
        if((int)(data[i].size()) +1 > 256) tl++;

        // Stores the size of the length container.
        tls.push_back(tl);

        
        // Used to fix the ASN.1 Encoding for sizes divisible by 8.
        char adjust = 0;

        if(vec[i] != 0)
        {
            // Detects if divisible by 8.
            adjust = (msb(vec[i]) + 1) % 8;

            // Adjusts the value
            if(adjust) adjust = 0;
            else adjust = 1;
        }
        else
        {
            // Ensures that there is a single byte to use.
            data[i].push_back(0);
        }

        // Adjust the overall length.
        length += tl + (int)(data[i].size() + adjust);
    }



    // Increases the sequence "length" container size.
    if(length > 127)
    {
        tl1++;
        length++;
    }

    // Increases the sequence "length" container size.
    if(length > 256)
    {
        tl1++; 
        length++;
    }

    // Offsets are computed by the number of bytes that are supposed to trail behind the section, you'll see them in the "O"s written below.
    // The "T"s are just tags that indicate what type of value they are in ASN1 Encoding.
    // The "TL"s are apparently the containers that describe the length of the section.
    // The "V"s are the actual size of the section. 
    // The "A"s describes what type it is.

    // Write the Sequence Header Tag
    file << "<C " << quotes("O", "0") << " " << quotes("T", "[UNIVERSAL 16]") << " " 
         << quotes("TL", tl1) << " " << quotes("V", length - tl1) << " " << quotes("A", "SEQUENCE") << ">" << endl; 

    int offset = tl1;

    // adds each element to the sequence.
    for(int i = 0; i < vec.size(); i++)
    {
        // Used to fix the ASN.1 Encoding for sizes divisible by 8.
        char adjust = 0;

        if(vec[i] != 0)
        {
            // Detects if divisible by 8.
            adjust = (msb(vec[i]) + 1) % 8;

            // Adjusts the value
            if(adjust) adjust = 0;
            else adjust = 1;
        }

         // Write the Integer Header Tag (for the Modulus).
        file << "<P " << quotes("O", offset) << " " << quotes("T", "[UNIVERSAL 2]") << " " 
            << quotes("TL", tls[i]) << " " << quotes("V", (int)(data[i].size() + adjust)) << " " << quotes("A", "INTEGER") << ">"; 

        if(adjust) file << "&#x00;";

        // Write the bytes for the Modulus.
        printByteArray(data[i], file);
        
        // Close the Integer Tag
        file << "</P>" << endl;

        // Adjusts the offset for the next item by the current item's length.
        offset += tls[i] + (int)(data[i].size()) + adjust;
    }

    // Close the sequence.
    file << "</C " << quotes("O", length) << " " << quotes("T", "[UNIVERSAL 16]") << " " 
        << quotes("A", "SEQUENCE") << " " << quotes("L", length) << ">" << endl;
         
}

// Displays the flags possible.
// Needs more verbosity.
void displayHelp(char* exeName, bool verbose = false)
{
    // Basic Help
    cout << exeName << " [<bitsize>|-help|-verify <filename>] [-out <fileName>] [args]" << endl;
    
    // Attempts to give some documentation of the options.
    // Todo: Clarify -tolerate, and include the shortened tags.
    if(verbose)
    {
        cout << "-randomize" << "\tRandomizes the generator" << endl;
        cout << "-quadratic" << "\tExports a Quadratic Residue Generator in the out file. (Wei Dai Recommendation)" << endl;
        cout << "-smallest"  << "\tExports the smaller of the two generators (between Quadratic Residue and Primitive Root)" << endl;
        cout << "-multiple"  << "\tSearches for all associated NSPs associated with the Pohlig-Hellman Prime." << endl;
        cout << "-max"  << "\t\tSpecifies the maximum number of NSPs to find associated with the Pohlig-Hellman prime. (-multiple required)." << endl;
        cout << "-n <value>" << "\tSets the restricted range for finding NSPs." << endl;
        cout << "-bn <value>" << "\tSets the restricted range (in bits) for finding NSPs." << endl;
        cout << "-8" << "\t\tForces the Modulus length to be divisible by 8." << endl;
        cout << "-verify <file>" << "\tVerifies the parameters in a file. OpenSSL Formats are not supported (yet)." << endl;
        cout << "-tuple <size>" << "\tFinds a smaller prime of the specified size that will be related to the Pohlig-Hellman prime, creating a Tri-Prime Tuple. For use with signatures." << endl;
        cout << "-out <fileName>" << "\tExports the parameters to a file." << endl;
        #ifdef LINKASN1C
        cout << "-der" << "\t\tExports an OpenSSL-compatible .der file." << endl;
        cout << "-pem" << "\t\tExports an OpenSSL-compatible .pem file. (Experimental)" << endl;
        #endif
        cout << flush;
    }
}


static const std::string base64_chars = 
             "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
             "abcdefghijklmnopqrstuvwxyz"
             "0123456789+/";

std::string base64_encode(unsigned char const* bytes_to_encode, unsigned int in_len) {
  std::string ret;
  int i = 0;
  int cnt = 0;
  unsigned char char_array_3[3];
  unsigned char char_array_4[4];

  while (in_len--) {
    char_array_3[i++] = *(bytes_to_encode++);
    if (i == 3) {
      char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
      char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
      char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
      char_array_4[3] = char_array_3[2] & 0x3f;

      for(i = 0; (i <4) ; i++)
      {
        ret += base64_chars[char_array_4[i]];
        if(!(++cnt % 64)) ret += '\n';
      }
      i = 0;
    }
  }

  if (i)
  {
    int j; 
    for(j = i; j < 3; j++)
      char_array_3[j] = '\0';

    char_array_4[0] = (char_array_3[0] & 0xfc) >> 2;
    char_array_4[1] = ((char_array_3[0] & 0x03) << 4) + ((char_array_3[1] & 0xf0) >> 4);
    char_array_4[2] = ((char_array_3[1] & 0x0f) << 2) + ((char_array_3[2] & 0xc0) >> 6);
    char_array_4[3] = char_array_3[2] & 0x3f;

    for (j = 0; (j < i + 1); j++)
    {
      ret += base64_chars[char_array_4[j]];
      if(!(++cnt % 64)) ret += '\n';
    }

    while((i++ < 3))
      ret += '=';
  }

  return ret;

}

// Exports the parameters.
void exportParameters(const string& outFile, vector<cpp_int>& params, char convert)
{
    // Creates a string stream for an easy memory buffer.
    ostringstream conversion;

    
    // creates the xml version of the parameters.
    createXML(params, conversion);

    #ifdef LINKASN1C
    // If ASN1C is linked, it checks if the conversion flag is enabled.
    if(convert)
    {
        // If so, it converts the parameters and outputs a .der
        char outputBuffer[8192*2];
        int size; 
        
        // Processes the buffer using the utility from enber. 
        processBuffer(conversion.str().c_str(), conversion.str().length(), outputBuffer, &size);

        if(convert == 1)
        {
            // Outputs the returned buffer.
            ofstream ofs((outFile + ".der").c_str(), ios::binary);
            ofs.write(outputBuffer, size);
            ofs.close();
        }
        else if((convert & 2) == 2)
        {
            ofstream ofs((outFile + ".pem").c_str());
            string outputString = base64_encode((unsigned char*)outputBuffer, size);
            
            // Used for the correct PEM Header.
            string pemType;

            // Erases the first two bits, for easy comparisons.
            convert &= ~3; 

            // Sets the correct PEM Header
            if(convert == NSPDH_DHPARAM)
            {
                pemType = "DH PARAMETERS";
            }
            else if (convert == NSPDH_DSAPARAM)
            {
                pemType = "DSA PARAMETERS";
            }

            // Outputs the header, body, and footer.
            ofs << "-----BEGIN " << pemType << "-----" << endl;
            ofs << outputString << endl;
            ofs << "-----END " << pemType << "-----" << endl;

            ofs.close(); 
        }
        else
        {
            // Nothing yet.
        }
    }
    else
    {
    // Outputs the xml. 
    #endif
        ofstream ofs((outFile + ".xml").c_str());
        ofs << conversion.str();
        ofs.close();
    #ifdef LINKASN1C
    }
    #endif
}


// Exports the parameters. (Wrapper function)
inline void exportParameters(const string& outFile, cpp_int& modulusPrime, cpp_int& generator, char convert)
{
    // Creates a vector for the parameters and pushes them in order.
    vector<cpp_int> params;
    params.push_back(modulusPrime);
    params.push_back(generator);

    // Exports the parameters.
    exportParameters(outFile, params, convert);
}

// This is used as a hack to allow me to decide whether I want to return something or not.
static cpp_int hack; 

// Verifier function.
// Needs refractoring and some expansion. 
// It works, but I need to make it so that it's a pure function without output.
// Note that this is not proper XML Parsing, but it is not completely necessary.
// I am relying on the user to pass in valid input. 
char verifier(char* p, int range = 4096, int trials = 10, cpp_int& rmodulus = hack, cpp_int& rgenerator = hack)
{
    char results[3];
    int i(0), N(1);

    // Used for hex conversion
    char b(0);

    // Used for validation.
    cpp_int modulus, pohlig, generator;

    for(int j = 0; j < 2; j++) 
    {
        // Skips to the first character after the first '>'
        while(p[i++] != '>'); 
        // Skips to the first character after the second '>'
        while(p[i++] != '>');
        
        // Checks if we're currently on a hex value.
        while(p[i] == '&')
        {
            // Skips to the first character after the 'x'.
            while(p[i++] != 'x');

            // Convert two digit hex to decimal.
            // Top 4 Bits 
            if(p[i] >= '0' && p[i] <= '9')
            b = (p[i++] - '0')*16;
            else
            b = (p[i++] - 'a' + 10)*16;

            // Bottom 4 bits.
            if(p[i] >= '0' && p[i] <= '9')
            b += (p[i++] - '0');
            else
            b += (p[i++] - 'a' + 10);
            // End hex conversion.
            
            // Consume the hex value. 
            generator <<= 8;
            generator |= (unsigned char)b; 

            // Go onto the next hex value or end.
            i++;    
        }

        // This is a shortcut done to shorten the code. 
        if(!j)
        {
            modulus = generator;
            generator = 0;
        }
    }

    // Prints out the supposed Modulus Prime.
    cout << "Modulus Prime (" << msb(modulus)+1 << " bits): " << modulus << endl;
    cout << endl;
    
    // Used to compute the actual Pohlig.
    pohlig = (modulus - 1)/2;
    int c(0);
    
    // Factor out every prime within range.
    while (c < 10000 && N <= range && prime(c) <= range)
    {
        while(!(pohlig % prime(c))) 
        {
            pohlig /= prime(c);
            N *= prime(c);
        }
        c++;
    }

    // Prints out the computed Pohlig-Hellman prime.
    cout << "Pohlig-Hellman Prime (" << msb(pohlig)+1 << " bits): " << pohlig << endl;
    cout << endl << "N: " << N << endl;

    // Used to parallelize the tests.
    Mutex mut;
    int branchCount(0);

    // A silly (but simple) way to parallelize this process.
    #pragma omp parallel
    {
        do
        {
            // Grabs the current branch.
            mut.Lock();
            int branch(branchCount++);
            mut.Unlock();
            switch(branch)
            {
                case 0:
                    // Check the first prime value.
                    results[0] = miller_rabin_test(modulus, trials, gen2);
                    break;
                case 1:
                    // Check the second prime value.
                    results[1] = miller_rabin_test(pohlig, trials, gen2);
                    break;
                case 2: 
                    // Check the Generator value.
                    results[2] = checkGeneratorInclusive(generator, modulus, pohlig, N);
                    break;
                default:
                    break;
            }
        } while(branchCount < 3);
    }

    // Sets the generator to zero in case of some false positive.
    if(!(results[0] && results[1])) results[2] = 0;
    
    // Prints out the generator.
    cout << "Generator: " << generator << " ";
    
    // Outputs if it was quadratic.
    if(results[2] & 2) cout << "(Quadratic Residue)"; 
    
    // Print out verification results.
    cout << endl << endl;
    cout << ((results[0]) ? "Modulus Verified" : "Modulus Invalid") << endl;
    cout << ((results[1]) ? "Pohlig Verified" : "Pohlig Invalid (Consider expanding the tolerated range?)") << endl;
    
    // We can't verify the generator without knowing the Pohlig-Hellman prime, so we create a requirement for it.
    if(results[1])
    cout << ((results[2]) ? "Generator Verified" : "Generator Not Verified") << endl;

    // Used to return these parameters.
    rgenerator = generator;
    rmodulus = modulus;

    char res = 0 | (results[0]) | (results[1] << 1) | (results[2] << 2);
    return res;
}

// Generates a prime value of the requested size.
// While not preferable, the volatile char* is included to allow the function to sleep, yielding to other threads. 
cpp_int generatePrime(int size, volatile char* completionStatus)
{            
    // Generates a random value.
    cpp_int primeVal = generateIt(size);
            
    // Enforces size restrictions (and that it is odd)
    primeVal |= (((cpp_int)1) << (size-1)) | 1;
                
    while(!fastPrimeC(primeVal))
    {
        // If a Pohlig-Hellman Prime has been found by another thread, sleep. 
        while(*completionStatus == NSPDH_POHLIG_FOUND) Sleep(2000);
        
        // If a Modulus Prime was found, just quit.
        if(*completionStatus == NSPDH_MODULUS_FOUND) break;
        
        // If the value wasn't prime, skip ahead to the next odd integer.
        primeVal += 2;
    }

    return primeVal;
}


// Generates a Prime Tuple, and intentionally attempts to "maximize" its Pohlig-Hellman strength.
// While this is not technically necessary (not a requirement of DSA), I think that it is probably worth it.
// One thing to consider is that since the added Prime Factors are unknown, the GFNS (for factorization) must be applied.
cpp_int generatePrimeTuple(int size, cpp_int base, volatile char* completionStatus)
{
    // Generates a random value.
    cpp_int offset, enhancer;
    int count(65536), retries(0);

    while(true)
    {
        // Used in case there are no results to be had. 
        if(count++ & 65536)
        {
            if(++retries == 10) return 0;
            
            offset = gen2();
            enhancer = 2;

            // Attempts to boost the known strength of the tuple.
            if((msb(base)+1) < 1024 && (size - (int)(msb(base)+1) - 1024) >= 256)
            {
                enhancer = 2*generatePrime(1024, completionStatus);
            }
            else    
            // Attempts to boost it again, just smaller.
            if((msb(base)+1) < 512 && (size - (int)(msb(base)+1) - 512) >= 256)
            {
                enhancer = 2*generatePrime(512, completionStatus);
            }
            else 
            // Attempts to boost it again, just smaller.
            if((msb(base)+1) < 256 && (size - (int)(msb(base)+1) - 256) >= 256)
            {
                enhancer = 2*generatePrime(256, completionStatus);
            } 
 
            count = 0;
        }

        // If a Pohlig-Hellman Prime has been found by another thread, sleep. 
        while(*completionStatus == NSPDH_POHLIG_FOUND) Sleep(2000);
        
        // If a Modulus Prime was found, just quit.
        if(*completionStatus == NSPDH_MODULUS_FOUND) break;
    
        // This is used to place the randomly generated offset into a range where it can create the correct bit size. 
        cpp_int shiftedOffset = offset;
        
        // Shifts the offset.
        while(msb(shiftedOffset*base*enhancer + 1) + 1 < size) 
        {
            shiftedOffset <<= 1;

            //Gives it some extra entropy.
            shiftedOffset |= (gen2() & 1);
        }

        // If it created a value of a correct size,
         if(msb(shiftedOffset*base*enhancer + 1) + 1 == size)
    
        // check if it is prime and return it if it is. 
        if(fastPrimeC(shiftedOffset*base*enhancer + 1)) return base*shiftedOffset*enhancer + 1;
        offset--;
    }

    // Used just in case.
    return 0;
}


int main(int argc, char **args)
{
    // Todo: Create a C-Style struct for these parameters.
    // This is more of a collection than a class, 
    // so a struct is more reasonable.

    // Flags for the Program //
    int size(2048);
    int maxModuli(0); 
    int tuple(0);
    char convert(0);
    bool multiple(false);
    bool divisibleBy8(false);
    bool randomizeGenerator(false);
    char generatorMode(0);
    char *outFile = nullptr;
    char *verify = nullptr;
    long long maxN(4096ll | (1ll << 40));
    // End Flags //

    // Mutex for adding the parameters.
    Mutex mut; 

    // Simple command parser.
    for(int i = 1; i < argc; i++)
    {
        #ifdef LINKASN1C
        // Converts the parameters to the .der format OpenSSL understands.
        if(!strcmp(args[i], "-der")) convert = 1;
        else

        // Converts the parameters to the .pem format OpenSSL understands. (Experimental)
        if(!strcmp(args[i], "-pem")) convert = 2;
        else
        #endif

        // Displays the verbose version of the help.
        if(!strcmp(args[i], "-h") || !strcmp(args[i], "-help")) 
        {
            displayHelp(args[0], true);
            return 0;
        }
        else

        // Enforces a divisible by 8 restriction, which is required by some Certificate Authorities.
        // Some places enforce divisible by 32 and 64, but I feel those requirements are too harsh.
        // I have to put my foot down somewhere.  
        if(!strcmp(args[i], "-8")) divisibleBy8 = true;
        else

        // Used to generate a prime tuple structure, which is a set of parameters that can support both encryption and signatures efficiently.
        if(!strcmp(args[i], "-tuple") || !strcmp(args[i], "-t")) tuple = atoi(args[i++ + 1]);
        else
        
        // Sets the maximum size of the offset (n). (for 2np + 1) 
        if(!strcmp(args[i], "-n") || !strcmp(args[i], "-N")) maxN = atoi(args[i++ + 1]);
        else

        // Sets the maximum size of the offset in bits (n). (for 2np + 1) 
        if(!strcmp(args[i], "-bn") || !strcmp(args[i], "-bN")) maxN = 1 << atoi(args[i++ + 1]);
        else

        // Randomize tries to randomize the generator. 
        if(!strcmp(args[i], "-r") || !strcmp(args[i], "-randomize")) randomizeGenerator = true;    
        else

        // Attempts to use the same Pohlig-Hellman Prime to find multiple Moduli. 
        if(!strcmp(args[i], "-m") || !strcmp(args[i], "-multiple")) multiple = true;    
        else

        // Sets a max number of Moduli to find from a Pohlig-Hellman prime.
        if(!strcmp(args[i], "-max")) maxModuli = atoi(args[i++ + 1]);    
        else

        // Quadratic exports a quadratic residue generator in the output file.
        if(!strcmp(args[i], "-q") || !strcmp(args[i], "-quadratic")) generatorMode = 1;
        else

        // Guarantees that the generator exported is the smallest (between Quadratic and Primitive Root).
        if(!strcmp(args[i], "-s") || !strcmp(args[i], "-smallest")) generatorMode = 2;
        else

        // Attempts to verify the parameters within the given file.
        if(!strcmp(args[i], "-v") || !strcmp(args[i], "-verify")) verify = args[i++ + 1];
        else
        
        // Allows you to specify the output file.
        if(!strcmp(args[i], "-o") || !strcmp(args[i], "-out")) outFile = args[i++ + 1];
        else

        // Assumes it is the size parameter and gets it.
        {
            size = atoi(args[i]); 
        }    
    }

    // Checks if divisible by 8 flag enabled and not modified.
    if(divisibleBy8 && (maxN & (1ll << 40)))
    {
        cout << "Warning: Divisible by 8 restrictions enforced. We would recommend modifying your tolerated N range for increased speed." << endl << flush;
    }

    // Checks if the multiple flag is enabled to set default max number of Moduli.    
    if(multiple)
    {
        // if the max is currently set to zero, set it to be as many moduli as possible :).
        if(!maxModuli) maxModuli = maxN;
    }
    // otherwise set it to one.
    else maxModuli = 1;

    // Checks if the size isn't divisible by 8 while the divisibleBy8 flag is enabled.
    // The algorithm will find compatible primes more quickly if your Pohlig-Hellman is also divisible by 8.
    if(divisibleBy8 && size % 8)
    {
        cout << "Warning: Your requested Pohlig-Hellman size is not divisible by 8. The algorithm will work more quickly for these parameters if it is. (Yes, even if you increase the size). ";
        cout << "Could we recommend size " << size + (8 - (size % 8)) << "?" << endl << flush;
    }

    // Removes the "Not Modified" flag if it currently exists.
    if(maxN & (1ll << 40))
        maxN ^= (1ll << 40);

    // If there isn't an outFile, then don't convert.
    if(outFile == nullptr) convert = false;

    // Checks if a verification was requested.
    if(verify != nullptr)
    {
        char *buf = new char[65536];
        
        // Gets the file and its size.
        ifstream ifs(verify, ios::binary | ios::ate);
        ifstream::pos_type pos = ifs.tellg();
        
        // Reads in the entire file.
        ifs.seekg(0, ios::beg);
        ifs.read(buf, pos);
        
        // Used to return the parameters.
        cpp_int modulus, generator;

        // Verifies the data within the file.
        char result = verifier(buf, (int)maxN, 10, modulus, generator);

        // Allows the person to export the parameters upon verification.
        if((result & 7) == 7 && outFile != nullptr)
        {
            // Used to tell it to flag it as a DH Parameter.
            if(convert & 2) convert |= NSPDH_DHPARAM;
            exportParameters(string(outFile), modulus, generator, convert);
        }

        delete[] buf;
        return 0;
    }

    // Storage for the parameters.
    vector<cpp_int> phPrimes;
    vector<cpp_int> modulusPrimes;
    vector<cpp_int> tupleBases;
    vector<int> offsets;

    // hacks for the algorithm. 
    // technically there is a potential race condition but the odds of it happening are truly astronomical.
    // so I won't use a mutex.
    volatile char completionStatus = 0;

    // Allows nested parallelization, which will work well since I am making the other threads sleep.
    omp_set_nested(1);
    #pragma omp parallel
    {
        cout << "t";
        bool internal = false;
        cpp_int primeVal, tupleBase; 

        // Searching for Prime values.
        while(completionStatus != 2)
        {
            // If a Prime-Tuple is requested, generate one.
            // Tuple Parameters will compute a requested bit-size factor for the Pohlig-Hellman prime.
            // This can be used for signatures (DSA-like parameters). 
            if(tuple)
            {
                // Generates a Prime "Tuple".
                tupleBase = generatePrime(tuple, &completionStatus);
                primeVal = generatePrimeTuple(size, tupleBase, &completionStatus);
                if(primeVal == 0) continue;
            }
            else
            {
                // Generates a Prime Value.
                primeVal = generatePrime(size, &completionStatus);
            }
                   
            // Sleeps if another Pohlig was already found. (This exists here just in case two are found at the same time.)
            while(completionStatus == NSPDH_POHLIG_FOUND) Sleep(2000);
                
            // If a Modulus prime was found then break. 
            if(completionStatus == NSPDH_MODULUS_FOUND) break;

            // Tell the other threads to sleep.
            completionStatus = NSPDH_POHLIG_FOUND;

            // Output a dash to tell the user it found a Pohlig-Hellman Prime.
            cout << "-";

            // 2*p 
            cpp_int temp = primeVal << 1; 

            // Computes the n
            #pragma omp parallel for
            for(long long i = 1; i <= maxN; i++)
            { 
                unsigned long long mul = i;

                // Checks if a result has already been found/if we're supposed to be computing
                // more values
                if(phPrimes.size() != maxModuli)
                {

                    // Checks if divisible by 8 flag is enabled, if it is,
                    if(divisibleBy8)
                    {
                        // it attempts to hurry it up.
                        while((msb(temp*mul + 1)+1) % 8)
                        {
                            mul <<= 1;                             
			}

                        // Checks if the current multiplier is greater than the current max N.
                        if(mul > maxN) continue;
                    }

                    // is Prime(2*p*n + 1)
                    if(fastPrimeC(temp*mul + 1))
                    {                     
                        // Tell other threads it found an NSP, and wake them up.
                        completionStatus = 2; 
                                    
                        // Mark that this was the thread that found it.
                        internal = true;

                        // Add the parameters and break.
                        mut.Lock();
                        if(phPrimes.size() != maxModuli) 
                        {
                            phPrimes.push_back(primeVal);
                            modulusPrimes.push_back(temp*mul + 1); 
                            tupleBases.push_back(tupleBase);
                            offsets.push_back(mul);
                        }                         
                        mut.Unlock();

                        cout << "+" << flush;
                    }
                }
            }
            
            // If it didn't find an NSP within bounds, print an x and flag for randomization (it seems to find them faster with this behavior). 
            if(!internal)
            {
                cout << "x";
                
                // Wakes up the other threads.
                completionStatus = NSPDH_SEARCH;
            }
        }
    }
    
    for(int i = 0; i < phPrimes.size(); i++)
    {
        // Grab the computed parameters.
        cpp_int phPrime = phPrimes[i];
        cpp_int modulusPrime = modulusPrimes[i];
        cpp_int tupleBase = tupleBases[i];
        int offset = offsets[i];

        // Print out the parameters
        cout << endl << endl; 
        cout << "Pohlig-Hellman Prime: " << phPrime << endl << endl;
        if(tuple) cout << "Tuple Base: " << tupleBase << endl << endl;
        cout << "Modulus Prime: " << modulusPrime << endl << endl;
        cout << "Offset N: " << offset << endl; 
        
        // Print out the bit sizes.
        cout << logit(phPrime) << " "; 
        if(tuple) cout << logit(tupleBase) << " ";
        cout << logit(modulusPrime) << " " << logit(offset) << endl;

        // Find a valid primitive root.
        cpp_int c;

        if(randomizeGenerator)
        {
            // Finds a random primitive root.
            c = (gen() % modulusPrime).convert_to<unsigned int>();
        }
        else
        {
            // Finds the smallest possible primitive root.
            c = 2; 
        }
        
        // used to find a quadratic-residue generator
        // Note: Potential optimization, when c = 2 test quadratic = 3, if that fails, quadratic is immediately 4 without a test.
        cpp_int quadraticGenerator(c);
        
        // While the value isn't a generator, increment.
        while(!checkGenerator(c, modulusPrime, phPrime, offset))
        {
            c++;
        } 

        // Tries to find a quadratic generator, which will have g^(totient/2) == 1,
        while(powm(quadraticGenerator, (modulusPrime-1)/2, modulusPrime) != 1)
        {
            quadraticGenerator++;
        }
        
        // Prints out the generators.
        cout << "Primitive Root Generator: " << c << endl;
        cout << "Quadratic Residue Generator: " << quadraticGenerator << endl; 
        

        // If there was a specified output file, it exports the parameters.
        if(outFile != nullptr)
        {
            // Only export this once.
            if(tuple && i == 0)
            {
                cpp_int g = 2;     
                cpp_int tot = phPrime-1;
                tot = tot/tupleBase;

                // Finds a generator of Multiplicative Order q (tupleBase).
                while(powm(g, tot, phPrime) == 1) g++;
                g = powm(g, tot, phPrime);
                
                // Pushes the parameters (in the correct order) to the sequence.
                vector<cpp_int> vec;
                vec.push_back(phPrime);
                vec.push_back(tupleBase);
                vec.push_back(g);
                
                // Specifies that this is a DSA Parameter.
                char dsaconvert = convert; 
                if(dsaconvert & 2) dsaconvert |= NSPDH_DSAPARAM;

                /*
                cpp_int k = 3777;
                cpp_int x = 234;
                cpp_int y = powm(g, x, phPrime);
                cpp_int k1 = powm(k, tupleBase-2, tupleBase);
                cpp_int h = 7113;

                cpp_int r = powm(g, k, phPrime);
                r = r % tupleBase;

                cpp_int s = k1 * (h + r*x);
                s = s % tupleBase;

                cout << "r: " << r << endl;
                cout << "s: " << s << endl;

                cpp_int w = powm(s, tupleBase-2, tupleBase);
                cpp_int u1 = w*h;
                u1 = u1 % tupleBase;

                cpp_int u2 = w*r;
                u2 = u2 % tupleBase;

                cpp_int v1 = powm(g, u1, phPrime);
                cpp_int v2 = powm(y, u2, phPrime);

                cpp_int v = (v1*v2) % phPrime;
                v = v % tupleBase;

                cout << "v: " << v << endl;
                */

                // Exports the parameters.
                exportParameters(outFile + string("_dsa"), vec, dsaconvert);
            }

            // Used for outputting "multiple" files.
            string variant = "-" + lexical_cast<string>(i);
            if(!multiple) variant = "";

            //if quadratic mode is toggled, switch the generator into a quadratic residue.
            if(generatorMode == 1)
            {
                c = quadraticGenerator;
            }
            else
            // if "smallest" mode is enabled, switch the generator into the smallest between
            // the quadratic residue and the primitive root.
            if(generatorMode == 2)
            {
                c = min(c, quadraticGenerator);
            }
            
            // Specifies that these are Diffie-Hellman Parameters.
            char dhconvert = convert; 
            if(dhconvert & 2) dhconvert |= NSPDH_DHPARAM;

            // Exports the parameters.
            exportParameters(string(outFile) + variant, modulusPrime, c, dhconvert);
           
        } 
    }
    return 0;
}