Last active
December 21, 2015 05:38
-
-
Save stla/6258179 to your computer and use it in GitHub Desktop.
Knitr output with a too long character string
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| <!DOCTYPE html> | |
| <!-- saved from url=(0014)about:internet --> | |
| <html> | |
| <head> | |
| <meta http-equiv="Content-Type" content="text/html; charset=utf-8"/> | |
| <title>Introductory example: Euler's approximation of \( \pi \)</title> | |
| <style type="text/css"> | |
| body, td { | |
| font-family: sans-serif; | |
| background-color: white; | |
| font-size: 12px; | |
| margin: 8px; | |
| } | |
| tt, code, pre { | |
| font-family: 'DejaVu Sans Mono', 'Droid Sans Mono', 'Lucida Console', Consolas, Monaco, monospace; | |
| } | |
| h1 { | |
| font-size:2.2em; | |
| } | |
| h2 { | |
| font-size:1.8em; | |
| } | |
| h3 { | |
| font-size:1.4em; | |
| } | |
| h4 { | |
| font-size:1.0em; | |
| } | |
| h5 { | |
| font-size:0.9em; | |
| } | |
| h6 { | |
| font-size:0.8em; | |
| } | |
| a:visited { | |
| color: rgb(50%, 0%, 50%); | |
| } | |
| pre { | |
| margin-top: 0; | |
| max-width: 95%; | |
| border: 1px solid #ccc; | |
| white-space: pre-wrap; | |
| } | |
| pre code { | |
| display: block; padding: 0.5em; | |
| word-wrap: break-word; | |
| } | |
| code.r, code.cpp { | |
| background-color: #F8F8F8; | |
| } | |
| table, td, th { | |
| border: none; | |
| } | |
| blockquote { | |
| color:#666666; | |
| margin:0; | |
| padding-left: 1em; | |
| border-left: 0.5em #EEE solid; | |
| } | |
| hr { | |
| height: 0px; | |
| border-bottom: none; | |
| border-top-width: thin; | |
| border-top-style: dotted; | |
| border-top-color: #999999; | |
| } | |
| @media print { | |
| * { | |
| background: transparent !important; | |
| color: black !important; | |
| filter:none !important; | |
| -ms-filter: none !important; | |
| } | |
| body { | |
| font-size:12pt; | |
| max-width:100%; | |
| } | |
| a, a:visited { | |
| text-decoration: underline; | |
| } | |
| hr { | |
| visibility: hidden; | |
| page-break-before: always; | |
| } | |
| pre, blockquote { | |
| padding-right: 1em; | |
| page-break-inside: avoid; | |
| } | |
| tr, img { | |
| page-break-inside: avoid; | |
| } | |
| img { | |
| max-width: 100% !important; | |
| } | |
| @page :left { | |
| margin: 15mm 20mm 15mm 10mm; | |
| } | |
| @page :right { | |
| margin: 15mm 10mm 15mm 20mm; | |
| } | |
| p, h2, h3 { | |
| orphans: 3; widows: 3; | |
| } | |
| h2, h3 { | |
| page-break-after: avoid; | |
| } | |
| } | |
| </style> | |
| <!-- Styles for R syntax highlighter --> | |
| <style type="text/css"> | |
| pre .operator, | |
| pre .paren { | |
| color: rgb(104, 118, 135) | |
| } | |
| pre .literal { | |
| color: rgb(88, 72, 246) | |
| } | |
| pre .number { | |
| color: rgb(0, 0, 205); | |
| } | |
| pre .comment { | |
| color: rgb(76, 136, 107); | |
| } | |
| pre .keyword { | |
| color: rgb(0, 0, 255); | |
| } | |
| pre .identifier { | |
| color: rgb(0, 0, 0); | |
| } | |
| pre .string { | |
| color: rgb(3, 106, 7); | |
| } | |
| </style> | |
| <!-- R syntax highlighter --> | |
| <script type="text/javascript"> | |
| var hljs=new function(){function m(p){return p.replace(/&/gm,"&").replace(/</gm,"<")}function f(r,q,p){return RegExp(q,"m"+(r.cI?"i":"")+(p?"g":""))}function b(r){for(var p=0;p<r.childNodes.length;p++){var q=r.childNodes[p];if(q.nodeName=="CODE"){return q}if(!(q.nodeType==3&&q.nodeValue.match(/\s+/))){break}}}function h(t,s){var p="";for(var r=0;r<t.childNodes.length;r++){if(t.childNodes[r].nodeType==3){var q=t.childNodes[r].nodeValue;if(s){q=q.replace(/\n/g,"")}p+=q}else{if(t.childNodes[r].nodeName=="BR"){p+="\n"}else{p+=h(t.childNodes[r])}}}if(/MSIE [678]/.test(navigator.userAgent)){p=p.replace(/\r/g,"\n")}return p}function a(s){var r=s.className.split(/\s+/);r=r.concat(s.parentNode.className.split(/\s+/));for(var q=0;q<r.length;q++){var p=r[q].replace(/^language-/,"");if(e[p]){return p}}}function c(q){var p=[];(function(s,t){for(var r=0;r<s.childNodes.length;r++){if(s.childNodes[r].nodeType==3){t+=s.childNodes[r].nodeValue.length}else{if(s.childNodes[r].nodeName=="BR"){t+=1}else{if(s.childNodes[r].nodeType==1){p.push({event:"start",offset:t,node:s.childNodes[r]});t=arguments.callee(s.childNodes[r],t);p.push({event:"stop",offset:t,node:s.childNodes[r]})}}}}return t})(q,0);return p}function k(y,w,x){var q=0;var z="";var s=[];function u(){if(y.length&&w.length){if(y[0].offset!=w[0].offset){return(y[0].offset<w[0].offset)?y:w}else{return w[0].event=="start"?y:w}}else{return y.length?y:w}}function t(D){var A="<"+D.nodeName.toLowerCase();for(var B=0;B<D.attributes.length;B++){var C=D.attributes[B];A+=" "+C.nodeName.toLowerCase();if(C.value!==undefined&&C.value!==false&&C.value!==null){A+='="'+m(C.value)+'"'}}return A+">"}while(y.length||w.length){var v=u().splice(0,1)[0];z+=m(x.substr(q,v.offset-q));q=v.offset;if(v.event=="start"){z+=t(v.node);s.push(v.node)}else{if(v.event=="stop"){var p,r=s.length;do{r--;p=s[r];z+=("</"+p.nodeName.toLowerCase()+">")}while(p!=v.node);s.splice(r,1);while(r<s.length){z+=t(s[r]);r++}}}}return z+m(x.substr(q))}function j(){function q(x,y,v){if(x.compiled){return}var u;var s=[];if(x.k){x.lR=f(y,x.l||hljs.IR,true);for(var w in x.k){if(!x.k.hasOwnProperty(w)){continue}if(x.k[w] instanceof Object){u=x.k[w]}else{u=x.k;w="keyword"}for(var r in u){if(!u.hasOwnProperty(r)){continue}x.k[r]=[w,u[r]];s.push(r)}}}if(!v){if(x.bWK){x.b="\\b("+s.join("|")+")\\s"}x.bR=f(y,x.b?x.b:"\\B|\\b");if(!x.e&&!x.eW){x.e="\\B|\\b"}if(x.e){x.eR=f(y,x.e)}}if(x.i){x.iR=f(y,x.i)}if(x.r===undefined){x.r=1}if(!x.c){x.c=[]}x.compiled=true;for(var t=0;t<x.c.length;t++){if(x.c[t]=="self"){x.c[t]=x}q(x.c[t],y,false)}if(x.starts){q(x.starts,y,false)}}for(var p in e){if(!e.hasOwnProperty(p)){continue}q(e[p].dM,e[p],true)}}function d(B,C){if(!j.called){j();j.called=true}function q(r,M){for(var L=0;L<M.c.length;L++){if((M.c[L].bR.exec(r)||[null])[0]==r){return M.c[L]}}}function v(L,r){if(D[L].e&&D[L].eR.test(r)){return 1}if(D[L].eW){var M=v(L-1,r);return M?M+1:0}return 0}function w(r,L){return L.i&&L.iR.test(r)}function K(N,O){var M=[];for(var L=0;L<N.c.length;L++){M.push(N.c[L].b)}var r=D.length-1;do{if(D[r].e){M.push(D[r].e)}r--}while(D[r+1].eW);if(N.i){M.push(N.i)}return f(O,M.join("|"),true)}function p(M,L){var N=D[D.length-1];if(!N.t){N.t=K(N,E)}N.t.lastIndex=L;var r=N.t.exec(M);return r?[M.substr(L,r.index-L),r[0],false]:[M.substr(L),"",true]}function z(N,r){var L=E.cI?r[0].toLowerCase():r[0];var M=N.k[L];if(M&&M instanceof Array){return M}return false}function F(L,P){L=m(L);if(!P.k){return L}var r="";var O=0;P.lR.lastIndex=0;var M=P.lR.exec(L);while(M){r+=L.substr(O,M.index-O);var N=z(P,M);if(N){x+=N[1];r+='<span class="'+N[0]+'">'+M[0]+"</span>"}else{r+=M[0]}O=P.lR.lastIndex;M=P.lR.exec(L)}return r+L.substr(O,L.length-O)}function J(L,M){if(M.sL&&e[M.sL]){var r=d(M.sL,L);x+=r.keyword_count;return r.value}else{return F(L,M)}}function I(M,r){var L=M.cN?'<span class="'+M.cN+'">':"";if(M.rB){y+=L;M.buffer=""}else{if(M.eB){y+=m(r)+L;M.buffer=""}else{y+=L;M.buffer=r}}D.push(M);A+=M.r}function G(N,M,Q){var R=D[D.length-1];if(Q){y+=J(R.buffer+N,R);return false}var P=q(M,R);if(P){y+=J(R.buffer+N,R);I(P,M);return P.rB}var L=v(D.length-1,M);if(L){var O=R.cN?"</span>":"";if(R.rE){y+=J(R.buffer+N,R)+O}else{if(R.eE){y+=J(R.buffer+N,R)+O+m(M)}else{y+=J(R.buffer+N+M,R)+O}}while(L>1){O=D[D.length-2].cN?"</span>":"";y+=O;L--;D.length--}var r=D[D.length-1];D.length--;D[D.length-1].buffer="";if(r.starts){I(r.starts,"")}return R.rE}if(w(M,R)){throw"Illegal"}}var E=e[B];var D=[E.dM];var A=0;var x=0;var y="";try{var s,u=0;E.dM.buffer="";do{s=p(C,u);var t=G(s[0],s[1],s[2]);u+=s[0].length;if(!t){u+=s[1].length}}while(!s[2]);if(D.length>1){throw"Illegal"}return{r:A,keyword_count:x,value:y}}catch(H){if(H=="Illegal"){return{r:0,keyword_count:0,value:m(C)}}else{throw H}}}function g(t){var p={keyword_count:0,r:0,value:m(t)};var r=p;for(var q in e){if(!e.hasOwnProperty(q)){continue}var s=d(q,t);s.language=q;if(s.keyword_count+s.r>r.keyword_count+r.r){r=s}if(s.keyword_count+s.r>p.keyword_count+p.r){r=p;p=s}}if(r.language){p.second_best=r}return p}function i(r,q,p){if(q){r=r.replace(/^((<[^>]+>|\t)+)/gm,function(t,w,v,u){return w.replace(/\t/g,q)})}if(p){r=r.replace(/\n/g,"<br>")}return r}function n(t,w,r){var x=h(t,r);var v=a(t);var y,s;if(v){y=d(v,x)}else{return}var q=c(t);if(q.length){s=document.createElement("pre");s.innerHTML=y.value;y.value=k(q,c(s),x)}y.value=i(y.value,w,r);var u=t.className;if(!u.match("(\\s|^)(language-)?"+v+"(\\s|$)")){u=u?(u+" "+v):v}if(/MSIE [678]/.test(navigator.userAgent)&&t.tagName=="CODE"&&t.parentNode.tagName=="PRE"){s=t.parentNode;var p=document.createElement("div");p.innerHTML="<pre><code>"+y.value+"</code></pre>";t=p.firstChild.firstChild;p.firstChild.cN=s.cN;s.parentNode.replaceChild(p.firstChild,s)}else{t.innerHTML=y.value}t.className=u;t.result={language:v,kw:y.keyword_count,re:y.r};if(y.second_best){t.second_best={language:y.second_best.language,kw:y.second_best.keyword_count,re:y.second_best.r}}}function o(){if(o.called){return}o.called=true;var r=document.getElementsByTagName("pre");for(var p=0;p<r.length;p++){var q=b(r[p]);if(q){n(q,hljs.tabReplace)}}}function l(){if(window.addEventListener){window.addEventListener("DOMContentLoaded",o,false);window.addEventListener("load",o,false)}else{if(window.attachEvent){window.attachEvent("onload",o)}else{window.onload=o}}}var e={};this.LANGUAGES=e;this.highlight=d;this.highlightAuto=g;this.fixMarkup=i;this.highlightBlock=n;this.initHighlighting=o;this.initHighlightingOnLoad=l;this.IR="[a-zA-Z][a-zA-Z0-9_]*";this.UIR="[a-zA-Z_][a-zA-Z0-9_]*";this.NR="\\b\\d+(\\.\\d+)?";this.CNR="\\b(0[xX][a-fA-F0-9]+|(\\d+(\\.\\d*)?|\\.\\d+)([eE][-+]?\\d+)?)";this.BNR="\\b(0b[01]+)";this.RSR="!|!=|!==|%|%=|&|&&|&=|\\*|\\*=|\\+|\\+=|,|\\.|-|-=|/|/=|:|;|<|<<|<<=|<=|=|==|===|>|>=|>>|>>=|>>>|>>>=|\\?|\\[|\\{|\\(|\\^|\\^=|\\||\\|=|\\|\\||~";this.ER="(?![\\s\\S])";this.BE={b:"\\\\.",r:0};this.ASM={cN:"string",b:"'",e:"'",i:"\\n",c:[this.BE],r:0};this.QSM={cN:"string",b:'"',e:'"',i:"\\n",c:[this.BE],r:0};this.CLCM={cN:"comment",b:"//",e:"$"};this.CBLCLM={cN:"comment",b:"/\\*",e:"\\*/"};this.HCM={cN:"comment",b:"#",e:"$"};this.NM={cN:"number",b:this.NR,r:0};this.CNM={cN:"number",b:this.CNR,r:0};this.BNM={cN:"number",b:this.BNR,r:0};this.inherit=function(r,s){var p={};for(var q in r){p[q]=r[q]}if(s){for(var q in s){p[q]=s[q]}}return p}}();hljs.LANGUAGES.cpp=function(){var a={keyword:{"false":1,"int":1,"float":1,"while":1,"private":1,"char":1,"catch":1,"export":1,virtual:1,operator:2,sizeof:2,dynamic_cast:2,typedef:2,const_cast:2,"const":1,struct:1,"for":1,static_cast:2,union:1,namespace:1,unsigned:1,"long":1,"throw":1,"volatile":2,"static":1,"protected":1,bool:1,template:1,mutable:1,"if":1,"public":1,friend:2,"do":1,"return":1,"goto":1,auto:1,"void":2,"enum":1,"else":1,"break":1,"new":1,extern:1,using:1,"true":1,"class":1,asm:1,"case":1,typeid:1,"short":1,reinterpret_cast:2,"default":1,"double":1,register:1,explicit:1,signed:1,typename:1,"try":1,"this":1,"switch":1,"continue":1,wchar_t:1,inline:1,"delete":1,alignof:1,char16_t:1,char32_t:1,constexpr:1,decltype:1,noexcept:1,nullptr:1,static_assert:1,thread_local:1,restrict:1,_Bool:1,complex:1},built_in:{std:1,string:1,cin:1,cout:1,cerr:1,clog:1,stringstream:1,istringstream:1,ostringstream:1,auto_ptr:1,deque:1,list:1,queue:1,stack:1,vector:1,map:1,set:1,bitset:1,multiset:1,multimap:1,unordered_set:1,unordered_map:1,unordered_multiset:1,unordered_multimap:1,array:1,shared_ptr:1}};return{dM:{k:a,i:"</",c:[hljs.CLCM,hljs.CBLCLM,hljs.QSM,{cN:"string",b:"'\\\\?.",e:"'",i:"."},{cN:"number",b:"\\b(\\d+(\\.\\d*)?|\\.\\d+)(u|U|l|L|ul|UL|f|F)"},hljs.CNM,{cN:"preprocessor",b:"#",e:"$"},{cN:"stl_container",b:"\\b(deque|list|queue|stack|vector|map|set|bitset|multiset|multimap|unordered_map|unordered_set|unordered_multiset|unordered_multimap|array)\\s*<",e:">",k:a,r:10,c:["self"]}]}}}();hljs.LANGUAGES.r={dM:{c:[hljs.HCM,{cN:"number",b:"\\b0[xX][0-9a-fA-F]+[Li]?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+(?:[eE][+\\-]?\\d*)?L\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\b\\d+\\.(?!\\d)(?:i\\b)?",e:hljs.IMMEDIATE_RE,r:1},{cN:"number",b:"\\b\\d+(?:\\.\\d*)?(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"number",b:"\\.\\d+(?:[eE][+\\-]?\\d*)?i?\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"keyword",b:"(?:tryCatch|library|setGeneric|setGroupGeneric)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\.",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\.\\.\\d+(?![\\w.])",e:hljs.IMMEDIATE_RE,r:10},{cN:"keyword",b:"\\b(?:function)",e:hljs.IMMEDIATE_RE,r:2},{cN:"keyword",b:"(?:if|in|break|next|repeat|else|for|return|switch|while|try|stop|warning|require|attach|detach|source|setMethod|setClass)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"literal",b:"(?:NA|NA_integer_|NA_real_|NA_character_|NA_complex_)\\b",e:hljs.IMMEDIATE_RE,r:10},{cN:"literal",b:"(?:NULL|TRUE|FALSE|T|F|Inf|NaN)\\b",e:hljs.IMMEDIATE_RE,r:1},{cN:"identifier",b:"[a-zA-Z.][a-zA-Z0-9._]*\\b",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"<\\-(?!\\s*\\d)",e:hljs.IMMEDIATE_RE,r:2},{cN:"operator",b:"\\->|<\\-",e:hljs.IMMEDIATE_RE,r:1},{cN:"operator",b:"%%|~",e:hljs.IMMEDIATE_RE},{cN:"operator",b:">=|<=|==|!=|\\|\\||&&|=|\\+|\\-|\\*|/|\\^|>|<|!|&|\\||\\$|:",e:hljs.IMMEDIATE_RE,r:0},{cN:"operator",b:"%",e:"%",i:"\\n",r:1},{cN:"identifier",b:"`",e:"`",r:0},{cN:"string",b:'"',e:'"',c:[hljs.BE],r:0},{cN:"string",b:"'",e:"'",c:[hljs.BE],r:0},{cN:"paren",b:"[[({\\])}]",e:hljs.IMMEDIATE_RE,r:0}]}}; | |
| hljs.initHighlightingOnLoad(); | |
| </script> | |
| <!-- MathJax scripts --> | |
| <script type="text/javascript" src="https://c328740.ssl.cf1.rackcdn.com/mathjax/2.0-latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"> | |
| </script> | |
| </head> | |
| <body> | |
| <p>In this article you will firstly see how to get rational numbers arbitrary close to \( \pi \) | |
| by performing the <em>binary splitting algorithm</em> with the <code>gmp</code> package. </p> | |
| <p>The <em>binary splitting algorithm</em> fastly calculates the partial sums of a rational hypergeometric series by manipulating only integer numbers. But these integer numbers are generally gigantic hence they cannot be handled by ordinary arithmetic computing. After describing the binary splitting algorithm we will show how to implement it in R with the <code>gmp</code> package which allows <em>arithmetic without limitation</em>. Our main application is the evaluation of the Gauss hypergeometric function.</p> | |
| <h2>Introductory example: Euler's approximation of \( \pi \)</h2> | |
| <p>The following formula is due to Euler: | |
| \[ | |
| \frac{\pi}{2} = 1 + \frac{1}{3} + \frac{1\times 2}{3\times 5} + | |
| \frac{1\times 2 \times 3}{3\times 5 \times 7} + \cdots | |
| + \frac{n!}{3\times 5 \times 7 \times \cdots \times (2n+1)} + \cdots, | |
| \] | |
| that is, \( \pi = \lim S_n \) where | |
| \[ | |
| \begin{aligned} | |
| S_n & = 1 + \frac{u_1}{v_1} + \frac{u_1 u_2}{v_1v_2} + | |
| \frac{u_1u_2 u_3}{v_1v_2v_3} + \cdots + | |
| \frac{u_1u_2\ldots u_{n-1}u_n}{v_1v_2\ldots v_{n-1}v_n} \\ | |
| & = 1 + \sum_{k=1}^n \prod_{i=1}^k\frac{u_i}{v_i} \\ | |
| \end{aligned} | |
| \] | |
| with \( u_i=i \) and \( v_i=2i+1 \).</p> | |
| <p>Using new notations | |
| \( \alpha_i = \delta_i = u_i \) and \( \beta_i=v_i \) and writing | |
| \[ | |
| S_n -1 = \frac{\alpha_1}{\beta_1} + \frac{\delta_1 \alpha_2}{\beta_1\beta_2} + | |
| \frac{\delta_1\delta_2 \alpha_3}{\beta_1\beta_2\beta_3} + \cdots + | |
| \frac{\delta_1\delta_2\ldots\delta_{n-1}\alpha_n}{\beta_1\beta_2\ldots\beta_{n-1}\beta_n} | |
| \] | |
| could sound silly at first glance. But now assume \( \boxed{n=2^m} \). Then, by summing each \( (2i-1) \)-st term with the \( (2i) \)-th term, we can write \( S_n-1 \) as a sum of \( n/2 \) terms with a similar expression: | |
| \[ | |
| S_n - 1 = \frac{\alpha'_1}{\beta'_1} + \frac{\delta'_1 \alpha'_2}{\beta'_1\beta'_2} + | |
| \frac{\delta'_1\delta'_2 \alpha'_3}{\beta'_1\beta'_2\beta'_3} + \cdots + | |
| \frac{\delta'_1\delta'_2\ldots\delta'_{\frac{n}{2}-1}\alpha'_\frac{n}{2}}{\beta'_1\beta'_2\ldots\beta'_{\frac{n}{2}-1}\beta'_{\frac{n}{2}}} | |
| \] | |
| where \( \alpha'_i \), \( \delta'_i \) and \( \beta'_i \) are given by | |
| \[ | |
| \begin{aligned} | |
| \alpha'_i = \alpha_{2i-1}\beta_{2_i} + \alpha_{2i}\delta_{2i-1}, \quad | |
| \delta'_i = \delta_{2i-1}\delta_{2i} | |
| \qquad \text{and } \quad | |
| \beta'_i = \beta_{2i-1}\beta_{2i} | |
| \end{aligned} | |
| \] | |
| for all \( i \in \{1, \ldots, n/2\} \). </p> | |
| <p>Continuing so on, after \( m \) steps we obtain | |
| \[ | |
| S_n - 1 = \frac{\alpha^{(m)}}{\beta^{(m)}} | |
| \] | |
| where \( \alpha^{(m)} \) and \( \beta^{(m)} \) are integer numbers obtained by applying above formulas</p> | |
| <p>The above method is the <em>binary splitting algorithm</em> for evaluating \( S_n \) with \( n=2^m \), summarized as follows:</p> | |
| <ol> | |
| <li><p>Initialization: put \( \alpha^{(0)}_i = \delta^{(0)}_i = u_i \) and \( \beta^{(0)}_i=v_i \) for \( i \in \{1,n\} \);</p></li> | |
| <li><p>Compute recursively for \( k \) going from \( 1 \) to \( m \) | |
| \[ | |
| \begin{aligned} | |
| \alpha^{(k)}_i = \alpha^{(k-1)}_{2i-1}\beta^{(k-1)}_{2_i} + \alpha^{(k-1)}_{2i}\delta^{(k-1)}_{2i-1}, \quad | |
| \delta^{(k)}_i = \delta^{(k-1)}_{2i-1}\delta^{(k-1)}_{2i} | |
| \qquad \text{and } \quad | |
| \beta^{(k)}_i = \beta^{(k-1)}_{2i-1}\beta^{(k-1)}_{2i} | |
| \end{aligned} | |
| \] | |
| for \( i \in \{1,n/2^k\} \);</p></li> | |
| <li><p>Evaluate \( S_n = 1 + \frac{\alpha^{(m)}}{\beta^{(m)}} \).</p></li> | |
| </ol> | |
| <p>The advantage of the binary splitting as compared to a direct evaluation of \( S_n \) by summing its \( 2^m \) terms is twofold:</p> | |
| <ul> | |
| <li>the binary splitting only performs operations on integer numbers;</li> | |
| <li>it returns an exact expression of \( S_n \) as a ratio of two integer numbers.</li> | |
| </ul> | |
| <pre><code class="r">## example: rational approximation of pi ## | |
| bs.pi <- function(m) { | |
| u <- function(i) as.numeric(i) | |
| v <- function(i) 2 * i + 1 | |
| n <- 2^m | |
| indexes <- c(1:n) | |
| delta <- alpha <- u(indexes) | |
| beta <- v(indexes) | |
| j <- 1 | |
| l <- n | |
| while (j < n) { | |
| l <- l/2 | |
| odd <- 2 * c(1:l) | |
| even <- odd - 1 | |
| alpha <- beta[odd] * alpha[even] + delta[even] * alpha[odd] | |
| j <- 2 * j | |
| beta <- beta[odd] * beta[even] | |
| delta <- delta[even] * delta[odd] | |
| } | |
| Sn <- alpha/beta + 1 | |
| out <- list(alpha = alpha, beta = beta, Sn = Sn) | |
| return(out) | |
| } | |
| </code></pre> | |
| <p>The method very well performs while \( m\leq 7 \) :</p> | |
| <pre><code class="r">print(bs.pi(7), digits = 22) | |
| </code></pre> | |
| <pre><code>## $alpha | |
| ## [1] 9.589805429639700552931e+254 | |
| ## | |
| ## $beta | |
| ## [1] 1.680074832206408008955e+255 | |
| ## | |
| ## $Sn | |
| ## [1] 1.570796326794896557999 | |
| </code></pre> | |
| <pre><code class="r">print(pi/2, digits = 22) | |
| </code></pre> | |
| <pre><code>## [1] 1.570796326794896557999 | |
| </code></pre> | |
| <p>But the numerator and the denominator become too gigantic when \( m=8 \):</p> | |
| <pre><code class="r">bs.pi(8) | |
| </code></pre> | |
| <pre><code>## $alpha | |
| ## [1] Inf | |
| ## | |
| ## $beta | |
| ## [1] Inf | |
| ## | |
| ## $Sn | |
| ## [1] NaN | |
| </code></pre> | |
| <h2>Second example: exponential of a rational number</h2> | |
| <p>It is well known that \( \exp(x)=\lim S_n(x) \) where \( S_n(x)=\sum_{k=0}^n\frac{x^n}{n!} \). | |
| Thus, when \( x=p/q \) for some integers \( p \) and \( q \), we can write as before | |
| \[ | |
| S_n(x) = 1 + \sum_{k=1}^n \prod_{i=1}^k\frac{u_i}{v_i} | |
| \] | |
| where \( u_i \equiv p \) and \( v_i= i q \) are integer numbers. Thus, we can use the binary splitting algorithm to compute \( S_{2^m} \):</p> | |
| <pre><code class="r">## example: rational approximation of exp(p/q) ## | |
| bs.exp <- function(p, q, m) { | |
| v <- function(i) i * q | |
| n <- 2^m | |
| indexes <- 1:n | |
| delta <- alpha <- rep(p, n) | |
| beta <- v(indexes) | |
| j <- 1 | |
| l <- n | |
| while (j < n) { | |
| l <- l/2 | |
| odd <- 2 * c(1:l) | |
| even <- odd - 1 | |
| alpha <- beta[odd] * alpha[even] + delta[even] * alpha[odd] | |
| j <- 2 * j | |
| beta <- beta[odd] * beta[even] | |
| delta <- delta[even] * delta[odd] | |
| } | |
| Sn <- alpha/beta + 1 | |
| out <- list(alpha = alpha, beta = beta, Sn = Sn) | |
| return(out) | |
| } | |
| </code></pre> | |
| <p>Let us try to evaluate \( \exp(1) \). For \( m=7 \), the approximation is not entirely satisfactory:</p> | |
| <pre><code class="r">print(bs.exp(1, 1, 7), digits = 22) | |
| </code></pre> | |
| <pre><code>## $alpha | |
| ## [1] 6.626046675252336548741e+215 | |
| ## | |
| ## $beta | |
| ## [1] 3.85620482362580407204e+215 | |
| ## | |
| ## $Sn | |
| ## [1] 2.718281828459045534885 | |
| </code></pre> | |
| <pre><code class="r">print(exp(1), digits = 22) | |
| </code></pre> | |
| <pre><code>## [1] 2.718281828459045090796 | |
| </code></pre> | |
| <p>And for \( m=8 \), it crashes:</p> | |
| <pre><code class="r">bs.exp(1, 1, 8) | |
| </code></pre> | |
| <pre><code>## $alpha | |
| ## [1] Inf | |
| ## | |
| ## $beta | |
| ## [1] Inf | |
| ## | |
| ## $Sn | |
| ## [1] NaN | |
| </code></pre> | |
| <h2>The <code>gmp</code> package comes to our rescue</h2> | |
| <p>As we noted above, the binary splitting manipulates only <em>integer</em> numbers. | |
| The evaluation of \( \exp(1) \) has crashed because the numerator and the denominator were too big integers. | |
| The crantastic <a href="http://www.inside-r.org/packages/cran/gmp"><code>gmp</code></a> package overcomes this problem because it allows “arithmetic without limitations'' using the <a href="http://gmplib.org/">C library GMP (GNU Multiple Precision Arithmetic)</a>.</p> | |
| <p>Let us show how the <code>gmp</code> works on the \( \pi \) example. This is very easy: we only have to convert the two input sequences of integers \( (u_i) \) and \( (v_i) \) to sequences of <code>bigz</code> integers:</p> | |
| <pre><code class="r">library(gmp) | |
| ## rational approximation of pi with gmp ## | |
| bs.pi.gmp <- function(m) { | |
| u <- function(i) as.numeric(i) | |
| v <- function(i) 2 * i + 1 | |
| n <- 2^m | |
| indexes <- 1:n | |
| delta <- alpha <- as.bigz(u(indexes)) | |
| beta <- as.bigz(v(indexes)) | |
| j <- 1 | |
| l <- n | |
| while (j < n) { | |
| l <- l/2 | |
| odd <- 2 * c(1:l) | |
| even <- odd - 1 | |
| alpha <- beta[odd] * alpha[even] + delta[even] * alpha[odd] | |
| j <- 2 * j | |
| beta <- beta[odd] * beta[even] | |
| delta <- delta[even] * delta[odd] | |
| } | |
| Sn <- alpha/beta + 1 | |
| out <- list(Sn = Sn, eval.Sn = format(as.numeric(Sn), digits = 22)) | |
| return(out) | |
| } | |
| </code></pre> | |
| <p>The evaluation of \( S_n \) with \( n=2^3 \) illustrates the first advantage of the <code>gmp</code> package:</p> | |
| <pre><code class="r">bs.pi.gmp(3) | |
| </code></pre> | |
| <pre><code>## $Sn | |
| ## Big Rational ('bigq') : | |
| ## [1] 1202048/765765 | |
| ## | |
| ## $eval.Sn | |
| ## [1] "1.569734840323075530932" | |
| </code></pre> | |
| <pre><code class="r">bs.pi(3) | |
| </code></pre> | |
| <pre><code>## $alpha | |
| ## [1] 19632735 | |
| ## | |
| ## $beta | |
| ## [1] 34459425 | |
| ## | |
| ## $Sn | |
| ## [1] 1.57 | |
| </code></pre> | |
| <p>As you can see, \( S_n \) is written as an irreducible fraction with the <code>gmp</code> approach. | |
| But this is not the main strength of the <code>gmp</code> package. Now we have (almost) no limitation on \( m \) for evaluating \( S_{2^m} \):</p> | |
| <pre><code class="r">bs.pi.gmp(8) | |
| </code></pre> | |
| <pre><code>## $Sn | |
| ## Big Rational ('bigq') : | |
| ## [1] 115056663317199981372832786803399641133848259535718238578854114440177847232763528127119686643465544336537363974090559640151844992619459739337642897335661405374200830442503779326745081494631228217510085926896107230240702464/73247346810369298651903071099557979072216039642432949710389234675732768750102001285974817825809831148661290123993641325086924401900965008305646606428886048721946203288377842830920059623434101646117412656625454480462852875 | |
| ## | |
| ## $eval.Sn | |
| ## [1] "1.570796326794896557999" | |
| </code></pre> | |
| <p>Obviously the first limitation is the width of your screen. The more serious limitations of the <code>gmp</code> package are beyond the scope of this article. </p> | |
| <p>Let us come back to the exponential example:</p> | |
| <pre><code class="r">## rational approximation of exp(p/q) with gmp ## | |
| bs.exp.gmp <- function(p, q, m) { | |
| v <- function(i) i * q | |
| n <- 2^m | |
| indexes <- 1:n | |
| delta <- alpha <- as.bigz(rep(p, n)) | |
| beta <- as.bigz(v(indexes)) | |
| j <- 1 | |
| l <- n | |
| while (j < n) { | |
| l <- l/2 | |
| odd <- 2 * c(1:l) | |
| even <- odd - 1 | |
| alpha <- beta[odd] * alpha[even] + delta[even] * alpha[odd] | |
| j <- 2 * j | |
| beta <- beta[odd] * beta[even] | |
| delta <- delta[even] * delta[odd] | |
| } | |
| Sn <- alpha/beta + 1 | |
| out <- list(Sn = Sn, eval.Sn = format(as.numeric(Sn), digits = 22)) | |
| return(out) | |
| } | |
| </code></pre> | |
| <pre><code class="r">bs.exp.gmp(1, 1, 8) | |
| </code></pre> | |
| <pre><code>## $Sn | |
| ## Big Rational ('bigq') : | |
| ## [1] 63021364076854400517126597190157042974914655085470311494152999074896589361987361775329179623527760806690590676400388872831695705790559736341994225392293021235691155101792729596391087505487119686065032680426816409018591609682896947897581062232056198801713371950662092427153111247485380584396839593243205795931189046725531379112787311119506517584752693953099433873873085939642331053890371322719954788883613838912023544946108979472116077229049863887551154910123100635718060217444974605564852221865532212127661/23184264198455206868083304640033314193453554602148259996206909469655931150085069983174061928660848877037186090333421197463708022559289093927629440229660162856206414393604561795747978584507961086161320755987057927235191284503958147694842900705427915576370346458939828967066328925689811313743116731571304256245141968042147553432082017992236165926654195533967789698937870367867112218743295876678624370999142239502871990876622238944437605633097728000000000000000000000000000000000000000000000000000000000000000 | |
| ## | |
| ## $eval.Sn | |
| ## [1] "2.718281828459045090796" | |
| </code></pre> | |
| <p>Very well.</p> | |
| <h2>A general function for the binary splitting algorithm </h2> | |
| <p>Before turning to the Gauss hypergeometric function we write a general function for the binary splitting taking as arguments the two sequences \( (u_i) \) and \( (v_i) \):</p> | |
| <pre><code class="r">bs.gmp <- function(u, v, m = 7, value = "eval") { | |
| n <- 2^m | |
| indexes <- 1:n | |
| delta <- alpha <- as.bigz(u(indexes)) | |
| beta <- as.bigz(v(indexes)) | |
| j <- 1 | |
| l <- n | |
| while (j < n) { | |
| l <- l/2 | |
| odd <- 2 * c(1:l) | |
| even <- odd - 1 | |
| alpha <- beta[odd] * alpha[even] + delta[even] * alpha[odd] | |
| j <- 2 * j | |
| beta <- beta[odd] * beta[even] | |
| delta <- delta[even] * delta[odd] | |
| } | |
| Sn <- alpha/beta + 1 | |
| eval.Sn <- format(as.numeric(Sn), digits = 22) | |
| out <- switch(value, eval = eval.Sn, exact = Sn, both = list(Sn = Sn, | |
| eval.Sn = eval.Sn)) | |
| return(out) | |
| } | |
| </code></pre> | |
| <h2>The Gauss hypergeometric function </h2> | |
| <p>Now consider the <em>Gauss hypergeometric function</em> \( {}_2\!F_1 \). | |
| This is the function \( {}_2\!F_1(\alpha,\beta,\gamma; \cdot) \)<br/> | |
| with complex parameters \( \alpha \), \( \beta \), \( \gamma \not\in \mathbb{Z}^- \) and complex variable \( z \) defined for \( |z|<1 \) as the sum of an absolute convergent series: | |
| \[ | |
| {}_2\!F_1(\alpha,\beta,\gamma; z) = \sum_{n=0}^{\infty}\frac{{(\alpha)}_{n}{(\beta)}_n}{{(\gamma)}_{n}}\frac{z^n}{n!}, | |
| \] | |
| and extended by analytical continuation in the complex plane with the cut | |
| along \( (1,+\infty) \). Here \( {(a)}_n:=a(a+1)\cdots(a+n-1) \) denotes Pochhammer's symbol used | |
| to represent the \( n \)-th ascending factorial of \( a \). . </p> | |
| <p>The binary splitting allows to evaluate \( {}_2\!F_1(\alpha,\beta,\gamma; z) \) for rational values of | |
| \( \alpha,\beta,\gamma, z \) by manipulating only integer numbers. | |
| This is performed by the R function below</p> | |
| <pre><code class="r">## rational approximation of 2F1(a1/a2, b1/b2, c1/c2; p/q) | |
| ## with gmp ## | |
| hypergeo_bs <- function(a1, a2, b1, b2, c1, c2, p, q, m) { | |
| u <- function(i) c2 * (a1 + (i - 1) * a2) * (b1 + (i - 1) * | |
| b2) * p | |
| v <- function(i) a2 * b2 * i * (c1 + (i - 1) * c2) * q | |
| bs.gmp(u, v, m) | |
| } | |
| </code></pre> | |
| <p>For more convenience I have firstly written the function below which returns the irreducible rational | |
| notation of a given number \( x \). The user can also specify a rounding order for \( x \). </p> | |
| <pre><code class="r">n.decimals <- function(x, tol = .Machine$double.eps) { | |
| sapply(x, function(x) { | |
| i <- 0 | |
| while (abs(x - round(x, i)) > tol) { | |
| i <- i + 1 | |
| } | |
| return(i) | |
| }) | |
| } | |
| irred.frac <- function(x, rnd = n.decimals(x)) { | |
| b <- 10^rnd | |
| a <- as.bigz(b * round(x, rnd)) | |
| num <- a/gcd.bigz(a, b) | |
| den <- b/gcd.bigz(a, b) | |
| list(num = num, den = den) | |
| } | |
| </code></pre> | |
| <p>For example:</p> | |
| <pre><code class="r">irred.frac(pi) | |
| </code></pre> | |
| <pre><code>## $num | |
| ## Big Rational ('bigq') : | |
| ## [1] 3141592653589793 | |
| ## | |
| ## $den | |
| ## Big Rational ('bigq') : | |
| ## [1] 1000000000000000 | |
| </code></pre> | |
| <pre><code class="r">irred.frac(pi, rnd = 7) | |
| </code></pre> | |
| <pre><code>## $num | |
| ## Big Rational ('bigq') : | |
| ## [1] 31415927 | |
| ## | |
| ## $den | |
| ## Big Rational ('bigq') : | |
| ## [1] 10000000 | |
| </code></pre> | |
| <p>Finally, here is a user-friendly function for evaluating \( {}_2\!F_1 \) with the binary splitting:</p> | |
| <pre><code class="r">Hypergeometric2F1 <- function(a, b, c, z, m = 7, rnd.params = max(n.decimals(c(a, | |
| b, c))), rnd.z = n.decimals(z), check.cv = FALSE) { | |
| frac.a <- irred.frac(a, rnd.params) | |
| frac.b <- irred.frac(b, rnd.params) | |
| frac.c <- irred.frac(c, rnd.params) | |
| a1 <- frac.a$num | |
| a2 <- frac.a$den | |
| b1 <- frac.b$num | |
| b2 <- frac.b$den | |
| c1 <- frac.c$num | |
| c2 <- frac.c$den | |
| frac.z <- irred.frac(z, rnd.z) | |
| p <- frac.z$num | |
| q <- frac.z$den | |
| out <- hypergeo_bs(a1, a2, b1, b2, c1, c2, p, q, m) | |
| if (check.cv) { | |
| x <- hypergeo_bs(a1, a2, b1, b2, c1, c2, p, q, m + 1) | |
| cv <- x == out | |
| out <- list(result = out, convergence = cv) | |
| if (!cv) { | |
| out$convergence <- paste(out$convergence, " - m=", | |
| m, " need to be increased", sep = "") | |
| } | |
| } | |
| return(out) | |
| return(a) | |
| } | |
| </code></pre> | |
| <p>For example:</p> | |
| <pre><code class="r">a <- 20.5 | |
| b <- 11.92 | |
| c <- 19 | |
| z <- 0.5 | |
| Hypergeometric2F1(a, b, c, z) | |
| </code></pre> | |
| <pre><code>## [1] "8057.994139606238604756" | |
| </code></pre> | |
| <pre><code class="r">Hypergeometric2F1(a, b, c, z, m = 3, check.cv = TRUE) | |
| </code></pre> | |
| <pre><code>## $result | |
| ## [1] "1522.06880440136683319" | |
| ## | |
| ## $convergence | |
| ## [1] "FALSE - m=3 need to be increased" | |
| </code></pre> | |
| <pre><code class="r">Hypergeometric2F1(a, b, c, z, m = 7, check.cv = TRUE) | |
| </code></pre> | |
| <pre><code>## $result | |
| ## [1] "8057.994139606238604756" | |
| ## | |
| ## $convergence | |
| ## [1] TRUE | |
| </code></pre> | |
| <p>Note that Robin Hankin's <code>gsl</code> package does an excellent job:</p> | |
| <pre><code class="r">library(gsl) | |
| hyperg_2F1(a, b, c, z) | |
| </code></pre> | |
| <pre><code>## [1] 8058 | |
| </code></pre> | |
| </body> | |
| </html> | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment