### Author Topic: Praying Pythagorean triples  (Read 5346 times)

#### B+

• Guest
##### Praying Pythagorean triples
« on: January 21, 2016, 03:24:11 pm »
Hi Tomaaz (and/or all),

If you have some time, could you show me this in Ruby. JB is the only Basic I have with extended integers and even that won't do SQR of arbitrary long integers. sbmeTW has found cool workaround that seems to work in JB, I wonder if it could be verified in another language, (or a better Basic would be even better!)

Here is sbmeTW code in JB that is confusing to me but seems to work
Code: [Select]
`'sbmeTW Pray Pyth Tri 100.txt copied and modified 2016-01-20print" gives prayers"print "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^"   'Input "input  number list ";xa=6 :b=5:c=1:s=5: D=1: z =8for p=1 to 100    if p=2 then D=7    gosub [legs]    s=(a*b)-c    print p;")  ";La;",  ";Lb;",  ";s;    if La*La+Lb*Lb<>s*s then        print "  Error !!!"; s*s-La*La-Lb*Lb        ecount=ecount+1    else        print "  Checks OK"    end if    print La*La;" + ";Lb*Lb;" ?=? ";s*s    print    c=b:b=snextend[legs]    La=int(((z*b)+D)/2)    Lb=La+1    TEMPB=La+Lb    D=TEMP    TEMP=TEMPBreturn`
Here is my code in JB that bombs after the 20th round that I am wondering if would work in Ruby or better BASIC
Code: [Select]
`'Pray Pyth Tri.txt for JB 2016-01-19 [B+MGA]' 2016-01-19 started from copy of my post at JB,' old business I wanted to check up on' from copy of my post at JB,' Have you heard of praying Pythagorean triangles?' There is only a difference of 1 between the two legs (hence they are always primitive).' The first is the first 3-4-5' The secound is 20-21-29 but they are rare.' They can be constructed from sequence that starts at 1' Next term is INT(Last term times (SQR(2)+1)squared)' 1, 5, 29, 169, 985,..' Those are hypotenuse c, the legs are:' a is INT ( SQR( (c)squared /2) )' b is a+1' 2016-01-19 show code of the program outlinedc=1:ecount=0:pptcount=0while ecount<3    c= INT(c*(SQR(2)+1)^2)    a= INT( SQR((c^2)/2) )    b=a+1    pptcount=pptcount+1    print pptcount;") ";a;",  ";b;",  ";c    print a*a;" + ";b*b;" = ";c*c;    if a*a+b*b<>c*c then        print "  Error !!!"        ecount=ecount+1    else        print "  Checks OK"    end if    printwend `

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #1 on: January 22, 2016, 12:49:21 am »
I'm not sure what is that you wanted to achieve, but here are two programs in Python:

This one
Code: [Select]
`from math import sqrtfor x in range (1, 1000000001): y = sqrt(x * x + (x + 1) * (x + 1)) if y == int(y): print(x, x + 1, int(y))`
finds all Pythagorean Triples where the legs are smaller than billion and there is only a difference of 1 between them. The result is:
Quote
3 4 5
20 21 29
119 120 169
696 697 985
4059 4060 5741
23660 23661 33461
137903 137904 195025
803760 803761 1136689
4684659 4684660 6625109
27304196 27304197 38613965
65918161 65918162 93222358
104532126 104532127 147830751
120526554 120526555 170450288
159140519 159140520 225058681
197754484 197754485 279667074
213748912 213748913 302286611
229743340 229743341 324906148
236368449 236368450 334275467
252362877 252362878 356895004
274982414 274982415 388883860
290976842 290976843 411503397
306971270 306971271 434122934
313596379 313596380 443492253
345585235 345585236 488731327
368204772 368204773 520720183
384199200 384199201 543339720
390824309 390824310 552709039
400193628 400193629 565959257
406818737 406818738 575328576
416188056 416188057 588578794
422813165 422813166 597948113
429438274 429438275 607317432
438807593 438807594 620567650
445432702 445432703 629936969
454802021 454802022 643187187
461427130 461427131 652556506
477421558 477421559 675176043
484046667 484046668 684545362
493415986 493415987 697795580
500041095 500041096 707164899
509410414 509410415 720415117
516035523 516035524 729784436
532029951 532029952 752403973
554649488 554649489 784392829
570643916 570643917 807012366
586638344 586638345 829631903
593263453 593263454 839001222
609257881 609257882 861620759
625252309 625252310 884240296
631877418 631877419 893609615
641246737 641246738 906859833
647871846 647871847 916229152
654496955 654496956 925598471
657241165 657241166 929479370
663866274 663866275 938848689
670491383 670491384 948218008
679860702 679860703 961468226
686485811 686485812 970837545
695855130 695855131 984087763
702480239 702480240 993457082
709105348 709105349 1002826401
718474667 718474668 1016076619
725099776 725099777 1025445938
741094204 741094205 1048065475
750463523 750463524 1061315693
757088632 757088633 1070685012
763713741 763713742 1080054331
770338850 770338851 1089423650
773083060 773083061 1093304549
779708169 779708170 1102673868
782452379 782452380 1106554767
786333278 786333279 1112043187
789077488 789077489 1115924086
795702597 795702598 1125293405
798446807 798446808 1129174304
802327706 802327707 1134662724
808952815 808952816 1144032043
811697025 811697026 1147912942
818322134 818322135 1157282261
824947243 824947244 1166651580
827691453 827691454 1170532479
834316562 834316563 1179901798
840941671 840941672 1189271117
843685881 843685882 1193152016
850310990 850310991 1202521335
856936099 856936100 1211890654
863561208 863561209 1221259973
866305418 866305419 1225140872
872930527 872930528 1234510191
875674737 875674738 1238391090
879555636 879555637 1243879510
882299846 882299847 1247760409
888924955 888924956 1257129728
895550064 895550065 1266499047
898294274 898294275 1270379946
904919383 904919384 1279749265
908800282 908800283 1285237685
911544492 911544493 1289118584
914288702 914288703 1292999483
918169601 918169602 1298487903
920913811 920913812 1302368802
924794710 924794711 1307857222
927538920 927538921 1311738121
934164029 934164030 1321107440
936908239 936908240 1324988339
943533348 943533349 1334357658
950158457 950158458 1343726977
952902667 952902668 1347607876
956783566 956783567 1353096296
959527776 959527777 1356977195
966152885 966152886 1366346514
968897095 968897096 1370227413
972777994 972777995 1375715833
975522204 975522205 1379596732
979403103 979403104 1385085152
982147313 982147314 1388966051
988772422 988772423 1398335370
991516632 991516633 1402216269
995397531 995397532 1407704689
998141741 998141742 1411585588

and this one
Code: [Select]
`from math import sqrtplik = open('pythagorean-triples.txt', 'w')for x in range(3, 10001): for y in range(1, x - 1): plik.write(str(x * x - y * y) + " " + str(2 * x * y) + " " + str(x * x + y * y) + "\n") z += 1plik.close()`
finds and writes to the file 49985001 (almost 50 million) Pythagorean Triples. The size of the file is 1.29 GB.

EDIT I've chosen Python, cause it's faster than Ruby.
« Last Edit: January 22, 2016, 01:01:22 am by Tomaaz »

#### B+

• Guest
##### Re: Praying Pythagorean triples
« Reply #2 on: January 22, 2016, 03:31:51 am »
Hi Tomaaz,

Thanks for your prompt post. Our tables match exactly item for item up through triple #10 then diverge rapidly. I suspect something is wrong after 10th when the numbers cease to expand but remain of almost of constant length for remainder of table. The table should expand wider all the way down like one side of pyramid. Your code is certainly to most elegant.

I just spent afternoon and evening getting a SQR funtion working for arbitrary long integers and was able to confirm sbmeTW's results from a different method as shown in code above.

Here is our 100th c of a, b, c st a^2+b^2=c^2 and b=a+1:
30645573943232956180057972969833245887630954508753693529117371074705767728665

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #3 on: January 22, 2016, 11:03:24 am »
The table should expand wider all the way down like one side of pyramid.

Yes. One condition was missing from my code, but with this condition added the simple method I was using is very slow. To confirm your number I need to find something different.

#### B+

• Guest
##### Re: Praying Pythagorean triples
« Reply #4 on: January 22, 2016, 02:38:20 pm »
Thank you Tomaaz, I hope the following is more clear, I know I can be difficult to follow.

Yes, I will give you formula for c (the hypotenuse of a,b,c of PRAY Pythagorean triangle) from which all else is determined quite quickly hold on...

Code: [Select]
`'Pray Pyth Tri 100 hypotenuse.txt for JB 2016-01-21 [B+MGA]'just the hypotenuse please, from which legs can be calculated and then a^2+b^2=c^2 can be confirmed'                                             TO AVOID MASSIVE CALCULATIONS AND TESTS'  THIS IS THE SIMPLE FORMULA FOR THE SEQUENCE OF HYPOTENUSE FOR PRAYING PYTHAGOREAN TRIANGLES' C = Hypotenuse of PRAY Pythagorean Triangle, starting C=1 next C formula for loop:' C = C * [SquareRoot(2) + 1] ^ 2'AS YOU MIGHT NOTICE IF ONE C IS WRONG, ALL THAT FOLLOW ARE USELESS' THE TRICK IS TO USE AN APPROPRIATELY LONG ENOUGH EXPANSION OF THE SQUARE ROOT OF 2 TO MAINTAIN ACCURACY OF C'for Pray Pyth Tri Leg A = INT((Hypotenuse/2)^2)  and leg B is just Leg A+1 then confirm A^2+B^2=C^2 that you have everything right' AGAIN WITH LEG A YOU NEED THE CALCULATION EXPANDED ENOUGH THAT ALL THE INTEGER PLACES ARE CORRECT'100th c from this program then sbmeTW's 100th c and now with Newtons Method for SQR 'now with abbreviated code for hypotenuse only'30645573943232956180057972969833245887630954508753693529117371074705767728665 (my first run with expanded SQR(2) const from NASA)'30645573943232956180057972969833245887630954508753693529117371074705767728665  Checks OK (sbmeTW)'30645573943232956180057972969833245887630954508753693529117371074705767728665 using bSQR with Newton's Method for SQR'30645573943232956180057972969833245887630954508753693529117371074705767728665 confirmed again with this abreviated version for Hypotenuse only'I have to look up square root of 2 from table on-line as float precision for JB is 14 places max'OK lets see if we can get the first  200+ listed with the help of SQR(2) listing from http://apod.nasa.gov/htmltest/gifcity/sqrt2.1mil'square root of 2 + 1 with decimal removedaconst\$="241421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851"'print len(aconst\$) :end'use 200 digits to right of decimal = (sqr(2) + 1)*10^200aconst\$=left\$(aconst\$,201)bconst=int(val(aconst\$))cconst=bconst*bconst              'avoid float math by doing cconst/dconst at the last calculation for new hypotenusedconst=10^200*10^200            ' cconst/dconst is just (SQR(2)+1)^2 to 200 places past decimalc=1:ecount=0:pptcount=0while pptcount<100              ' the new c is calculated from the last c times const (SQR(2)+1)^2 and then rounded down as INT function does    c= int(c*cconst/dconst)  'cconst/dconst = (SQR(2)+1)^2 but you need square root of 2 expanded way past 100 places closer to 200 places    pptcount=pptcount+1    print pptcount;") ";c    scanwend`
Hypotenuse Table for PRAY PYTH TRI
Code: [Select]
`1) 52) 293) 1694) 9855) 57416) 334617) 1950258) 11366899) 662510910) 3861396511) 22505868112) 131173812113) 764537004514) 4456048214915) 25971752284916) 151374465494517) 882275040682118) 5142275778598119) 29971379630906520) 174686002006840921) 1018144632410138922) 5934181792453992523) 34586946122313816124) 201587494941428904125) 1174938023526259608526) 6848040646216128746927) 39913305853770512872928) 232631794476406948490529) 1355877461004671178070130) 7902632971551620119930131) 46059920368305049541510532) 268456889238278677129132933) 1564681415061367013233286934) 9119631601129923402270588535) 53153108191718173400390244136) 309799017549179117000070876137) 1805640997103356528600035012538) 10524046965070960054600139198939) 61338640793322403799000800180940) 357507797794863462739404661886541) 2083708145975858372637427171138142) 12144741078060286773085158364942143) 70784738322385862265873523018514544) 412563688856254886822155979746144945) 2404597394815143458667062355458354946) 14015020680034605865180218153003984547) 81685526685392491732414246562565552148) 476098139432320344529305261222389328149) 2774903309908529575443417320771770416550) 16173321720018857108131198663408233170951) 94265027010204613073343774659677628608952) 549416840341208821331931449294657538482553) 3202236015037048314918244921108267602286154) 18663999249881081068177538077354948075234155) 108781759484249438094146983543021420849118556) 634026557655615547496704363180773577019476957) 3695377586449443846886079195541620041267742958) 21538238961041047533819770810068946670586980559) 125534056179796841356032545664872059982254140160) 731666098117740000602375503179163413222937860161) 4264462532526643162258220473410108419355373020562) 24855109097042118972946947337281487102909300262963) 144866192049726070675423463550278814198100428556964) 844342043201314305079593833964391398085693271078565) 4921186067158159759802139540236069574316059197914166) 28682774359747644253733243407452026047810661916406167) 167175460091327705762597320904476086712547912300522568) 974369986188218590321850682019404494227476811886728969) 5679044457037983836168506771211950878652312959019850970) 33099896756039684426689189945252300777686400942232376571) 192920336079200122723966632900301853787466092694374408172) 1124422119719161051917110607456558821947110155224014072173) 6553612382235766188778697011839051077895194838649710024574) 38197252173695436080755071463577747645424058876674246074975) 222629900659936850295751731769627434794649158421395766424976) 1297582151785925665693755319154186861122470891651700352474577) 7562863010055617143866780183155493731940176191488806348422178) 44079595908547777197506925779778775530518586257281137738058179) 256914712441231046041174774495517159451171341352198020079926580) 1497408678738838499049541721193324181176509461855906982741500981) 8727537359991799948256075552664427927607885429783243876369078982) 50867815481211961190486911594793243384470803116843556275472972583) 296479355527279967194665394016095032379216933271278093776468756184) 1728008317682467841977505452501776950890830796510825006383339564185) 10071570550567527084670367320994566672965767845793671944523568628586) 58701414985722694666044698473465623086903776278251206660758072206987) 342136919363768640911597823519799171848456889823713568020024864612988) 1994120101196889150803542242645329408003837562664030201459391115470589) 11622583687817566263909655632352177276174568486160467640736321828210190) 67741382025708508432654391551467734249043573354298775642958539853790191) 394825708466433484332016693676454228218086871639632186217014917294530592) 2301212868772892397559445770507257635059477656483494341659130963913392993) 13412451504170920901024657929367091582138779067261333863737770866185826994) 78173496156252633008588501805695291857773196747084508840767494233201568595) 455628525433344877150506352904804659564500401415245719180867194533023584196) 2655597656443816629894449615623132665529229211744389806244435672964939936197) 15477957413229554902216191340833991333610874869051093118285746843256616032598) 90212146822933512783402698429380815336136020002562168903470045386574756258999) 5257949235243715217981999992354509006832052451463219203025345254761919215209100) 30645573943232956180057972969833245887630954508753693529117371074705767728665`
Although I have now pretty much confirmed to myself this is accurate, the perspective from another dialect or PL can catch things I might miss, specially errors using float math which I will make over and over because I use that same dialect or PL.
« Last Edit: January 22, 2016, 03:55:42 pm by B+ »

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #5 on: January 22, 2016, 05:05:27 pm »
OK, I see.

Now I can confirm that 30645573943232956180057972969833245887630954508753693529117371074705767728665 is a correct number. And this is a 1000th number 3038095123235207225722532398795271888915972182692220979441709989450932553745204304366970394001744336363561330945028635859982168423497065504307505845123315541212047129732091243719906387145725142746680845253661481550583558682474536256339697014842860398275579017532206595897411636873353687605060717900641424711342658860465041261456442244068576488739643985259772988571962412743970628410345853384053301908371649136570973582212620809195940648180427240766040576866115744453515052131389778963569464801938728739281184420051637609468823775061145715348032393933387142361323849988656096551850301152046220836131098950815100998117189916885430865859461271720719771561381398781642295115105584950514614814370923444030465233817437069543426599429189080433481691668956913306952960375665. And the one you get after 10000 rounds is 7656 digits long (attached). To calculate it I used sqr(2) that is 100000 digits long. And the code looks like this:
Code: [Select]
`from decimal import *getcontext().prec = 100000s = 2 * (Decimal(2).sqrt())c = 1for x in range(10000): c = int(c * (3 + s)) print(c)`
Is this formula mathematically proven or do you need additionally to check the numbers?
« Last Edit: January 22, 2016, 05:08:26 pm by Tomaaz »

#### jbk

• Guest
##### Re: Praying Pythagorean triples
« Reply #6 on: January 22, 2016, 06:38:41 pm »
Hi Tomaaz
I timed your program for x in range(100) 141.53 seconds
Code: [Select]
`from decimal import *import timestart = time.time()getcontext().prec = 100000s = 2 * (Decimal(2).sqrt())c = 1prec=10for x in range(100): getcontext().prec = prec c = int(c * (3 + s)) print(c) prec +=10 end = time.time()-startprint end`time 77.44 seconds
using cdecimal ( http://www.bytereef.org/mpdecimal/ ) instead of decimal the times are 0.67 and 0.14 seconds respectively
« Last Edit: January 22, 2016, 06:42:59 pm by jbk »

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #7 on: January 22, 2016, 06:53:30 pm »
Hi Tomaaz
I timed your program for x in range(100) 141.53 seconds

For 100 numbers, change the precision of sqrt to 1000 (or even 100). The whole program takes less than a second. What you did, doesn't change the precision, cause variable s (that is used in calculations) still has 1000000 digits. Also, I think that getcontex().prec in a loop makes the whole program much slower. Try this:

Code: [Select]
`from decimal import *import timestart = time.time()getcontext().prec = 1000s = 2 * (Decimal(2).sqrt())c = 1for x in range(100): c = int(c * (3 + s)) print(c) end = time.time()-startprint end`
Shouldn't take more than a second.

#### B+

• Guest
##### Re: Praying Pythagorean triples
« Reply #8 on: January 22, 2016, 06:56:29 pm »
Hi Tomaaz,

I am very impressed! So with Python you can set the level of precision for variables it seems,
Code: [Select]
`getcontext().prec = 100000 `, I am guessing. Nice!

I realized this was a sort of test for float math with the dialect/PL one is using.

Tomaaz
Quote
Is this formula mathematically proven or do you need additionally to check the numbers?

No, no additional check needed I've the confirmation I need, to go farther (and/or is it further?) I might have to extend my own Square root function precision into decimal places (well that might be useful (unless I get Python and it's already done for me)).

Proven mathematically? here is my source: Wonder of Numbers by Clifford A. Pickover CR 2002 Chapter "Further Exploring" section on Chapter 99 "Everything you wanted to know about triangles but were afraid to ask", good credentials but I saw no mathematical proof. But if you are taking bets...

The curious thing about this whole thing is SbmeTW's code (at JB forum), who started the thread there on Pythagorean Triples which started me on the PRAYing type. That code bypasses the square root calculations altogether, the code though confusing might be arranged into even more elegance than your Python code because it is worry free of float precision though you do still need big integer maths.

#### B+

• Guest
##### Re: Praying Pythagorean triples
« Reply #9 on: January 22, 2016, 06:58:24 pm »
To Tomaaz and jbk (hi!),

For precision setting for 100 Pray Pyth Tri, you may need more than 100 but less than 200, at least I did for SQR(2) to get 100 PPT's. With 100 places past decimal my code failed checks just below 70 Pray PTs.

If you alter the code allot trying to get more speed, I recommend doing proper checks that you are actually getting Pray Pyth Triangles and haven't mutated to some other thing. (OTOH if getting same c for 100, you are probably good, for 100).

EDIT: this post was overhauled because I had to leave when first posting.
« Last Edit: January 22, 2016, 07:22:28 pm by B+ »

#### jbk

• Guest
##### Re: Praying Pythagorean triples
« Reply #10 on: January 22, 2016, 07:16:09 pm »
Also, I think that getcontex().prec in a loop makes the whole program much slower. Try this:
actually getcontex().prec in a loop works quite well, 141 seconds without changing the prec vs 77 seconds with the change.
I realize that the precision of s does not change but the calculations are faster when lowering the prec.
I also realize that for x=100 the starting precision for s could be set as you suggested.
« Last Edit: January 22, 2016, 07:19:56 pm by jbk »

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #11 on: January 22, 2016, 07:26:46 pm »
actually getcontex().prec in a loop works quite well, 141 seconds without changing the prec vs 77 seconds with the change.

Well, don't get me wrong, but... With precision set to 100 digits the whole program takes less than a second for 100 numbers. Even with precision set to 100000 digits, it takes no more than 10 sec. How did you manage to get 141 and 77 sec.? Also, the sqrt(2) is calculated only once and then is used in a main loop. Does changing precision in a loop affects calculation that were already made? I'm not 100% sure, but, looking at your 141 and 77 sec. (comparing to mine less than 1 sec.), I don't think so.

#### jbk

• Guest
##### Re: Praying Pythagorean triples
« Reply #12 on: January 22, 2016, 07:38:18 pm »
the point I was in vain trying to make is that setting the precision in the loop speeds things up
Code: [Select]
`from decimal import *import timestart = time.time()getcontext().prec = 10000s = 2 * (Decimal(2).sqrt())c = 1for x in range(1000): c = int(c * (3 + s)) print(c) end = time.time()-startprint end`15.78 seconds
Code: [Select]
`from decimal import *import timestart = time.time()getcontext().prec = 10000s = 2 * (Decimal(2).sqrt())c = 1prec=10for x in range(1000): getcontext().prec = prec c = int(c * (3 + s)) print(c) prec +=10 end = time.time()-startprint end`12.04 seconds

#### Tomaaz

• Guest
##### Re: Praying Pythagorean triples
« Reply #13 on: January 22, 2016, 08:07:20 pm »
I agree, but still the most important thing is setting a proper precision for s. Now, your results are OK, but you must agree that 77 sec. for 100 numbers was a bit to long.

B+, how fast is your code?

#### B+

• Guest
##### Re: Praying Pythagorean triples
« Reply #14 on: January 22, 2016, 09:09:58 pm »
I agree, but still the most important thing is setting a proper precision for s. Now, your results are OK, but you must agree that 77 sec. for 100 numbers was a bit to long.

B+, how fast is your code?

B+ laughed loudly when he read Tomaaz question.

OK maybe you need a laugh too, I will time it.

EDIT: averages under 10.5 secs with scans and printing, cut the scans and print only last c, JB is down to just over 9 secs so scans and print added less than 1.5 secs. Oh so 77 secs is a bit long.

It is a good idea to adjust precision but the calls will suck up more time than the precision will save. In JB I could use standard float up to get first 20 Pray PTs right, that was 14 digit precision max, that flies, but of course smaller numbers as well.

Still if you want to see kick ass times, try sbmeTW's method of producing c's. I will test that now since this little game is started.

OK the results are in: 16 ms is longest time in 5 tests, it is either 15 or 16 or 0, the correct 100th c is found (no printing or checking in between).

Here is code stripped to bare bones (no legs calculations or checks or in loop printing)
Code: [Select]
`'100 Pray Pyth Tri using sbmeTW method.txt  for JB very modified for timed tests 2106-01-22t=time\$("ms")a=6 :b=5:c=1:s=5for p=1 to 100    s=(a*b)-c    c=b:b=snextprint time\$("ms")-tprint cend`
Can you see why I wanted to gets these results checked?
« Last Edit: January 22, 2016, 10:09:27 pm by B+ »