def radical(n):
"""
Return the radical of n which is the product of all of the primes that
divide n unless n is 0, 1 or -1. Then it will return 0, 1, and 1
respectively.
EXAMPLES:
sage: radical(44)
22
Here we will also need to check the behavior with respect to negative
numbers.
sage: radical(-1)
1
sage: radical(0)
0
AUTHORS:
- Chris Hurlburt through total imitation of
William Stein and Alex Clemesha (2006-01-10): some examples
"""
if n==0:
return 0
return misc.mul([p for p, r in factor(n)])
def integer_sqrt_floor(n):
"""
Given n, return the floor(sqrt(n)) as an integer data type. This is a hack to
deal with the fact that python is a strongly typed language and that the python
types and the sage types get mixed. This function therefore will return a python
integer which is the floor (ie greatest integer less than or equal to a rational
number) of an imputed integer n. Obviously since square roots are involved, this
function does not work for negative integers.
EXAMPLES:
sage: integer_sqrt_floor(16)
4
sage: integer_sqrt_floor(516)
22
sage: integer_sqrt_floor(519)
22
sage: integer_sqrt_floor(1000)
31
AUTHORS:
-Chris Hurlburt
DATE:
May 10, 2006
"""
return int(sqrt(int(n)*log(2,2)).floor())
def radlist(n):
"""
This returns an array s such that the ith entry of s, s[i] = radical(i). This
will be more clear by looking at the following examples:
EXAMPLES:
sage: radlist(1)
[0, 1]
sage: radlist(2)
[0, 1, 2]
sage: radlist(9)
[0, 1, 2, 3, 2, 5, 6, 7, 2, 3]
sage: radlist(15)
[0, 1, 2, 3, 2, 5, 6, 7, 2, 3, 10, 11, 6, 13, 14, 15]
sage: radlist(22)
[0, 1, 2, 3, 2, 5, 6, 7, 2, 3, 10, 11, 6, 13, 14, 15, 2, 17, 6, 19, 10, 21, 22]
AUTHORS:
-Chris Hurlburt
"""
s = [0] #Generate a list where s[0] = 0
for k in range(n): # and s[i] = 1 for i>0.
s.append(1)
i = 2
while i <= n: #If an entry still has a 1 in it when
if s[i] == 1: # this routine gets to it, then it is
k=i # a prime and so its radical must contain
while k<=n: # a single occurance of this prime as must
s[k] = i*s[k] # all multiples of it. Therefore we go
k = k+i # through the list of all multiples and multiply
# by a single occurance of the prime which is
i = i+1 # the index.
return s
def radtuple(n):
"""
This returns a python dictionary, fry, where the key is x and the entry for
x is radical(x), ie
{x:radical(x)}
for all 01.
NOTE: The entire reason for creating this version as opposed to radlist is
to conserve memory. For example if n = 10^6, then radlist(n) has 10^6 + 1
entries but radtuple(n) has 392074 entries, less than half the number of entries
in radlist(n).
EXAMPLES:
sage: radtuple(10)
{8: 2, 9: 3, 4: 2}
sage: radtuple(20)
{4: 2, 8: 2, 9: 3, 12: 6, 16: 2, 18: 6, 20: 10}
sage: radtuple(30)
{4: 2, 8: 2, 9: 3, 12: 6, 16: 2, 18: 6, 20: 10, 24: 6, 25: 5, 27: 3, 28: 14}
AUTHORS:
-Chris Hurlburt
"""
fry = { } #Note this is the notation for a dictionary list, not just a set
i = 2
while i<= n:
if i not in fry:
k = 2*i
while k<=n:
if k in fry:
fry[k] = i*fry[k]
if k not in fry:
fry[k] = i
k = k+i
if i in fry:
if fry[i] == i:
del fry[i]
i = i+1
return fry
def lizardList(m,n):
"""
This is a test function to generate a python dictionary
where the key is x and the entry for x is 1, ie {x:1}
for all m<=x<=n where x/radical(x)>1.
NOTES: The function in fact works. We do add the caveat though
that it is necessary for m-n to be less than 10^8 or 10^9. Namely, this
function will die when the limit of xrange/prange in python is exceeded.
This function also requires the integer_sqrt_floor(n) function.
The name is lizardList because it is not meant to be used as a stand alone
routine. The idea behind this function is the fact that if x>rad(x), then
x is a multiple of p^2 for some prime p that is <= sqrt(n).
EXAMPLES:
sage: lizardList(1,20)
{4: 1, 8: 1, 9: 1, 12: 1, 16: 1, 18: 1, 20: 1}
sage: lizardList(2050,2065)
{2050: 1, 2052: 1, 2056: 1, 2057: 1, 2058: 1, 2060: 1, 2061: 1, 2064: 1}
lizardList(20,48)
{32: 1, 36: 1, 40: 1, 44: 1, 45: 1, 48: 1, 20: 1, 24: 1, 25: 1, 27: 1, 28: 1}
AUTHOR:
-Chris Hurlburt
"""
ending = integer_sqrt_floor(n)+1
coverage = prange(ending)
fryer = { }
for p in coverage:
mplier = p**2
start = m//mplier + 1
stop = n//mplier + 1
distance = stop - start #This is necessary to not exceed the maximum values
# in range/xrange of somewhere just under 10^9.
for k in xrange(distance):
fryer[mplier*(start + k)] = 1
if mod(m,mplier) == 0:
fryer[m] = 1
return fryer
def withradlist(n,radnum):
"""
This returns a list of numbers whose values are less than n and whose radical is the
same as the radical of radnum.
AUTHOR:
-Chris Hurlburt
DATE:
-March 20, 2006
"""
if radnum == 1:
return [1]
gecko = factor(radnum)
curnum = gecko[0][0]
stop = int(log(n*1.0,curnum).floor())
salamander = [curnum**i for i in range(1,stop+1)]
next = radnum/(gecko[0][0]**(gecko[0][1]))
rlist = [x*y for x in salamander for y in withradlist(n,next) if x*y <= n]
return rlist
def root_rad_list(n):
"""
This generates the set of all positive numbers x where 1 radical(abc).
NOTE: This function overwrites whatever file name is given each time it is called.
AUTHORS:
-Chris Hurlburt
DATE: May 13, 2006
"""
checklist = root_rad_list(n) #Generate a dictionary indexed by the rad value
# for all rad values less than or equal to sqrt(n)
# where each item is a list of elements <=n with
# that rad value. This will be used to generate
# the set of possible a/b values that need to
# be checked for a given possible c value.
f = open(outfile, 'w') #Open the file given by the outfile parameter which
# has a default value of DIPLOD.txt for writing.
notdone = 1 #This is because if we check more than 10^5 numbers
start = m # at once the process slows significantly do to
stop = min(start + 10**5,n) # exhaustion of memory. Hence we break the
# checking process into segments.
while notdone:
gecko = lizardList(start,stop) #Generate a dictionary {x:1} of numbers with x>rad x
for x in gecko: #Generate testlist from gecko where is a dictionary
gecko[x] = x/radical(x) # such that 1) testlist[x] = x/radical(x) and
# testlist[x]>2 for all x in testlist.
testlist = { }
tail = gecko.keys()
for k in tail:
if gecko[k]>2:
testlist[k] = gecko[k]
tail = testlist.keys()
tail.sort()
gecko.clear() #This is to clear up memory
for c in tail:
if testlist[c] > radical(c-1): #Check the 1 case.
stuff = str(1) + '\t' + str(c-1) + '\t' + str(c) + '\n'
f.write(stuff)
if testlist[c] not in [3,4,5,6,8,9,10,12,14,18,20,24,30,42,60]: #For c/rad(c) in this set one
top = integer_sqrt_floor(testlist[c]) # case is only possible case
abPossibleRad = [x for x in checklist if x<=top if gcd(x,c) == ZZ(1)]
for radval in abPossibleRad:
tester=[x for x in checklist[radval] if x radval*radical(c-i):
if i < c-i:
stuff = str(i) + '\t' + str(c-i) + '\t' + str(c) + '\n'
f.write(stuff)
else:
if radical(c-i) not in abPossibleRad + [1]: #Must check that we didn't already write this triple in checking
stuff = str(c-i) + '\t' + str(i) + '\t' + str(c) + '\n' #c-i. If c-i were checked it would either
f.write(stuff) #have to be 1 or have a radical in abPossibleRad
if stop == n: #If we haven't hit n yet get the next 10^5 numbers or
notdone = 0 # the rest of the numbers to n.
else:
start = stop +1
stop = min(start + 10**5,n)
f.close()
return gecko