import numpy as np

def prime_list(bound):
  '''Make a list of prime numbers
  up to bound. Requires numpy package.'''
  primes=[]
  for i in range(2,bound+1):
    proceed=True
    for p in primes:
      if np.sqrt(i)<p:
        break
      if i%p==0:
        proceed=False
        break
      else:
        continue
    if proceed==True:
      primes.append(i)
  return primes

def sieve(bound):
  '''Sieves out primes
  up to bound.'''
  bound +=1
  primes=range(2,bound) #this will only contain primes at the end
  for i in range(2,bound):
    if i**2>bound: #more optimized halting condition
      break
    if i in primes:
      mult_i=[i*j for j in range(i,(bound//i)+1)] #list of prime mults
      primes=[n for n in primes if n not in mult_i] #sieve out mult_i
  return primes

import numpy as np

def fermat(n):
  '''Tests a positive integer
  n>3 for probable primality.'''
  if n <= 3 or n%2==0:
    raise Exception('test only valid for odd n>3')
  a=0
  while not gcd(a,n)==1: #finds random a coprime to n
    a=np.random.randint(1,n)
  test=pow(a,n-1,n)
  if test==1:
    print('{} is a probable prime with witness {}.'.format(n,a))
  else:
    print('{} is composite.'.format(n))

from numpy import gcd

def selfridge(n):
  '''Selfridge primality heuristic.'''
  if n < 3 or n%2==0:
    raise Exception('test only valid for odd n>1')
  if not (n%5==2 or n%5==3):
    raise Exception('not cong 2,3 mod 5')
  if not pow(2,n-1,n)==1:
    return False
  if fd_fib(n+1)%n==0: 
    print("It's possible...")
    return True
  else:
    return False

import numpy as np

def mr_test(n):
  '''Miller-Rabin test for
  probable primality.'''
  if n<=2 or n%2==0:
    raise Exception('test only valid for odd n>2')
  # write n=(2**s)*d+1
  b=n-1; s=0
  while b%2==0:
    b=b//2
    s+=1
  d=b
  a=np.random.randint(1,n)
  test=pow(a,d,n)
  if test==1 or test==(n-1):
    print('{} is a probable prime with witness {}.'.format(n,a))
    return
  for r in range(1,s): # r=0 case above
    if pow(a,(2**r)*d,n)==(n-1):
      print('{} is a probable prime with witness {}.'.format(n,a))
      return
  print('{} is composite.'.format(n))

import numpy as np

def det_mr_test(n):
  '''Deterministic Miller-Rabin
  primality test.'''
  if n<=2 or n%2==0:
    raise Exception('test only valid for odd n>2')
  b=n-1; s=0
  while b%2==0:
    b=b//2
    s+=1
  d=b
  for a in range(2,min(n-2,int(2*(np.log(n)**2)))+1):
    x=pow(a,d,n)
    if x==1 or x==(n-1):
      continue
    flag=False
    for r in range(1,s):
      if pow(a,(2**r)*d,n)==(n-1):
        flag=True
        break
    if flag==False:
      print('{} is composite.'.format(n))
      return
  print('{} is prime.'.format(n))

def collatz(n,max_step):
  '''
  Runs the Collatz function
  on n; returns number of steps to 1.
  Will time out after max_step.
  '''
  #Hardcode n=0 case so it doesn't
  #have to loop max_step times to get step=0
  if n==0:
    return 0
  #general logic
  step=0
  while n!=1:
    #timeout condition
    if step>max_step:
      return 0
    #collatz function
    if n%2==0:
      n=n//2
    else:
      n=3*n+1
    step+=1
  return step

def egcd_w_print(a,b):
  '''Computes GCD(a,b) using the Euclidean algorithm.
  Prints formatted division steps.'''
  r=max(a,b)%min(a,b)
  c=min(a,b)
  if r==0:
    print('GCD({},{})={}'.format(a,b,c))
    return
  print('{}={}({})+{}'.format(max(a,b),max(a,b)//c,c,r))
  while c%r != 0:
    d=c
    s=r
    r=c%r
    c=s
    print('{}={}({})+{}'.format(d,d//s,s,r))
  print('GCD({},{})={}'.format(a,b,r))

def gcd(a,b):
  '''Computes GCD(a,b) using the Euclidean algorithm.'''
  if a==0:
    return b
  if b==0:
    return a
  c=min(a,b)
  r=max(a,b)%c
  if r==0:
    return c
  while c%r != 0: #while remainder is nonzero
    s=r
    r=c%r
    c=s
  return r

# adapted from https://stackoverflow.com/a/9758173

def modinv(a, m):
  def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)
  ###
  g, x, y = egcd(a, m)
  if g != 1:
      raise Exception('modular inverse does not exist')
  else:
      return x % m

# encoding dictionary
rsa_alphabet={' ':'00','a':'01','b':'02','c':'03','d':'04','e':'05',
              'f':'06','g':'07','h':'08','i':'09','j':'10','k':'11',
              'l':'12','m':'13','n':'14','o':'15','p':'16','q':'17',
              'r':'18','s':'19','t':'20','u':'21','v':'22','w':'23',
              'x':'24','y':'25','z':'26'}
# decoding dictionary
inv_dict={num:let for let,num in rsa_alphabet.items()}

def mess_to_int(mess):
  '''Translates string to int via RSA dictionary.
  Must declare rsa_alphabet first.'''
  m=''
  for let in mess:
    if let not in rsa_alphabet:
      raise Exception('invalid character; lower case english and space allowed')
    m+=rsa_alphabet[let]
  return int(m)
  
def rsa_encrypt(message,modulus,exponent):
  '''Encrypts a string with RSA encryption.
  Only lower case characters and spaces allowed.
  Limited by size of modulus. 
  Note: public key=(modulus,exponent).'''
  c=pow(mess_to_int(message),exponent,modulus)
  return c

def rsa_test(p,q,e):
  '''This is just to test'''
  N=p*q
  phi_N=(p-1)*(q-1)
  d=modinv(e,phi_N)
  m=int(input("Message:"))
  print(pow(m,d,N))

def init_rsa(p,q):
  '''Decryption setup. Generates a public key
  given the private keys. Returns (N,e,d),
  where (N,e) is the public key and
  d is the decryption key.'''
  fine=False
  while not fine:
    e=int(input('Suggest exponent:'))
    if gcd(e,(p-1)*(q-1))==1:
      fine=True
  d=modinv(e,(p-1)*(q-1))
  return (p*q,e,d)

def rsa_decrypt(cipher):
  '''Decrypt RSA cypher.'''
  m=str(pow(cipher,d,N))
  if len(m)%2==1:
    m='0'+m
  split_list=[m[i:i+2] for i in range(0,len(m),2)]
  mess=''
  for num in split_list:
    mess+=inv_dict[num]
  return mess

def fac(n):
  '''The factorial function
  for positive integers.'''
  if n==1:
    return 1
  else:
    return n*fac(n-1)

from sympy import gamma

def efac(x):
  '''Extended factorial function.
  Takes positive real-valued inputs.'''
  return gamma(x+1)

import numpy as np

# Terrible warning: For some implementation or resource
# reason, this function stops working after n=92. See below.

def fib_mat_exp(n):
  '''Calculates the nth Fibonacci
  number using matrix eponentiation.'''
  f=np.linalg.matrix_power([[1,1],[1,0]],n-1)
  return f[0][0]

# The trouble it seems, is that 
# the matrix exponentiation method eventually
# stops returning the correct answer.
# I suspect it's an issue with resources or the 
# implementation of numpy.linalg.matrix_power.
def lazy_fib(n):
  '''Classic Fibonacci.'''
  a=0;b=1
  for i in range(n-1):
    c=a+b
    a=b
    b=c
  return c

# This dictionary increases the working range of fd_fib()
fib_dict={0:0,1:1}
plist=prime_list(100000)
for i in range(2,plist[-1]+1):
  fib_dict[i]=fib_dict[i-1]+fib_dict[i-2]
for i in range(2,plist[-1]+1):
  if i not in plist:
    fib_dict.pop(i)

def fd_fib(n):
  '''Another Fibonacci algorithm;
  this one uses fast doubling.'''
  if n in fib_dict:
    return fib_dict[n]
  if n==1:
    return 1
  if n%2==0:
    f=fd_fib(n//2)
    return f*(2*fd_fib((n//2)+1)-f)
  else:
    k=(n-1)//2
    return fd_fib(k+1)**2+fd_fib(k)**2