124-Invoer van parameters.sagews

KIES EERST HIERBOVEN "Action" OPTIE "Evaluate All"

 

Uitleg ABC algoritme

De opzet is overgenomen van Reken mee met ABC:

• een groot getal NN wordt genomen
• in 5 stappen wordt getracht alle ABC-triples te vinden met radiaal kleiner dan NN
• men gebruikt het feit dat radiaal(A), radiaal(B) en radiaal(C) verschillende getallen zijn!
• dat maakt dat je deze drie radiaal-getallen in volgorde kunt plaatsen: r(x)<r(y)<r(z)
• waaruit direct volgt dat 1 < r(y) < sqrt(NN) en r(x) < NN/r(y)^2.

De 5 stappen zijn:

1. bij het getal NN bepaal alle mogelijke waarden voor r(y)
2. bij de waarden r(y) bepaal alle mogelijke waarden voor y
3. bij de waarden r(y) bepaal alle mogelijke waarden voor r(x)
4. bij de waarden r(x) bepaal alle mogelijke waarden voor x
5. bij elke waarde van x en y zijn er twee mogelijke uitkomsten: {x,y,c=x+y} of {|x-y|,min(x,y),max(x,y)}

Tenslotte gaan we na of de gevonden drietallen ABC-triples zijn.

Gebruik deze link hier om abc-triples te zoeken

 

Uitbreiding

Een slimme uitbreiding is om elk relatie A+B=C te vermenigvuldigen met respectievelijk (C+A), (C+B) of (B-A). Dat geeft de volgende kansrijke drietallen:

• A²+B*(C+A)=C²
• B²+A*(C+B)=C²
• A²+C*(B-A)=B²

Als je wilt, worden ook deze drietallen worden doorgerekend.

# Default bovengrens tot waar we abc-triples zoeken (o.a. handig voor testen)
NN = 400

# Om het programma sneller te laten zijn wordt eenmalig de lijst van mogelijke radialen bepaald
radiaal_range=[k for k in range(1,1000+1) if is_squarefree(k)]     #608 stuks

print radiaal_range[:80]

[1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38, 39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73, 74, 77, 78, 79, 82, 83, 85, 86, 87, 89, 91, 93, 94, 95, 97, 101, 102, 103, 105, 106, 107, 109, 110, 111, 113, 114, 115, 118, 119, 122, 123, 127, 129, 130]

#REKENPROGRAMMA's

#    Hulpfunctie om maximale exponent snel te bepalen
var('x,y')
f_exp=fast_callable(floor(log_b(x,y)+1),vars=[x,y],domain=float)

#----------------------------------------------------------------------------------------------------------

def mogelijk(reeks,maxum=400):

#    COMMENTAAR
#    We bepalen nu alle y =< maxum met y opgebouwd uit machten uit reeks

 if reeks in [[],[1]] or 1 in reeks: return [1]
 uitkomst=[1]
 for p in reeks: uitkomst=uitkomst+[k*p^e for k in uitkomst for e in xrange(1,int(f_exp(maxum//k,p)))]
 if maxum in uitkomst: return sorted(uitkomst[1:])[:-1]
 else: return sorted(uitkomst[1:])

#----------------------------------------------------------------------------------------------------------

def uitbreiding(triples):

#    COMMENTAAR
#    Een gevonden abc-triple a+b=c wordt vermenigvuldigd met (c+a),(c+b) of (b-a) om te zien of dat een #    nieuw abc drietal geeft. We gebruiken hier meer rekenwerk dan strikt noodzakelijk.

 uitkomst=[]
 for (p,q) in triples:
  b=max(p,q)
  a=min(p,q) 
  c=a+b
  if prod(prime_divisors(a*b*c*(b-a)))<b^2: uitkomst.append((a^2,(b-a)*c))
  if prod(prime_divisors(a*b*c*(c+a)))<c^2: uitkomst.append((a^2,(c+a)*b))
  if prod(prime_divisors(a*b*c*(c+b)))<c^2: uitkomst.append((b^2,(c+b)*a))
 return uitkomst

#----------------------------------------------------------------------------------------------------------

#HOOFDPROGRAMMA REKENEN

def abc_vinden(nmax):

#    COMMENTAAR
#    We bepalen nu in 5 stappen alle abc-triples. Getracht is effici gebruik te maken van wat we al
#    weten (i.h.b. r(y) en r(x)) en niet zaken opnieuw door te rekenen.

 abc_triples=[]
 ry=[k for k in radiaal_range if 1<k<isqrt(nmax)]                                 #stap 1
 for ny in ry: 
  yw=mogelijk(prime_divisors(ny),nmax)                                            #stap 2
  rx=[n for n in radiaal_range if 0<n<nmax/ny^2 and n<ny and gcd(n,ny)==1]        #stap 3
  for nx in rx:
   xw=mogelijk(prime_divisors(nx),nmax)                                           #stap 4
   for y in yw:
    for x in xw:
     #Geval x+y=c                                                                 #stap 5
     if nx*ny*prod(prime_divisors(x+y))<x+y: abc_triples.append((min(x,y),max(x,y)))
     #Geval |x-y|=a of b
     if nx*ny*prod(prime_divisors(abs(x-y)))<max(x,y): abc_triples.append((min(x,y),abs(x-y)))
 return uniq(abc_triples)

# TEST mogelijk

print mogelijk([2,3,5])

[2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160, 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384]

# TEST uitbreiding

print uitbreiding([(1, 8), (1, 48), (1, 63), (1, 80), (5, 27), (8, 1), (27, 5), (32, 49)])

[(1, 63), (1, 80), (1, 2303), (1, 2400), (1, 3968), (1, 4095), (1, 6399), (1, 6560), (6400, 161), (25, 704), (1, 63), (1, 80), (25, 704), (1024, 1377), (1024, 5537), (2401, 4160)]

# TEST abc_vinden

abc_vinden(100)

[(1, 8), (1, 48), (1, 63), (1, 80), (5, 27), (8, 1), (27, 5), (32, 49)]

#PROGRAMMA's VOOR UITVOER

"""  COMMENTAAR 
     Deze routine is bijna gelijk aan het hoofdprogramma, en is vanwege de overzichtelijkheid apart
     gehouden. Anders zijn:
     - een begrenzing voor de zoekopdracht (anders teveel uitvoer)
     - printopdrachten voor elke stap
     - uitbreiding met de truc van Ismael Jimz Calvo
"""

#----------------------------------------------------------------------------------------------------------

def abc_stappen(nmax,stap):

#    COMMENTAAR
#    We bepalen nu (desgewenst) in 5 stappen alle abc-triples en geven extra analyse informatie

 [teller_ry,teller_yw,teller_rx,teller_xw,teller_abc]=[0]*5
 if stap>1: nmax=min(400,nmax)
 if nmax>400: nmax=400                                                            #Maximum gesteld
 print "r(y)<sqrt(%d) ofwel 1 < r(y) < %d"%(nmax,isqrt(nmax))                     #Extra informatie
 print ""
 abc_triples=[]
 ry=[k for k in radiaal_range if 1<k<isqrt(nmax)]                                 #stap 1
 teller_ry=len(ry)
 if stap>=1: print "1: Mogelijke radialen voor y:",ry," [r(y)]"
 for ny in ry:
  yw=mogelijk(prime_divisors(ny),nmax)                                            #stap 2
  teller_yw=teller_yw+len(yw)
  if stap>=2: print "2: Mogelijke y bij r(y)=%d:"%ny,yw
  rx=[n for n in radiaal_range if 0<n<nmax/ny^2 and n<ny and gcd(n,ny)==1]        #stap 3
  teller_rx=teller_rx+len(rx)
  if stap>=3: print "3: Mogelijke radialen voor x bij r(y)=%d:"%ny,rx
  for nx in rx:
   xw=mogelijk(prime_divisors(nx),nmax)                                           #stap 4
   teller_xw=teller_xw+len(xw)
   if stap>=4: print "4: Mogelijke x bij r(x)=%d:"%nx,xw
   for y in yw:
    for x in xw:
     #Geval x+y=c                                                                 #stap 5
     if nx*ny*prod(prime_divisors(x+y))<x+y: abc_triples.append((min(x,y),max(x,y)))
     #Geval |x-y|=a of b
     if nx*ny*prod(prime_divisors(abs(x-y)))<max(x,y): abc_triples.append((min(x,y),abs(x-y)))
     teller_abc=teller_abc+2
 print "" 
 return uniq(abc_triples),[teller_ry,teller_yw,teller_rx,teller_xw,teller_abc]

#----------------------------------------------------------------------------------------------------------

# TEST abc_stappen

(drietallen,[teller_ry,teller_yw,teller_rx,teller_xw,teller_abc])=abc_stappen(60,1)

print '#ry:',teller_ry,'#y:',teller_yw,'#rx:',teller_rx,'#x:',teller_xw,'#3-tallen:',teller_abc
print ''
print drietallen

r(y)<sqrt(60) ofwel 1 < r(y) < 7 1: Mogelijke radialen voor y: [2, 3, 5, 6] [r(y)] #ry: 4 #y: 25 #rx: 6 #x: 14 #3-tallen: 100 [(1, 8), (1, 48), (5, 27), (8, 1), (27, 5)]

# HOOFPRGRAMMA INVOER

html('<A name="start"></A>')

html('<h1>Invoer van parameters</h1>')


@interact

def intake(stap=input_box(default=1,label='Weergeven tot stap', width=1),nmax=input_box(default=NN,label='Invoer maximum (vanaf stap 2 beperkt tot 400)', width=2),alle_resultaten=("Alle resultaten weergeven?",True),uitbreiding_keuze=("Uitbreiding meenemen?",False),pythagoras_keuze=("Check of Pythagoras drietal?",False)):
      
   html('</>')
   (drietallen,[teller_ry,teller_yw,teller_rx,teller_xw,teller_abc])=abc_stappen(nmax,stap)
   print "Aantal r(y)      :",teller_ry
   print "Aantal y         :",teller_yw
   print "Aantal r(x)      :",teller_rx
   print "Aantal x         :",teller_xw
   print "Aantal drietallen:",teller_abc
   print "--------------"
   print drietallen
   print ""
   if uitbreiding_keuze: 
    print uitbreiding(drietallen)
    print ""
   if pythagoras_keuze: print "Pythagoras drietallen",[k for k in drietallen if k[0]^2+k[1]^2==(k[0]+k[1])^2]
   html('</>')

###########