Ulam

Computer Science Level 4

Image Credit: Wolfram Image Credit: Wolfram What percent of the first three million positive integers yield a prime when the following polynomial is evaluated by them:

n 2 + n + 41 n^2 + n + 41


The answer is 24.

View wiki
Thaddeus Abiy by Thaddeus Abiy , 25, Ethiopia

This section requires Javascript.
You are seeing this because something didn't load right. We suggest you, (a) try refreshing the page, (b) enabling javascript if it is disabled on your browser and, finally, (c) loading the non-javascript version of this page . We're sorry about the hassle.

2 solutions

Rushikesh Jogdand
Jun 29, 2016

Cheating! 

1
2
3
4
5
 
from gmpy import is_prime
count=0
for i in range(1,3000001):
    if is_prime(i**2+i+41):count+=1
print(count*100/3000000)

1 Helpful 0 Interesting 0 Brilliant 0 Confused
Piotr Idzik
May 9, 2020

Here is my parallel C++ code. I was using Fermat primality test (after all, there are only few pseudo primes more than real primes). The code is iterative. I start with the list [ 1 , 2 , … , 3 ⋅ 1 0 6 ] [1, 2, \ldots, 3\cdot10^6] . In each iteration cf. function s i n g l e I t e r a t i o n P a r \mathtt{singleIterationPar} ), I remove from the list the numbers n n , for which the value n 2 + n + 41 n^2 + n + 41 does not pass the Fermat primality test (with single base/witness).

The Fermat primality test is implemented in the function i s P r o b P r i m e \mathtt{isProbPrime} . 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
 
//solution of: https://brilliant.org/problems/ulam/

#include <iostream>
#include <vector>
#include <thread>
#include <random>
#include <ctime>
#include <chrono>

template <typename ValueType>
ValueType getPolynomialValue(const ValueType& inVal)
{
    //returns the value of the polynomial from the problem/puzzle
    return (inVal+ValueType(1))*inVal+ValueType(41);
}

template <typename ValueType>
ValueType binPowMod(const ValueType& aVal, const ValueType& expVal, const ValueType& modVal)
{
    //returns (aVal**expVal)%modVal
    ValueType resVal = ValueType(1);
    ValueType curExp = expVal;
    ValueType curMultip = aVal;
    while(curExp > 0)
    {
        if(curExp%ValueType(2) == ValueType(1))
        {
            resVal *= curMultip;
            resVal %= modVal;
        }
        curMultip *= curMultip;
        curMultip %= modVal;
        curExp /= ValueType(2);
    }
    return resVal;
}

template<typename UnsignedIntType>
UnsignedIntType gcd(const UnsignedIntType _a, UnsignedIntType _b)
{
    UnsignedIntType a, b;
    if(_a > _b)
    {
        a = _a;
        b = _b;
    }
    else
    {
        a = _b;
        b = _a;
    }
    while(b > 0)
    {
        const UnsignedIntType div = a%b;
        a = b;
        b = div;
    }
    return a;
}

template<typename ValueType, typename BigValueType, typename RandomEngineType>
bool isProbPrime(const ValueType& inVal, RandomEngineType& randomEngine)
{
    bool res;
    if(inVal <= ValueType(3))
    {
        if(inVal <= ValueType(1))
        {
            res = false;
        }
        else //inVal == 2 || inVal == 3
        {
            res = true;
        }
    }
    else if(inVal%ValueType(2) == ValueType(0) || inVal%ValueType(3) == ValueType(0))
    {
        res = false;
    }
    else
    {
        std::uniform_int_distribution<ValueType> numGenerator(2, inVal-2);
        auto aVal = numGenerator(randomEngine);
        if(gcd(aVal, inVal) == ValueType(1))
        {
            res = binPowMod<BigValueType>(aVal, inVal-1, inVal) == 1;
        }
        else
        {
            res = false;
        }
    }
    return res;
}

template <typename ValueType, typename BigValueType, typename RandomEngineType>
class Worker
{
    private:
        std::vector<ValueType> inData;
        std::vector<ValueType> outData;
        std::size_t startPos, endLimit;
        RandomEngineType randomEngine;
    public:
        Worker(const std::vector<ValueType>& _inData,
               const std::size_t _startPos,
               const std::size_t _endLimit,
               const int inSeed) :
            inData(_inData),
            outData(),
            startPos(_startPos),
            endLimit(_endLimit),
            randomEngine(inSeed)
        {}
        void operator()()
        {
            this->outData.reserve(this->inData.size());
            for(std::size_t i = this->startPos; i < this->endLimit; ++i)
            {
                const auto curNum = this->inData.at(i);
                const auto curVal = getPolynomialValue(curNum);
                if(isProbPrime<ValueType, BigValueType, RandomEngineType>(curVal, this->randomEngine))
                {
                    this->outData.push_back(curNum);
                }
            }
        }
        size_t getOutputSize() const
        {
            return this->outData.size();
        }
        std::vector<ValueType> getOutputData() const
        {
            return this->outData;
        }
};

void updatePosData(const std::size_t dataSize,
                   const std::size_t pieceSize,
                   std::size_t& startPos,
                   std::size_t& endLimit)
{
    if(endLimit < dataSize)
    {
        startPos = endLimit;
        endLimit = startPos+pieceSize;
        if(endLimit < dataSize)
        {
            if(dataSize-endLimit < pieceSize)
            {
                endLimit = dataSize;
            }
        }
        else
        {
            endLimit = dataSize;
        }
    }
    else
    {
        endLimit = dataSize+1;
        startPos = endLimit;
    }

}

template <typename ValueType, typename BigValueType, typename RandomEngineType>
std::vector<ValueType> singleIterationPar(const std::vector<ValueType>& inData)
{
    const unsigned numberOfThreads = std::thread::hardware_concurrency();
    const std::size_t pieceSize = std::max<std::size_t>(inData.size()/numberOfThreads, 1);
    std::vector<Worker<ValueType, BigValueType, RandomEngineType> > workerVector;
    workerVector.reserve(numberOfThreads);
    std::vector<std::thread> threadVector;
    threadVector.reserve(numberOfThreads);
    std::size_t startPos = 0;
    std::size_t endLimit = pieceSize;
    unsigned procNum = 0;
    while(endLimit <= inData.size())
    {
        workerVector.push_back(Worker<ValueType, BigValueType, RandomEngineType>(inData,
                                                                                 startPos,
                                                                                 endLimit,
                                                                                 time(NULL)+startPos));
        threadVector.push_back(std::thread(std::ref(workerVector[procNum])));
        updatePosData(inData.size(), pieceSize, startPos, endLimit);
        procNum++;
    }
    for(auto& curThread : threadVector)
    {
        curThread.join();
    }
    std::size_t outputSize = 0;
    for(const auto& curWorker : workerVector)
    {
        outputSize += curWorker.getOutputSize();
    }
    std::vector<ValueType> outData;
    outData.reserve(outputSize);
    for(const auto& curWorker : workerVector)
    {
        for(const auto curNum : curWorker.getOutputData())
        {
            outData.push_back(curNum);
        }
    }
    return outData;
}

template<typename ValueType, typename BigValueType, typename RandomEngineType>
void runSimulation(const ValueType maxN, const std::size_t iterMax)
{
    std::vector<ValueType> numData;
    numData.reserve(maxN);
    for(ValueType i = 1; i <= maxN; ++i)
    {
        numData.push_back(i);
    }
    for(std::size_t iterNum = 0; iterNum < iterMax; ++iterNum)
    {
        auto startTime = std::chrono::steady_clock::now();
        numData = singleIterationPar<ValueType, BigValueType, RandomEngineType>(numData);
        auto endTime = std::chrono::steady_clock::now();
        std::chrono::duration<double> runTime = endTime-startTime;
        ValueType gcdVal = gcd<ValueType>(100*numData.size(), maxN);
        double resVal = static_cast<double>((100*numData.size())/gcdVal)/static_cast<double>(maxN/gcdVal);
        std::cout<<"Iteration: "<<iterNum<<":\t";
        std::cout<<"resVal: "<<resVal<<"%\t";
        std::cout<<"runtime: "<<runTime.count()<<" [s]"<<std::endl;
    }
}

int main()
{
    runSimulation<std::uintmax_t, __uint128_t, std::mt19937_64>(3000000, 10);
    return 0;
}

0 Helpful 0 Interesting 0 Brilliant 0 Confused

0 pending reports

×

Problem Loading...

Note Loading...

Set Loading...