#HW2+HW3, Jan. 26, 2024, OK TO SHARE   
#For reasons of class pride I found it more fitting to use the homemade versions of functions   #rather than the maple versions, so I have copy and pasted them here below    
#Also, I was having some trouble converting to a txt file so this all was copy and pasted
#and then manually edited, so there might be some weird formatting at times
#I think the code is mostly intact but let me know if you want me to send the actual 
#worksheet file since that might be easier



#NextPrime1(n): inputs a pos. integer n  and outputs the first prime >=n  
#NextPrime1:=proc(n) local i:    
for i from n while not isprime(i) do od:  
i:  
end:    

#MakeRSAkey: The key [n,e,d]: a->a^e mod n  n must be a product of two primes  inputs D1 and outputs #an RSA key  [n,e,d], where n is a product of two primes  with D1 digits. Try:  MakeRSAkey(100);  MakeRSAkey:=proc(D1) local n,d,S,m,p,q,e:  

p:=NextPrime1(rand(10^(D1-1)..10^D1-1)()):  
q:=NextPrime1(rand(10^(D1-1)..10^D1-1)()):  
if p=q then    
     RETURN(FAIL):  
fi:  

n:=p*q:  
m:=(p-1)*(q-1):  
S:=rand(m/2..m-1)():  
for e from S to m-1 while gcd(e,m)>1 do od: 
d:=e&^(-1) mod m:    

[n,e,d]:  

end:   
#m>n m=n*q+r (r= m mod n)  
#gcd(m,n)=gcd(n,r)    
#gcd(n,0)=n  
#EA(m,n): inputs m>n and outputs the gcd(m,n)  
EA:=proc(m,n) local q,r:  

if n=0 then   RETURN(m):  fi:  
r:=m mod n:  

EA(n,r): 
 
end:    

#EAnew(m,n): inputs m>n and outputs the gcd(m,n) using iteration rather than recursion  EAnew:=proc(m,n) local m1,n1,m1new:  

m1:=m:  
n1:=n:    

while n1>0 do    
     m1new:=n1:  
     n1:=m1 mod n1:  
     m1:=m1new:  
od:    
end:  





#MyIFactor(n): inputs integers n and returns the prime factorization in ascending order 
#For example, MyIFactor(100) should output [2,2,5,5] 
MyIFactor:=proc(n) local L,p:  
 
if n=1 then return []: end if:

p := 2;
while n mod p != 0 do
    p := NextPrime1(p + 1):
end do:

L := MyIFactor(n/p):
L := [p, op(L)]:
L:

end:

time(MyIFactor(21*10^2000));
time(ifactor(21*10^20000));
                             0.049

                               0.


#EAnew is better because recursion is generally slower due to all of the repeated function calls    a:=rand(10^50000..10^50001)():   
b:=rand(10^50000..10^50001)():   
time(gcd(a,b));  
time(EAnew(a,b));
                             0.007

                             1.358





#EstimateProbGCD(N,K): estimates the chance that two random natural numbers less than or equal to N  #are prime running K times and tallying the results  
#EstimateProbGCD:=proc(N,K) local hits, total, i:    

hits:=0:  
total:=0:  

for i from 1 to K do       
     if EAnew(rand(1..N)(),rand(1..N)())=1 then ++hits: end if:       
     ++total;  
end do:   

(hits)/(total):    
end:



 
#I started trying to find the limit by thinking of ways to describe the probability exactly for some #given n. My first thought was that saying two numbers are coprime is the same as saying that neither #is divisible by  p for all primes p <= n. The probability of a number being divisible by p is just #1/(p), so the probability  that neither number is divisible by p is just (1- (1)/(p^(2))), so we can #take the product of this for all primes p <= n.  I honestly have no clue how to characterize the #product of such an irregular index set though so I thought further.     

#My second thought was that probability in discrete cases is just the ratio of two counts, so I could #just count the   number of pairs that are coprime and the total number of pairs. The latter is #clearly n^(2), but the former is a little more   tricky. Counting the complement with Inclusion-#Exclusion was the only thing that came to mind. There are   floor(n/(p))integers that are divisible #by p, so we can take the square of this so we get all the pairs and deduct it from n^(2), the total #number of pairs, for all   prime numbers p <= n, but this of course deducts numbers that are #divisible by more than one prime number more  than once, so we would have to add back in the numbers #that are divisible by two primes p1p2, but then re-remove  all the primes divisble by three primes #p1p2p3 for the same reason as before. Compared to my last approach this  makes it so that I don't #have to deal with an index set of primes, but it's just shifted the issue of somehow  manipulating #only the primes, which I honestly have no clue how to do, to the terms rather than the indices. At this   point I didn't think I could make much more progress at my current level of math knowledge as #I'm currently  a 1 st year undergraduate.     

#In the spirit of this class however I thought I'd type up a function for both methods and see what #the probability is  at very high n and compare the two. The summation method is faster than brute-#forcing all pairs but will   undoubtedly be terribly slow relative to the product method. The   nice thing about both methods listed above is that they can be generalized to any m-tuple of random integers by   just raising to the power of m rather than two, which I will also include in my functions since it's a trivial addition.  Implementing the first method is dead simple, but the second is a little trickier. I know how to select combinations  of n size using for loops but I don't know how to do it in a manner that's scalable for arbitrary n. I considered using   some kind of array of indices and using a while loop to replicate the nested for loops I'd normally use for selecting  combinations of a small known size n but I realized it would just be easier to just iterate through all combinations   of any size and then add or remove those pairs depending on whether the number of prime factors is even or odd.  Iterating through these combinations is really simply due to a trick I learned several eons ago when I was a   highschool freshman programming competitively. There's a bijection between all of the non-empty subsets of a set    of size n and the integers from 0 to 2^(n)-1. You simply convert the integer to binary and whether the ith digit   is 1 or 0 to decide whether the ith element is in this particular subset. A lot of the variables are going to be overly long for the sake of readability.        

#I start by just extending our previous functions to check my work for n-tuples.    #checkCoprime(nums,N): This function checks if a set of integers are all relatively prime to each #other by just   checking if any of the primes <= N divide all of the integers.  checkCoprime:=proc(nums,N)local i,p,isFactor:    
p:=2:    
while p <=N/(2) do       
     isFactor:=true:     
  
     for i from 1 to nops(nums) do            
         if nums[i] mod p <>0 then isFactor:=false: end if:       
     end do:              

     if isFactor = true then return false: end if:                   
     p:= NextPrime1(p+1):  
end do:    

true:    

end:    

#EstimateProbGCDGeneral(N,K,M): estimates the chance that M random natural numbers less than or equal #to N  are prime running K times and tallying the results  
EstimateProbGCDGeneral:=proc(N,K,M) local hits, total,i,j:    
hits:=0:  
total:=0:  
for i from 1 to K do       
     if checkCoprime([seq(rand(1..N)(),i=1.. M )],N) then ++hits: end if:       
     ++total;  
end do:   

evalf((hits)/(total)):    

end:  

#Finds the probability that a set of M random integers between 1 and N are coprime using the product #idea  laid out previously. I use NextPrime to iterate through all primes though just iterating #through 1..N and then  checking primality for each would be similar.  
findProbProd:=proc(N,M)local p, runningProb, i:    
p:=2:  
runningProb :=1:    

while p <=N/(2) do       
     runningProb:=runningProb*(1-1/(p^(M))):       
     p:= NextPrime1(p+1):  
end do:    

evalf(runningProb):    

end:
 
#Finds the probability that a set of M random integers between 1 and N are coprime using the #Inclusion-Exclusion  idea laid out previously. This is a ton of local variable but I wanted to #follow the convention in class of declaring  these at the beginning of a function.  findProbSum:=proc(N,M)local p, primeList,currSubset, runningTotal,numOfFactors,digitList,i,currFactor,length:    

primeList:=[]:  
p:=2:  
runningTotal:=N^(M);    

while p<=N/(2) do       
     primeList:=[op(primeList),p]:       
     p:=NextPrime1(p+1):  
end do: 
#The variable length shouldn't be necessary but for some reason it wouldn't let me raise 2 to the #power of nops(primeList)  

length:=nops(primeList):  
for currSubset from 1 to (2^(length)-2) do        
#This looks odd but it's really not, the integer is being converted to binary and then its        #digits are turned into an list.       Extra zeroes are then added to the beginning of this list to #even out the size.       
     digitList :=       
     map(parse,StringTools:-Explode(convert(       
     convert(currSubset,binary),       
     string))):         
     digitList:=[seq(0,i=1..(nops(primeList)-nops(digitList))),op(digitList)]: 
             
     #Now we're turning digitList into the factor that we'll be checking for       
     currFactor :=1:       
     numOfFactors:=0:       
     for i from 1 to nops(digitList)do            
          if digitList[i]=1 then                  
               currFactor:=currFactor*primeList[i]:                 
               ++numOfFactors:             
          end if:        
     end do:         

     runningTotal := runningTotal +((-1)^(numOfFactors)* floor(N/(currFactor))^(M)):    
end do:         

evalf(runningTotal/(N^(M))):    

end:

#I compared both methods against EstimateProbGCDGeneral to test for correctness.
#Below is a time comparison between the two methods, as expect sum is much slower.

time(findProbProd(100, 3));
time(findProbSum(100, 3));
                               0.

                             2.620

#At this point I just googled the answer to see how far off I was. Here's the percent error between #the real value, (6)/(pi^(2)),  and the highest value of N I could feasibly run on my Macbook Air 
#M1 run for the product metho for M=2 (pairs of integers). It took maybe ten-ish seconds. Not too #shabby.    

expand(evalf((findProbProd(1000000,2)-( 6)/(Pi^(2)))/((6/(Pi^(2))))));
                                      -7
                        1.426157836 10  





#MakeRSAkeyG: The key [n,e,d]: a->a^e mod n  n must be a product of G primes  
#inputs D1 and outputs an RSA key  [n,e,d], where n is a product of G primes  
#with D1 digits.   
MakeRSAkeyG:=proc(D1,G) local n,d,S,m,p,q,e,primeList,i,j:    
primeList:=[]:  

for i from 1 to G do       
     primeList:=[op(primeList),NextPrime1(rand(10^(D1-1)..10^D1-1)())]:  
end do:    

for i from 1 to G-1 do       
     for j from i+1 to G do            
          if primeList[i]=primeList[j] then RETURN(FAIL): end if:       
     end do:  
end do:    

n:=1:  
m:=1:  

for i from 1 to G do      
     n:=n*primeList[i];      
     m:=m*(primeList[i]-1);  
end do:  

S:=rand(m/2..m-1)():  
for e from S to m-1 while gcd(e,m)>1 do od:  
d:=e&^(-1) mod m:    
[n,e,d]:  
end:  

#Keys with a higher G are less secure as it means that for any given key length, someone needs to #search through less   values to factor it since each factor is at most key/(G)long.





#I'm unsure of what it means to implement the signatures, here's a decrypt function and an encrypt #function  that includes the signature system. The public and private keys have been named with #either an A or B   at the end for Alice or Bob to help demonstrate what's occuring. Both functions #are practically the same but   this is a bit more illustrative.  

#Encrypt(eB, dA, M): Encodes a message M signed by Alice to Bob using Alice's private key and Bob's public key  
Encrypt:=proc(eB,dA,M,n):
    
#Signing  
M:=M^(dA)mod n:  

#Encrypting  
M:=M^(eB)mod n:  

M:    
end:      

#Decrypt(eA, dB, C,n): Encodes a ciphertext C signed by Bob to Alice using Alice's public key and #Bob's private key `  
Decrypt:=proc(eA,dB,C,n):
    
#Decrypting  
C:=C^(dB)mod n:  

#`Checking Signature`  
C:=C^(eA)mod n:  

C:    
end: