### Experimental Math Spring 2012 Project [This code is completely open to the public; please credit when appropriate.] ### Pat Devlin --- 1-May-2012 (much of the code was originally written 28-March-2012) ## Maple supplement for "Primes Appearing in Prime Tower Factorization" ## --- Automatically addressing (and unfortunately disproving) a beautiful conjecture due to Edinah Gnang --- ## This script is to estimate a certain number theory sort of probability question posed by Edinah Gnang. ## Any natural number n can be writen as a prime factorization n = p1 ^ e1 * p2 ^e2 * ... * pk ^ek, (no pi = 1) # After this, write each exponent that appears here (e1, e2, ... , ek) ALSO # as a prime factorization: # # n = p1 ^ (p_1,1 ^ e_1,1 * p_1,2 ^ e_1,2 * ... * p_1,j1 ^ e_1,j1) * # p2 ^ ( ...) * pk ^ (...). # # Then factor each exponent, e_1,1 [et cetera], into prime factors and # continue until you have written the number as # n = PRODUCT( p_k ^ (PRODUCT p_(k,j) ^ PRODUCT( ... ) ) ). # ## Examples! # 144 = 2^4 * 3^2 = 2^(2^2) * 3^2. # 2^48 * 3^20 * 5^144 = 2^(2^4 * 3) * 3^(2^2 * 5) * 5 ^(2^4 * 3^2) = 2^(2^(2^2) * 3) * 3^(2^2 * 5) * 5^(2^(2^2) * 3^2) # ## Now let q be any prime number ## Then if you pick a random natural number, what is the probability that its "complete tower factorization" contains the prime q? ## Let d(q) denote this number. ## ## In other words, let M(q) denote the set of natural numbers such that the "complete tower factorization" contains the prime q. ## Then d(q) := lim_{N to infinity} |M(q) \cap {1, 2, 3, ... , N}| / N. # # ## See the arxiv paper with myself and Edina Gnang titled ## Primes Appearing in Prime Tower Factorization (http://arxiv.org/abs/1204.5251) ## for more information and results on this problem. ## ## More recent developed papers will be available on the arxiv when I get the chance to clean them up. # Throughout, the set M^c (q) will be {z in {0, 1, 2, ...} : z not in M(q)}. # In particular, note {0,1} is a subset of M^c (q). Help:=proc(): print(`Maple supplement for "Primes Appearing in Prime Tower Factorization"`): print(`--- Automatically addressing (and unfortunately disproving) a beautiful conjecture due to Edinah Gnang ---`): print(` `): print(`This is the Experimental Math (Spring 2012) project of Pat Devlin!`): print(`1-May-2012`): print(` `): print(`Type HelpGeneral() to get a blurb about the project`): print(`Type HelpFunctions() to see a list of functions`): print(`Type HelpAbout([function name in quotes]) to find about a specific function`): print(`~~~For example: type HelpAbout("containsPrimeInTower") to find out about that function`): end proc: HelpFunctions:=proc(): print(`containsPrimeInTower(q,n)`): print(`firstElementsInM(q, n), firstElementsNotInM(q, n)`): print(`upperBoundGivenStuffNotInM(q, A, P, longOutput:=true), lowerBoundGivenStuffInM(q, S, P, longOutput:=true)`): print(`boundWrapper(q, a, s, p, longOutput:=true)`): print(`Help(), HelpGeneral(), HelpAbout(functionNameInQuotes), HelpFunctions()`): end proc: ## HelpGeneral() and HelpAbout(...) are at the end of this file for clarity ## containsPrimeInTower(q, n) # This inputs q and an integer n and returns true if and only if the number # contains the prime q in its "complete tower prime factorization" # [i.e., it returns true iff n is in M(q)] # # For conveneince, it returns false if n<=1 [that is 1 and 0 "do not" have # any primes in their tower factorization] containsPrimeInTower:=proc(q,n) option remember: local i, L, currentExponent: if ((n <= 1) or (n < q) or (q < 2)) then return false: fi: if (type(n/q, integer)) then return true: fi: L:=ifactors(n)[2]; for i from 1 to nops(L) do: currentExponent:=L[i][2]: if(containsPrimeInTower(q, currentExponent)) then return true: fi: od: return false: end proc: ## firstElementsInM(q, n) # This inputs q and n, and it finds and outputs the set {0, 1, 2, 3, ..., n} # intersect M(q) # [that is, it outputs the set of natural numbers <= n that ARE in M(q)] firstElementsInM:=proc(q,n) local M, i: M := {}: for i from q to n do: if(containsPrimeInTower(q, i)) then M := M union {i}: fi: od: return M: end proc: ## firstElementsNotInM(q, n) # This inputs q and n, and it finds and outputs the set {0, 1, 2, 3, ..., n} # minus M(q) # [that is, it outputs the set of natural numbers <= n that are NOT in M(q)] firstElementsNotInM:=proc(q,n) local M, i: M := {seq(i, i=0..q-1)}: for i from q to n do: if(not containsPrimeInTower(q, i)) then M := M union {i}: fi: od: return M: end proc: ## upperBoundGivenStuffNotInM(q, A, P) # If A is a subset of M^c (q) [and A might as well contain {0, 1, 2, ... , q-1}] # and P is any subset of primes with q not in P, # Then we have the rigorous bound # # d(q) <= 1 - (q^q - q^(q-1) )/( q^q - 1) * 1/Zeta(q) * prod( (1- 1/p)/(1- # 1/p^q) * sum(1/p^a, a in A), p in P) # # This outputs that upper bound on d(q) upperBoundGivenStuffNotInM:=proc(q, A, P, longOutput:=true) local A2, P2, i, p, a, theAnswer: A2 := A union {seq(i, i=0..(q-1))}: P2 := P minus {q}: theAnswer:=1 - (q^q - q^(q-1) )/( q^q - 1) * 1/Zeta(q) * mul( (1-1/p)/(1-1/p^q) * add(1/p^a, a in A2), p in P2): if(longOutput) then printf("d(%d) <= %g\n", q, convert(theAnswer, float)): fi: return theAnswer: end proc: ## lowerBoundGivenStuffInM(q, S, P) # If S is a subset of M(q), and P is any subset of primes with q not in P, # # Then we have the rigorous bound # # d(q) >= 1 - (1 - 1/q) * prod( 1 - (1-1/p) * sum(1/p^s, s in S), p in P) # # This outputs that lower bound on d(q) given S and P (as a rational) lowerBoundGivenStuffInM:=proc(q, S, P, longOutput:=true) local P2, p, s, theAnswer: P2:= P minus {q}: theAnswer:=1 - (1 - 1/q) * mul( 1 - (1-1/p) * add(1/p^s, s in S), p in P2): if(longOutput) then printf("d(%d) >= %g\n", q, convert(theAnswer, float)): fi: return theAnswer: end proc: ## boundWrapper(q, a, s, p) # This is a convenient wrapper function to find rigorous # upper and lower bounds for d(q) simultaneously. # It inputs q and three positive integers, a, s, and p, which determine # the size of A, the size of S, and how many primes to use (respectively) # # Larger values of 'a' tighten the upper bound, # and larger values of 's' tighten the lower bound. # # It seems that using larger values of p is a more effective way to tighten either bound # # This returns a list [lowerBound, upperBound], which are completely rigorous # (usually very large!) rational numbers. Using rationals prevents any rounding error. boundWrapper:=proc(q, a, s, p, longOutput:=true) local A, S, P, i, lowerBound, upperBound: A:=firstElementsNotInM(q, a): S:=firstElementsInM(q, s): P:={seq(ithprime(i), i=1..p)}: lowerBound:=lowerBoundGivenStuffInM(q, S, P, false): upperBound:=upperBoundGivenStuffNotInM(q, A, P, false): if(longOutput) then printf("%g <= d(%d) <= %g\n", convert(lowerBound,float), q, convert(upperBound, float)): fi: return [lowerBound, upperBound]: end proc: ## HelpAbout(functionName:=""): # Displays detailed information about the function name given HelpAbout:=proc(functionName:=""): print(` `): if(functionName = "containsPrimeInTower") then print(`containsPrimeInTower(q,n):\n ~~ This inputs q and an integer n and returns true if and only if n contains the prime q in its "complete tower prime factorization".\n ~ [i.e., it returns true iff n is in M(q)]\n ~~ For conveneince, it returns false if n<=1 [that is 1 and 0 "do not" have any primes in their tower factorization].`): return: elif(functionName = "firstElementsInM") then print(`firstElementsInM(q, n):\n ~~ This inputs q and n, and [using containsPrimeInTower] it finds and outputs the set {0, 1, 2, 3, ..., n} intersect M(q).\n ~ [that is, it outputs the set of natural numbers (including 0) <= n that are NOT in M(q)].`): return: elif(functionName = "firstElementsNotInM") then print(`firstElementsNotInM(q, n):\n ~~ This inputs q and n, and [using containsPrimeInTower] it finds and outputs the set {0, 1, 2, 3, ..., n} minus M(q).\n ~ [that is, it outputs the set of natural numbers (including 0) <= n that are NOT in M(q)].`): return: elif(functionName = "upperBoundGivenStuffNotInM") then print(`upperBoundGivenStuffNotInM(q, A, P, longOutput:=true):\n ~~ The Devlin-Gnang paper proves that if A is a subset of M^c (q) [WLOG, assume {0, 1, 2, ..., q-1} in A],\n ~ and P is any subset of primes with q not in P, then we have the rigorous upper bound:\n ~ d(q) <= 1 - (q^q - q^(q-1) )/( q^q - 1) * 1/Zeta(q) * prod((1-1/p)/(1-1/p^q) * sum(1/p^a, a in A), p in P)\n ~~ This function inputs these parameters q, A, and P, and outputs this upper bound on d(q) (as a rational comabination a+b/Zeta(q)).\n ~ There is also an optional fourth boolean argument, which determines if the long output format should be used (default is true).`): return: elif(functionName = "lowerBoundGivenStuffInM") then print(`lowerBoundGivenStuffInM(q, S, P, longOutput:=true):\n ~~ The Devlin-Gnang paper proves that if S is a subset of M(q), and P is any subset of primes with q not in P,\n ~ then we have the rigorous lower bound:\n ~ d(q) >= 1 - (1 - 1/q) * prod( 1 - (1-1/p) * sum(1/p^s, s in S), p in P)\n ~~ This function inputs these parameters q, S, and P, and outputs this lower bound on d(q) (as a rational).\n ~ There is also an optional fourth boolean argument, which determines if the long output format should be used (default is true).`): return: elif(functionName = "boundWrapper") then print(`boundWrapper(q, a, s, p, longOutput:=true):\n ~~ This is a convenient wrapper function to find rigorous upper and lower bounds for d(q) simultaneously.\n ~~ Its inputs are the prime of interest, q, and three positive integers, a, s, and p, which determine\n ~ the size of A, the size of S, and how many primes to use (respectively).\n ~~ Larger values of 'a' tighten the upper bound, and larger values of 's' tighten the lower bound.\n ~~ It seems that using larger values of 'p' is a more effective way to tighten either bound.\n ~~ This returns a list [lowerBound, upperBound], which are completely rigorous\n ~ (usually very large!) rational numbers. Using rationals prevents any rounding error.\n ~ There is also an optional fifth boolean argument, which determines if the long output format should be used (default is true).`): return: elif(functionName = "Help") then print(`Help():\n ~~ Lists the various Help functions`): return: elif(functionName = "HelpGeneral") then print(`HelpGeneral():\n ~~ Provides a blurb about this script as a whole`): return: elif(functionName = "HelpAbout") then print(`HelpAbout(functionNameInQuotes):\n ~~ Displays detailed information about the function name given`): elif(functionName = "HelpFunctions") then print(`HelpFunctions():\n ~~ Lists all functions that can be called including their parameters.`): return: elif(not (functionName = "")) then print(`Sorry, that's not a valid function name`): fi: print(` `): print(`Valid options for HelpAbout(...) are:\n HelpAbout("containsPrimeInTower"); HelpAbout("firstElementsInM"); HelpAbout("firstElementsNotInM");\n HelpAbout("upperBoundGivenStuffNotInM"); HelpAbout("lowerBoundGivenStuffInM"); HelpAbout("boundWrapper");\n HelpAbout("Help"); HelpAbout("HelpGeneral"); HelpAbout("HelpAbout"); HelpAbout("HelpFunctions");`): end proc: ## HelpGeneral() # Displays a general blurb about this (tiny) Maple package HelpGeneral:=proc(): print(`Maple supplement for "Primes Appearing in Prime Tower Factorization"`): print(`--- Automatically addressing (and unfortunately disproving) a beautiful conjecture due to Edinah Gnang ---`): print(` `): print(`This is the Experimental Math (Spring 2012) project of Pat Devlin!`): print(`1-May-2012`): print(` `): print(`\n ~~ This script is to estimate a certain number theory sort of probability question posed by Edinah Gnang.\n\n ~ Any natural number n can be writen as a prime factorization n = p1 ^e1 * p2 ^e2 * ... * pk ^ek, (no pi = 1).\n ~ After this, write each exponent that appears here (e1, e2, ... , ek) ALSO as a prime factorization:\n ~ n = p1 ^ (p_1,1 ^ e_1,1 * p_1,2 ^ e_1,2 * ... * p_1,j1 ^ e_1,j1) * p2 ^ ( ...) * pk ^ (...).\n ~ Then factor each exponent, e_1,1 [et cetera], into prime factors and continue until you have written the number as\n ~ n = PRODUCT( p_k ^ (PRODUCT p_(k,j) ^ PRODUCT( ... ) ) ).\n\n ~~ Examples!\n ~ 144 = 2^4 * 3^2 = 2^(2^2) * 3^2.\n ~ 2^48 * 3^20 * 5^144 = 2^(2^4 * 3) * 3^(2^2 * 5) * 5 ^(2^4 * 3^2) = 2^(2^(2^2) * 3) * 3^(2^2 * 5) * 5^(2^(2^2) * 3^2)\n\n ~~ Now let q be any prime number. Then if you pick a random natural number,\n ~ what is the probability that its "complete tower factorization" contains the prime q?`): print(`~~ Let d(q) denote this number. ~~`): print(`\n ~ In other words, let M(q) denote the set of natural numbers such that the "complete tower factorization" contains the prime q.\n ~ Then what is d(q) := lim_{N to infinity} |M(q) \cap {1, 2, 3, ... , N}| / N ??`): print(`The (beautiful) conjecture due to Edinah Gnang was that d(2) = 1/sqrt(3).`): print(`Using some human mathematics, we have the algorithms of this script to rigorously get upper and lower bounds d(q).`): print(`Unfortunately, a call to boundWrapper(2, 20, 20, 25000) disproves Gnang's conjecture.`): print(`See the arxiv paper by Edinah Gnang and myself:\n "Primes Appearing in Prime Tower Factorization" (http://arxiv.org/abs/1204.5251), for more details\n and results on this problem, as well as human math that explains why on earth we're running these scripts!`): end proc: