Praying Pythagorean triples

[B+]

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-20
print" gives prayers"
print "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^"
   'Input "input  number list ";x
a=6 :b=5:c=1:s=5: D=1: z =8
for 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=s
next
end

[legs]
    La=int(((z*b)+D)/2)
    Lb=La+1
    TEMPB=La+Lb
    D=TEMP
    TEMP=TEMPB
return

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 outlined
c=1:ecount=0:pptcount=0
while 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
    print
wend
 

[Tomaaz]

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 sqrt
for 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:

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 sqrt
plik = 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 += 1
plik.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.

[B+]

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]

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+]

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 removed
aconst$="241421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206057147010955997160597027453459686201472851"
'print len(aconst$) :end

'use 200 digits to right of decimal = (sqr(2) + 1)*10^200
aconst$=left$(aconst$,201)
bconst=int(val(aconst$))
cconst=bconst*bconst              'avoid float math by doing cconst/dconst at the last calculation for new hypotenuse
dconst=10^200*10^200            ' cconst/dconst is just (SQR(2)+1)^2 to 200 places past decimal

c=1:ecount=0:pptcount=0
while 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
    scan
wend

Hypotenuse Table for PRAY PYTH TRI

Code: [Select]

1) 5
2) 29
3) 169
4) 985
5) 5741
6) 33461
7) 195025
8) 1136689
9) 6625109
10) 38613965
11) 225058681
12) 1311738121
13) 7645370045
14) 44560482149
15) 259717522849
16) 1513744654945
17) 8822750406821
18) 51422757785981
19) 299713796309065
20) 1746860020068409
21) 10181446324101389
22) 59341817924539925
23) 345869461223138161
24) 2015874949414289041
25) 11749380235262596085
26) 68480406462161287469
27) 399133058537705128729
28) 2326317944764069484905
29) 13558774610046711780701
30) 79026329715516201199301
31) 460599203683050495415105
32) 2684568892382786771291329
33) 15646814150613670132332869
34) 91196316011299234022705885
35) 531531081917181734003902441
36) 3097990175491791170000708761
37) 18056409971033565286000350125
38) 105240469650709600546001391989
39) 613386407933224037990008001809
40) 3575077977948634627394046618865
41) 20837081459758583726374271711381
42) 121447410780602867730851583649421
43) 707847383223858622658735230185145
44) 4125636888562548868221559797461449
45) 24045973948151434586670623554583549
46) 140150206800346058651802181530039845
47) 816855266853924917324142465625655521
48) 4760981394323203445293052612223893281
49) 27749033099085295754434173207717704165
50) 161733217200188571081311986634082331709
51) 942650270102046130733437746596776286089
52) 5494168403412088213319314492946575384825
53) 32022360150370483149182449211082676022861
54) 186639992498810810681775380773549480752341
55) 1087817594842494380941469835430214208491185
56) 6340265576556155474967043631807735770194769
57) 36953775864494438468860791955416200412677429
58) 215382389610410475338197708100689466705869805
59) 1255340561797968413560325456648720599822541401
60) 7316660981177400006023755031791634132229378601
61) 42644625325266431622582204734101084193553730205
62) 248551090970421189729469473372814871029093002629
63) 1448661920497260706754234635502788141981004285569
64) 8443420432013143050795938339643913980856932710785
65) 49211860671581597598021395402360695743160591979141
66) 286827743597476442537332434074520260478106619164061
67) 1671754600913277057625973209044760867125479123005225
68) 9743699861882185903218506820194044942274768118867289
69) 56790444570379838361685067712119508786523129590198509
70) 330998967560396844266891899452523007776864009422323765
71) 1929203360792001227239666329003018537874660926943744081
72) 11244221197191610519171106074565588219471101552240140721
73) 65536123822357661887786970118390510778951948386497100245
74) 381972521736954360807550714635777476454240588766742460749
75) 2226299006599368502957517317696274347946491584213957664249
76) 12975821517859256656937553191541868611224708916517003524745
77) 75628630100556171438667801831554937319401761914888063484221
78) 440795959085477771975069257797787755305185862572811377380581
79) 2569147124412310460411747744955171594511713413521980200799265
80) 14974086787388384990495417211933241811765094618559069827415009
81) 87275373599917999482560755526644279276078854297832438763690789
82) 508678154812119611904869115947932433844708031168435562754729725
83) 2964793555272799671946653940160950323792169332712780937764687561
84) 17280083176824678419775054525017769508908307965108250063833395641
85) 100715705505675270846703673209945666729657678457936719445235686285
86) 587014149857226946660446984734656230869037762782512066607580722069
87) 3421369193637686409115978235197991718484568898237135680200248646129
88) 19941201011968891508035422426453294080038375626640302014593911154705
89) 116225836878175662639096556323521772761745684861604676407363218282101
90) 677413820257085084326543915514677342490435733542987756429585398537901
91) 3948257084664334843320166936764542282180868716396321862170149172945305
92) 23012128687728923975594457705072576350594776564834943416591309639133929
93) 134124515041709209010246579293670915821387790672613338637377708661858269
94) 781734961562526330085885018056952918577731967470845088407674942332015685
95) 4556285254333448771505063529048046595645004014152457191808671945330235841
96) 26555976564438166298944496156231326655292292117443898062444356729649399361
97) 154779574132295549022161913408339913336108748690510931182857468432566160325
98) 902121468229335127834026984293808153361360200025621689034700453865747562589
99) 5257949235243715217981999992354509006832052451463219203025345254761919215209
100) 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.

[Tomaaz]

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 = 100000
s = 2 * (Decimal(2).sqrt())
c = 1
for 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?

[jbk]

Hi Tomaaz
I timed your program for x in range(100) 141.53 seconds
then made a small change to your program

Code: [Select]

from decimal import *
import time
start = time.time()
getcontext().prec = 100000
s = 2 * (Decimal(2).sqrt())
c = 1
prec=10
for x in range(100):
getcontext().prec = prec
c = int(c * (3 + s))
print(c)
prec +=10

end = time.time()-start
print 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

[Tomaaz]

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 time
start = time.time()
getcontext().prec = 1000
s = 2 * (Decimal(2).sqrt())
c = 1
for x in range(100):
c = int(c * (3 + s))
print(c)
end = time.time()-start
print end

Shouldn't take more than a second.

[B+]

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

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+]

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.

[jbk]

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.

[Tomaaz]

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]

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 time
start = time.time()
getcontext().prec = 10000
s = 2 * (Decimal(2).sqrt())
c = 1
for x in range(1000):
c = int(c * (3 + s))
print(c)

end = time.time()-start
print end
15.78 seconds

Code: [Select]

from decimal import *
import time
start = time.time()
getcontext().prec = 10000
s = 2 * (Decimal(2).sqrt())
c = 1
prec=10
for x in range(1000):
getcontext().prec = prec
c = int(c * (3 + s))
print(c)
prec +=10

end = time.time()-start
print end
12.04 seconds

[Tomaaz]

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+]

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-22
t=time$("ms")
a=6 :b=5:c=1:s=5
for p=1 to 100
    s=(a*b)-c
    c=b:b=s
next
print time$("ms")-t
print c
end

Can you see why I wanted to gets these results checked?