-
Notifications
You must be signed in to change notification settings - Fork 27
/
main.go
116 lines (92 loc) · 2.56 KB
/
main.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
package main
import (
"flag"
"fmt"
"slices"
"strings"
"github.com/itsubaki/q"
"github.com/itsubaki/q/math/number"
"github.com/itsubaki/q/math/rand"
)
// go run main.go --N 15
func main() {
var N, t, a int
var seed uint64
flag.IntVar(&N, "N", 15, "positive integer")
flag.IntVar(&t, "t", 3, "precision bits")
flag.IntVar(&a, "a", -1, "coprime number of N")
flag.Uint64Var(&seed, "seed", 0, "PRNG seed for measurements")
flag.Parse()
if N < 2 {
fmt.Printf("N=%d. N must be greater than 1.\n", N)
return
}
if number.IsEven(N) {
fmt.Printf("N=%d is even. p=%d, q=%d.\n", N, 2, N/2)
return
}
if a, b, ok := number.BaseExp(N); ok {
fmt.Printf("N=%d. N is exponentiation. %d^%d.\n", N, a, b)
return
}
if number.IsPrime(N) {
fmt.Printf("N=%d is prime.\n", N)
return
}
if a < 0 {
a = rand.Coprime(N)
}
if a < 2 || a > N-1 {
fmt.Printf("N=%d, a=%d. a must be 1 < a < N.\n", N, a)
return
}
if number.GCD(N, a) != 1 {
fmt.Printf("N=%d, a=%d. a is not coprime. a is non-trivial factor.\n", N, a)
return
}
fmt.Printf("N=%d, a=%d, t=%d, seed=%d.\n\n", N, a, t, seed)
qsim := q.New()
if seed > 0 {
qsim.Rand = rand.Const(seed)
}
r0 := qsim.ZeroWith(t)
r1 := qsim.ZeroLog2(N)
qsim.X(r1[len(r1)-1])
print("initial state", qsim, r0, r1)
qsim.H(r0...)
print("create superposition", qsim, r0, r1)
// qsim.CModExp2(a, N, r0, r1)
// print("apply controlled-U", qsim, r0, r1)
for i := 0; i < len(r0); i++ {
qsim.ControlledModExp2(a, i, N, r0[i], r1)
print(fmt.Sprintf("apply controlled-U[%d]", i), qsim, r0, r1)
}
qsim.InvQFT(r0...)
print("apply inverse QFT", qsim, r0, r1)
qsim.Measure(r1...)
print("measure reg1", qsim, r0, r1)
for _, st := range qsim.State(r0) {
i, m := st.Int(), st.BinaryString()
s, r, d, ok := number.FindOrder(a, N, fmt.Sprintf("0.%s", m))
if !ok || number.IsOdd(r) {
fmt.Printf(" i=%3d: N=%d, a=%d, t=%d; s/r=%2d/%2d ([0.%v]~%.4f);\n", i, N, a, t, s, r, m, d)
continue
}
p0 := number.GCD(number.Pow(a, r/2)-1, N)
p1 := number.GCD(number.Pow(a, r/2)+1, N)
if number.IsTrivial(N, p0, p1) {
fmt.Printf(" i=%3d: N=%d, a=%d, t=%d; s/r=%2d/%2d ([0.%v]~%.4f); p=%v, q=%v.\n", i, N, a, t, s, r, m, d, p0, p1)
continue
}
fmt.Printf("* i=%3d: N=%d, a=%d, t=%d; s/r=%2d/%2d ([0.%v]~%.4f); p=%v, q=%v.\n", i, N, a, t, s, r, m, d, p0, p1)
}
}
func print(desc string, qsim *q.Q, reg ...any) {
fmt.Println(desc)
max := slices.Max(qsim.Probability())
for _, st := range qsim.State(reg...) {
p := strings.Repeat("*", int(st.Probability()/max*32))
fmt.Printf("%s: %s\n", st, p)
}
fmt.Println()
}