Question: High Precision Computing in Python or R - Data Science Central2020-09-18T19:54:24Zhttps://www.datasciencecentral.com/forum/topics/question-how-precision-computing-in-python?utm_content=bufferf2696&utm_medium=social&utm_source=twitter.com&utm_campaign=buffer&feed=yes&xn_auth=noHere is code in R. I understa…tag:www.datasciencecentral.com,2020-04-14:6448529:Comment:9447912020-04-14T06:32:05.546Zvictor zurkowskihttps://www.datasciencecentral.com/profile/victorzurkowski
<p>Here is code in R. I understand that the initial value of X is X[1], and the code outputs the first few digits of X[10,000], starting with initial value 0.3. </p>
<p></p>
<p>According to a previous post, the values obtained form a (non-independent) sample drawn from the density</p>
<p>f(u) = 1/ ( 2* sqrt(1-u) ) .</p>
<p>The histogram of {X[1], ..., X[10001]} looks like a sample drawn from f.</p>
<p></p>
<p>#------------------------------</p>
<p>require("Rmpfr")</p>
<p>#--- precision in bits…</p>
<p>Here is code in R. I understand that the initial value of X is X[1], and the code outputs the first few digits of X[10,000], starting with initial value 0.3. </p>
<p></p>
<p>According to a previous post, the values obtained form a (non-independent) sample drawn from the density</p>
<p>f(u) = 1/ ( 2* sqrt(1-u) ) .</p>
<p>The histogram of {X[1], ..., X[10001]} looks like a sample drawn from f.</p>
<p></p>
<p>#------------------------------</p>
<p>require("Rmpfr")</p>
<p>#--- precision in bits <br/>D = 20000</p>
<p>#--- number of iterations<br/>N = 10001</p>
<p>#---- CHAOTIC process<br/>X=mpfr(rep(0,N), D )<br/>X[1]=.3</p>
<p>for(i in 2:N ){ X[i] = sqrt( 4*X[i-1]*(1-X[i-1]) ) }</p>
<p></p>
<p>#------------------ Terminal points<br/>Z = floor( .5 + 10^50*X[9999:10001] )/10^50</p>
<p>format(Z, digits= 10)<br/>roundMpfr( Z, precBits = 40)<br/># [1] "0.2929369280" "0.9102194993" "0.5717340724"<br/># 3 'mpfr' number of precision 40 bits <br/># [1] 0.2929369280459 0.9102194993147 0.5717340723877</p>
<p>#------------------------- Histogram and theoretical density</p>
<p>u=0:100/101</p>
<p>windows(8,8)</p>
<p>hist(as.numeric(X), col='blue', <br/> breaks=50, probability = TRUE,<br/> xlab = 'u' , ylab = 'y (Density)' ,<br/> main = 'Histogram of values and theoretical density')<br/>points(u, (1-u)^(-.5)/2, type='l', col='red', lwd=3)</p>
<p>text(.2, 1.5, expression(y == frac(1, 2*sqrt(1-u))))</p>
<p></p>
<p></p> R user here.
Browsing through…tag:www.datasciencecentral.com,2020-04-13:6448529:Comment:9450172020-04-13T23:42:32.391Zvictor zurkowskihttps://www.datasciencecentral.com/profile/victorzurkowski
<p>R user here.</p>
<p>Browsing through an old paper lead my to look at the chaotic dynamic of $x_i = 4 * x_{i-1}*( 1 - x_{i-1})$. The paper suggests looking at the cumulative sums of $x - .5$, i.e.: $P_i = P_{i-1} + x_{i-1} - .5$, so I did, but I computed $x-.5$ in two ways. One, by using the recursion above, two, by writing the recursion of $y_i = x_i - .5$ as $y_i = -.5 + 4*(.25 - y_i^2)$. The results are different, and by plotting the cumulated processes it is clear that after a few…</p>
<p>R user here.</p>
<p>Browsing through an old paper lead my to look at the chaotic dynamic of $x_i = 4 * x_{i-1}*( 1 - x_{i-1})$. The paper suggests looking at the cumulative sums of $x - .5$, i.e.: $P_i = P_{i-1} + x_{i-1} - .5$, so I did, but I computed $x-.5$ in two ways. One, by using the recursion above, two, by writing the recursion of $y_i = x_i - .5$ as $y_i = -.5 + 4*(.25 - y_i^2)$. The results are different, and by plotting the cumulated processes it is clear that after a few iterations the differences are noticeable. </p>
<p>R has the package Rmpfr to run computations using the MPFR C library.</p>
<p>Below is my R snippet. Varying the precision D one see the effect of precision on the computations. Taking D=1000 gives un-discernible results for both ways of computing the dynamics for 1000 iterations.</p>
<p></p>
<p>#--------------------------------------</p>
<p>require("Rmpfr")</p>
<p>#----------- precision in bits <br/>D = 100</p>
<p>#----------- number of iterations<br/>N=1000</p>
<p>#----------- quadratic dynamics parameter<br/>a=4</p>
<p>#----------- CHAOTIC process, and cumulative process<br/>X=mpfr(rep(0,N), D )<br/>X[1]=.3</p>
<p>P = mpfr(rep(0,N), D )<br/>P[1] = 1000</p>
<p>for(i in 2:N ){ X[i] = a*X[i-1]*(1-X[i-1]) <br/> P[i] = P[i-1] + X[i-1] - .5}</p>
<p><br/>#----------- CHAOTIC process centered, and cumulative process <br/>dP = mpfr(rep(0,N), D )<br/>dP[1] = X[1] - .5</p>
<p>P2 = mpfr(rep(0,N), D )<br/>P2[1] = 1000</p>
<p>for(i in 2:N ){ dP[i] = a*(.25 - dP[i-1]^2) -.5 <br/> P2[i] = P2[i-1] + dP[i-1]}</p>
<p><br/>#----------- PLOT<br/>windows(8,8)<br/>plot(1:1000, P, type='l', col='blue')<br/>points(1:1000, P2, type='l', col='red')</p>
<p>#--------------------------------------------</p> Using Python's mpmath (arbitr…tag:www.datasciencecentral.com,2017-12-30:6448529:Comment:6744182017-12-30T17:03:31.231ZEmanuel Woiskihttps://www.datasciencecentral.com/profile/EmanuelWoiski
<p>Using Python's mpmath (arbitrary precision) library, with the seed pi/11 and 10000,20000,30000,45000 and 50000 iterations, for 1000,10000 and 100000 decimal places precision. Notice that the decimal places of the difference between 10000 and 100000 decimal-place precision results depend almost linearly on the number of iterations. As the number of iterations grow, precision is quickly destroyed.</p>
<p>Using Python's mpmath (arbitrary precision) library, with the seed pi/11 and 10000,20000,30000,45000 and 50000 iterations, for 1000,10000 and 100000 decimal places precision. Notice that the decimal places of the difference between 10000 and 100000 decimal-place precision results depend almost linearly on the number of iterations. As the number of iterations grow, precision is quickly destroyed.</p> Using the Apfloat Java JAR li…tag:www.datasciencecentral.com,2017-11-20:6448529:Comment:6514172017-11-20T04:26:51.169ZMatthew Boernerhttps://www.datasciencecentral.com/profile/MatthewBoerner
<p>Using the Apfloat Java JAR library available at <a href="http://www.apfloat.org/apfloat_java/" target="_blank">http://www.apfloat.org/apfloat_java/</a>, you can generate very high precision numbers in Java code (hundreds or thousands of digits or any arbitrary number of digits, computation speed will degrade as you need more precision obviously). Below is my Java code generating numbers using the well known iterative function x(k+1) = 4*x(k)*(1-x(k)). You can compare numbers using this…</p>
<p>Using the Apfloat Java JAR library available at <a href="http://www.apfloat.org/apfloat_java/" target="_blank">http://www.apfloat.org/apfloat_java/</a>, you can generate very high precision numbers in Java code (hundreds or thousands of digits or any arbitrary number of digits, computation speed will degrade as you need more precision obviously). Below is my Java code generating numbers using the well known iterative function x(k+1) = 4*x(k)*(1-x(k)). You can compare numbers using this iteration vs. the actual known kth-element solution of pseudo-code formula of sin((2^k)*init*PI)^2. The caret represents exponentiation. Interesting that the sine function degrades our precision in the estimated and actual computation since we are taking the sine of a larger and larger number for each iteration (due to the 2^k part in the formula). Here's my Java code:</p>
<pre><span>public class </span>Main {<br/> <span>public static void </span>main(String[] args) {<br/> <span>long </span>k;<br/> String zp;<br/> <span>long </span>prec = <span>30</span>;<br/> <span>long </span>highprec = <span>1000</span>;<br/> Apfloat seed = ApfloatMath.<span>pi</span>(prec).divide(<span>new </span>Apfloat(<span>11</span>,prec));<br/> Apfloat estimate = seed;<br/> Apfloat init = <span>new </span>Apfloat(<span>1</span>).divide(ApfloatMath.<span>pi</span>(highprec)).multiply(ApfloatMath.<span>asin</span>(ApfloatMath.<span>sqrt</span>(ApfloatMath.<span>pi</span>(highprec).divide(<span>new </span>Apfloat(<span>11</span>)))));<br/><br/> <span>for </span>(k=<span>1</span>; k<<span>100</span>; k++) {<br/> estimate = <span>new </span>Apfloat(<span>4</span>).multiply(estimate).multiply(<span>new </span>Apfloat(<span>1</span>).subtract(estimate));<br/> Apfloat actual = ApfloatMath.<span>pow</span>(ApfloatMath.<span>sin</span>(init.multiply(ApfloatMath.<span>pi</span>(highprec).multiply(ApfloatMath.<span>pow</span>(<span>new </span>Apfloat(<span>2</span>),<span>new </span>Apfloat(k,highprec))))),<span>2</span>);<br/><br/> System.<span>out</span>.println(k);<br/> System.<span>out</span>.println(String.<span>format</span>(<span>"%#s"</span>, estimate));<br/> System.<span>out</span>.println(String.<span>format</span>(<span>"%#s"</span>, actual));<br/> }<br/> }<br/>}</pre> Here's a Java version. You s…tag:www.datasciencecentral.com,2017-11-15:6448529:Comment:6493852017-11-15T03:19:42.499ZMatthew Boernerhttps://www.datasciencecentral.com/profile/MatthewBoerner
<p>Here's a Java version. You should find your own implementation of bigSqrt. This Java version diverges from the correct answer at iteration #30. I was hoping it would be better. Unfortunately Java does not have a built-in version of square root for BigDecimal.</p>
<p></p>
<pre><span>public static void </span>main(String[] args) {<br></br> BigDecimal pi = <span>new </span>BigDecimal(<span>4 </span>* <span>atan2</span>(<span>1</span>,<span>1</span>));<br></br> BigInteger bigint = <span>new…</span></pre>
<p>Here's a Java version. You should find your own implementation of bigSqrt. This Java version diverges from the correct answer at iteration #30. I was hoping it would be better. Unfortunately Java does not have a built-in version of square root for BigDecimal.</p>
<p></p>
<pre><span>public static void </span>main(String[] args) {<br/> BigDecimal pi = <span>new </span>BigDecimal(<span>4 </span>* <span>atan2</span>(<span>1</span>,<span>1</span>));<br/> BigInteger bigint = <span>new </span>BigInteger(<span>"10"</span>).pow(<span>50</span>);<br/> BigDecimal seed = pi.divide(<span>new </span>BigDecimal(<span>11</span>), MathContext.<span>DECIMAL128</span>);<br/> BigDecimal z = seed;<br/> <span>long </span>k;<br/> String zp;<br/><br/> <span>for </span>(k=<span>1</span>; k<<span>10000</span>; k++) {<br/> z = z.multiply(<span>new </span>BigDecimal(<span>"4"</span>)).multiply(<span>new </span>BigDecimal(<span>"1"</span>).subtract(z));<br/> z = <span>bigSqrt</span>(z, MathContext.<span>DECIMAL128</span>);<br/> z = <span>new </span>BigDecimal(bigint).multiply(z).add(<span>new </span>BigDecimal(<span>".5"</span>)).setScale(<span>0</span>, RoundingMode.<span>CEILING</span>).divide(<span>new </span>BigDecimal(bigint));<br/> zp = String.<span>format</span>(<span>"%.10f"</span>, z);<br/> System.<span>out</span>.println(k);<br/> System.<span>out</span>.println(zp);<br/> }<br/>}</pre> As for Python, someone on Red…tag:www.datasciencecentral.com,2017-11-13:6448529:Comment:6492122017-11-13T21:51:58.413ZVincent Granvillehttps://www.datasciencecentral.com/profile/VincentGranville
<p>As for Python, <a href="https://www.reddit.com/r/Python/comments/7ccdye/question_high_precision_computing_in_python/" target="_blank">someone on Reddit posted the following</a>, with seed = 3/10:</p>
<p><a href="https://storage.ning.com/topology/rest/1.0/file/get/2773309020?profile=original" target="_self"><img class="align-center" src="https://storage.ning.com/topology/rest/1.0/file/get/2773309020?profile=RESIZE_320x320" width="311"></img></a></p>
<p>Now, if I run that I get:</p>
<p><a href="https://storage.ning.com/topology/rest/1.0/file/get/2773309041?profile=original" target="_self"><img class="align-center" src="https://storage.ning.com/topology/rest/1.0/file/get/2773309041?profile=original" width="185"></img></a></p>
<p>With…</p>
<p>As for Python, <a href="https://www.reddit.com/r/Python/comments/7ccdye/question_high_precision_computing_in_python/" target="_blank">someone on Reddit posted the following</a>, with seed = 3/10:</p>
<p><a href="https://storage.ning.com/topology/rest/1.0/file/get/2773309020?profile=original" target="_self"><img width="311" src="https://storage.ning.com/topology/rest/1.0/file/get/2773309020?profile=RESIZE_320x320" width="311" class="align-center"/></a></p>
<p>Now, if I run that I get:</p>
<p><a href="https://storage.ning.com/topology/rest/1.0/file/get/2773309041?profile=original" target="_self"><img src="https://storage.ning.com/topology/rest/1.0/file/get/2773309041?profile=original" width="185" class="align-center"/></a></p>
<p>With fixed being the result using float, arbitrary being decimal.Decimal, and diff being what I believe represents the accumulated rounding error. That said I can't really run your Perl, and I'm not sure I haven't missed something transcribing your code, and since we're starting with different seeds the values in your spreadsheet bear no relation.</p>
<p><em>My comment is that over the 10,000 iterations, the difference at any given iteration, is much bigger than 0.1229, on average. The difference was computed only at iteration 10,000 in the above Python code. </em></p> Hi David,
My goal here is to…tag:www.datasciencecentral.com,2017-11-13:6448529:Comment:6488992017-11-13T18:51:56.313ZVincent Granvillehttps://www.datasciencecentral.com/profile/VincentGranville
<p>Hi David,</p>
<p>My goal here is to illustrate high precision computing, in a context where many scientists are not aware that their simulations are wrong after just 50 iterations (though, due to the nature of the process, it does not really matter.) </p>
<p>I am aware of the Mersenne twister, and yes it is very good for random number generation. Would love to see how R and Python compare in terms of high precision computing, regardless of the application. Yep, C has nice libraries too for…</p>
<p>Hi David,</p>
<p>My goal here is to illustrate high precision computing, in a context where many scientists are not aware that their simulations are wrong after just 50 iterations (though, due to the nature of the process, it does not really matter.) </p>
<p>I am aware of the Mersenne twister, and yes it is very good for random number generation. Would love to see how R and Python compare in terms of high precision computing, regardless of the application. Yep, C has nice libraries too for that purpose.</p>
<p>Thanks,</p>
<p>Vincent</p> I don't see any way to do it…tag:www.datasciencecentral.com,2017-11-13:6448529:Comment:6487482017-11-13T02:18:02.492ZDavid Reinkehttps://www.datasciencecentral.com/profile/DavidReinke826
<p>I don't see any way to do it in R or Python in a way that would execute quickly enough. Two possible alternatives:</p>
<ol>
<li>Write your routine in C# using the BigInteger class, which has pretty much everything you need to do math on integers of unlimited size. I haven't tried it, but I understand that it is now possible to link C# routines in .NET to R.</li>
<li>Write from scratch in C and link to R. There is a mechanism for doing this already in R. Knuth's TAOCP ch 4 goes through the…</li>
</ol>
<p>I don't see any way to do it in R or Python in a way that would execute quickly enough. Two possible alternatives:</p>
<ol>
<li>Write your routine in C# using the BigInteger class, which has pretty much everything you need to do math on integers of unlimited size. I haven't tried it, but I understand that it is now possible to link C# routines in .NET to R.</li>
<li>Write from scratch in C and link to R. There is a mechanism for doing this already in R. Knuth's TAOCP ch 4 goes through the classical algorithms for arithmetic on big integers.</li>
</ol>
<p>Why do you need such big integers? I find that R's standard RNG -- a Mersenne twister -- or L'Ecuyer's WELL package in R work for me.</p>