Your skill will accomplishwhat the force of many cannot PerlMonks

### Fun Question: Square Root of 5 to 10000 digits

by zentara (Archbishop)
 on Feb 27, 2007 at 18:36 UTC Need Help??

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 :-)

[download]

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,$/' [download] 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,$/"
[download]

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";
[download]

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
[download]

___________
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 :-)

[download]

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

Create A New User
Node Status?
node history
Node Type: perlquestion [id://602366]
Approved by Zaxo
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others chilling in the Monastery: (3)
As of 2019-12-15 21:36 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?

No recent polls found

Notices?