loom/test.lum at main · raaymax/loom · GitHub

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
// Euclid GCD
fn gcd(a, b) {
  while (b != 0) {
    t = b;
    b = a % b;
    a = t;
  };
  a
};

// fast (a^e) mod m
fn pow_mod(a, e, m) {
  res = 1;
  base = a % m;
  while (e > 0) {
    if (e % 2 == 1) { res = (res * base) % m; };
    base = (base * base) % m;
    e = e / 2;
  };
  res
};

// multiplicative order of a mod n (assumes gcd(a,n)==1)
fn findOrder(a, n) {
  y = 1;
  t = pow_mod(a, y, n);
  while (t != 1) {
    y = y + 1;
    t = pow_mod(a, y, n);
  };
  y
};

// classical Shor post-processing (trial a)
fn shor(x) {
  if (x % 2 == 0) { return "[ " + 2 + ", " + (x/2) + " ]"; };

  a = 2;
  while (a < x) {
    d = gcd(a, x);
    if (d > 1 && d < x) { return "[ " + d + ", " + (x/d) + " ]"; };

    if (d == 1) {
      r = findOrder(a, x);
      if (r % 2 == 0) {
        ar2 = pow_mod(a, r/2, x);         // a^(r/2) mod x
        if (ar2 != x - 1) {               // not -1 mod x
          p = gcd(ar2 + 1, x);
          q = gcd(ar2 - 1, x);
          if (p*q == x && p > 1 && q > 1) { return "[ " + p + ", " + q + " ]"; };
        };
      };
    };
    a = a + 1;
  };
  ""
};

// smoke tests
x = 21; y = 49;
print("gcd("+x+","+y+") = " + gcd(x,y));
print("shor(15) = " + shor(15));
assert(gcd(15, 21) == 3);
assert(shor(15) == "[ 5, 3 ]");
0