Shor’s Algorithm - Period Finding Co-routine

Presentation
Import QC framework
<<Quantum`
In[]:=
Set variables:
  • N1: Number to be factored.
    • ◼
  • a: Initial guess (a<N1).
    • ◼
  • q: Counting qubits
    (
    2
    N1
    ≤
    q
    2
    <2
    2
    N1
    )
    .
    • ◼
  • n: Ancilla qubits (n ≥
    log
    2
    N1
    ).
    • ◼
    N1=6​​a=7;​​N1^2
    In[]:=
    6
    Out[]=
    36
    Out[]=
    q=3;​​2^q
    In[]:=
    8
    Out[]=
    n=2;​​Ceiling[Log[2,N1]]
    In[]:=
    3
    Out[]=
    Prepare initial circuit state (counting qubits already in plus state, ancilla qubits in 0 state)
    initialstates=QuantumTensorProduct[Catenate[{Table[QuantumDiscreteState["Plus"],{i,1,q}],Table[QuantumDiscreteState[{0,1},2],{j,1,n-1}],{QuantumDiscreteState[{1,0},2]}}]]
    In[]:=
    QuantumDiscreteState
    State Type: Pure
    Qudits: 5
    Dimensions: 2
    ​
    
    Out[]=
    Naming basis elements
    map=AssociationThread[Flatten[Table[Ket[i,j],{i,0,2^q-1},{j,0,2^n-1}]]Range[2^(n+q)]]
    In[]:=
    |0,0〉1,|0,1〉2,|0,2〉3,|0,3〉4,|1,0〉5,|1,1〉6,|1,2〉7,|1,3〉8,|2,0〉9,|2,1〉10,|2,2〉11,|2,3〉12,|3,0〉13,|3,1〉14,|3,2〉15,|3,3〉16,|4,0〉17,|4,1〉18,|4,2〉19,|4,3〉20,|5,0〉21,|5,1〉22,|5,2〉23,|5,3〉24,|6,0〉25,|6,1〉26,|6,2〉27,|6,3〉28,|7,0〉29,|7,1〉30,|7,2〉31,|7,3〉32
    Out[]=
    Defining required operator action
    action=Association[Table[Ket[i,0]Ket[i,Mod[a^i,N1]],{i,0,2^q-1}]]
    In[]:=
    |0,0〉|0,1〉,|1,0〉|1,1〉,|2,0〉|2,1〉,|3,0〉|3,1〉,|4,0〉|4,1〉,|5,0〉|5,1〉,|6,0〉|6,1〉,|7,0〉|7,1〉
    Out[]=
    Constructing matrix form
    inputoutput=KeyMap[Replace[map]]@(action/.map);
    In[]:=
    replacements=Flatten[KeyValueMap[{{#1,#2}1,{#2,#1}1,{#1,#1}0,{#2,#2}0}&,inputoutput]];
    In[]:=
    matrixForm=ReplacePart[IdentityMatrix[2^(q+n)],replacements];
    In[]:=
    Constructing and applying CondModExp operator (may take a while)
    condModExp=QuantumDiscreteOperator[N[matrixForm],Range[q+n]]
    In[]:=
    QuantumDiscreteOperator
    Arity: 5
    Dimensions: 2
    
    Out[]=
    moddedstates=condModExp[initialstates]
    In[]:=
    QuantumDiscreteState
    State Type: Pure
    Qudits: 5
    Dimensions: 2
    ​
    
    Out[]=
    Constructing and applying inverse QFT (May take a while)
    QFT=QuantumDiscreteOperator[{"Fourier",1,2^q},Range[q]]
    In[]:=
    QuantumDiscreteOperator
    Arity: 3
    Dimensions: 2
    
    Out[]=
    iQFT=QuantumDiscreteOperator[Inverse[N[QFT["MatrixRepresentation"]]],Range[q]]
    In[]:=
    QuantumDiscreteOperator
    Arity: 3
    Dimensions: 2
    
    Out[]=
    finalstates=iQFT[moddedstates]
    In[]:=
    QuantumDiscreteState
    State Type: Pure
    Qudits: 5
    Dimensions: 2
    ​
    
    Out[]=
    Evaluating final state
    KeyValueMap[If[Abs[#2]>0.001,First[First[Position[#1,1]]]#2,##&[]]&,KeyMap[Replace[finalstates["Basis"]["BasisElementAssociation"]]]@finalstates["Amplitudes"]]/.Association[KeyValueMap[#2#1&,map]]​​
    In[]:=
    {|7,1〉1.-1.2326×
    -32
    10
    }
    Out[]=
    QMS=ResourceFunction["QuantumToMultiwaySystem"]
    In[]:=
    [◼]
    QuantumToMultiwaySystem
    Out[]=
    Converting superpositions to proportional initial states
    toProp[stateV_]:=(probs=Table[If[#[[i]]>0,#[[i]]^2ReplacePart[ConstantArray[0,Length[#]],{i}1],##&[]],{i,1,Length[#]}]&@stateV;​​out=Catenate[Table[Table[#[[i,2]],{j,1,#[[i,1]]/Min[probs[[All,1]]]}],{i,1,Length[#]}]]&@probs;​​Return[out];)
    In[]:=
    QMS[Floor[condModExp["MatrixRepresentation"]],toProp[initialstates["StateVector"]],1,"EvolutionGraph","IncludeStatePathWeights"True,VertexLabels"VertexWeight"]
    In[]:=
    Out[]=
    QMS[Floor[iQFT["MatrixRepresentation"]],toProp[moddedstates["StateVector"]],2]
    In[]:=
    {{{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0},{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0}},{},{}}
    Out[]=