# Setup básico
!pip install qiskit -q
!pip install qiskit[visualization] -q
import qiskit as qk

!pip install qiskit-aer -q
import qiskit_aer as qk_aer

import numpy as np
np.set_printoptions(precision=3, suppress=True)

from matplotlib import pyplot as plt
%matplotlib inline

import pandas as pd
import sklearn as sk

import qiskit as qk

qrx = qk.QuantumRegister(2,'x')
qry = qk.QuantumRegister(2,'y')
qc = qk.QuantumCircuit(qrx,qry)

qc.barrier()
qc.cx(0,2)
qc.cx(0,3)
qc.cx(1,2)
qc.cx(1,3)
qc.barrier()

qc.draw('mpl')

def histograma(qc):
  simulador = qk_aer.Aer.get_backend('aer_simulator')
  resultado = simulator.run(qk.transpile(qc,simulator),shots=10000).result()
  contagens = resultado.get_counts()
  grafico = qk.visualization.plot_histogram(contagens)
  return(grafico)

def simon(initial_y, initial_x):
  import qiskit as qk
  qrx = qk.QuantumRegister(2,'x')
  crx = qk.ClassicalRegister(4,'c')
  qry = qk.QuantumRegister(2,'y')

  qc = qk.QuantumCircuit(qrx,qry,crx)

  qc.initialize(initial_x,qrx)
  qc.initialize(initial_y,qry)

  qc.barrier()
  qc.cx(0,2)
  qc.cx(0,3)
  qc.cx(1,2)
  qc.cx(1,3)
  qc.barrier()

  qc.measure([0,1,2,3],[0,1,2,3])

  display(qc.draw('mpl'))
  grafico = histograma(qc)
  return(grafico)

simon([1,0,0,0],[1,0,0,0])

simon([1,0,0,0],[0,0,0,1])

simon([1,0,0,0],[0,1,0,0])

simon([1,0,0,0],[0,0,1,0])

import qiskit as qk
qrx = qk.QuantumRegister(2,'x')
crx = qk.ClassicalRegister(2,'c')
qry = qk.QuantumRegister(2,'y')

qc = qk.QuantumCircuit(qrx,qry,crx)

qc.initialize([1,0,0,0],qrx)
qc.initialize([1,0,0,0],qry)
qc.h([0,1])

qc.barrier()
qc.cx(0,2)
qc.cx(0,3)
qc.cx(1,2)
qc.cx(1,3)
qc.barrier()

qc.h([0,1])
qc.measure([0,1],[0,1])

display(qc.draw('mpl'))
histograma(qc)

import numpy as np
x0 = np.array([[1,0]]).T
x1 = np.array([[1,0]]).T
x1x0 = np.kron(x1,x0)

y0 = np.array([[1,0]]).T
y1 = np.array([[1,0]]).T
y1y0 = np.kron(y1,y0)

psi0 = np.kron(y1y0,x1x0)
psi0.T

H = np.array([[1,1],[1,-1]],dtype='int')
H2 = np.kron(H,H)
I2 = np.eye(4)
I2H2 = np.kron(I2,H2)
psi1 = I2H2.dot(psi0)/(np.sqrt(2)*np.sqrt(2))
I2H2.dot(psi0).T

Tf1_2 = np.eye(16)
Tf1_2 = Tf1_2[:,[0,5,2,7,4,1,6,3,8,9,10,11,12,13,14,15]]
psi2 = Tf1_2.dot(psi1); psi2.T*2

Tf2_3 = np.eye(16)
Tf2_3 = Tf2_3[:,[0,1,2,3,4,13,6,15,8,9,10,11,12,5,14,7]]
psi3 = Tf2_3.dot(psi2); psi3.T*2

Tf3_4 = np.eye(16)
Tf3_4 = Tf3_4[:,[0,1,6,3,4,5,2,7,8,9,10,15,12,13,14,11]]
psi4 = Tf3_4.dot(psi3); psi4.T*2

Tf4_5 = np.eye(16)
Tf4_5 = Tf4_5[:,[0,1,2,11,4,5,14,7,8,9,10,3,12,13,6,15]]
psi5 = Tf4_5.dot(psi4); psi5.T*2

I2 = np.eye(4)
H = np.array([[1,1],[1,-1]])/np.sqrt(2)
H2 = np.kron(H,H)
I2H2 = np.kron(I2,H2)
psi6 = I2H2.dot(psi5)
print(psi6.T*2)

import numpy as np
I2 = np.eye(4, dtype='int');
H2 = np.array([[1, 1, 1, 1],
               [1,-1, 1,-1],
               [1, 1,-1,-1],
               [1,-1,-1, 1]])/2;
I2H2 = np.kron(I2,H2);
print(I2H2*2)

import numpy as np
x0 = np.array([[1,0]]).T
x1 = np.array([[1,0]]).T
x2 = np.array([[1,0]]).T
x2x1x0 = np.kron(x2,np.kron(x1,x0))
x2x1x0.T

y0 = np.array([[1,0]]).T
y1 = np.array([[1,0]]).T
y2 = np.array([[1,0]]).T
y2y1y0 = np.kron(y2,np.kron(y1,y0))
y2y1y0.T

psi0 = np.kron(y2y1y0,x2x1x0)
psi0.T

I3 = np.eye(2**3)
H = np.array([[1,1],[1,-1]])/np.sqrt(2)
H2 = np.kron(H,H)
H3 = np.kron(H,H2)
I3H3 = np.kron(I3,H3)
psi1 = I3H3.dot(psi0)
psi1.T.round(2)

n = 2*3 # 2 números de 3 digitos cada
Uf = np.eye(2**n) # 2^6 = 64 posições
Uf = Uf[:,[ 0, 9,18,27,12, 5,30,23,
            8, 1,10,11, 4,13,14,15,
           16,17, 2,19,20,21,22, 7,
           24,25,26, 3,28,29, 6,31,
           32,33,34,35,36,37,38,39,
           40,41,42,43,44,45,46,47,
           48,49,50,51,52,53,54,55,
           56,57,58,59,60,61,62,63]]

psi5 = Uf.dot(psi1)
psi5.T.round(2)

psi6 = I3H3.dot(psi5)
psi6.T

np.where(psi6 > 0)[0]

[format(i,'#08b')[2:] for i in np.where(psi6 > 0)[0]]

import numpy as np
CCNOT = np.eye(8)
CCNOT = CCNOT[:,[0,1,2,7,4,5,6,3]]
CCNOT.dot(CCNOT)

import qiskit as qk

def qNAND(y0,x0):
  x = qk.QuantumRegister(1,"x")
  y = qk.QuantumRegister(1,"y")
  z = qk.QuantumRegister(1,"z")
  u = qk.QuantumRegister(1,"u")
  c = qk.ClassicalRegister(3,'c')
  qc = qk.QuantumCircuit(x,y,z,u,c)

  # Put qubits u and z in state |0>
  qc.initialize([1,0],z)
  qc.initialize([1,0],u)
  qc.initialize(y0,y)
  qc.initialize(x0,x)

  # Perform computation
  qc.barrier()
  qc.ccx(x,y,z)
  qc.x(z)
  # Copy CLASSICALY state of z to u
  qc.cx(z,u)
  qc.barrier()

  # Reverse computation
  qc.x(z)
  qc.ccx(x,y,z)
  qc.barrier()

  # Measure results
  qc.measure(x,c[0])
  qc.measure(y,c[1])
  qc.measure(u,c[2])
  display(qc.draw('mpl'))

  simulator = qk_aer.Aer.get_backend("aer_simulator")
  results = simulator.run(qk.transpile(qc,simulator), shots=10000).result().get_counts()
  grafico = qk.visualization.plot_histogram(results)
  return(grafico)

print("0 NAND 0 = 1")
qNAND([1,0],[1,0])

print("0 NAND 1 = 1")
qNAND([1,0],[0,1])

print("1 NAND 0 = 1")
qNAND([0,1],[1,0])

print("1 NAND 1 = 0")
qNAND([0,1],[0,1])

import qiskit as qk

def qOR(y0,x0):
  x = qk.QuantumRegister(1,"x")
  y = qk.QuantumRegister(1,"y")
  z = qk.QuantumRegister(1,"z")
  u = qk.QuantumRegister(1,"u")
  c = qk.ClassicalRegister(3,'c')

  # Initialize qubits
  qc = qk.QuantumCircuit(x,y,z,u,c)
  qc.initialize(x0,x)
  qc.initialize(y0,y)
  qc.initialize([1,0],z)
  qc.initialize([1,0],u)

  # Compute function
  qc.barrier()
  qc.x(x)
  qc.x(y)
  qc.ccx(x,y,z)
  qc.x(z)
  qc.cx(z,u)
  qc.barrier()

  # Reverse computation
  qc.x(z)
  qc.ccx(x,y,z)
  qc.x(y)
  qc.x(x)
  qc.barrier()

  # Measure results
  qc.measure(x,c[0])
  qc.measure(y,c[1])
  qc.measure(u,c[2])
  display(qc.draw('mpl'))

  # Simulate circuit
  simulator = qk_aer.Aer.get_backend('aer_simulator')
  results = simulator.run(qk.transpile(qc,simulator),shots=10000).result().get_counts()
  grafico = qk.visualization.plot_histogram(results)
  return(grafico)

qOR([1,0],[1,0])

qOR([1,0],[0,1])

qOR([0,1],[1,0])

qOR([0,1],[0,1])

import qiskit as qk

def qNXOR(y0,x0,calc=True,show=False):
  ##################################################
  # Create registers                               #
  ##################################################
  # Base registers
  x = qk.QuantumRegister(1,"x")
  y = qk.QuantumRegister(1,"y")

  # Auxiliary register (for x and y)
  z = qk.QuantumRegister(1,"z")

  # Auxiliary register (to store x AND y)
  a1 = qk.QuantumRegister(1,"a1")

  # Auxiliary register (to store NOT(x) AND NOT(y))
  a2 = qk.QuantumRegister(1,"a2")

  # Auxiliary register (for a1 and a2)
  b1 = qk.QuantumRegister(1,"b1")

  # Auxiliary register (to store a1 OR a2)
  b2 = qk.QuantumRegister(1,"b2")

  # Classical Registers to store x,y and final measurement
  c = qk.ClassicalRegister(3,"c")
  ##################################################
  # Create Circuit                                 #
  ##################################################
  qc = qk.QuantumCircuit(x,y,z,a1,a2,b1,b2,c)

  ##################################################
  # Initialize registers                           #
  ##################################################
  qc.initialize(x0,x)
  qc.initialize(y0,y)
  qc.initialize([1,0],z)
  qc.initialize([1,0],a1)
  qc.initialize([1,0],a2)
  qc.initialize([1,0],b1)
  qc.initialize([1,0],b2)

  ###################################################
  # Calculate x AND y. Copy result to a1. Reverse z #
  ###################################################
  qc.barrier()
  qc.ccx(x,y,z)
  qc.cx(z,a1)
  qc.ccx(x,y,z)
  qc.barrier()

  #########################################################
  # Calc. NOT(x) AND NOT(y). Copy result to a2. Reverse z #
  #########################################################
  qc.barrier()
  qc.x(x)
  qc.x(y)
  qc.ccx(x,y,z)
  qc.cx(z,a2)
  qc.ccx(x,y,z)
  qc.x(y)
  qc.x(x)
  qc.barrier()

  #################################################
  # Calc. a1 OR a2. Copy result to b2. Reverse b1 #
  #################################################
  qc.barrier()
  qc.x(a1)
  qc.x(a2)
  qc.ccx(a1,a2,b1)
  qc.x(b1)
  qc.cx(b1,b2)
  qc.x(b1)
  qc.ccx(a1,a2,b1)
  qc.barrier()

  #################################################
  # Measure b2                                    #
  #################################################
  qc.measure(x,c[0])
  qc.measure(y,c[1])
  qc.measure(b2,c[2])

  #################################################
  # Draw circuit                                  #
  #################################################
  if show:
    display(qc.draw("mpl"))

  #################################################
  # Run circuit. Collect results                  #
  #################################################
  if calc:
    simulator = qk_aer.Aer.get_backend("aer_simulator")
    results = simulator.run(qk.transpile(qc,simulator),shots=10000).result().get_counts()
    grafico = qk.visualization.plot_histogram(results)
    return(grafico)
  else:
    return()


qNXOR([1,0],[1,0],False,True)

qNXOR([1,0],[1,0])

qNXOR([1,0],[0,1])

qNXOR([0,1],[1,0])

qNXOR([0,1],[0,1])

x0 = np.array([1,1])/np.sqrt(2)
y0 = np.array([1,1])/np.sqrt(2)
qNXOR(y0,x0)

import qiskit as qk
from qiskit.quantum_info.operators import Operator

qr = qk.QuantumRegister(6,"q")
cr = qk.ClassicalRegister(6,"c")
qc = qk.QuantumCircuit(qr,cr)

for i in range(6):
  qc.initialize([1,0],i)

for i in range(3):
  qc.h(i)

oracle = np.eye(2**6)
oracle[:,[ 0,  0]] = oracle[:,[  0, 0]]
oracle[:,[ 1,  9]] = oracle[:,[  9, 1]]
oracle[:,[ 2, 18]] = oracle[:,[ 18, 2]]
oracle[:,[ 3, 27]] = oracle[:,[ 27, 3]]
oracle[:,[ 4, 12]] = oracle[:,[ 12, 4]]
oracle[:,[ 5,  5]] = oracle[:,[  5, 5]]
oracle[:,[ 6, 30]] = oracle[:,[ 30, 6]]
oracle[:,[ 7, 23]] = oracle[:,[ 23, 7]]
oracle = Operator(oracle)

qc.append(oracle,qr)

for i in range(3):
  qc.h(i)

qc.barrier()

for i in range(3):
  qc.measure(i,i)

display(qc.draw("mpl"))

simulator = qk_aer.Aer.get_backend("aer_simulator")
results = simulator.run(qk.transpile(qc,simulator),shots=10000).result().get_counts()
qk.visualization.plot_histogram(results)

import numpy as np
oracle = np.eye(8)
oracle[:,[2,6]] = oracle[:,[6,2]]
oracle

import numpy as np
x0 = np.array([[1,0]]).T
x1 = np.array([[1,0]]).T
psi0 = np.kron(x1,x0)

H = np.array([[1,1],[1,-1]])/np.sqrt(2)
H2 = np.kron(H,H)
psi1 = H2.dot(psi0)
psi1.T

# Initialize a qubit y in state |1>
y  = np.array([[0,1]]).T

# Pass it through a Hadamard gate and obtain state psi2a
H = np.array([[1,1],[1,-1]])/np.sqrt(2)
psi2a = H.dot(y)
psi2a.T

psi2 = np.kron(psi2a,psi1)
psi2.T

oracle = np.eye(8)
oracle[:,[2,6]] = oracle[:,[6,2]]
oracle

psi3 = oracle.dot(psi2)
psi3.T

v0 = np.array([[10,20,-30,40,50]]).T
mu = np.mean(v0)
v0.T, mu

D = v0 - mu
D.T

v1 = 2*np.mean(v0)-v0
v1.T

import numpy as np
A = np.ones(5**2).reshape(5,5)
A = A/5
I = np.eye(5)
R = 2*A-I
R

R.dot(R)

n = 4
I = np.eye(n)
A = np.ones(n**2).reshape(n,n)/n
R = 2*A-I
R

I2 = np.eye(2)
I2R = np.kron(I2,R)
I2R

psi4 = I2R.dot(psi3)
psi4.T

import qiskit as qk
import numpy as np

from qiskit.quantum_info.operators import Operator

q = qk.QuantumRegister(3,'q')
c = qk.ClassicalRegister(3,'c')
qc = qk.QuantumCircuit(q,c)

qc.initialize([1,0],0)
qc.initialize([1,0],1)
qc.h(0)
qc.h(1)

qc.initialize([0,1],2)
qc.h(2)

qc.barrier()
oracle = np.eye(8)
oracle[:,[2,6]] = oracle[:,[6,2]]
oracle = Operator(oracle)
qc.append(oracle,q)

# n = 2 (two digits number to look for)
A = np.ones(4*4).reshape(4,4)/4
I4 = np.eye(4)
R = 2*A - I4
I2 = np.eye(2)
boost = np.kron(I2,R)
boost = Operator(boost)
qc.append(boost,q)

qc.barrier()
qc.measure(q[:2],c[:2])

display(qc.draw('mpl'))

simulator = qk_aer.Aer.get_backend('aer_simulator')
results = simulator.run(qk.transpile(qc,simulator),shots=1000).result().get_counts()
qk.visualization.plot_histogram(results)

import numpy as np
x0 = np.array([[1,0]]).T
x1 = np.array([[1,0]]).T
x2 = np.array([[1,0]]).T
x2x1x0 = np.kron(x2,np.kron(x1,x0))
x2x1x0.T

H = np.array([[1,1],[1,-1]])/np.sqrt(2)
H2 = np.kron(H,H)
H3 = np.kron(H,H2)
x = H3.dot(x2x1x0)
x.T

y0 = np.array([[0,1]]).T
y0.T

H = np.array([[1,1],[1,-1]])/np.sqrt(2)
y = H.dot(y0)
y.T

psi0 = np.kron(y,x)
psi0.T

oracle = np.eye(16)
oracle[:,[5,13]] = oracle[:,[13,5]]
psi1 = oracle.dot(psi0)
psi1.T

# n = 8
A = np.ones(8*8).reshape(8,8)/8
I8 = np.eye(8)
R = 2*A-I8
R

I2 = np.eye(2)
I2R = np.kron(I2,R)
psi2 = I2R.dot(psi1)
psi2.T

0.625**2 + 0.625**2

psi3 = I2R.dot(oracle).dot(psi2)
psi3.T

0.687**2 + 0.687**2

I2R.dot(oracle).dot(psi3).T

T = I2R.dot(oracle)

psi = psi1
prob = []
for i in range(25):
  prob.append(psi[5,0]**2 + psi[13,0]**2)
  psi = T.dot(psi)

from matplotlib import pyplot as plt
plt.plot(list(range(25)), prob);

import qiskit as qk
import numpy as np

from qiskit.quantum_info.operators import Operator

q = qk.QuantumRegister(4,'q')
c = qk.ClassicalRegister(4,'c')
qc = qk.QuantumCircuit(q,c)

qc.initialize([1,0],0)
qc.initialize([1,0],1)
qc.initialize([1,0],2)
qc.h(0)
qc.h(1)
qc.h(2)

qc.initialize([0,1],3)
qc.h(3)

qc.barrier()

# Create oracle matrix
oracle = np.eye(16)
oracle[:,[5,13]] = oracle[:,[13,5]]

# Create rotation about the mean matrix
# n = 3 (three digits number to look for)
A = np.ones(8*8).reshape(8,8)/8
I8 = np.eye(8)
R = 2*A - I8

# Identity matrix to leave fourth qubit undisturbed
I2 = np.eye(2)

# Create second transformation matrix
boost = np.kron(I2,R)

# Combine oracle with second transformation matrix in
#  an unique operator
T = boost.dot(oracle)
T = Operator(T)

# Apply operator twice (2 phase inversions and
#  rotations about the mean)
qc.append(T,q)
qc.append(T,q)

qc.barrier()
qc.measure(q[:3],c[:3])

display(qc.draw('mpl'))

simulator = qk_aer.Aer.get_backend('aer_simulator')
results = simulator.run(qk.transpile(qc,simulator),shots=1000).result().get_counts()
qk.visualization.plot_histogram(results)

import qiskit as qk
import numpy as np

from qiskit.quantum_info.operators import Operator

q = qk.QuantumRegister(5,'q')
c = qk.ClassicalRegister(5,'c')
qc = qk.QuantumCircuit(q,c)

qc.initialize([1,0],0)
qc.initialize([1,0],1)
qc.initialize([1,0],2)
qc.initialize([1,0],3)
qc.h(0)
qc.h(1)
qc.h(2)
qc.h(3)

qc.initialize([0,1],4)
qc.h(4)

qc.barrier()

# Create oracle matrix
oracle = np.eye(32)
oracle[:,[13,29]] = oracle[:,[29,13]]

# Create rotation about the mean matrix
# n = 4 (four digits number to look for)
A = np.ones(16*16).reshape(16,16)/16
I16 = np.eye(16)
R = 2*A - I16

# Identity matrix to leave fifth qubit undisturbed
I2 = np.eye(2)

# Create second transformation matrix
boost = np.kron(I2,R)

# Combine oracle with second transformation matrix in
#  an unique operator
T = boost.dot(oracle)

# n = 4. Therefore it will be necessary 4 T operations
#  to get maximum probability.
Tn = T.dot(T).dot(T)
Tn = Operator(Tn)

# Apply operator Tn once (n phase inversions and
#  rotations about the mean)
qc.append(Tn,q)

qc.barrier()
qc.measure(q[:4],c[:4])

display(qc.draw('mpl'))

simulator = qk_aer.Aer.get_backend('aer_simulator')
results = simulator.run(qk.transpile(qc,simulator),shots=1000).result().get_counts()
qk.visualization.plot_histogram(results)



import numpy as np
import pandas as pd
N = 35; a = 6

GCD = pd.DataFrame({'N':[N], 'a':[a], 'N-a':[N-a]})

while GCD.iloc[-1,2] != 0:
  temp = GCD.iloc[-1].values
  temp = np.sort(temp); temp[-1] = temp[-2] - temp[-3]
  temp = pd.DataFrame(data = temp.reshape(-1,3), columns = ['N','a','N-a'])
  GCD = pd.concat([GCD,temp], ignore_index=True)

GCD

import numpy as np
import pandas as pd
N = 15; a = 9

GCD = pd.DataFrame({'N':[N], 'a':[a], 'N-a':[N-a]})

while GCD.iloc[-1,2] != 0:
  temp = GCD.iloc[-1].values
  temp = np.sort(temp); temp[-1] = temp[-2] - temp[-3]
  temp = pd.DataFrame(data = temp.reshape(-1,3), columns = ['N','a','N-a'])
  GCD = pd.concat([GCD,temp], ignore_index=True)

GCD

import numpy as np
np.gcd(6,35), np.gcd(9,15)

import numpy as np
np.gcd(215,35), np.gcd(217,35)

20 % 7

import numpy as np
f = lambda x,a,N: a**x % N

[f(x,2,15) for x in range(10)]

import numpy as np
(13*35)%11, ((2%11)*(24%11))%11

import numpy as np

def f(x,a,N):
  a0 = 1
  for i in range(x):
    a0 = ((a%N)*(a0))%N
  return(a0)

[f(x,2,371) for x in range(8)]

[f(x,2,371) for x in range(154,159)]

[f(x,24,371) for x in range(77,81)]

import numpy as np

def fr(a,N):
  ax = ((a%N)*(1))%N
  x = 1
  while ax != 1:
    ax = ((a%N)*(ax))%N
    x = x+1
  return(x)

fr(2,371), fr(6,371), fr(24,371)

import numpy as np

def fr(a,N):
  a0 = 1
  a0 = ((a%N)*a0)%N
  x = 2
  find = False
  while find==False:
    a0 = ((a%N)*(a0))%N
    if a0 == 1:
      find = True
    x = x+1
  return(x-1)

N = 35
a = 2
fr(a,N)

2**(12/2)-1, 2**(12/2)+1

np.gcd(63,35), np.gcd(65,35)

import numpy as np
# Function to calculate the period of f(x) = a^x MOD N
def fr(a,N):
  a1 = ((a%N)*(1))%N
  x = 1
  while a1 != 1:
    a1 = ((a%N)*a1)%N
    x = x+1
  return(x)

# Function to factor N based on a
def factor(N, a=2):
  r = fr(a,N)
  g1 = a**(int(r/2)) - 1
  g2 = a**(int(r/2)) + 1
  f1 = np.gcd(g1,N)
  f2 = np.gcd(g2,N)
  return(f1,f2, f1*f2)

factor(247)

factor(1045)