T PExp = (BigInt prime, Int exp)

F isqrt(self)
   V b = self
   L
      V a = b
      b = (self I/ a + a) I/ 2
      I b >= a
         R a

F factor(BigInt n)
   [PExp] pf
   V nn = n
   V b = 0
   L ((nn % 2) == 0)
      nn I/= 2
      b++

   I b > 0
      pf [+]= PExp(BigInt(2), b)

   V s = isqrt(nn)
   V d = BigInt(3)
   L nn > 1
      I d > s
         d = nn
      V e = 0
      L
         V (div, rem) = divmod(nn, d)
         I bits:length(rem) > 0
            L.break
         nn = div
         e++

      I e > 0
         pf [+]= PExp(d, e)
         s = isqrt(nn)

      d += 2

   R pf

F moBachShallit58(BigInt a, BigInt n; pf)
   V n1 = n - 1
   V mo = BigInt(1)
   L(pe) pf
      V y = n1 I/ pow(pe.prime, BigInt(pe.exp))
      V o = 0
      V x = pow(a, y, n)
      L x > 1
         x = pow(x, pe.prime, n)
         o++
      V o1 = pow(pe.prime, BigInt(o))
      o1 I/= gcd(mo, o1)
      mo *= o1
   R mo

F moTest(a, n)
   I bits:length(a) < 100
      print(‘ord(’a‘)’, end' ‘’)
   E
      print(‘ord([big])’, end' ‘’)
   print(‘ mod ’n‘ = ’moBachShallit58(a, n, factor(n - 1)))

moTest(37, 3343)

moTest(pow(BigInt(10), 100) + 1, 7919)
moTest(pow(BigInt(10), 1000) + 1, 15485863)
moTest(pow(BigInt(10), 10000) - 1, BigInt(22801763489))

moTest(1511678068, 7379191741)
moTest(BigInt(‘3047753288’), BigInt(‘2257683301’))