Classical.cs

using System;
using Microsoft.Quantum.Simulation.Core;
using Microsoft.Quantum.Simulation.Simulators;

using Integer = System.Numerics.BigInteger;

namespace ShorAlgorithm {
    enum FactoringType {
        Classical,
        Quantum,
    }

    class ClassicalPart {
        private static Integer SelectRandomInt(Integer n) {
            var rng = new Random();
            var nBytes = n.ToByteArray();
            Integer result;
            do {
                rng.NextBytes(nBytes);
                result = new Integer(nBytes);
            } while (result <= 1 || result >= n);
            return result;
        }

        // Uses Miller-Rabin to check if the number is prime. If so, no point in trying to factor it.
        private static bool IsPrime(Integer n) {
            if (n == 2) {
                return true;
            }

            Integer powerOf2 = 0;
            Integer rest = n - 1;
            while (rest % 2 == 0) {
                ++powerOf2;
                rest /= 2;
            }

            for (var witnessIter = 0; witnessIter < 20; ++witnessIter) {
                var a = SelectRandomInt(n);
                var x = Integer.ModPow(a, rest, n);
                var foundWitness = x != 1 && x != n - 1;
                for (var power = 0; foundWitness && power < powerOf2; ++power) {
                    x = (x * x) % n;
                    if (x == n - 1) {
                        foundWitness = false;
                    }
                }
                if (foundWitness) {
                    return false;
                }
            }

            return true;
        }

        // Checks if the number is a perfect power of a prime. If so, we immediately know how to factor it.
        private static Integer? FindPowerBase(Integer n) {
            var exponent = 2;
            while (Integer.Pow(2, exponent) <= n) {
                Integer min = 1;
                Integer max = n;
                while (max - min >= 2) {
                    var candidateBase = (min + max) / 2;
                    var power = Integer.Min(n + 1, Integer.Pow(candidateBase, exponent));
                    if (power == n) {
                        return candidateBase;
                    } else if (power < n) {
                        min = candidateBase;
                    } else {
                        max = candidateBase;
                    }
                }
                ++exponent;
            }
            return null;
        }

        private static Integer? ExtractPeriodFromBits(IQArray<Result> bits, Integer a, Integer n) {
            var len = bits.Count;
            var numer = Integer.Zero;
            for (var idx = 0; idx < len; ++idx) {
                if (bits[idx] == Result.One) {
                    numer += Integer.Pow(2, len - 1 - idx);
                }
            }
            var denom = Integer.Pow(2, len);

            Integer? prevCand = null;
            Integer? prevPrevCand = null;
            Integer? result = null;
            while (denom != 0) {
                var nextValue = numer / denom;
                var tmp = denom;
                denom = numer - nextValue * denom;
                numer = tmp;

                Integer cand = prevCand.HasValue
                    ? (prevPrevCand.HasValue ? (nextValue * prevCand.Value + prevPrevCand.Value) : nextValue)
                    : 1;

                if (cand < n && Integer.ModPow(a, cand, n) == 1) {
                    result = cand;
                }

                prevPrevCand = prevCand;
                prevCand = cand;
            }
            return result;
        }

        public static (Integer, Integer, FactoringType)? Factor(Integer n, QuantumSimulator sim) {
            if (n <= 1) {
                throw new Exception("can't factor an integer that is less than 2");
            }

            if (IsPrime(n)) {
                return null;
            }

            if (n % 2 == 0) {
                return (2, n / 2, FactoringType.Classical);
            }

            var powerBase = FindPowerBase(n);
            if (powerBase.HasValue) {
                return (powerBase.Value, n / powerBase.Value, FactoringType.Classical);
            }

            for (var iter = 0; iter < 100; ++iter) {
                var candidate = SelectRandomInt(n);

                var g1 = Integer.GreatestCommonDivisor(candidate, n);
                if (g1 != 1) {
                    // If gcd is not 1, quantum factoring won't work. In fact, if gcd is not 1,
                    // we have already factored the number, but let's just pretend this candidate
                    // doesn't exist to make things more interesting.
                    continue;
                }

                var periodBits = GetPeriodBits.Run(sim, candidate, n).Result;
                var period = ExtractPeriodFromBits(periodBits, candidate, n);
                
                var periodRatio = 0.0;
                for (var idx = 0; idx < periodBits.Count; ++idx) {
                    if (periodBits[idx] == Result.One) {
                        periodRatio += 1.0 / Math.Pow(2, idx + 1);
                    }
                }
                Console.WriteLine($"a = {candidate}, n = {n}, periodBits = {periodRatio}, period = {period}");

                if (period.HasValue && period.Value % 2 == 0) {
                    var root = Integer.ModPow(candidate, period.Value / 2, n);
                    var mul1 = root - 1;
                    var mul2 = root + 1;
                    if (mul1 % n != 0 && mul2 % n != 0) {
                        var g2 = Integer.GreatestCommonDivisor(mul1, n);
                        if (n % g2 != 0) {
                            throw new Exception($"{g2} doesn't divide {n}");
                        }

                        return (g2, n / g2, FactoringType.Quantum);
                    }
                }
            }

            return null;
        }
    }
}

Main.cs

using System;
using System.Linq;
using Microsoft.Quantum.Simulation.Simulators;

using Integer = System.Numerics.BigInteger;

namespace ShorAlgorithm {
    class Program {
        static (Integer, Integer, FactoringType)? Factor(Integer n, QuantumSimulator sim, bool shouldBePrime) {
            var factors = ClassicalPart.Factor(n, sim);
            if (shouldBePrime && factors.HasValue) {
                throw new Exception($"{n} should be prime, got factors {factors}");
            }
            if (!shouldBePrime) {
                if (!factors.HasValue) {
                    throw new Exception($"failed to find factors for non-prime {n}");
                }

                var (a, b, factoredQuantumly) = factors.Value;
                var good = a > 1 && a < n && b > 1 && b < n && a * b == n;
                if (!good) {
                    throw new Exception($"{factors} are not valid factors for {n}");
                }
            }
            return factors;
        }

        static string FactorsToString(Integer n, (Integer, Integer, FactoringType)? factors) {
            var factorsStr = factors.HasValue
                ? $"{factors.Value.Item1} * {factors.Value.Item2}"
                : "N/A";
            var result = $"{n} = {factorsStr}";
            if (factors.HasValue && factors.Value.Item3 == FactoringType.Quantum) {
                result = $"𝚿 {result}";
            }
            return result;
        }

        static void Main(string[] args) {
            using (var sim = new QuantumSimulator()) {
                var cmdLineN = args.Length > 0 ? int.Parse(args[0]) : (int?)null;
                var lowerBound = cmdLineN ?? 2;
                var upperBound = cmdLineN ?? 100;
                for (var n = lowerBound; n <= upperBound; ++n) {
                    bool isPrime = Enumerable.Range(2, n - 2).All(x => n % x != 0);
                    Console.WriteLine(FactorsToString(n, Factor(n, sim, isPrime)));
                }
            }
        }
    }
}

Quantum.cs

namespace ShorAlgorithm {
    open Microsoft.Quantum.Arrays;
    open Microsoft.Quantum.Arithmetic;
    open Microsoft.Quantum.Canon;
    open Microsoft.Quantum.Convert;
    open Microsoft.Quantum.Diagnostics;
    open Microsoft.Quantum.Intrinsic;
    open Microsoft.Quantum.Math;

    // The implementation is a slightly modified version of
    // "Circuit for Shor's algorithm using 2n+3 qubits" by Beauregard

    // For extra fun, we roll our own quantum Fourier transform
    // instead of re-using the one from the standard library.
    operation QFT(register : Qubit[]) : Unit is Adj + Ctl {
        let len = Length(register);
        for (idx in 0..len - 1) {
            let q = register[idx];
            H(q);
            for (otherIdx in idx + 1..len - 1) {
                let power = otherIdx - idx + 1;
                Controlled R1Frac([register[otherIdx]], (2, power, q));
            }
        }
        for (idx in 0..len / 2 - 1) {
            SWAP(register[idx], register[len - 1 - idx]);
        }
    }

    // φADD(addend) from the paper. The register should be in the Fourier space.
    operation AddToRegister(register : Qubit[], addend : BigInt) : Unit is Adj + Ctl {
        let len = Length(register);
        for (idx in 0..len - 1) {
            let q = register[len - 1 - idx];
            for (otherIdx in idx..len - 1) {
                if ((addend &&& (1L <<< (len - 1 - otherIdx))) != 0L) {
                    R1Frac(2, otherIdx - idx + 1, q);
                }
            }
        }
    }

    // φADD(addend)MOD(n) from the paper. The register should be in the Fourier space.
    // The ancilla qubit should be set to |0⟩.
    operation AddModNToRegister(register: Qubit[], addend : BigInt, n : BigInt, anc : Qubit) : Unit is Adj + Ctl {
        AddToRegister(register, addend);
        Adjoint AddToRegister(register, n);

        let overflowQubit = Head(register);

        Adjoint QFT(register);
        CX(overflowQubit, anc);
        QFT(register);

        Controlled AddToRegister([anc], (register, n));

        Adjoint AddToRegister(register, addend);

        Adjoint QFT(register);
        X(overflowQubit);
        CX(overflowQubit, anc);
        X(overflowQubit);
        QFT(register);

        AddToRegister(register, addend);
    }

    // CMULT(multiplier)MOD(n) from the paper.
    // The register should NOT be in the Fourier space.
    // The ancilla register should have the same size as the register and be set to |0...0⟩.
    // The ancilla qubit is passed to AddModNToRegister() internally, so it should be set to |0⟩.
    operation AddMulModNToRegister(
        register : Qubit[], multiplier : BigInt, n : BigInt,
        ancRegister : Qubit[], anc : Qubit, readout : Qubit
    ) : Unit is Adj + Ctl {
        let len = Length(register);
        QFT(ancRegister);
        // Ignore the overflow prevention qubit since it is always zero.
        for (idx in 0..len - 2) {
            let mulPower = multiplier * (2L ^ idx) % n;
            let control = [register[len - 1 - idx], readout];
            Controlled AddModNToRegister(control, (ancRegister, mulPower, n, anc));
        }
        Adjoint QFT(ancRegister);
    }

    // The Uₐ gate from the paper, controlled on readout.
    operation ControlledMulRegisterModN(
        register : Qubit[], multiplier : BigInt, n : BigInt,
        ancRegister : Qubit[], anc : Qubit, readout : Qubit
    ) : Unit is Adj+Ctl {
        let len = Length(register);
        Fact(len == Length(ancRegister), "The register and the ancilla are expected to have the same number of qubits");

        AddMulModNToRegister(register, multiplier, n, ancRegister, anc, readout);
        let inverseMul = InverseModL(multiplier, n);

        // Don't swap the overflow prevention qubits since they are always zero.
        for (idx in 0..Length(register) - 2) {
            Controlled SWAP([readout], (register[len - 1 - idx], ancRegister[len - 1 - idx]));
        }

        Adjoint AddMulModNToRegister(register, inverseMul, n, ancRegister, anc, readout);
    }

    operation GetPeriodBits(a : BigInt, n : BigInt) : Result[] {
        let len = BitSizeL(n);
        let iterCount = 2 * len;
        mutable result = new Result[iterCount];

        using ((register, ancRegister, anc1, anc2, ancReadout) = (Qubit[len], Qubit[len], Qubit(), Qubit(), Qubit())) {
            let fullRegister = [anc1] + register;
            let fullAncRegister = [anc1] + ancRegister;

            X(register[len - 1]);

            for (iter in 0..iterCount - 1) {
                let phaseBitIdx = iterCount - 1 - iter;
                let multiplier = ExpModL(a, 2L ^ phaseBitIdx, n);

                H(ancReadout);
                ControlledMulRegisterModN(fullRegister, multiplier, n, fullAncRegister, anc2, ancReadout);
                // Rotate the qubit appropriately before the Hadamard transform (as in https://arxiv.org/pdf/quant-ph/0001066.pdf)
                for (mIdx in phaseBitIdx + 1..iterCount - 1) {
                    if (result[mIdx] == One) {
                        Adjoint R1Frac(2, mIdx - phaseBitIdx + 1, ancReadout);
                    }
                }
                H(ancReadout);

                let readout = M(ancReadout);
                set result = result w/ phaseBitIdx <- readout;

                Reset(ancReadout);
            }

            ResetAll(register);
        }

        return result;
    }
}