Support N=33 in Shor's algorithm

[schrodinteq]

Issue #, if available:

Description of changes:
This PR adds support for modular exponentiation circuits for N = 33 in Shor's algorithm.
Specifically:

• A new subroutine modular_exponentiation_amod33() is added.

• shors_algorithm() now calls modular_exponentiation_amod33() when N == 33.

Although the current implementation does not appear to produce the expected nontrivial factors (i.e., 3 and 11) in practice, these additions allow users to run meaningful Shor factorization examples with N = 33.

Testing done:
Yes

Merge Checklist

Put an x in the boxes that apply. You can also fill these out after creating the PR. If you're unsure about any of them, don't hesitate to ask. We're here to help! This is simply a reminder of what we are going to look for before merging your pull request.

General

• I have read the CONTRIBUTING doc
• I used the commit message format described in CONTRIBUTING
• I have updated any necessary documentation, including READMEs and API docs (if appropriate)

Tests

• I have added tests that prove my fix is effective or that my feature works (if appropriate)
• I have checked that my tests are not configured for a specific region or account (if appropriate)

By submitting this pull request, I confirm that my contribution is made under the terms of the Apache 2.0 license.

[licedric]

Thank you for the pull request! Could you tell us a bit why you think this change is necessary, i.e. why you would like to see the N=33 case?

[schrodinteq]

Thank you for your comment!

Shor’s algorithm is fundamentally designed to work with a variety of composite numbers N.
The notebook at notebooks/textbook/Shors_Algorithm.ipynb currently mentions support for N = 15, 21, and 35.
To improve the generality of the implementation, I thought it would be meaningful to support N = 33 — another composite number between 21 and 35, and a product of two distinct prime numbers.

While the current library technically allows users to set N = 33 and returns a result, it internally calls modular_exponentiation_amod15, which is tailored specifically for N = 15.
As a result, the algorithm runs under incorrect theoretical assumptions for N = 33.

To address this, I added a new modular_exponentiation_amod33() subroutine that correctly handles this case.
With this fix, the algorithm is now able to return the correct nontrivial factors (3 and 11), making the N = 33 case a valid and educationally useful example of Shor’s algorithm in action.

[licedric]

Hi schrodinteq,

The notebooks in the amazon-braket-algorithm-library are not meant to be comprehensive; rather they serve as an example and/or initial building block from which a user can construct their algorithms. As such, I think the N=15 case is sufficient for a simple demonstration.

If you're interested in making this example more general, could I suggest implementing a more general modular exponentiation routine, e.g. along the lines in this paper https://arxiv.org/abs/2311.08555 ? That would indeed be a significant improvement to this notebook's applicability.

[schrodinteq]

Hi licedric,

Understood — I’ll avoid adding a dedicated function for N = 33.
If I manage to implement a more general modular exponentiation approach, as suggested in the paper you shared, I’d be happy to submit a new PR in the future.

That said, I noticed that the notebook (notebooks/textbook/Shors_Algorithm.ipynb) mentions that the library works for N = 21 and N = 35 as well.

However, in the current implementation, when setting N = 21 or 35, the function modular_exponentiation_amod15() is still called internally. This leads to incorrect theoretical behavior. This issue is also mentioned in Issue #156.

For example, when using a = 8 with N = 21, the current implementation adds the following gates:

            if integer_a in [7, 8]:
                mod_exp_amod15.cswap(x, aux_qubits[2], aux_qubits[3])
                mod_exp_amod15.cswap(x, aux_qubits[1], aux_qubits[2])
                mod_exp_amod15.cswap(x, aux_qubits[0], aux_qubits[1])

However, the correct circuit for modular exponentiation in this case should look like:

            if integer_a == 8:
                mod_exp_amod15.cnot(x, aux_qubits[3])

If it would be helpful, I’d be happy to work on adding dedicated subroutines for N = 21 and N = 35 to ensure theoretical correctness.

Please let me know your thoughts.