my algorithm

_DEFINE A-Z AS _UNSIGNED _INTEGER64

N = 13537 * 21019 'Semiprime to factorise

prs = 16 'prime parity bits in power vector

'generate small primes
DIM pr(255)
i = 0: a = 2
WHILE i < prs: ok = 1
  FOR j = 1 TO i
    IF a MOD pr(j) = 0 THEN ok = 0
  NEXT
  IF ok = 1 THEN i = i + 1: pr(i) = a
  a = a + 1
WEND

'look for squares composed solely of these prs primes
DIM hi(65535), hs(65535)
i0 = INT(SQR(N) + 1): s = i0 * i0 - N: a = 2 * i0 + 1
PRINT i0; "->"; N, i0; "^ 2 ="; s: SLEEP
i = i0: WHILE i < N
  s0 = s
  pv = 0: p = 1 'set bits of pv equal to the parity of prime occurance
  FOR j = 1 TO prs
    WHILE s0 MOD pr(j) = 0: s0 = s0 \ pr(j): pv = pv XOR p: WEND
    p = p * 2
  NEXT

  IF s0 = 1 THEN 's is prs-smooth, pv is it's power vector
    PRINT i; "^ 2 ="; s; TAB(50); bin$(pv, prs)
    IF pv = 0 THEN 'direct root
      im = i: sm = SQR(s): GOSUB roots
    END IF
    'NB: looking for matching pv only is naive.
    'the great power of this method is realised when we xor elements of our
    'table together to find matching combos. ie, where xor(combo)=0
    IF hi(pv) > 0 THEN 'matching root
      im = (i * hi(pv)) MOD N: sm = SQR(s * hs(pv)): GOSUB roots
    ELSE
      hi(pv) = i: hs(pv) = s
    END IF
  END IF
  s = s + a: IF s >= N THEN s = s - N
  a = a + 2: IF a >= N THEN a = a - N
  i = i + 1
WEND
END

roots:
PRINT "after"; i - i0; "iterations:";
s1 = im + sm: GOSUB gcd: PRINT s2;
s1 = im - sm: GOSUB gcd: PRINT "*"; s2; TAB(50); im; "+/-"; sm
IF s2 > 1 AND s2 < N THEN END
RETURN

gcd:
s2 = N
WHILE s1 > 0
  s0 = s1: s1 = s2 - (s2 \ s1) * s1: s2 = s0
WEND
RETURN

FUNCTION bin$ (pv, bits)
  x = pv: bin$ = ""
  FOR k = 1 TO bits
    IF x MOD 2 = 0 THEN ch$ = "0" ELSE ch$ = "1"
    bin$ = ch$ + bin$: x = x \ 2
  NEXT
END FUNCTION