Hello, Anton
Here is my solution of the SUM(1/N!) problem, as for the low accuracy of a 32bit double value, there is no need to calculate 1/N! when N >= 178.
Sample Code
1: //QUESTION
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum2"> 2:</span> <span style="color:#008000;">//Calculate the SUM(1/N!), this problem is from Anton's blog, see:</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum3"> 3:</span> <span style="color:#008000;">//http://jetcracker.wordpress.com/2012/03/07/mpi-calculating-sum/</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum4"> 4:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum5"> 5:</span> <span style="color:#cc6633;">#include</span> <stdio.h></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum6"> 6:</span> <span style="color:#cc6633;">#include</span> <omp.h> </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum7"> 7:</span> <span style="color:#cc6633;">#include</span> <ctime></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum8"> 8:</span> <span style="color:#cc6633;">#include</span> <cmath></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum9"> 9:</span> <span style="color:#cc6633;">#include</span> <cstdlib></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum10"> 10:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum11"> 11:</span> <span style="color:#008000;">//threads' number set to openmp.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum12"> 12:</span> <span style="color:#008000;">//#define NUMBER_THREADS 4</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum13"> 13:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum14"> 14:</span> <span style="color:#008000;">//*FACTORIAL*</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum15"> 15:</span> <span style="color:#008000;">//Theory 1: While N >= 178, the value of 1/N! would be nearly to 0 and won't be </span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum16"> 16:</span> <span style="color:#008000;">// saved as 'double' in 32bit machine, so it is actually a 1ms task.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum17"> 17:</span> <span style="color:#008000;">//Theory 2: 1/(N-1)! can be used to calculate 1/N! by multiply 1/N.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum18"> 18:</span> <span style="color:#0000ff;">double</span> factorial()</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum19"> 19:</span> {</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum20"> 20:</span> <span style="color:#0000ff;">int</span> n = 177;</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum21"> 21:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum22"> 22:</span> <span style="color:#0000ff;">double</span> factor = 1;</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum23"> 23:</span> <span style="color:#0000ff;">double</span> retval = 1;</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum24"> 24:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum25"> 25:</span> <span style="color:#008000;">//#pragma omp parallel for</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum26"> 26:</span> <span style="color:#0000ff;">for</span> (<span style="color:#0000ff;">int</span> i = 2; i <= n; i++)</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum27"> 27:</span> {</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum28"> 28:</span> factor *= (1.0 / (<span style="color:#0000ff;">double</span>)i);</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum29"> 29:</span> retval += factor;</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum30"> 30:</span> }</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum31"> 31:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum32"> 32:</span> <span style="color:#0000ff;">return</span> retval;</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum33"> 33:</span> }</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum34"> 34:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum35"> 35:</span> <span style="color:#008000;">//*TEST MAIN*</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum36"> 36:</span> <span style="color:#008000;">//Caculate the result, and then print out to console.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum37"> 37:</span> <span style="color:#0000ff;">void</span> main()</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum38"> 38:</span> { </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum39"> 39:</span> clock_t start = clock();</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum40"> 40:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum41"> 41:</span> <span style="color:#008000;">//setup the threads number.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum42"> 42:</span> <span style="color:#008000;">//omp_set_num_threads(NUMBER_THREADS);</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum43"> 43:</span> </pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum44"> 44:</span> <span style="color:#008000;">//calculate the result of sum.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum45"> 45:</span> printf(<span style="color:#006080;">"SUM(1/N!): %.15fn"</span>, factorial());</pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum46"> 46:</span>  </pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum47"> 47:</span> <span style="color:#008000;">//print the executing period.</span></pre>
<pre style="text-align:left;line-height:12pt;background-color:#f4f4f4;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum48"> 48:</span> printf(<span style="color:#006080;">"Total execute: %ldmsn"</span>, clock() - start);</pre>
<pre style="text-align:left;line-height:12pt;background-color:white;width:100%;font-family:'direction:ltr;color:black;font-size:8pt;overflow:visible;border-style:none;margin:0;padding:0;"><span style="color:#606060;" id="lnum49"> 49:</span> }</pre>
And here is the result of execution:
SUM(1/N!): 1.718281828459046
Total execute: 1ms
If you use a hi-accuracy float type, then you might need to parallelize the code, but as for ‘double’ type here, such algorithm should be fast enough.
Thanks!
近期评论