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 0<x<=n where x/radical(x)>1.
	
	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<x<=n and rad(x) <= sqrt(n). 
	The format of the data structure returned is that of a dictionary where the keys
	are in image of radical(x) for x in [2,sqrt(n)] and the item for each key, k, is 
	the set [x for x in range(1,n+1) if radical(x) = k].
	
	NOTES:  This is a modification of a previous function, piclist, in bigabc.py.  This 
	function depends on the function, withradlist, which given numbers n and radnum 
	returns a list of all numbers less than or equal to n whose radicals are equal 
	to the radical of radnum.  It also depends on radtuple.
	
	
	EXAMPLES:
	
	sage: root_rad_list(10)
	{2: [2, 4, 8], 3: [3, 9]}
	
	sage: root_rad_list(100)
	{2: [2, 4, 8, 16, 32, 64],
 	 3: [3, 9, 27, 81],
 	 5: [5, 25],
	 6: [6, 18, 54, 12, 36, 24, 72, 48, 96],
	 7: [7, 49],
	 10: [10, 50, 20, 100, 40, 80]}

	sage: root_rad_list(375)
	{2: [2, 4, 8, 16, 32, 64, 128, 256],
	 3: [3, 9, 27, 81, 243],
	 5: [5, 25, 125],
	 6: [6, 18, 54, 162, 12, 36, 108, 324, 24, 72, 216, 48, 144, 96, 288, 192],
	 7: [7, 49, 343],
	 10: [10, 50, 250, 20, 100, 40, 200, 80, 160, 320],
	 11: [11, 121],
	 13: [13, 169],
	 14: [14, 98, 28, 196, 56, 112, 224],
	 15: [15, 75, 375, 45, 225, 135],
	 17: [17, 289],
	 19: [19, 361]}
	
	AUTHOR:
		-Chris Hurlburt
		
	DATE:
		-May 11, 2006
	
	"""
	
	maxrad = integer_sqrt_floor(n)
	
	rad_image = radlist(maxrad)									#Generate the image of 
	odd = [x for x in range(2,maxrad+1) if x in rad_image]  	# the radical function.
	
	critic = { }
	
	for x in odd:
		critic[x] = withradlist(n,x)
		
	return critic
	
	
def abc_triple_list(m, n, outfile="DIPLOD.txt"):
	"""
	This function writes all good abc triples, (a,b,c), with m <= c <= n to a file 
	that is by default named DIPLOD.txt where each line in the file contains one triple.
	A good abc triple here is an order triple of three integers (a,b,c) with a<b<c such that 
		1) a,b,c are positive integers;
		2) a+b=c;
		3) gcd(a,b,c) = 1;
		4) c> 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<c]
					for i in tester:
						if testlist[c] > 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