http://qs321.pair.com?node_id=602366

zentara has asked for the wisdom of the Perl Monks concerning the following question:

UPDATE: See followup node below. FIXED as per syphilis 's correction. I leave this node untouched, so you can appreciate the correction and source of error

A recent Chatterbox comment, got me thinking about all the various ways to do this. Take the square root of 5 to 10000 digits, print it, and return back a 5 in the inverse operation.

My submission:

#!/usr/bin/perl use warnings; use strict; use Math::MPFR qw(:mpfr); #use FileHandle; #STDOUT->autoflush(1); $|++; Rmpfr_set_default_prec(10000); #print Rmpfr_get_default_prec(),"\n"; my $rop = Rmpfr_init2(10000); # return my $op = Rmpfr_init2(10000); # operand my $flip = Rmpfr_init2(10000); # return test Rmpfr_set_d ($op, 5.0 , GMP_RNDN); + #Set $rop to the square root of the 2nd arg rounded in the #direction $rnd. Set $rop to NaN if 2nd arg is negative. #Return 0 if the operation is exact, a non-zero value otherwise. my $bool = Rmpfr_sqrt($rop, $op, GMP_RNDN); print "$bool\n"; if($bool == 0){print "Exact\n";} print "$rop\n\n"; # dosn't print 10000 digits unless asked for print "length ",length $rop,"\n\n"; #my $s3 = Rmpfr_get_str($rop,10,0, GMP_RNDD); #the 0 means print as many digits needed to be exact #in the reverse operation # actually ask for 10000 my $s3 = Rmpfr_get_str($rop,10,10000, GMP_RNDN); print "$s3\n\n"; print "length ", length $s3,"\n\n\n"; # see if return is accurate #Set $flip to the square of $op, rounded in direction $rnd. my $si = Rmpfr_sqr($flip, $rop, GMP_RNDN ); print "go back\n"; print "$flip\n"; # pretty close :-)

I'm not really a human, but I play one on earth. Cogito ergo sum a bum

Replies are listed 'Best First'.
Re: Fun Question: Square Root of 5 to 10000 digits
by Zaxo (Archbishop) on Feb 27, 2007 at 18:42 UTC

    My entry was,

    perl -MMath::Pari -e'Math::Pari::setprecision(10000); my $foo = PARI(5 +); my $bar = $foo->sqrt; print $bar,$/,$bar**2,$/'
    which is cheating I suppose. But it's really fast.

    After Compline,
    Zaxo

      You can do essentially the same with Math::MPFR:
      perl -MMath::MPFR=:mpfr -e "Rmpfr_set_default_prec(10000); $foo = Math +::MPFR->new(5); $bar = $foo ** 0.5; print $bar, $/, $bar ** 2, $/"

      With Math::MPFR, setting a default precision of 10000, means you have 10000 bits of precision, not 10000 decimal digits of precision. To be guaranteed 10000 decimal digits of precision I think you need to ask for 33220 bits - ie ceil(10000 / log(2)). I haven't checked whether Math::PARI sets its precision in bits or decimal digits.

      Cheers,
      Rob
        With Math::MPFR, setting a default precision of 10000, means you have 10000 bits of precision, not 10000 decimal digits of precision.

        Doh! You are right, I was wondering last night why I was getting a 1/3 factor of digits. 10000 gave me ~ 3000 decimal places, while 100000 precision gave me ~30000 decimal places. I thought it might be the print accuracy needed to do the inverse, but I was wrong. Thanks for clearly that up. :-)


        I'm not really a human, but I play one on earth. Cogito ergo sum a bum
Re: Fun Question: Square Root of 5 to 10000 digits
by eric256 (Parson) on Feb 27, 2007 at 19:48 UTC

    Computing 10k digits the hard slow way (at this time it is still computing, climbing through digit 2000 as I type.

    use strict; use warnings; =pod Algorithm http://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Digit_b +y_digit_calculation Write the decimal expansion of the square. The numbers are laid out in + a fashion similar to the long division algorithm, and the root will +appear above its square. Divide the digits of the square into pairs, +starting from the decimal point and going both left and right. The de +cimal point of the root should be placed above the decimal point of t +he square. One digit of the root will appear above each pair of digit +s of the square. Begin at the left-most (non-zero) pair of digits of the square. Initia +lize the remainder to zero. Initialize p to zero. For each iteration (pair of digits) do: 1. Bring down the most significant (leftmost) pair of digits of the + square not yet used (if all the digits have been used, use the digit +s 00) and append them to the remainder from the previous step, i.e. m +ultiply the remainder by one hundred and add the two digits. This wil +l be the current value c. 2. Let p be the part of the root found so far, ignoring any decimal + point (for the first step, p = 0). Determine the greatest digit x su +ch that y = (20 \cdot p + x) \cdot x does not exceed c (notice 20p + +x is simply twice p, with the digit x appended to the right). You can + find x by guessing what c/(20·p) is and doing a trial calculation of + y, then adjusting x upward or downward as necessary. Place the digit + x as the next digit of the root, i.e above the two digits of the squ +are which you just brought down. Thus the next p will be the old p ti +mes 10 plus x. 3. Subtract y from c to form a new remainder. 4. If the remainder is zero and there are no more digits to bring d +own, then the algorithm has terminated. Otherwise go back to step 1 f +or another iteration. =cut $|++; use Data::Dumper; use bignum; my $number = "5.0"; my ($front, $back) = split '\.', $number; $front = '0' . $front if length $front % 2; $back .= '0' if length $back % 2; print "Front $front, Back $back\n"; my @digits; push @digits, grep {length} split(/(..)/, $front), ".", grep {length} +split(/(..)/, $back); print "Digits ", join("-", @digits), "\n"; my $i = 0; my $remainder = 0; my $p = 0; # root without decimal my $root = 0; # root with decimal my $x = 0; my $y = 0; my $c = 0; while ($i < 10_000) { my $pair = shift @digits; $pair ||= '00'; if ($pair eq '.') { $root .= '.'; print "$i\t.\n"; next; } $i++; $c = (100*$remainder) + $pair; $x = int($c / (20 - $p)); $x = 0 if $x < 0; #print "Guessing X is $x when C = $c and P = $p\n"; $y = ((20 * $p) + $x) * $x; if ($y <= $c) { while (1) { $x++; $y = ((20 * $p) + $x) * $x; #print "X = $x, Y = $y, C = $c\n" if $i == 308; if ($y > $c) { $x--; $y = ((20 * $p) + $x) * $x; last; } } } elsif ($y > $c) { while (1) { #print "X = $x, Y = $y, C = $c\n" if $i == 308; $x--; die if $x < 0; $y = ((20 * $p) + $x) * $x; last if ($y < $c) } } $p = $p . $x; $root = $root . $x; print "$i\t$x"; #print "X = $x, Y = $y, C = $c"; print "\n"; $remainder = $c -$y; last if ($remainder == 0 and scalar @digits == 0); } print "\nAfter $i iterations:\nRoot=$root\n";

    It completed about 5 hours later and it returns 4.999999999999999.....9999999......9999 ;)

    2.23606797749978969640917366873127623544061835961152572427089724541052 +092563780489941441440837878227496950817615077378350425326772444707386 +358636012153345270886677817319187916581127664532263985658053576135041 +753378500342339241406444208643253909725259262722887629951740244068161 +177590890949849237139072972889848208864154268989409913169357701974867 +888442508975413295618317692149997742480153043411503595766833251249881 +517813940800056242085524354223555610630634282023409333198293395974635 +227120134174961420263590473788550438968706113566004575713995659556695 +691756457822195250006053923123400500928676487552972205676625366607448 +585350526233067849463342224231763727702663240768010444331582573350589 +309813622634319868647194698997018081895242644596203452214119223291259 +819632581110417049580704812040345599494350685555185557251238864165501 +026243631257102444961878942468290340447471611545572320173767659046091 +852957560357798439805415538077906439363972302875606299948221385217734 +859245351512104634555504070722787242153477875291121212118433178933519 +103800801111817900459061884624964710424424830888012940681131469595327 +944789899893169157746079246180750067987712420484738050277360829155991 +396244891494356068346252906440832794464268088898974604630835353787504 +206137475760688340187908819255911797357446419024853787114619409019191 +368803511039763843604128105811037869895185201469704564202176389289088 +444637782638589379244004602887540539846015606170522361509038577541004 +219368498725427185037521555769331672300477826986666244621067846427248 +638527457821341006798564530527112418059597284945519545131017230975087 +149652943628290254001204778032415546448998870617799819003360656224388 +640963928775351726629597143822795630795614952301544423501653891727864 +091304197939711135628213936745768117492206756210888781887367167162762 +262337987711153950968298289068301825908140100389550972326150845283458 +789360734639611723667836657198260792144028911900899558424152249571291 +832321674118997572013940378819772801528872341866834541838286730027431 +532022960762861252476102864234696302011180269122023601581012762843054 +186171761857514069010156162909176398126722596559628234906785462416185 +794558444265961285893756485497480349011081355751416647462195183023552 +595688656949581635303619557453683223526500772242258287366875340470074 +223266145173976651742067264447621961802422039798353682983502466268030 +546768767446900186957209958589198316440251620919646185105744248274087 +229820410943710992236175285315302212109176295120886356959716907946257 +260325089752229704043412880822332153390119551566514079022175646165421 +295787804223138207855367690772666643131659319546206872064645091487274 +408248812817765347516867907359186246442687464199149977893991312947201 +459199967825762063948526250359428286402462255910378955634538283178235 +598391296251160036910131265905719718200181724360595512757851998329989 +285638604458710469334951865390330842804218272603638944541578024417457 +472341469729999631251094562274695974331390549780162887681065496756275 +649338348884592698294163140147050914141795453509386876452390937230662 +419067158476029218547020420238380436721350194617915057915493628459086 +788770986310679260761458338351692202921990110129607358608294473144079 +720147101521804634625003226409687167296354096963621983204885046543344 +380378669192757217575057403478718606026718022474204783425318094052698 +805661533753487277302654212560646348138634668964687129063701162706217 +099466701519933557424898116727350826578172481264912790714425048522340 +556057312086469885674603451148811674556535992063478728026575255402487 +359662289287389534106254498482094334002764956625731301298686836078008 +203561067901175449173311510458783164794168354596674564623051385218599 +188448000112125335734871584794490816963530394687253053788977710544054 +955749467196707345552281518342410265386896750598329996187204923568853 +514555358003838381407610440922464964782652208654383369024612047255787 +090864923539951507378083527300509570276492629316672766752047155798534 +597726432371679180727996367691655289824919618740861111192275946865226 +966098989937362179071392696563562577250729216840678930763888389142853 +336474367898366474181714970053313607979488132421072061280052163422533 +199083987374632189144577621841557645544072733689630651234568235381958 +533331044769376622743705983843263810403137262445641431199752936847104 +118570743515615312100735462619503824479477744493651611431948914809685 +975614704431968533532515615412403886080108510031662500603506818823438 +203859780768945006659760490028735936883389591909060918206276232437409 +135995732732349211914000689194022705036269201313107040695776234829988 +254965283042711355278807814207741763646761360670007609360034961644118 +219368840528928043754106802006360576332883061827878963128063856455709 +648290210063776503799401497245758843117914856404333141243761811561761 +006493539825945744207741473948128713349178405173131479571217191330682 +140719956640892672692970978995327770702091054596484581399697707393656 +092919491525302868718101876642487486667741033314298011814211340497759 +716087436302522008807629760814504881232858044956454305448224170131577 +677424987270213612730333486444655535511594798540752463829409464791024 +121411007984176885207417581686668523676827194156329659107428643922379 +007595429260015111950759140710454289863826434511288025661836100900179 +843741024237213867146307791870158060147345404662833064084680310748288 +537430811023295922286646049708808188138229122797460520790365633606896 +505086534771518011208640490745438582497291626668833970598782714957397 +915972878996046094233934314724567824036254625833179905519838440636744 +713654558771274662530959971824926550060121185134909958870176238590113 +709865187106374582836022728243749415052562137396602715210494388911864 +391071922090566062976782353860239317166862884978979713116850166821885 +900554395166704488582514729876150834227478487520287013659756986542599 +502457376392099671550317543560821394263933506954389584527303803267954 +256947815867222238281799661120672212197434356611087080712179058581636 +928287427858875627120964077895825149015415115020600484145325800361808 +458684988518121332282664574453961380291989023990956032798302825225051 +456561328662523314938776390212884334774360002200843696605161833086767 +498472823677771293702863001274638085902962938848629217905094144074811 +133826138441981609638905950221300928562108355105181903742637767182953 +199208263592041883061717106647754507604654552659547442862559364334324 +688423664036057628254948863376944369187855628709481999981444664061185 +259532224766559666339765078625240130074057689565733388089461589420952 +251173167505972472501999646467194310144676766648816305155638672852526 +086605317916341600902557746231871175494434512989400103273345154307841 +968190065490224307374601824399259045531826327418793721454268538524630 +950660875986633162214739286284343958868112783102421621627252537771394 +961361202338378835005445974831739835829069989248883880243957172027473 +216573814447302954278253748419330275751241183708657776683485841803126 +266566387151244179427531261957003126309964912891730849585871445657501 +216962906702704363459175865982342006495244410438929021072490102597686 +174268887901448853470292572359836646729196739275265445151383194479087 +661041732949484763022158546989673904792958537987396649835999005579000 +120919322626926726049899902961610658035805950365031750098014870375967 +023672065445545203434809071143317711156594582123916387034211096515861 +418201152717398038594435990337462351126288971296200440028509081108585 +469176742320419895891441647560873743788961127378365160488999263756684 +054982030671582145467250657813866948247604444023252554238617089700590 +838264008019997311333035513281907312395795636760902060713020263178917 +805743722173811787894273602969140036732991294406588668748597892854825 +102871811686968183909740304722806347827807232880396910102098242339584 +002403999210139899328060701727385807882014038901064032469745526465464 +898879260961781108502759446629503704141820501273719633590609636201478 +849063400477609519668646900828516862812722544219204564846756456180559 +531921554216987830349774633755427044780182342347018372013092401980499 +517055585085563319407669901160212523106673821875693195421059500446346 +148243556688378823691931722059603755748548912773393225544900769172105 +283020608179651555508948230664152815176335502995107609423259335542011 +753292319099355385410992478797141851014054813995628168624993072614373 +306743612117484485196330614105147669083154108584325996229835017222623 +531546344191231312957390948978542641216127091558924829062133967484227 +596337927647066608955766338679457457836207328166539713976508877033351 +724579861392869369795029681758079295208407220412043434788940526975267 +308786390458154767233477962356248496731156210068338903127252086007331 +486216953309755560257155847290837044894472342748458511683186271225732 +743340656144343106785292651461346127821708217736171485677176561204606 +682817100781947077452269023925852831990425578622708862920305461805107 +654208651932453487807491127224572278156638867141180076297401797322630 +796391714884660883941713393444586285461482769765577951177721599477408 +940406333669713883981930960596498639635315853659711259446021365554470 +325481567614863755654636423839390560103217583144257652675936462545125 +740003000365859515459987158189839281526885723151427088855796766080909 +405420389160085164042404689161260690067316294437098407359979945870707 +839362426639033075949079822298848890366067717682580803635637642752013 +318569882735508634903210818774220737430423280811643868942408965551921 +083389729079756652539096027830038077991862613406373233413274392805138 +573427774293262378385371365383990552915995436557518921923234437736218 +909303157738244821219628394537238510957579852630718945845650161085033 +813628021563592043770661524611127632632803844909006515313478519953035 +061602854321428617437725719672074930114268409340166865508460554095586 +622367333806465770613747759814271801480609814919779027295375217356886 +476496437861235140639127606461639438727134548392877452517412308661459 +274076255030340812010115189765447712690312781053154208529189520811139 +0191968177780752415991327760357237118318882234501846265595422760

    ___________
    Eric Hodges
Re: Fun Question: Square Root of 5 to 10000 digits
by zentara (Archbishop) on Feb 28, 2007 at 14:01 UTC
    Here is a better script, which includes a generic method for determining the bits/digit needed.
    #!/usr/bin/perl use warnings; use strict; use POSIX qw(ceil); use Math::MPFR qw(:mpfr); sub log10{ my $n = shift; return log($n)/log(10); } my $req_digits = 10000; my $acc = ceil( $req_digits/log10(2) ); print "$acc\n"; $|++; # 33220 is bits needed for 10000 digits Rmpfr_set_default_prec($acc); my $rop = Rmpfr_init2($acc); # return + my $op = Rmpfr_init2($acc); # operand + my $flip = Rmpfr_init2($acc); # return test + Rmpfr_set_d ($op, 5.0 , GMP_RNDN); + #Set $rop to the square root of the 2nd arg rounded in the #direction $rnd. Set $rop to NaN if 2nd arg is negative. #Return 0 if the operation is exact, a non-zero value otherwise. my $bool = Rmpfr_sqrt($rop, $op, GMP_RNDN); print "$bool\n"; if($bool == 0){print "Exact\n";} print "$rop\n\n"; print "length ",length $rop,"\n\n"; # see if return is accurate #Set $flip to the square of $op, rounded in direction $rnd. my $si = Rmpfr_sqr($flip, $rop, GMP_RNDN ); print "go back\n"; print "$flip\n"; # pretty close :-)

    I'm not really a human, but I play one on earth. Cogito ergo sum a bum